text:
module PSO

using Printf, CUDA
using PV

export initialize, evolve!, pso

"""
    initialize(x_lb::AbstractVector{T}, x_ub::AbstractVector{T},
        d::Integer, np::Integer) where T <: AbstractFloat

Initialize a population of `np` particles with random positions `X` and velocities `V`.
`d` is the dimension (i.e., number of variables). `x_lb` and `x_ub` are lower bound and upper bound
of variables.

Two matrices `X` and `V` of size `d × np` are returned, each column of which denotes a position
or velocity.
"""
function initialize(x_lb::AbstractVector{T}, x_ub::AbstractVector{T},
        d::Integer, np::Integer) where T <: AbstractFloat
    @assert d >= 1 && np >= 1 && length(x_lb) == d && length(x_ub) == d
    dx = x_ub .- x_lb
    X = rand(T, d, np) .* dx .+ x_lb
    V = rand(T, d, np) .* 2 .* dx .- dx
    return X, V
end

@inline function _clamp!(x::AbstractVector, l::AbstractVector, u::AbstractVector)
    @inbounds for i in eachindex(x)
        x[i] = clamp(x[i], l[i], u[i])
    end
    nothing
end

@inline function bounce_back!(x, v, l, u)
    T = eltype(x)
    @inbounds for i in eachindex(x)
        t = x[i] + v[i]
        if t < l[i]
            t = l[i] + rand(T) * (x[i] - l[i])
        elseif t > u[i]
            t = x[i] + rand(T) * (u[i] - x[i])
        end
        x[i] = t
    end
end

"""
    evolve!(f::Function, X::Matrix{T}, V::Matrix{T}, x_lb::AbstractVector{T},
        x_ub::AbstractVector{T}, niter::Integer=100; w, c1, c2, log=false) where T <: AbstractFloat

Evolve particles with initial positions `X` and velocities `V` for at most `niter` iterations
to minimize the function `f`. `x_lb` and `x_ub` are lower bound and upper bound of variables.

`w, c1, c2` are parameters of PSO. If `log` is `true`, then progress information is printed in each
iteration.

The best particle's position, its fitness, and the best fitness history for `niter` iterations are
returned.

This method is designed to work on CPU to validate the correctness of PSO.
"""
function evolve!(f::Function, X::Matrix{T}, V::Matrix{T}, x_lb::AbstractVector{T},
        x_ub::AbstractVector{T}, niter::Integer=100;
        w, c1, c2, log=false) where T <: AbstractFloat
    @assert size(X) == size(V)
    @assert size(X, 1) == length(x_lb) == length(x_ub)
    @assert w_min <= w_max
    np = size(X, 2)
    fitness = zeros(T, np)
    pbest = copy(X)
    pbest_fitness = fill(T(Inf), np)
    gbest = X[:, 1]
    gbest_fitness = Inf
    gbest_fitness_hist = zeros(T, niter)

    if log
        @printf("%-10s%-10s\n", "iter", "sse")
    end

    ws = range(w_max, w_min; length=niter_w_decay)  # inertia weight decay

    w = w_max
    for iter = 1:niter
        # w = iter <= niter_w_decay ? ws[iter] : w_min
        # w *= 0.998
        w = max(w, w_min)
        # evaluate fitness and update pbest
        @inbounds for i = 1:np
            fitness[i] = f(@view X[:, i])
            if fitness[i] <= pbest_fitness[i]
                pbest_fitness[i] = fitness[i]
                pbest[:, i] .= @view X[:, i]
            end
        end
        # update gbest
        pbest_idx = argmin(pbest_fitness)
        @inbounds if pbest_fitness[pbest_idx] <= gbest_fitness
            gbest_fitness = pbest_fitness[pbest_idx]
            gbest .= @view pbest[:, pbest_idx]
        end
        gbest_fitness_hist[iter] = gbest_fitness
        if log
            @printf("%-10d%-10.5e%10f\n", iter, gbest_fitness,w)
        end
        # terminate?
        if iter == niter
            break
        end
        # update velocity/position: parallel
        @inbounds @views for i = 1:np
            V[:, i] .= w .* V[:, i] .+ c1 .* rand() .* (pbest[:, i] .- X[:, i]) .+
                        c2 .* rand() .* (gbest .- X[:, i])
            # X[:, i] .+= V[:, i]
            # _clamp!(X[:, i], x_lb, x_ub)
            bounce_back!(X[:, i], V[:, i], x_lb, x_ub)
        end
    end
    return gbest, gbest_fitness, gbest_fitness_hist
end

@inline function _copy!(dest::AbstractVector, src::AbstractVector)
    @inbounds for i = 1:length(src)
        dest[i] = src[i]
    end
end

"""
    evolve!(evaluator::Function, data::CuMatrix{T}, temperature::T,
        X::CuMatrix{T}, V::CuMatrix{T}, x_lb::CuVecOrMat{T},
        x_ub::CuVecOrMat{T}, niter::Integer=100;
        w, c1, c2, log=false, threads::Integer=160,
        min_delta=0.0f0, patience=niter) where T <: AbstractFloat

This method is designed to work on GPU.
Refer to the CPU-based [`evolve!`](@ref) for details. Note the two parameters for early termination:
- `min_delta`: a threshold to decide whether the best fitness (i.e., fitness of *gbest*) has reduced
- `patience`: if the best fitness has not reduced for continuous `patience` iterations, terminate.

A `NamedTuple` is returned with three fields
- `sse`: SSE of the best solution
- `solution`: best solution (*gbest*)
- `stopped_iter`: the iteration at which PSO is terminated
"""
function evolve!(evaluator::Function, data::CuMatrix{T}, temperature::T,
        X::CuMatrix{T}, V::CuMatrix{T}, x_lb::CuVecOrMat{T},
        x_ub::CuVecOrMat{T}, niter::Integer=100;
        w, c1, c2, log=false, threads::Integer=160,
        min_delta=0.0f0, patience=niter) where T <: AbstractFloat
    np = size(X, 1)
    nv = size(X, 2)
    fitness = CUDA.zeros(T, np)
    pbest = copy(X)
    pbest_fitness = CUDA.fill(T(Inf), np)
    gbest = X[1, :]  # copied
    gbest_fitness = T(Inf)
    gbest_fitness_hist = CUDA.zeros(T, niter)
    R = CUDA.rand(T, 4 * np, niter)   # pre-generate random numbers

    function evaluate!(iter::Integer)
        i = threadIdx().x + (blockIdx().x - 1) * blockDim().x
        if i > np
            return nothing
        end
        @inbounds fitness[i] = fit = evaluator(data, temperature, @view X[i, :])
        @inbounds if fit <= pbest_fitness[i]
            pbest_fitness[i] = fit
            @views _copy!(pbest[i, :], X[i, :])
        end
        nothing
    end

    function update_vel_pos!(iter)
        i = threadIdx().x + (blockIdx().x - 1) * blockDim().x
        if i > np
            return nothing
        end
        # update velocity and position
        @inbounds r1 = R[i, iter]
        @inbounds r2 = R[i + np, iter]
        @inbounds for j = 1:nv  # update each component individually
            V[i, j] = new_v = w * V[i, j] + c1 * r1 * (pbest[i, j] - X[i, j]) + c2 * r2 * (gbest[j] - X[i, j])
            old_x = X[i, j]
            t = old_x  + new_v
            l = x_lb[j]
            u = x_ub[j]
            if t < l
                t = l + R[i + 2 * np, iter] * (old_x - l)
            elseif t > u
                t = old_x + R[i + 3 * np, iter] * (u - old_x)
            end
            X[i, j] = t
        end
        nothing
    end

    # kernel config
    blocks = np ÷ threads + (np % threads > 0)
    if log
        @show threads, blocks
        @printf("%-10s%-10s\n", "iter", "sse")
    end

    wait = 0
    stopped_iter = 0

    for iter = 1:niter
        # fitness
        @cuda threads = threads blocks = blocks evaluate!(iter)
        # update gbest： parallel reduction
        pbest_idx = argmin(pbest_fitness)
        min_pbest_fitness = pbest_fitness[pbest_idx]
        if min_pbest_fitness <= gbest_fitness
            if gbest_fitness - min_pbest_fitness >= min_delta
                wait = 0
            else
                wait += 1
            end
            gbest_fitness = min_pbest_fitness
            @inbounds _copy!(gbest, @view pbest[pbest_idx, :])
        end
        if log
            @printf("%-10d%-10.5e\n", iter, gbest_fitness)
        end
        # terminate?
        if wait >= patience
            stopped_iter = iter
            if log
                @printf("Early stopped at iteration %d\n", iter)
            end
            break
        end
        if iter == niter
            stopped_iter = niter
            break
        end

        # update velocity/position: parallel
        @cuda threads = threads blocks = blocks update_vel_pos!(iter)
    end
    return (sse = gbest_fitness, solution = Array(gbest), stopped_iter = stopped_iter)
end

"""
    pso(data::AbstractMatrix{T}, temperature::T, model::Symbol, x_lb::AbstractVector{T},
        x_ub::AbstractVector{T}; niter::Integer=400, np::Integer=1600,
        w::T, c1::T, c2::T, log=false, threads::Integer=160,
        min_delta::T=0.0, patience=niter + 1) where T <: AbstractFloat

Parallel PSO for PV parameter identification. Each column of `data` is a V-I point. `model` can be either
`:SDM` or `:DDM`. `x_lb` and `x_ub` denote the lower bound and upper bound of each variable to
confine the search range.
- `niter`: number of iterations
- `np`: number of particles
- `w, c1, c2`: PSO velocity update parameters
- `log`: print progress information if `true`
- `threads`: number of threads per block in CUDA kernel
- `min_delta, patience`: early stopping. Similar to *keras* API:
    https://www.tensorflow.org/api_docs/python/tf/keras/callbacks/EarlyStopping
"""
function pso(data::AbstractMatrix{T}, temperature::T, model::Symbol, x_lb::AbstractVector{T},
        x_ub::AbstractVector{T}; niter::Integer=400, np::Integer=1600,
        w::T, c1::T, c2::T, log=false, threads::Integer=160,
        min_delta::T=0.0, patience=niter + 1) where T <: AbstractFloat
    @assert model in (:SDM, :DDM)
    nv = model === :SDM ? 5 : 7
    @assert length(x_lb) == length(x_ub) == nv
    @assert size(data, 1) == 2
    @assert threads % 32 == 0

    # initialize: each row is a particle (for GPU)
    dx = reshape(x_ub .- x_lb, 1, :) |> CuArray
    x_lb = reshape(x_lb, 1, :) |> CuArray
    x_ub = reshape(x_ub, 1, :) |> CuArray
    X = CUDA.rand(T, np, nv) .* dx .+ x_lb
    V = CUDA.rand(T, np, nv) .* 2 .* dx .- dx

    # evolve particles
    if model === :SDM
        return evolve!(PV.sse_sdm, data |> CuArray, temperature,
            X, V, x_lb |> CuArray, x_ub |> CuArray, niter;
            w, c1, c2, log, threads, min_delta, patience)
    else
        return evolve!(PV.sse_ddm, data |> CuArray, temperature,
            X, V, x_lb |> CuArray, x_ub |> CuArray, niter;
            w, c1, c2, log, threads, min_delta, patience)
    end
end

end # module

max_stars_repo_name: ShuhuaGao/PV-PPSO

max_issues_repo_name: ShuhuaGao/PV-PPSO

max_forks_repo_name: ShuhuaGao/PV-PPSO


text:
#!$JULIA_HOME/julia
using PyConversions
using Documenter

makedocs()

fname = joinpath("build", "README.md")
if isfile(fname)
    info("Updating README...")
    cp(fname, joinpath("..", "README.md"), remove_destination=true)
end

ext: jl


text:
using DataFrames, CSV
import DataFrames.DataFrame
#bpath  = "./benchmark/out/"
bpath  = "./benchmark/out/64/"

function collate_benchmarks(bpath, out)
    bfiles = readdir(bpath)
    bfiles = setdiff(bfiles,["64", "old"])

    @time bdf = CSV.read.(bpath .* bfiles)

    hehe = r"(?<algo>sp|mrs|srg|srs)(?<tries>[0-9]+)"

    hehe1  = match.(hehe, bfiles)

    hehe2  = [hehe2.captures for hehe2 in hehe1]

    hehe3 = [(hehe3[1], eval(parse(hehe3[2])) , bdf1[1,1]) for (hehe3, bdf1) in zip(hehe2,bdf)]

    algo = [h[1] for h in hehe3]
    n = [h[2] for h in hehe3]
    elapsed = [h[3] for h in hehe3]

    hehe4 = DataFrame(algo = algo, n = n, elapsed =  elapsed)

    hehe5 = DataFrames.by(hehe4, [:algo, :n], df -> mean(df[:elapsed]))

    sort!(hehe5, cols = [:algo,:n])

    CSV.write(out,hehe5)

    hehe5
end

res64 = collate_benchmarks(bpath, "hehe6.csv")

using Gadfly
plot(res64, x = :n, y = :x1, color = :algo, Geom.point)


text:
# Does not test leafsize
@testset "inrange" begin
    @testset "metric" for metric in [Euclidean()]
        @testset "tree type" for TreeType in trees_with_brute
            function test(data)
                tree = TreeType(data, metric; leafsize=2)
                dosort = true

                idxs = inrange(tree, [1.1, 1.1, 1.1], 0.2, dosort)
                @test idxs == [8] # Only corner 8 at least 0.2 distance away from [1.1, 1.1, 1.1]

                idxs = inrange(tree, [0.0, 0.0, 0.5], 0.6, dosort)
                @test idxs == [1, 2] # Corner 1 and 2 at least 0.6 distance away from [0.0, 0.0, 0.5]

                idxs = inrange(tree, [0, 0, 0], 0.6, dosort)
                @test idxs == [1]

                X = [0.0 0.0; 0.0 0.0; 0.5 0.0]
                idxs1 = inrange(tree, X, 0.6, dosort)
                idxs2 = inrange(tree, view(X,:,1:2), 0.6, dosort)
                @test idxs1 == idxs2
                @test idxs1[1] == [1,2]
                @test idxs1[2] == [1]

                idxs = inrange(tree, [SVector{3,Float64}(0.0, 0.0, 0.5), SVector{3,Float64}(0.0, 0.0, 0.0)], 0.6, dosort)
                @test idxs[1] == [1,2]
                @test idxs[2] == [1]

                idxs = inrange(tree, [0.33333333333, 0.33333333333, 0.33333333333], 1, dosort)
                @test idxs == [1, 2, 3, 5]

                idxs = inrange(tree, [0.5, 0.5, 0.5], 0.2, dosort)
                @test idxs == []

                idxs = inrange(tree, [0.5, 0.5, 0.5], 1.0, dosort)
                @test idxs == [1, 2, 3, 4, 5, 6, 7, 8]

                @test_throws ArgumentError inrange(tree, rand(3), -0.1)
                @test_throws ArgumentError inrange(tree, rand(5), 1.0)

                empty_tree = TreeType(rand(3,0), metric)
                idxs = inrange(empty_tree, [0.5, 0.5, 0.5], 1.0)
                @test idxs == []

                one_point_tree = TreeType([0.5, 0.5, 0.5], metric)
                idxs = inrange(one_point_tree, data, 1.0)
                @test idxs == repeat([[1]], size(data, 2))
            end
            data = [0.0 0.0 0.0 0.0 1.0 1.0 1.0 1.0;
                    0.0 0.0 1.0 1.0 0.0 0.0 1.0 1.0;
                    0.0 1.0 0.0 1.0 0.0 1.0 0.0 1.0] # 8 node cube
            test(data)
            test(view(data, :, :))
        end
    end
end


text:
"""
    Diagnostics

Accumulate mean fields and covariance statistics on the computational grid.

"""
module Diagnostics

using Dates
using FileIO
using JLD2
using MPI
using StaticArrays

using ..DGmethods
using ..DGmethods: num_state, vars_state, num_aux, vars_aux, vars_diffusive, num_diffusive
using ..Mesh.Topologies
using ..Mesh.Grids
using ..MPIStateArrays
using ..VariableTemplates

Base.@kwdef mutable struct Diagnostic_Settings
    mpicomm::MPI.Comm = MPI.COMM_WORLD
    dg::Union{Nothing,DGModel} = nothing
    Q::Union{Nothing,MPIStateArray} = nothing
    starttime::Union{Nothing,String} = nothing
    outdir::Union{Nothing,String} = nothing
end

const Settings = Diagnostic_Settings()

include("diagnostic_vars.jl")

"""
    visitQ()

Helper macro to iterate over the DG grid. Generates the needed loops
and indices: `eh`, `ev`, `e`, `k,`, `j`, `i`, `ijk`.
"""
macro visitQ(nhorzelem, nvertelem, Nqk, Nq, expr)
    return esc(
        quote
            for eh in 1:nhorzelem
                for ev in 1:nvertelem
                    e = ev + (eh - 1) * nvertelem
                    for k in 1:Nqk
                        for j in 1:Nq
                            for i in 1:Nq
                                ijk = i + Nq * ((j-1) + Nq * (k-1))
                                $expr
                            end
                        end
                    end
                end
        end
    end)
end

# Helpers to extract data from `Q`, etc.
function extract_state(dg, localQ, ijk, e)
    bl = dg.balancelaw
    FT = eltype(localQ)
    nstate = num_state(bl, FT)
    l_Q = MArray{Tuple{nstate},FT}(undef)
    for s in 1:nstate
        l_Q[s] = localQ[ijk,s,e]
    end
    return Vars{vars_state(bl, FT)}(l_Q)
end
function extract_aux(dg, auxstate, ijk, e)
    bl = dg.balancelaw
    FT = eltype(auxstate)
    nauxstate = num_aux(bl, FT)
    l_aux = MArray{Tuple{nauxstate},FT}(undef)
    for s in 1:nauxstate
        l_aux[s] = auxstate[ijk,s,e]
    end
    return Vars{vars_aux(bl, FT)}(l_aux)
end
function extract_diffusion(dg, localdiff, ijk, e)
    bl = dg.balancelaw
    FT = eltype(localdiff)
    ndiff = num_diffusive(bl, FT)
    l_diff = MArray{Tuple{ndiff},FT}(undef)
    for s in 1:ndiff
        l_diff[s] = localdiff[ijk,s,e]
    end
    return Vars{vars_diffusive(bl, FT)}(l_diff)
end

include("AtmosDiagnostics.jl")

"""
    init(mpicomm, dg, Q, starttime, outdir)

Initialize the diagnostics collection module.
"""
function init(mpicomm::MPI.Comm,
              dg::DGModel,
              Q::MPIStateArray,
              starttime::String,
              outdir::String)
    Settings.mpicomm = mpicomm
    Settings.dg = dg
    Settings.Q = Q
    Settings.starttime = starttime
    Settings.outdir = outdir

    init(dg.balancelaw, mpicomm, dg, Q, starttime)
end

"""
    collect(currtime)

Collect various diagnostics specific to the balance law and write them to a
JLD2 file, indexed by `currtime`, in the output directory specified in
`init()`.
"""
function collect(currtime)
    mpicomm = Settings.mpicomm
    mpirank = MPI.Comm_rank(mpicomm)

    current_time = string(currtime)

    # make sure this time step is not already recorded
    docollect = [false]
    if mpirank == 0
        try
            jldopen(joinpath(Settings.outdir,
                    "diagnostics-$(Settings.starttime).jld2"), "r") do file
                diagnostics = file[current_time]
                @warn "diagnostics for time step $(current_time) already collected"
        end
        catch e
            docollect[1] = true
        end
    end
    MPI.Bcast!(docollect, 0, mpicomm)
    if !docollect[1]
        return nothing
    end

    # collect diagnostics
    diagnostics = collect(Settings.dg.balancelaw, currtime)

    # write results
    if mpirank == 0 && diagnostics !== nothing
        jldopen(joinpath(Settings.outdir,
                "diagnostics-$(Settings.starttime).jld2"), "a+") do file
            file[current_time] = diagnostics
        end
    end

    return nothing
end # function collect

end # module Diagnostics

max_stars_repo_name: favba/CLIMA

max_issues_repo_name: favba/CLIMA

max_forks_repo_name: favba/CLIMA


text:
using Combinatorics
using ArgCheck

"""
    Tensor{T<:Number,N}

Wrapper type for native Julia arrays. Provides most of the standard array functionality and
extends it by arithmetic and algebraic tensor operations.

# Examples
```jldoctest
julia> t1 = Tensor([1 2; 3 4])
2×2 Tensor{Int64, 2}:
 1  2
 3  4

julia> 3.14 * t1
2×2 Tensor{Float64, 2}:
 3.14   6.28
 9.42  12.56

julia> t2 = Tensor(Float32[1.1, -2.2, 3.3])
3-element Tensor{Float32, 1}:
  1.1
 -2.2
  3.3

julia> t1 ⊗ t2
2×2×3 Tensor{Float32, 3}:
[:, :, 1] =
 1.1  2.2
 3.3  4.4

[:, :, 2] =
 -2.2  -4.4
 -6.6  -8.8

[:, :, 3] =
 3.3   6.6
 9.9  13.2
```
"""
struct Tensor{T<:Number,N}
    data::AbstractArray{T,N}
    Tensor(arr::AbstractArray{T,N}) where {T,N} = new{T,N}(arr)
    Tensor(val::T) where {T} = new{T,0}(fill(val))
end

# extend array property methods to tensors
Base.size(t::Tensor) = size(t.data)
Base.size(t::Tensor, dim::Integer) = size(t.data, dim)
Base.ndims(t::Tensor) = ndims(t.data)
Base.length(t::Tensor) = length(t.data)
Base.eltype(t::Tensor) = eltype(t.data)

# extent array equality to tensors
Base.:(==)(t1::Tensor, t2::Tensor) = t1.data == t2.data
function Base.isapprox(t1::Tensor, t2::Tensor; kwargs...)
    return isapprox(t1.data, t2.data; kwargs...)
end

# type conversions
(::Type{Tensor{T}})(t::Tensor) where {T} = convert(Tensor{T}, t)
Base.convert(::Type{Tensor{T}}, t::Tensor) where {T} = Tensor(convert(Array{T}, t.data))
(::Type{T})(t::Tensor{T,0}) where {T} = t.data[1]
Base.convert(::Type{T}, t::Tensor{T,0}) where {T} = t.data[1]

# extend array size manipulation methods to tensors
Base.reshape(t::Tensor, dims...) = Tensor(reshape(t.data, dims...))
Base.permutedims(t::Tensor{T,1}) where {T} = Tensor(permutedims(t.data))
Base.permutedims(t::Tensor{T,2}) where {T} = Tensor(permutedims(t.data))
Base.permutedims(t::Tensor, perm) = Tensor(permutedims(t.data, perm))

# allow indexing tensors like arrays
@inline function Base.getindex(t::Tensor, inds...)
    @argcheck ndims(t) > 0 BoundsError(t)
    @boundscheck checkbounds(t, inds...)
    val = getindex(t.data, inds...)
    return Tensor(typeof(val) <: AbstractArray ? val : fill(val))
end

@inline function Base.setindex!(t1::Tensor, t2::Tensor, inds...)
    @boundscheck checkbounds(t1, inds...)
    val = ndims(t2) == 0 ? t2.data[1] : t2.data
    return setindex!(t1.data, val, inds...)
end

@inline function Base.setindex!(t::Tensor, val::S, inds...) where {S<:Number}
    @boundscheck checkbounds(t, inds...)
    return setindex!(t.data, val, inds...)
end

function Base.checkbounds(::Type{Bool}, t::Tensor, inds...)
    return checkbounds(Bool, t.data, inds...)
end

function Base.checkbounds(t::Tensor, inds...)
    if !checkbounds(Bool, t::Tensor, inds...)
        throw(BoundsError(t, inds))
    end
end

Base.axes(t::Tensor, inds...) = axes(t.data, inds...)
Base.eachindex(t::Tensor) = eachindex(t.data)
Base.CartesianIndices(t::Tensor) = CartesianIndices(t.data)

# Base.show overloads
Base.show(io::IO, t::Tensor{T,0}) where {T} = print(io, "$(typeof(t))($(t.data[1]))")
# TODO: Add ellipsis handling for large arrays
function Base.show(io::IO, ::MIME"text/plain", t::Tensor)
    dims = ndims(t)
    shape = if dims == 0
        "scalar"
    elseif dims == 1
        "$(length(t))-element"
    else
        join(size(t), "×")
    end
    print(io, "$shape $(typeof(t)):\n")
    tmp_io = IOBuffer()
    show(tmp_io, "text/plain", t.data)
    data_str = String(take!(tmp_io))
    data_str = join(split(data_str, "\n")[2:end], "\n")
    return print(io, data_str)
end

# tensor hash should be slightly different than the hash of the underlying data
Base.hash(t::Tensor) = hash(hash(t.data) + hash("t"))

# TODO: probably relocate these definitions to another (new?) file
function δ(::Val{N}) where {N}
    id = zeros(Bool, N, N)
    for i in axes(id, 1)
        @inbounds id[i, i] = true
    end
    return Tensor(id)
end
δ(N::Integer) = δ(Val(N))

function ϵ(::Val{N}) where {N}
    lc = zeros(Int8, ntuple(d -> N, Val(N)))
    for perm in permutations(1:N)
        @inbounds lc[perm...] = levicivita(perm)
    end
    return Tensor(lc)
end
ϵ(N::Integer) = ϵ(Val(N))


text:
@testset "Geometry" begin
    tests = ["vector",
             "point",
             "plane",
             "axisalignedbox",
             "polytopes/polytope",
             "polytopes/interpolate",
             "polytopes/jacobian",
             "polytopes/faces",
             "polytopes/edges",
             "polytopes/measure",
             "polytopes/triangulate",
            ]
    for t in tests
      include("$(t).jl")
    end
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
"""
```
t = find_threshold(MinimumIntermodes(),  histogram, edges)
```

Under the assumption that the histogram is bimodal the histogram is smoothed
using a length-3 mean filter until two modes remain. The threshold is then set
to the minimum value between the two modes.

# Output

Returns `t`, a real number that specifies the threshold.

# Arguments

The function arguments are described in more detail below.

##  `histogram`

An `AbstractArray` storing the frequency distribution.

##  `edges`

An `AbstractRange` specifying how the intervals for the frequency distribution
are divided.

# Example

Compute the threshold for the "camerman" image in the `TestImages` package.

```julia
using TestImages, ImageContrastAdjustment, HistogramThresholding

img = testimage("cameraman")
edges, counts = build_histogram(img,256)
#=
  The `counts` array stores at index 0 the frequencies that were below the
  first bin edge. Since we are seeking a threshold over the interval
  partitioned by `edges` we need to discard the first bin in `counts`
  so that the dimensions of `edges` and `counts` match.
=#
t = find_threshold(MinimumIntermodes(),counts[1:end], edges)
```

# Reference

1. C. A. Glasbey, “An Analysis of Histogram-Based Thresholding Algorithms,” CVGIP: Graphical Models and Image Processing, vol. 55, no. 6, pp. 532–537, Nov. 1993. [doi:10.1006/cgip.1993.1040](https://doi.org/10.1006/cgip.1993.1040)
"""
function find_threshold(algorithm::MinimumIntermodes, histogram::AbstractArray, edges::AbstractRange)
    histogram_local = smooth_histogram(histogram, 1000)
    indices = find_maxima_indices(histogram_local)
    if length(indices) != 2
        return 0
    end
    t = 0
    min_value = typemax(Float64)
    for i = min(indices[1],indices[2]):max(indices[1],indices[2])
        if histogram_local[i] < min_value
            min_value = histogram_local[i]
            t = i
        end
    end
    index = round(Int, t)
    return edges[index]
end


text:
# Compute Wilcoxon signed rank tests
# Author: Thomas Minier

using HypothesisTests
using DataFrames
using CalliPlots

# load data, compute mean execution time and then convert them to arrays
time_ref = Array(computeMean(
    :time,
    "/home/stagiaire-gdd/experiments-ulysses/reference-tpf-1servers-run1/execution_times_ref_1servers.csv"
)[:time])

ulysses_1s = Array(computeMean(
    :time,
    "/home/stagiaire-gdd/experiments-ulysses/ulysses-1servers-run1/execution_times_ulysses_1servers.csv",
    "/home/stagiaire-gdd/experiments-ulysses/ulysses-1servers-run2/execution_times_ulysses_1servers.csv",
    "/home/stagiaire-gdd/experiments-ulysses/ulysses-1servers-run3/execution_times_ulysses_1servers.csv"
)[:time])

ulysses_2s = Array(computeMean(
    :time,
    "/home/stagiaire-gdd/experiments-ulysses/ulysses-2servers-run1/execution_times_ulysses_2servers.csv",
    "/home/stagiaire-gdd/experiments-ulysses/ulysses-2servers-run2/execution_times_ulysses_2servers.csv",
    "/home/stagiaire-gdd/experiments-ulysses/ulysses-2servers-run3/execution_times_ulysses_2servers.csv"
)[:time])

ulysses_3s = Array(computeMean(
    :time,
    "/home/stagiaire-gdd/experiments-ulysses/ulysses-3servers-run1/execution_times_ulysses_3servers.csv",
    "/home/stagiaire-gdd/experiments-ulysses/ulysses-3servers-run2/execution_times_ulysses_3servers.csv",
    "/home/stagiaire-gdd/experiments-ulysses/ulysses-3servers-run3/execution_times_ulysses_3servers.csv"
)[:time])

pval_1s = pvalue(SignedRankTest(time_ref, ulysses_1s), tail=:left)
pval_2s = pvalue(SignedRankTest(time_ref, ulysses_2s), tail=:left)
pval_3s = pvalue(SignedRankTest(time_ref, ulysses_3s), tail=:left)
results = DataFrame(servers = [1, 2, 3], pvalue = [pval_1s, pval_2s, pval_3s])
println(results)

max_stars_repo_name: Callidon/ulysses-tpf

max_issues_repo_name: Callidon/ulysses-tpf

max_forks_repo_name: Callidon/ulysses-tpf


text:
println("==============flag! solver.jl======================================")

#=
A solver for MINLP problems using a NLP solver and Branch and Bound
=#

function get_default_options()
    nl_solver                           = nothing
    log_levels                          = [:Options,:Table,:Info]
    silent                              = false
    atol                                = 1e-6
    num_resolve_root_relaxation         = 3
    branch_strategy                     = :StrongPseudoCost
    gain_mu                             = 0.167 # http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.92.7117&rep=rep1&type=pdf
    # Strong branching
    strong_branching_perc               = 100
    strong_branching_nsteps             = 1
    strong_branching_time_limit         = 100
    strong_restart                      = true
    # Reliability branching
    reliability_branching_threshold     = 5 # reliability param
    reliability_branching_perc          = 25
    # Obj cuts
    incumbent_constr                    = false
    obj_epsilon                         = 0
    # :UserLimit
    time_limit                          = Inf
    mip_gap                             = 1e-4
    best_obj_stop                       = NaN
    solution_limit                      = 0
    all_solutions                       = false
    list_of_solutions                   = false
    # Parallel
    processors                          = 1
    # Traversing
    traverse_strategy                   = :BFS
    # Feasibility Pump
    feasibility_pump                    = true # changes to false if mip_solver not provided
    feasibility_pump_time_limit         = 60.0
    feasibility_pump_tolerance_counter  = 5
    tabu_list_length                    = 30
    num_resolve_nlp_feasibility_pump    = 1
    mip_solver                          = nothing
    allow_almost_solved                 = true
    allow_almost_solved_integral        = true
    registered_functions                = nothing

    # Only for testing
    force_parallel                      = false
    debug                               = false
    debug_write                         = false
    debug_file_path                     = "debug.json"

    fixed_gain_mu                       = false

    return SolverOptions(nl_solver,log_levels,silent,atol,num_resolve_root_relaxation,branch_strategy,gain_mu,
        strong_branching_perc,strong_branching_nsteps,strong_branching_time_limit,strong_restart,
        reliability_branching_threshold,reliability_branching_perc,
        incumbent_constr,obj_epsilon,time_limit,mip_gap,best_obj_stop,solution_limit,all_solutions,
        list_of_solutions,processors,traverse_strategy,
        feasibility_pump,feasibility_pump_time_limit,feasibility_pump_tolerance_counter,
        tabu_list_length,num_resolve_nlp_feasibility_pump,
        mip_solver, allow_almost_solved, allow_almost_solved_integral, registered_functions,
        force_parallel, debug, debug_write, debug_file_path, fixed_gain_mu)
end

function combine_options(options)
    branch_strategies = Dict{Symbol,Bool}()
    for strat in [:StrongPseudoCost, :PseudoCost, :Reliability, :MostInfeasible]
        branch_strategies[strat] = true
    end

    traverse_strategies = Dict{Symbol,Bool}()
    traverse_strategies[:BFS] = true
    traverse_strategies[:DFS] = true
    traverse_strategies[:DBFS] = true

    options_dict = Dict{Symbol,Any}()
    for kv in options
        if !in(kv[1], fieldnames(SolverOptions))
            if kv[1] == :strong_branching_approx_time_limit
                @warn "The option `strong_branching_approx_time_limit` got replaced by `strong_branching_time_limit`"
                @info "Setting strong_branching_time_limit = $(kv[2])"
                options_dict[:strong_branching_time_limit] = kv[2]
            else
                @warn "Option $(string(kv[1])) is not available"
            end
        else
            options_dict[kv[1]] = kv[2]
        end
    end
    if haskey(options_dict, :log_levels)
        if length(options_dict[:log_levels]) == 0
            options_dict[:log_levels] = Symbol[]
        end
    end

    defaults = get_default_options()

    # if gain mu is specified we shouldn't change it later
    if haskey(options_dict, :gain_mu)
        options_dict[:fixed_gain_mu] = true
    end

    # if time limit is set make sure to reduce the time limit for the feasibility pump
    if haskey(options_dict, :time_limit) && !haskey(options_dict, :feasibility_pump_time_limit)
        if options_dict[:time_limit] < defaults.feasibility_pump_time_limit
            setfield!(defaults, :feasibility_pump_time_limit, convert(Float64, options_dict[:time_limit]))
        end
    end

    for fname in fieldnames(SolverOptions)
        if haskey(options_dict, fname)
            # check that mip_solver is defined if feasibile pump should be used
            if fname == :feasibility_pump && options_dict[:feasibility_pump] == true
                if !haskey(options_dict, :mip_solver) || options_dict[:mip_solver] === nothing
                    @warn "The feasibility pump can only be used if a mip solver is defined."
                    options_dict[:feasibility_pump] = false
                end
            end

            # check branch strategy
            if fname == :branch_strategy
                if !haskey(branch_strategies, options_dict[fname])
                    error("Branching strategy "*string(options_dict[fname])* " is not supported")
                end
            end

             # check traverse strategy
             if fname == :traverse_strategy
                if !haskey(traverse_strategies, options_dict[fname])
                    error("Traverse strategy "*string(options_dict[fname])* " is not supported")
                end
            end

            if fieldtype(SolverOptions, fname) != typeof(options_dict[fname])
                options_dict[fname] = convert(fieldtype(SolverOptions, fname), options_dict[fname])
            end
            setfield!(defaults, fname, options_dict[fname])
        end
    end

    if defaults.allow_almost_solved_integral && !defaults.allow_almost_solved
        defaults.allow_almost_solved_integral = false
    end

    return defaults
end


text:
"""
struct to hold the parameters for Medlyn et al. (2011) stomatal
conductance model for CO₂.

# Arguments

- `g0`: intercept.
- `g1`: slope.
- `gs_min = 0.001`: residual conductance. We consider the residual conductance being different
 from `g0` because in practice `g0` can be negative when fitting real-world data.

# Useage

Then used for example as follows:
Gs = Medlyn(0.03,0.1)
gs_mod = gs(Gs,(Cₛ = 400.0, VPD = 1.5))
Gₛ = Gs.g0 + gs_mod * A

# Examples

```julia
meteo = Atmosphere(T = 20.0, Wind = 1.0, P = 101.3, Rh = 0.65)

leaf = LeafModels(stomatal_conductance = Medlyn(0.03, 12.0),
            A = A, Cₛ = 380.0, Dₗ = meteo.VPD)
gs(leaf,meteo)
```

# References

Medlyn, Belinda E., Remko A. Duursma, Derek Eamus, David S. Ellsworth, I. Colin Prentice,
Craig V. M. Barton, Kristine Y. Crous, Paolo De Angelis, Michael Freeman, et Lisa Wingate.
2011. « Reconciling the optimal and empirical approaches to modelling stomatal conductance ».
Global Change Biology 17 (6): 2134‑44. https://doi.org/10.1111/j.1365-2486.2010.02375.x.

"""
struct Medlyn{T} <: AbstractGsModel
    g0::T
    g1::T
    gs_min::T
end

function Medlyn(g0,gs,gs_min)
    Medlyn(promote(g0,gs,gs_min))
end

Medlyn(g0,g1) = Medlyn(g0,g1,oftype(g0,0.001))

Medlyn(;g0,g1) = Medlyn(g0,g1,oftype(g0,0.001))

function inputs(::Medlyn)
    (:Dₗ,:Cₛ,:A)
end

function outputs(::Medlyn)
    (:Gₛ,)
end

Base.eltype(x::Medlyn) = typeof(x).parameters[1]

"""
    gs_closure(leaf::LeafModels{I,E,A,<:Medlyn,S},meteo)

Stomatal closure for CO₂ according to Medlyn et al. (2011). Carefull, this is just a part of
the computation of the stomatal conductance.

The result of this function is then used as:

    gs_mod = gs_closure(leaf,meteo)

    # And then stomatal conductance (μmol m-2 s-1) calling [`gs`](@ref):
    Gₛ = leaf.stomatal_conductance.g0 + gs_mod * leaf.status.A

# Arguments

- `leaf::LeafModels{.,.,<:Fvcb,<:Medlyn,.}`: A [`LeafModels`](@ref) struct holding the parameters for
the model.
- `meteo`: meteorology structure, see [`Atmosphere`](@ref). Is not used in this model.

# Details

Use `variables()` on Medlyn to get the variables that must be instantiated in the `LeafModels` struct.

# Notes

- `Cₛ` is used instead of `Cₐ` because Gₛ is between the surface and the intercellular space. The conductance
between the atmosphere and the surface is accounted for using the boundary layer conductance
(`Gbc` in [`Monteith`](@ref)). Medlyn et al. (2011) uses `Cₐ` in their paper because they relate their models
to the measurements made at leaf level, with a well-mixed chamber where`Cₛ ≈ Cₐ`.
- `Dₗ` is forced to be >= 1e-9 because it is used in a squared root. It is prefectly acceptable to
get a negative Dₗ when leaves are re-hydrating from air. Cloud forests are the perfect example.
See *e.g.*: Guzmán‐Delgado, P, Laca, E, Zwieniecki, MA. Unravelling foliar water uptake pathways:
The contribution of stomata and the cuticle. Plant Cell Environ. 2021; 1– 13.
https://doi.org/10.1111/pce.14041

# Examples

```julia
meteo = Atmosphere(T = 20.0, Wind = 1.0, P = 101.3, Rh = 0.65)

leaf = LeafModels(stomatal_conductance = Medlyn(0.03, 12.0),
            Cₛ = 380.0, Dₗ = meteo.VPD)

gs_mod = gs_closure(leaf, meteo)

A = 20 # example assimilation (μmol m-2 s-1)
Gs = leaf.stomatal_conductance.g0 + gs_mod * A

# Or more directly using `gs()`:

leaf = LeafModels(stomatal_conductance = Medlyn(0.03, 12.0),
            A = A, Cₛ = 380.0, Dₗ = meteo.VPD)
gs(leaf,meteo)
```

# References

Medlyn, Belinda E., Remko A. Duursma, Derek Eamus, David S. Ellsworth, I. Colin Prentice,
Craig V. M. Barton, Kristine Y. Crous, Paolo De Angelis, Michael Freeman, et Lisa Wingate.
2011. « Reconciling the optimal and empirical approaches to modelling stomatal conductance ».
Global Change Biology 17 (6): 2134‑44. https://doi.org/10.1111/j.1365-2486.2010.02375.x.
"""
function gs_closure(leaf::LeafModels{I,E,A,<:Medlyn,S},meteo) where {I,E,A,S}
    (1.0 + leaf.stomatal_conductance.g1 / sqrt(max(1e-9,leaf.status.Dₗ))) / leaf.status.Cₛ
end


text:
# @doc """
# x=hard_sphere(y,H,Nc)
#   Solve the problem argmin_x ||y-Hx||, where the elements of x are
#   integers.  The input vector y is of dimension N, with H of dimension N
#   by M, and the returned vector x of dimension M and comes from the
#   integers 0:Nc-1.
# X=hard_sphere(Y,H,Nc)
#   The input may be a matrix Y of dimension N by Ns, in which case
#   the problem is solved for each column of Y, with the solutions in
#   the columns of X.
#
# Examples:
#   X = hard_sphere([1 2]', [1 2; 3 4],2)
#   X = hard_sphere(rand(0:20,2,15), [1 2; 3 4],2)
# """ ->
function hard_sphere{Td}(Y::Array{Td,2},H::Array{Td,2},Nc)
    (N,M) = size(H);

    if length(Nc)==1;
        Nc = ones(M,1)*Nc;
    end
    (Q,R) = qr(H);
    (N,Ns) = size(Y);
    Xh = zeros(Int,M,Ns);
    for ns = 1:Ns
        yp = (Q'*Y[:,ns] + R*ones(M,1).*(Nc-1))/2;
        xp = algII_smart(yp,R,Nc);
        #    Xh[:,ns] = xp*2 -(Nc-1);
        Xh[:,ns] = xp
    end

    return Xh
end

#########################################################################%
function algII_smart(yp,R,Qc)
    # xh = algII_smart(yp,R,Qc)
    #
    # Find the closest xh for yp=R*xh + w, where w is some unknown noise and
    # x(i) is from the range [0,...,Qc(i)-1]. This implements the "Alg II-smart"
    # hard-decision sphere decoder as in Damen, El Gamal, and Caire, Trans It
    # 03.
    (N,M) = size(R);

    xh = zeros(M,1); # the best solution
    xp = zeros(M,1); # work vector (a possible solution)
    ksi = zeros(M,1);
    TT = zeros(M,1);
    delta = zeros(M,1);

    # Step 1: initialization
    ix = M;
    TT[ix] = 0;
    ksi[ix] = 0;
    dc = Inf;
    state = 2;
    terminate= false;

    # show(Qc)
    # println(" ")
    #    @printf("state=%d, ix=%d, xp[ix] = %6.3f, xh[ix] = %6.3f, "*
    #       "delta[ix] = %6.3f\n",state,ix,xp[ix],xh[ix],delta[ix])
    while ~terminate

        if state==2 # Step 2: DFE
            xp[ix] = round((yp[ix]-ksi[ix])/R[ix,ix]);
            delta[ix] = sign(yp[ix]-ksi[ix] - R[ix,ix]*xp[ix]);
            state = 3;
        elseif state== 3 # Step 3: Main loop
            ytmp = yp[ix]-ksi[ix] - R[ix,ix]*xp[ix];
            # @printf("ytmp'*ytmp = %6.3f, TT[ix] = %6.3f, "*
            #   "TT[ix] + ytmp'*ytmp = %6.3f, dc=%6.3f\n",
            #   ytmp'*ytmp,TT[ix],TT[ix] + ytmp'*ytmp,dc)
            if TT[ix] + ytmp'*ytmp >= dc
                # We're outside the sphere
                state = 4;
                # elseif !(0<=real(xp[ix])<Qc[ix]) || !(0<=imag(xp[ix])<Qc[ix])
            elseif !(0 <= xp[ix] < Qc[ix])
                # inside sphere but outside signal set boundary
                state = 6;
            else # inside sphere and inside signal set
                if ix>1
                    # ksi[ix-1] = (conj(R[ix-1,ix:M])*xp[ix:M])[1];
                    ksi[ix-1] = (R[ix-1,ix:M]*xp[ix:M])[1];
                    TT[ix-1] = TT[ix] + ytmp'*ytmp;
                    ix = ix-1;
                    state = 2;
                else
                    state = 5;
                end
            end
        elseif state==4 # Step 4: Check termination criteria
            if ix==M
                terminate=true;
            else
                ix = ix+1;
                state = 6;
            end
        elseif state==5 # Step 5: A valid point is found
            ytmp = yp[1] -ksi[1] - R[1,1]*xp[1];
            dc = TT[1] + ytmp'*ytmp;
            xh[:] = xp[:];
            ix = ix+1;
            state = 6;
        elseif state==6 # Step 6: Schnorr-Euchner
            xp[ix] = xp[ix] + delta[ix];
            delta[ix] = -delta[ix] -sign(delta[ix]);
            state = 3;
        else
            error("Unknown state");
        end
        #@printf("state=%d, ix=%d, xp = [%1.0f %1.0f %1.0f %1.0f], "*
        #    xh = [%1.0f %1.0f %1.0f %1.0f], δ[ix] = %2.0f\n",
        #    state,ix,xp[1],xp[2],xp[3],xp[4],xh[1],xh[2],xh[3],xh[4],delta[ix])
        #    @printf("state=%d, ix=%d, xp = [%1.0f %1.0f], xh = [%1.0f "*
        #      "%1.0f],  δ[ix] = %2.0f\n",state,ix,xp[1],xp[2],xh[1],xh[2],delta[ix])
    end
    #println("xh=$(xh)")
    return xh
end


text:
# check overall performance using randomly chosen points
using Polylogarithms
using DataFrames, CSV
using PyPlot
using PyCall
using Printf
using StatsBase
fs = (4,4)
fs2 = (5,4)
PyCall.PyDict(matplotlib["rcParams"])["legend.markerscale"] = 3.0
# https://matplotlib.org/3.1.0/api/legend_api.html

L = Symbol("Li_s(z)")

# input data from Mathematica and reparse into complex numbers
C = 2
filename = @sprintf("../data/polylog_test_data_rand_%d.csv", C)
data1 = CSV.read(filename; delim=",", type=String)

# has trouble reading in numbers like "2." so read all into strings, and parse
data1[!,:s] = parse.(Complex{Float64}, data1[!,:s] )
data1[!,:z] = parse.(Complex{Float64}, data1[!,:z] )
data1[!, L] = parse.(Complex{Float64}, data1[!, L] )

m = size(data1,1)
Li = data1[!, L]
s  = data1[!,:s]
z  = data1[!,:z]

# plot the locations of points
markerSize = 1

figure("z", figsize=fs)
title("distribution of z")
plot(real.(z), imag.(z), "."; markersize=markerSize)
ylabel("Im(z)")
xlabel("Re(z)")
axis("equal")

figure("s", figsize=fs)
title("distribution of s")
plot(real.(s), imag.(s), "."; markersize=markerSize)
ylabel("Im(s)")
xlabel("Re(s)")
axis("equal")

S1 = zeros(Complex{Float64}, m)
S2 = zeros(Complex{Float64}, m)
n1 = zeros(Int64, m)
series1 = zeros(Int64, m)
error1 = zeros(Float64, m)
rel_error1 = zeros(Float64, m)

accuracy = 1.0e-12 # this is the default anyway, but can change it for testing
# polylog(1.5, 0.4)
polylog( complex(1.5), complex(0.4), Diagnostics() )
val, t, bytes, gctime, memallocs = @timed begin
    for i=1:m
        # print(".")
        S1[i],  n1[i], series1[i] = polylog(s[i], z[i], Diagnostics(); min_iterations=0, accuracy=accuracy )
        # S1[i] = result1[1]
        # n1[i] = result1[2]
        # series1[i] = result1[3]
    end
end
@printf("  finished %d calculations in time t=%.2f seconds, i.e., %.1f microsec each\n", m, t, (t/m)*1.0e6)

error1 = abs.( S1 - Li  )
rel_error1 = error1 ./ abs.( Li )
println("   max abs. rel. error1 = $(maximum( abs.(rel_error1) ))")

fig = figure(@sprintf("../data/polylog_bench_rand_1.csv"), figsize=(6,4))
clf()
h = hist(log10.(rel_error1), collect(-16 : 0.5 : -8); align="mid" )
ylabel("number")
xlabel("log10(relative absolute error)")
plot( log10(maximum(rel_error1))*[1,1], [0,maximum(h[1])], "--")
savefig(@sprintf("plots/polylog_bench_rand_%d.pdf", C); bbox_inches="tight")

figure("series", figsize=fs2)
ms = maximum(series1)
x = hist( series1, collect(0.5:1:ms+0.5); align="mid" )
for i=1:length(x[1])
    if i<=3
        @printf("  series %d used %d times\n", i, x[1][i])
    elseif i<8
        @printf("  reciprocal + series %d used %d times\n", i-3, x[1][i])
    elseif i< 20
        @printf("  duplication + series %d used %d times\n", i-10, x[1][i])
    else
        @printf("  duplication used more than once + series %d used %d times\n", i-10, x[1][i])
    end
end

# look at bad points
k = findall( rel_error1 .> 1.0e-12 )
count_bad = length(k)
bad_guys = DataFrame(s=s[k], z=z[k], Li=Li[k], err=rel_error1[k], series=series1[k],  n_terms=n1[k])

figure("z")
plot(real.(z[k]), imag.(z[k]), "r."; markersize=markerSize+1)

figure("s")
plot(real.(s[k]), imag.(s[k]), "r."; markersize=markerSize+1)

s_min = floor( minimum(min.(real.(s), imag.(s)) ))
s_max = ceil(  maximum(max.(real.(s), imag.(s)) ))
s_bins = collect(  s_min : 0.5 : s_max )
m = length(s_bins)
s_b = zeros(m-1)
error_v_s = zeros(m-1)
for i=1:m-1
    s_b[i] = (s_bins[i] + s_bins[i+1])/2
    k_i = findall( s_bins[i] .<= imag.(s) .< s_bins[i+1] )
    error_v_s[i] = log10.( mean( rel_error1[k_i] ) )
end
figure("error vs s", figsize=fs2)
plot( s_b, error_v_s)
ylabel("number")
xlabel("series (values >4 are a sum of two components)")

# look at really bad points
k = findall( (rel_error1 .> 1.0e-8) .| isinf.(S1) )
count_really_bad = length(k)
really_bad_guys = DataFrame(s=s[k], z=z[k], Li=Li[k], err=rel_error1[k], series=series1[k],  n_terms=n1[k])

figure("z")
plot(real.(z[k]), imag.(z[k]), "r."; markersize=markerSize+3)
savefig( @sprintf("plots/polylog_bench_rand_errors_z_%d.pdf", C); bbox_inches="tight")

figure("s")
plot(real.(s[k]), imag.(s[k]), "r."; markersize=markerSize+3)
savefig( @sprintf("plots/polylog_bench_rand_errors_s_%d.pdf", C); bbox_inches="tight")

# show which series is used where
figure("z2", figsize=fs)
for i=1:ms
    k = findall( series1 .== i )
    if length(k) > 0
        if i<10
            plot(real.(z[k]), imag.(z[k]), "."; markersize=markerSize, label=@sprintf("Series %d",i))
        elseif i<50
            n_series = (i - floor(i,digits=-1)) / 2 # not right when we get to 50
            plot(real.(z[k]), imag.(z[k]), "."; markersize=markerSize, label=@sprintf("Series 2 used %d times", n_series))
        elseif i>=50
            i -= 50
            n_series = 5
            n_series += (i - floor(i,digits=-1)) / 2 # not right when we get to 50
            plot(real.(z[k]), imag.(z[k]), "."; markersize=markerSize, label=@sprintf("Series 2 used %d times", n_series))
        end
    end
end
legend()
ylabel("Im(z)")
xlabel("Re(z)")
axis("equal")
savefig( @sprintf("plots/polylog_bench_rand_domains_%d.pdf", C); bbox_inches="tight")


text:
# conv-histogram-picture-generator.jl

# This code will generate a bunch of easy, but random, convolution, and then two pictures apiece: One showing the x[n] and h[n] as two separate histograms; the next, showing x[n] and x[n]*h[n]=h[n]*x[n] as a new histogram.

using DSP, Plots

#---------------------------------------------------------#

# First, let's test to make sure this works out correctly.

print("\n\n\n")
println("Hello, Julia! Long time no see! :)")

#---------------------------------------------------------#

println("Okay, let's start easy with an example of conv().")

x0 = conv([1 1], [1 2])

println("x0 = $(x0)")

println("Works like a charm. :)")

#---------------------------------------------------------#

# Okay, so now we need to actually generate some histograms.
# Time to dive back into the Plots documentation.

LENGTH = 5
x = rand((-1,0,1), LENGTH)
h = rand((-1,0,1), LENGTH)
y = conv(x, h)

plot(1:LENGTH, [x,h], seriestype=:bar)

# Okay, that's starting to look like what we want.

#---------------------------------------------------------#

# Let's first figure out how to plot x as a histogram.

LENGTH = 12
x = rand((-1,0,1), LENGTH)
plot(1:LENGTH, x, seriestype=:bar)

# Okay, that seems to work. Now let's see if we can get x and h.

LENGTH = 5
x = rand((-1,0,1), LENGTH)
h = rand((-1,0,1), LENGTH)
y = conv(x, h)

p1 = plot(1:LENGTH, x,
          seriestype=:bar,
          title="x[n]",
          linealpha=0,
          color=:blue,
          alpha=0.7,
          legend=:none)
plot!(showaxis=true)
plot!(ylims=(-1,1))

p2 = plot(1:LENGTH, h,
          seriestype=:bar,
          title="h[n]",
          linealpha=0,
          color=:grey,
          alpha=0.7,
          legend=:none)
plot!(showaxis=true)
plot!(ylims=(-1,1))

p3 = plot(1:LENGTH*2-1,
          y,
          seriestype=:bar,
          linealpha=0,
          alpha=0.7,
          color=:red)
plot!(1:LENGTH*2-1,
      [x;zeros(LENGTH-1)],
      seriestype=:bar,
      linealpha=0,
      alpha=0.2,
      color=:blue,
      legend=:none)
plot!(1:LENGTH*2-1,
      [h;zeros(LENGTH-1)],
       seriestype=:bar,
       linealpha=0,
       alpha=0.1,
       color=:grey)
plot!(title="y[n] = x[n] * h[n] = h[n] * x[n]")
plot!(ylims=(-5,5))
plot!(showaxis=true)
plot!(minorgrid=true)

l = @layout [a b; c]
plot(p1, p2, p3, layout = l)

# Sweet! Looks like we're really playing with fire now.

#---------------------------------------------------------#

# Okay, so now, I have to think here. How do I want to actually
# generate these? What would be useful?

# I think that 4 or 5 data points is about the sweet spot for this.
# That's few enough that I can visualize the two going together well.
# Let's generate a bunch of cards with that.

question_layout = @layout[a;b]
plot(p1, p2, layout = question_layout)

#---------------------------------------------------------#

answer_layout = @layout[a b;c]
plot(p3)
plot!(legend=true)
plot!(label=["y[n]" "x[n]" "h[n]"])


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

include("bat_default.jl")
include("initval_algorithm.jl")
include("transform_algorithm.jl")
include("sampling_algorithm.jl")
include("autocor_len.jl")
include("eff_sample_size.jl")
include("mode_estimator.jl")
include("median_estimator.jl")
include("integration_algorithm.jl")


text:
# Arm with Variation V_t as mentioned in
#       Besbes, O., Gur, Y., & Zeevi, A. (2014). Optimal Exploration-Exploitation in a Multi-Armed-Bandit Problem with Non-stationary Rewards, 1–20.

type Variational <: BanditArmBase
    step::Int64         # Time step
    variation::Float64  # Variation of the arms
    period::Float64     # Period of the sinusoidal wave
    isRestless::Bool    # Is this arm restless - True by default
    offset::Float64     # Initial pahse, random by default

    function Variational( _variation::Float64, _period::Float64, )
        new(
            0,
            _variation,
            _period,
            true,
            2*π*rand()
        )
    end

    function Variational( _variation::Float64, _period::Float64, _isRestless::Bool )
        new(
            0,
            _variation,
            _period,
            _isRestless,
            2*π*rand()
        )
    end

    function Variational( _variation::Float64, _period::Float64, _offset::Real )
        new(
            0,
            _variation,
            _period,
            true,
            _offset
        )
    end

    function Variational( _variation::Float64, _period::Int64, _isRestless::Bool, _offset::Real )
        new(
            0,
            _variation,
            _period,
            _isRestless,
            _offset
        )
    end
end

function pull!( arm::Variational )
    arm.step = arm.step + 1
    return 1/2 + 1/2 * sin( arm.variation * π * (arm.step-1)/arm.period + arm.offset )  # To limit rewards between 0 and 1
end

function tick!( arm::Variational )
    if arm.isRestless
        arm.step = arm.step + 1
    end
    nothing
end

function reset!( arm::Variational )
    arm.step = 0
    nothing
end


text:
using CayleyDickson
using Base.Test: @test, @test_throws

@test begin
    a = CayleyDickson.randomBigFloat()
    isreal(SplitQuaternion(a))
end

@test begin
    a = 1
    !isreal(SplitQuaternion(a, a, a, a))
end

@test begin
    z = random(SplitQuaternion{Int})
    z == +(z)
end

@test begin
    a = 1
    real(SplitQuaternion(a)) == a
end

@test begin
    a = 1
    b = 2
    c = 3
    l = unreal(SplitQuaternion(0, a, b, c))
    r = [a, b, c]
    l == r
end

@test begin
    io = IOBuffer()
    show(io, SplitQuaternion(1, 2, 3, 4))
    l = String(take!(io))
    r = "[1: 1, i: 2, s: 3, is: 4]"
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    b = CayleyDickson.randomBigFloat()
    l = SplitQuaternion(Binion(a, b))
    r = SplitQuaternion(a, b)
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    b = CayleyDickson.randomBigFloat()
    c = CayleyDickson.randomBigInt()
    l = SplitQuaternion(a, b, c)
    r = SplitQuaternion(a, b, c, 0)
    l == r
end

@test begin
    x = random(SplitQuaternion{BigFloat})
    y = random(SplitQuaternion{BigFloat})
    x != y
end

@test begin
    a = CayleyDickson.randomBigInt()
    z = random(SplitQuaternion{BigInt})
    l = a + z
    r = z + a
    l == r
end

@test begin
    a = random(Binion{BigInt})
    z = random(SplitQuaternion{BigInt})
    l = a + z
    r = z + a
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    z = random(SplitQuaternion{BigInt})
    l = a - z
    r = -(z - a)
    l == r
end

@test begin
    a = random(Binion{BigInt})
    z = random(SplitQuaternion{BigInt})
    l = a - z
    r = -(z - a)
    l == r
end

@test begin
    x = random(SplitQuaternion{BigInt})
    y = random(SplitQuaternion{BigInt})
    l = x + y
    r = y + x
    l == r
end

@test begin
    x = random(SplitQuaternion{BigInt})
    y = random(SplitQuaternion{BigInt})
    l = x * y
    r = y * x
    l != r
end

@test begin
    x = random(SplitQuaternion{BigInt})
    y = random(SplitQuaternion{BigInt})
    z = random(SplitQuaternion{BigInt})
    l = associator(x, y, z)
    r = zero(SplitQuaternion{BigInt})
    l == r
end

@test begin
    x = random(SplitQuaternion{BigInt})
    y = random(SplitQuaternion{BigInt})
    l = x - y
    r = -(y - x)
    l == r
end

@test begin
    x = random(SplitQuaternion{BigInt})
    y = random(SplitQuaternion{BigInt})
    z = random(SplitQuaternion{BigInt})
    l = (x + y) + z
    r = x + (y + z)
    l == r
end

@test begin
    x = random(SplitQuaternion{BigInt})
    y = random(SplitQuaternion{BigInt})
    z = random(SplitQuaternion{BigInt})
    l = (x * y) * z
    r = x * (y * z)
    l == r
end

@test begin
    x = random(SplitQuaternion{BigInt})
    conj(conj(x)) == x
end

@test begin
    x = random(SplitQuaternion{Rational{BigInt}})
    inv(inv(x)) == x
end

@test begin
    x = random(SplitQuaternion{BigInt})
    y = random(SplitQuaternion{BigInt})
    l = conj(x * y)
    r = conj(y) * conj(x)
    l == r
end

@test begin
    x = random(SplitQuaternion{BigInt})
    y = random(SplitQuaternion{BigInt})
    l = conj(x + y)
    r = conj(x) + conj(y)
    l == r
end

@test begin
    x = random(SplitQuaternion{BigInt})
    y = random(SplitQuaternion{BigInt})
    l = conj(x - y)
    r = conj(x) - conj(y)
    l == r
end

@test begin
    x = random(SplitQuaternion{Rational{BigInt}})
    y = random(SplitQuaternion{Rational{BigInt}})
    l = inv(x * y)
    r = inv(y) * inv(x)
    l == r
end

@test begin
    x = random(SplitQuaternion{Rational{BigInt}})
    y = random(SplitQuaternion{Rational{BigInt}})
    l = inv(x / y)
    r = y / x
    l == r
end

@test begin
    x = random(SplitQuaternion{Rational{BigInt}})
    y = random(SplitQuaternion{Rational{BigInt}})
    l = inv(x \ y)
    r = y \ x
    l == r
end

@test begin
    x = random(SplitQuaternion{BigInt})
    y = random(SplitQuaternion{BigInt})
    z = random(SplitQuaternion{BigInt})
    l = x * (y + z)
    r = (x * y) + (x * z)
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    x = random(SplitQuaternion{BigInt})
    y = random(SplitQuaternion{BigInt})
    l = a * (x + y)
    r = (a * x) + (a * y)
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    x = random(SplitQuaternion{BigInt})
    l = a * x
    r = x * a
    l == r
end

@test begin
    a = random(Binion{BigInt})
    x = random(SplitQuaternion{BigInt})
    l = a * x
    r = x * a
    l != r
end

@test begin
    x = random(SplitQuaternion{BigInt})
    y = random(SplitQuaternion{BigInt})
    z = random(SplitQuaternion{BigInt})
    l = (x + y) * z
    r = (x * z) + (y * z)
    l == r
end

@test begin
    x = random(SplitQuaternion{BigInt})
    y = random(SplitQuaternion{BigInt})
    z = random(SplitQuaternion{BigInt})
    l = x * (y - z)
    r = (x * y) - (x * z)
    l == r
end

@test begin
    x = random(SplitQuaternion{BigInt})
    y = random(SplitQuaternion{BigInt})
    z = random(SplitQuaternion{BigInt})
    l = (x - y) * z
    r = (x * z) - (y * z)
    l == r
end

@test begin
    x = random(SplitQuaternion{BigInt})
    quadrance(x) > big(0) || quadrance(x) < big(0) || quadrance(x) == big(0)
end

@test begin
    x = random(SplitQuaternion{BigInt})
    y = random(SplitQuaternion{BigInt})
    l = quadrance(x * y)
    r = quadrance(x) * quadrance(y)
    l == r
end

@test_throws ErrorException begin
    inv(SplitQuaternion(0, 0.0))
end

@test_throws ErrorException begin
    inv(SplitQuaternion(1, 1, 1, 1))
end

@test_throws ErrorException begin
    inv(SplitQuaternion(1, 2, 2, 1))
end

@test_throws ErrorException begin
    random(SplitQuaternion{Int}) / 0
end

@test_throws ErrorException begin
    0 \ random(SplitQuaternion{Int})
end

@test begin
    a = CayleyDickson.randomBigRational()
    z = random(SplitQuaternion{Rational{BigInt}})
    l = z / a
    r = a \ z
    l == r
end

@test begin
    a = CayleyDickson.randomBigRational()
    z = random(SplitQuaternion{Rational{BigInt}})
    l = a / z
    r = z \ a
    l == r
end

@test begin
    z = random(SplitQuaternion{BigInt})
    z == z + zero(SplitQuaternion{BigInt})
end

@test begin
    z = random(SplitQuaternion{BigInt})
    z == z * one(SplitQuaternion{BigInt})
end

@test begin
    z = random(SplitQuaternion{BigInt})
    z == z + zero(z)
end

@test begin
    z = random(SplitQuaternion{BigInt})
    z == z * one(z)
end


text:
@inline signbit(a::DoubleFloat{T}) where {T} = signbit(HI(a))
@inline sign(a::DoubleFloat{T}) where {T} = sign(HI(a))

@inline function (-)(a::DoubleFloat{T}) where {T}
    if iszero(LO(a))
        DoubleFloat(-HI(a), LO(a))
    else
        DoubleFloat(-HI(a), -LO(a))
    end
end

@inline function abs(a::DoubleFloat{T}) where {T}
    if HI(a) >= 0.0
        a
    else # HI(a) < 0.0
        -a
    end
end

@inline function flipsign(x::DoubleFloat{T}, y::T) where {T<:AbstractFloat}
    signbit(y) ? -x : x
end
@inline function flipsign(x::DoubleFloat{T}, y::DoubleFloat{T}) where {T<:AbstractFloat}
    signbit(y) ? -x : x
end

@inline function copysign(x::DoubleFloat{T}, y::T) where {T<:AbstractFloat}
    signbit(y) ? -abs(x) : abs(x)
end
@inline function copysign(x::DoubleFloat{T}, y::DoubleFloat{T}) where {T<:AbstractFloat}
    signbit(y) ? -abs(x) : abs(x)
end

flipsign(x::DoubleFloat{T}, y::U) where {T<:AbstractFloat, U<:Unsigned} = +x
copysign(x::DoubleFloat{T}, y::U) where {T<:AbstractFloat, U<:Unsigned} = +x

flipsign(x::DoubleFloat{T}, y::S) where {T<:AbstractFloat, S<:Signed} =
    signbit(y) ? -x : x
copysign(x::DoubleFloat{T}, y::S) where {T<:AbstractFloat, S<:Signed} =
    signbit(y) ? -abs(x) : abs(x)

function frexp(x::DoubleFloat{T}) where {T<:AbstractFloat}
    frhi, exphi = frexp(HI(x))
    frlo, explo = frexp(LO(x))
    return (frhi, exphi), (frlo, explo)
end

function ldexp(x::DoubleFloat{T}, exponent::I) where {T, I<:Integer}
    return DoubleFloat(ldexp(HI(x), exponent), ldexp(LO(x), exponent))
end

function ldexp(dhi::Tuple{T,I}, dlo::Tuple{T,I}) where {T<:AbstractFloat, I<:Integer}
    return DoubleFloat(ldexp(dhi[1], dhi[2]), ldexp(dlo[1], dlo[2]))
end

function ldexp(dhilo::Tuple{Tuple{T,I}, Tuple{T,I}}) where {T<:AbstractFloat, I<:Integer}
    return ldexp(dhilo[1], dhilo[2])
end

function exponent(x::DoubleFloat{T}) where {T<:AbstractFloat}
    ehi = Base.exponent(HI(x))
    elo = Base.exponent(LO(x))
    return ehi, elo
end

function significand(x::DoubleFloat{T}) where {T<:AbstractFloat}
    shi = Base.significand(HI(x))
    slo = Base.significand(LO(x))
    return shi, slo
end

function signs(x::DoubleFloat{T}) where {T<:AbstractFloat}
    shi = sign(HI(x))
    slo = sign(LO(x))
    return shi, slo
end

ulp(x::T) where {T<:Base.IEEEFloat} = significand(x) !== -one(T) ? eps(x) : eps(x)/2

posulp(x::T) where {T<:Base.IEEEFloat} = significand(x) !== -one(T) ? eps(x) : eps(x)/2
negulp(x::T) where {T<:Base.IEEEFloat} = significand(x) !== one(T) ? -eps(x) : -eps(x)/2

ulp(::Type{T}) where {T<:Base.IEEEFloat} = ulp(one(T))
posulp(::Type{T}) where {T<:Base.IEEEFloat} = ulp(one(T))
negulp(::Type{T}) where {T<:Base.IEEEFloat} = ulp(one(T))

function eps(x::DoubleFloat{T}) where {T<:AbstractFloat}
    return LO(x) !== zero(T) ? eps(LO(x)) : eps(posulp(HI(x)))
end

function ulp(x::DoubleFloat{T}) where {T<:AbstractFloat}
    return LO(x) !== zero(T) ? posulp(LO(x)) : posulp(posulp(HI(x)))
end

eps(::Type{D}) where {T<:AbstractFloat, D<:DoubleFloat{T}} = D(eps(posulp(one(T))))
ulp(::Type{D}) where {T<:AbstractFloat, D<:DoubleFloat{T}} = D(posulp(poslulp(one(T))))

function nextfloat(x::DoubleFloat{T}) where {T<:AbstractFloat}
    !isfinite(x) && return(x)
    if !iszero(LO(x))
        DoubleFloat{T}(HI(x)) + nextfloat(LO(x))
    else
        DoubleFloat{T}(HI(x), ulp(ulp(HI(x))))
    end
end

function prevfloat(x::DoubleFloat{T}) where {T<:AbstractFloat}
    !isfinite(x) && return(x)
    if !iszero(LO(x))
        DoubleFloat{T}(HI(x)) + prevfloat(LO(x))
    else
        DoubleFloat{T}(HI(x), -ulp(ulp(HI(x))))
    end
end

function nextfloat(x::DoubleFloat{T}, n::Int) where {T<:AbstractFloat}
    !isfinite(x) && return(x)
    if !iszero(LO(x))
        x + n*ulp(LO(x))
    else
        DoubleFloat{T}(HI(x), n*ulp(ulp(HI(x))))
    end
end

function prevfloat(x::DoubleFloat{T}, n::Int) where {T<:AbstractFloat}
    !isfinite(x) && return(x)
    if !iszero(LO(x))
        x - n*ulp(LO(x))
    else
        DoubleFloat{T}(HI(x), -n*ulp(ulp(HI(x))))
    end
end

@inline function intpart(x::Tuple{T,T}) where {T<:AbstractFloat}
    hi, lo = x
    ihi = modf(hi)[2]
    ilo = modf(lo)[2]
    return add_2(ihi, ilo)
end

@inline function fracpart(x::Tuple{T,T}) where {T<:AbstractFloat}
    hi, lo = x
    fhi = modf(hi)[1]
    flo = modf(lo)[1]
    return add_2(fhi, flo)
end

function intpart(x::DoubleFloat{T}) where {T<:AbstractFloat}
    ihi, ilo = intpart(HILO(x))
    return DoubleFloat(ihi, ilo)
end

function fracpart(x::DoubleFloat{T}) where {T<:AbstractFloat}
    fhi, flo = fracpart(HILO(x))
    return DoubleFloat(fhi, flo)
end

function modf(x::DoubleFloat{T}) where {T<:AbstractFloat}
    hi, lo = HILO(x)
    fhi, ihi = modf(hi)
    flo, ilo = modf(lo)
    ihi, ilo = add_2(ihi, ilo)
    fhi, flo = add_2(fhi, flo)
    i = DoubleFloat(ihi, ilo)
    f = DoubleFloat(fhi, flo)
    return f, i
end

function fmod(fpart::DoubleFloat{T}, ipart::DoubleFloat{T}) where {T<:AbstractFloat}
   return ipart + fpart
end

function fmod(parts::Tuple{DoubleFloat{T}, DoubleFloat{T}}) where {T<:AbstractFloat}
   return parts[1] + parts[2]
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule UDUNITS_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("UDUNITS")
JLLWrappers.@generate_main_file("UDUNITS", UUID("2619efb9-15dd-5c8d-95a3-19a84f6456dc"))
end  # module UDUNITS_jll


text:
using RecursiveArrayTools: recursivecopy

struct SimpleATsit5 <: DiffEqBase.AbstractODEAlgorithm end
export SimpleATsit5

# PI-adaptive stepping parameters:
const beta1 = 7/50
const beta2 = 2/25
const qmax = 10.0
const qmin = 1/5
const gamma = 9/10
const qoldinit = 1e-4

mutable struct SimpleATsit5Integrator{IIP, S, T, P, F, N} <: DiffEqBase.AbstractODEIntegrator{SimpleATsit5, IIP, S, T}
    f::F                  # eom
    uprev::S              # previous state
    u::S                  # current state
    tmp::S                # dummy, same as state
    tprev::T              # previous time
    t::T                  # current time
    t0::T                 # initial time, only for reinit
    tf::T
    dt::T                 # step size
    dtnew::T
    tdir::T
    p::P                  # parameter container
    u_modified::Bool
    ks::Vector{S}         # interpolants of the algorithm
    cs::SVector{6, T}     # ci factors cache: time coefficients
    as::SVector{21, T}    # aij factors cache: solution coefficients
    btildes::SVector{7,T}
    rs::SVector{22, T}    # rij factors cache: interpolation coefficients
    qold::Float64
    abstol::Float64
    reltol::Float64
    internalnorm::N       # function that computes the error EEst based on state
end
const SAT5I = SimpleATsit5Integrator

DiffEqBase.isinplace(::SAT5I{IIP}) where {IIP} = IIP
DiffEqBase.u_modified!(i::SAT5I, bool) = (i.u_modified = bool)

#######################################################################################
# Initialization
#######################################################################################
function DiffEqBase.__init(prob::ODEProblem,alg::SimpleATsit5;
                         dt = 0.1,
                         abstol = 1e-6, reltol = 1e-3,
                         internalnorm = DiffEqBase.ODE_DEFAULT_NORM, kwargs...)
  simpleatsit5_init(prob.f,DiffEqBase.isinplace(prob),prob.u0,
                   prob.tspan[1], prob.tspan[2], dt, prob.p, abstol, reltol,
                   internalnorm)
end

function DiffEqBase.__solve(prob::ODEProblem,alg::SimpleATsit5;
                          dt = 0.1, saveat = nothing, save_everystep = true,
                          abstol = 1e-6, reltol = 1e-3,
                          internalnorm = DiffEqBase.ODE_DEFAULT_NORM)
  u0 = prob.u0
  tspan = prob.tspan
  ts = Vector{eltype(dt)}(undef,1)

  if saveat === nothing
    ts = Vector{eltype(dt)}(undef,1)
    ts[1] = prob.tspan[1]
    us = Vector{typeof(u0)}(undef,0)
    push!(us,recursivecopy(u0))
  else
    ts = saveat
    cur_t = 1
    us = MVector{length(ts),typeof(u0)}(undef)
    if prob.tspan[1] == ts[1]
      cur_t += 1
      us[1] = u0
    end
  end

  integ = simpleatsit5_init(prob.f,DiffEqBase.isinplace(prob),prob.u0,
                            tspan[1], tspan[2], dt, prob.p, abstol, reltol, internalnorm)
  # FSAL
  while integ.t < tspan[2]
    step!(integ)
    if saveat === nothing && save_everystep
        push!(us,recursivecopy(integ.u))
        push!(ts,integ.t)
    else saveat !== nothing
      while cur_t <= length(ts) && ts[cur_t] <= integ.t
        if !isinplace(prob)
            us[cur_t] = integ(ts[cur_t])
        else
            integ(us[cur_t],ts[cur_t])
        end
        cur_t += 1
      end
    end
  end

  if saveat === nothing && !save_everystep
      push!(us,recursivecopy(u))
      push!(ts,t)
  end

  sol = DiffEqBase.build_solution(prob,alg,ts,us,
                                  calculate_error = false)
  DiffEqBase.has_analytic(prob.f) && DiffEqBase.calculate_solution_errors!(sol;timeseries_errors=true,dense_errors=false)
  sol
end

@inline function simpleatsit5_init(f::F,
                         IIP::Bool, u0::S, t0::T, tf::T, dt::T, p::P,
                         abstol, reltol,
                         internalnorm::N) where {F, P, S, T, N}

    cs, as, btildes, rs = _build_atsit5_caches(T)
    ks = _initialize_ks(u0)

    !IIP && @assert S <: SArray

    integ = SAT5I{IIP, S, T, P, F, N}(
        f, recursivecopy(u0), recursivecopy(u0), recursivecopy(u0), t0, t0, t0, tf, dt, dt,
        sign(tf-t0), p, true, ks, cs, as, btildes, rs,
        qoldinit,abstol,reltol, internalnorm
    )
end

@inline _initialize_ks(u0::AbstractArray{T}) where {T<:Number} = [zero(u0) for i in 1:7]
@inline function _initialize_ks(u0::Vector{<:AbstractVector})
    return [[zero(u0[j]) for j in 1:length(u0)] for i in 1:7]
end

#######################################################################################
# IIP version for vectors and matrices
#######################################################################################
@inline @muladd function DiffEqBase.step!(integ::SAT5I{true, S, T}) where {S, T}

    L = length(integ.u)

    c1, c2, c3, c4, c5, c6 = integ.cs;
    dt = integ.dtnew; t = integ.t; p = integ.p; tf = integ.tf
    a21, a31, a32, a41, a42, a43, a51, a52, a53, a54,
    a61, a62, a63, a64, a65, a71, a72, a73, a74, a75, a76 = integ.as
    btilde1, btilde2, btilde3, btilde4, btilde5, btilde6, btilde7 = integ.btildes

    k1, k2, k3, k4, k5, k6, k7 = integ.ks
    tmp = integ.tmp; f! = integ.f

    integ.uprev .= integ.u; uprev = integ.uprev; u = integ.u

    qold = integ.qold
    abstol = integ.abstol
    reltol = integ.reltol

    if integ.u_modified
      f!(k1, uprev, p, t)
      integ.u_modified=false
    else
      k1 .= k7
    end

    EEst = Inf

    @inbounds while EEst>1
      dt < 1e-14 && error("dt<dtmin")

      for i in 1:L
          tmp[i] = uprev[i]+dt*a21*k1[i]
      end
      f!(k2, tmp, p, t+c1*dt)
      for i in 1:L
          tmp[i] = uprev[i]+dt*(a31*k1[i]+a32*k2[i])
      end
      f!(k3, tmp, p, t+c2*dt)
      for i in 1:L
          tmp[i] = uprev[i]+dt*(a41*k1[i]+a42*k2[i]+a43*k3[i])
      end
      f!(k4, tmp, p, t+c3*dt)
      for i in 1:L
          tmp[i] = uprev[i]+dt*(a51*k1[i]+a52*k2[i]+a53*k3[i]+a54*k4[i])
      end
      f!(k5, tmp, p, t+c4*dt)
      for i in 1:L
          tmp[i] = uprev[i]+dt*(a61*k1[i]+a62*k2[i]+a63*k3[i]+a64*k4[i]+a65*k5[i])
      end
      f!(k6, tmp, p, t+dt)
      for i in 1:L
          u[i] = uprev[i]+dt*(a71*k1[i]+a72*k2[i]+a73*k3[i]+a74*k4[i]+a75*k5[i]+a76*k6[i])
      end
      f!(k7, u, p, t+dt)

      for i in 1:L
        tmp[i] = dt*(btilde1*k1[i]+btilde2*k2[i]+btilde3*k3[i]+btilde4*k4[i]+
                     btilde5*k5[i]+btilde6*k6[i]+btilde7*k7[i])
        tmp[i] = tmp[i]/(abstol+max(abs(uprev[i]),abs(u[i]))*reltol)
      end

      EEst = integ.internalnorm(tmp, integ.t)

      if iszero(EEst)
        q = inv(qmax)
      else
        @fastmath q11 = EEst^beta1
        @fastmath q = q11/(qold^beta2)
      end

      if EEst > 1
        dt = dt/min(inv(qmin),q11/gamma)
      else # EEst <= 1
        @fastmath q = max(inv(qmax),min(inv(qmin),q/gamma))
        qold = max(EEst,qoldinit)
        dtnew = dt/q #dtnew
        dtnew = min(abs(dtnew),abs(tf-t-dt))

        integ.dt = dt
        integ.dtnew = dtnew
        integ.qold = qold
        integ.tprev = t

        if (tf - t - dt) < 1e-14
          integ.t = tf
        else
          integ.t += dt
        end
      end
    end

    return  nothing
end

#######################################################################################
# OOP version for vectors and matrices
#######################################################################################
@inline @muladd function DiffEqBase.step!(integ::SAT5I{false, S, T}) where {S, T}

    c1, c2, c3, c4, c5, c6 = integ.cs;
    dt = integ.dtnew; t = integ.t; p = integ.p; tf = integ.tf
    a21, a31, a32, a41, a42, a43, a51, a52, a53, a54,
    a61, a62, a63, a64, a65, a71, a72, a73, a74, a75, a76 = integ.as
    btilde1, btilde2, btilde3, btilde4, btilde5, btilde6, btilde7 = integ.btildes

    tmp = integ.tmp; f = integ.f
    integ.uprev = integ.u; uprev = integ.u

    qold = integ.qold
    abstol = integ.abstol
    reltol = integ.reltol

    if integ.u_modified
      k1 = f(uprev, p, t)
      integ.u_modified=false
    else
      @inbounds k1 = integ.ks[7];
    end

    EEst = Inf

    while EEst > 1
      dt < 1e-14 && error("dt<dtmin")

      tmp = uprev+dt*a21*k1
      k2 = f(tmp, p, t+c1*dt)
      tmp = uprev+dt*(a31*k1+a32*k2)
      k3 = f(tmp, p, t+c2*dt)
      tmp = uprev+dt*(a41*k1+a42*k2+a43*k3)
      k4 = f(tmp, p, t+c3*dt)
      tmp = uprev+dt*(a51*k1+a52*k2+a53*k3+a54*k4)
      k5 = f(tmp, p, t+c4*dt)
      tmp = uprev+dt*(a61*k1+a62*k2+a63*k3+a64*k4+a65*k5)
      k6 = f(tmp, p, t+dt)
      u = uprev+dt*(a71*k1+a72*k2+a73*k3+a74*k4+a75*k5+a76*k6)
      k7 = f(u, p, t+dt)

      tmp = dt*(btilde1*k1+btilde2*k2+btilde3*k3+btilde4*k4+
                   btilde5*k5+btilde6*k6+btilde7*k7)
      tmp = tmp./(abstol.+max.(abs.(uprev),abs.(u))*reltol)
      EEst = integ.internalnorm(tmp, integ.t)

      if iszero(EEst)
        q = inv(qmax)
      else
        @fastmath q11 = EEst^beta1
        @fastmath q = q11/(qold^beta2)
      end

      if EEst > 1
        dt = dt/min(inv(qmin),q11/gamma)
      else # EEst <= 1
        @fastmath q = max(inv(qmax),min(inv(qmin),q/gamma))
        qold = max(EEst,qoldinit)
        dtnew = dt/q #dtnew
        dtnew = min(abs(dtnew),abs(tf-t-dt))

        @inbounds begin # Necessary for interpolation
            integ.ks[1] = k1; integ.ks[2] = k2; integ.ks[3] = k3
            integ.ks[4] = k4; integ.ks[5] = k5; integ.ks[6] = k6
            integ.ks[7] = k7;
        end

        integ.dt = dt
        integ.dtnew = dtnew
        integ.qold = qold
        integ.tprev = t
        integ.u = u

        if (tf - t - dt) < 1e-14
          integ.t = tf
        else
          integ.t += dt
        end
      end
    end

    return  nothing
end

#######################################################################################
# Vector of Vector (always in-place) stepping
#######################################################################################
# Vector{Vector}
@inline @muladd function DiffEqBase.step!(integ::SAT5I{true, S, T}) where {S<:Vector{<:Array}, T}

    M = length(integ.u) # number of states
    L = length(integ.u[1])

    c1, c2, c3, c4, c5, c6 = integ.cs;
    dt = integ.dtnew; t = integ.t; p = integ.p; tf = integ.tf
    a21, a31, a32, a41, a42, a43, a51, a52, a53, a54,
    a61, a62, a63, a64, a65, a71, a72, a73, a74, a75, a76 = integ.as
    btilde1, btilde2, btilde3, btilde4, btilde5, btilde6, btilde7 = integ.btildes

    k1, k2, k3, k4, k5, k6, k7 = integ.ks
    tmp = integ.tmp; f! = integ.f

    @inbounds for j in 1:M
        integ.uprev[j] .= integ.u[j]
    end
    uprev = integ.uprev; u = integ.u

    qold = integ.qold; abstol = integ.abstol; reltol = integ.reltol

    # FSAL
    if integ.u_modified
        f!(k1, uprev, p, t)
        integ.u_modified=false
    else
        @inbounds for j in 1:M; k1[j] .= k7[j]; end
    end

    EEst = Inf

    @inbounds while EEst>1
        dt < 1e-14 && error("dt<dtmin")

        for j in 1:M
            for i in 1:length(integ.u[j])
                tmp[j][i] = uprev[j][i]+dt*a21*k1[j][i]
            end
        end

        f!(k2, tmp, p, t+c1*dt)
        for j in 1:M
            for i in 1:length(integ.u[j])
                tmp[j][i] = uprev[j][i]+dt*(a31*k1[j][i]+a32*k2[j][i])
            end
        end

        f!(k3, tmp, p, t+c2*dt)
        for j in 1:M
            for i in 1:length(integ.u[j])
                tmp[j][i] = uprev[j][i]+dt*(a41*k1[j][i]+a42*k2[j][i]+a43*k3[j][i])
            end
        end

        f!(k4, tmp, p, t+c3*dt)
        for j in 1:M
            for i in 1:length(integ.u[j])
                tmp[j][i] = uprev[j][i]+dt*(a51*k1[j][i]+a52*k2[j][i]+a53*k3[j][i]+a54*k4[j][i])
            end
        end

        f!(k5, tmp, p, t+c4*dt)
        for j in 1:M
            for i in 1:length(integ.u[j])
                tmp[j][i] = uprev[j][i]+dt*(a61*k1[j][i]+a62*k2[j][i]+a63*k3[j][i]+a64*k4[j][i]+a65*k5[j][i])
            end
        end

        f!(k6, tmp, p, t+dt)
        for j in 1:M
            for i in 1:length(integ.u[j])
                u[j][i] = uprev[j][i]+dt*(a71*k1[j][i]+a72*k2[j][i]+a73*k3[j][i]+a74*k4[j][i]+a75*k5[j][i]+a76*k6[j][i])
            end
        end

        f!(k7, u, p, t+dt)
        for j in 1:M
            for i in 1:length(integ.u[j])
                tmp[j][i] = dt*(btilde1*k1[j][i]+btilde2*k2[j][i]+btilde3*k3[j][i]+btilde4*k4[j][i]+
                            btilde5*k5[j][i]+btilde6*k6[j][i]+btilde7*k7[j][i])
                tmp[j][i] = tmp[j][i]/(abstol+max(abs(uprev[j][i]),abs(u[j][i]))*reltol)
            end
        end

        EEst = integ.internalnorm(tmp, integ.t)

        if iszero(EEst)
            q = inv(qmax)
        else
            @fastmath q11 = EEst^beta1
            @fastmath q = q11/(qold^beta2)
        end

        if EEst > 1
            dt = dt/min(inv(qmin),q11/gamma)
        else # EEst <= 1
            @fastmath q = max(inv(qmax),min(inv(qmin),q/gamma))
            qold = max(EEst,qoldinit)
            dtnew = dt/q #dtnew
            dtnew = min(abs(dtnew),abs(tf-t-dt))

            integ.dt = dt
            integ.dtnew = dtnew
            integ.qold = qold
            integ.tprev = t

            if (tf - t - dt) < 1e-14
                integ.t = tf
            else
                integ.t += dt
            end
        end
    end
    return  nothing
end

# Vector{SVector}
@inline @muladd function DiffEqBase.step!(integ::SAT5I{true, S, T}) where {S<:Vector{<:SVector}, T}

    M = length(integ.u)
    L = length(integ.u[1])

    c1, c2, c3, c4, c5, c6 = integ.cs;
    dt = integ.dtnew; t = integ.t; p = integ.p; tf = integ.tf
    a21, a31, a32, a41, a42, a43, a51, a52, a53, a54,
    a61, a62, a63, a64, a65, a71, a72, a73, a74, a75, a76 = integ.as
    btilde1, btilde2, btilde3, btilde4, btilde5, btilde6, btilde7 = integ.btildes

    k1, k2, k3, k4, k5, k6, k7 = integ.ks
    tmp = integ.tmp; f! = integ.f

    integ.uprev .= integ.u; uprev = integ.uprev; u = integ.u

    qold = integ.qold
    abstol = integ.abstol
    reltol = integ.reltol

    @inbounds if integ.u_modified
        f!(k1, uprev, p, t)
        integ.u_modified=false
    else
        for j in 1:M
            k1[j] = k7[j]
        end
    end

    EEst = Inf

    @inbounds while EEst>1
        dt < 1e-14 && error("dt<dtmin")

        for j in 1:M
            tmp[j] = uprev[j]+dt*a21*k1[j]
        end

        f!(k2, tmp, p, t+c1*dt)
        for j in 1:M
            tmp[j] = uprev[j]+dt*(a31*k1[j]+a32*k2[j])
        end

        f!(k3, tmp, p, t+c2*dt)
        for j in 1:M
            tmp[j] = uprev[j]+dt*(a41*k1[j]+a42*k2[j]+a43*k3[j])
        end

        f!(k4, tmp, p, t+c3*dt)
        for j in 1:M
            tmp[j] = uprev[j]+dt*(a51*k1[j]+a52*k2[j]+a53*k3[j]+a54*k4[j])
        end

        f!(k5, tmp, p, t+c4*dt)
        for j in 1:M
            tmp[j] = uprev[j]+dt*(a61*k1[j]+a62*k2[j]+a63*k3[j]+a64*k4[j]+a65*k5[j])
        end

        f!(k6, tmp, p, t+dt)
        for j in 1:M
            u[j] = uprev[j]+dt*(a71*k1[j]+a72*k2[j]+a73*k3[j]+a74*k4[j]+a75*k5[j]+a76*k6[j])
        end

        f!(k7, u, p, t+dt)

        for j in 1:M
            tmp[j] = dt*(btilde1*k1[j]+btilde2*k2[j]+btilde3*k3[j]+btilde4*k4[j]+
                             btilde5*k5[j]+btilde6*k6[j]+btilde7*k7[j])
            tmp[j] = tmp[j]./(abstol.+max.(abs.(uprev[j]),abs.(u[j]))*reltol)
        end

        EEst = integ.internalnorm(tmp, integ.t)

        if iszero(EEst)
            q = inv(qmax)
        else
            @fastmath q11 = EEst^beta1
            @fastmath q = q11/(qold^beta2)
        end

        if EEst > 1
            dt = dt/min(inv(qmin),q11/gamma)
        else # EEst <= 1
            @fastmath q = max(inv(qmax),min(inv(qmin),q/gamma))
            qold = max(EEst,qoldinit)
            dtnew = dt/q #dtnew
            dtnew = min(abs(dtnew),abs(tf-t-dt))

            integ.dt = dt
            integ.dtnew = dtnew
            integ.qold = qold
            integ.tprev = t

            if (tf - t - dt) < 1e-14
                integ.t = tf
            else
                integ.t += dt
            end
        end
    end
    return  nothing
end

#######################################################################################
# Interpolation
#######################################################################################
# Interpolation function, both OOP and IIP
@inline @muladd function (integ::SAT5I{IIP, S, T})(t::Real) where {IIP, S<:AbstractArray{<:Number}, T}
    tnext, tprev, dt = integ.t, integ.tprev, integ.dt

    θ = (t - tprev)/dt
    b1θ, b2θ, b3θ, b4θ, b5θ, b6θ, b7θ = bθs(integ.rs, θ)

    ks = integ.ks
    if !IIP
        u = @inbounds integ.uprev + dt*(b1θ*ks[1] + b2θ*ks[2] + b3θ*ks[3] + b4θ*ks[4] +
                      b5θ*ks[5] + b6θ*ks[6] + b7θ*ks[7])
        return u
    else
        u = similar(integ.u)
        @inbounds for i in 1:length(u)
            u[i] = integ.uprev[i] + dt*(b1θ*ks[1][i] + b2θ*ks[2][i] + b3θ*ks[3][i] +
                   b4θ*ks[4][i] + b5θ*ks[5][i] + b6θ*ks[6][i] + b7θ*ks[7][i])
        end
        return u
    end
end

# Interpolation function, IIP only
@inline @muladd function (integ::SAT5I{true, S, T})(u,t::Real) where {S<:AbstractArray, T}
    tnext, tprev, dt = integ.t, integ.tprev, integ.dt

    θ = (t - tprev)/dt
    b1θ, b2θ, b3θ, b4θ, b5θ, b6θ, b7θ = bθs(integ.rs, θ)

    ks = integ.ks
    @inbounds for i in 1:length(u)
        u[i] = integ.uprev[i] + dt*(b1θ*ks[1][i] + b2θ*ks[2][i] + b3θ*ks[3][i] +
               b4θ*ks[4][i] + b5θ*ks[5][i] + b6θ*ks[6][i] + b7θ*ks[7][i])
    end
end

#######################################################################################
# Multiple steps at once
#######################################################################################
@inline function DiffEqBase.step!(integ::SimpleATsit5Integrator, dt::Real, stop_at_tdt::Bool = false)
    t = integ.t
    next_t = t+dt
    while integ.t < next_t
        if stop_at_tdt && integ.t+integ.dt >= next_t
            integ.dtnew = next_t-integ.t
        end
        step!(integ)
    end
end

#######################################################################################
# reinit!
#######################################################################################
@inline function DiffEqBase.reinit!(integrator::SimpleATsit5Integrator, u0 = integrator.u;
  t0 = integrator.t0, dt = integrator.dt)

    # Is modifying the `uprev` necessary? We do `u_modified(i, true)`
    if isinplace(integrator)
        recursivecopy!(integrator.u,u0)
        recursivecopy!(integrator.uprev,integrator.u)
    else
        integrator.u = u0
        integrator.uprev = integrator.u
    end
    u_modified!(integrator, true)

    integrator.t = t0
    integrator.tprev = t0
    integrator.dt = dt
    integrator.dtnew = dt
    integrator.qold = qoldinit
end

@inline function DiffEqBase.set_t!(integrator::SimpleATsit5Integrator, t::Real)
    reinit!(integrator; t0 = t)
end


text:
# SAM Writer
# ==========

"""
    SAMWriter(output::IO, header::SAMHeader=SAMHeader())

Create a data writer of the SAM file format.

# Arguments
* `output`: data sink
* `header=SAMHeader()`: SAM header object
"""
type SAMWriter{T<:IO} <: Bio.IO.AbstractWriter
    stream::T
end

function Bio.IO.stream(writer::SAMWriter)
    return writer.stream
end

function SAMWriter(output::IO, header::SAMHeader=SAMHeader())
    writer = SAMWriter(output)
    write(writer, header)
    return writer
end

function Base.write(writer::SAMWriter, header::SAMHeader)
    n = 0
    for metainfo in header
        n += write(writer, metainfo)
    end
    return n
end

function Base.write(writer::SAMWriter, metainfo::SAMMetaInfo)
    checkfilled(metainfo)
    return write(writer.stream, metainfo, '\n')
end

function Base.write(writer::SAMWriter, record::SAMRecord)
    checkfilled(record)
    return write(writer.stream, record, '\n')
end


text:
_eltype(::Type{T}) where {T <: AbstractArray} = eltype(T)
_eltype(::Type{T}) where {T} = T
for (i, field) in enumerate((:vertextype, :facetype, :normaltype,
                       :texturecoordinatetype, :colortype))
    @eval begin
        $field(t::Type{T}) where {T <: HomogenousMesh} = _eltype(fieldtype(t, $i))
        $field(mesh::HomogenousMesh) = $field(typeof(mesh))
    end
end

hasvertices(msh) = vertextype(msh) != Nothing
hasfaces(msh) = facetype(msh) != Nothing
hasnormals(msh) = normaltype(msh) != Nothing
hastexturecoordinates(msh) = texturecoordinatetype(msh) != Nothing
hascolors(msh) = colortype(msh) != Nothing

vertices(msh) = msh.vertices
faces(msh) = msh.faces
normals(msh) = msh.normals
texturecoordinates(msh) = msh.texturecoordinates
colors(msh) = msh.color

# Bad, bad name! But it's a little tricky to filter out faces and verts from the attributes, after get_attribute
function attributes_noVF(m::T) where T <: AbstractMesh
    fielddict = Dict{Symbol, Any}(map(fieldnames(T)[3:end]) do field
        field => getfield(m, field)
    end)
    return filter(fielddict) do (key, val)
        val != nothing && val != Nothing[]
    end
end
#Gets all non Nothing attributes from a mesh in form of a Dict fieldname => value
function attributes(m::AbstractMesh)
    filter(((key,val),) -> (val != nothing && val != Nothing[]), all_attributes(m))
end
#Gets all non Nothing attributes types from a mesh type fieldname => ValueType
function attributes(m::Type{M}) where M <: HMesh
    filter(((key,val),) -> (val != Nothing && val != Vector{Nothing}), all_attributes(M))
end

function all_attributes(m::Type{M}) where M <: HMesh
    Dict{Symbol, Any}(map(field -> (field => fieldtype(M, field)), fieldnames(M)))
end
function all_attributes(m::M) where M <: HMesh
    Dict{Symbol, Any}(map(field -> (field => getfield(m, field)),  fieldnames(M)))
end

# Needed to not get into an stack overflow
convert(::Type{M}, mesh::Union{AbstractGeometry, AbstractMesh}) where {M <: AbstractMesh} = M(mesh)
convert(::Type{T}, mesh::T) where T <: AbstractMesh = mesh
# (::Type{HM1}){HM1 <: AbstractMesh}(mesh::HM1) = mesh

isvoid(::Type{T}) where {T} = false
isvoid(::Type{Nothing}) = true
isvoid(::Type{Vector{T}}) where {T} = isvoid(T)

@generated function (::Type{HM1})(primitive::HM2) where {HM1 <: HomogenousMesh, HM2 <: HomogenousMesh}
    fnames = fieldnames(HM1)
    expr = Expr(:call, HM1)
    for i in 1:fieldcount(HM1)
        field = fnames[i]
        target_type = fieldtype(HM1, i)
        source_type = fieldtype(HM2, i)
        if !isconcretetype(fieldtype(HM1, i))  # target is not defined
            push!(expr.args, :(getfield(primitive, $(QuoteNode(field)))))
        elseif !isvoid(target_type) && isvoid(source_type) # target not there yet, maybe we can decompose though (e.g. normals)
            push!(expr.args, :(decompose($(HM1.parameters[i]), primitive)))
        elseif isvoid(target_type)
            push!(expr.args, :($(target_type())))
        else
            push!(expr.args, :(convert($target_type, getfield(primitive, $(QuoteNode(field))))))
        end
    end
    expr
end

#Should be:
#function call{M <: HMesh, VT <: Point, FT <: Face}(::Type{M}, vertices::Vector{VT}, faces::Vector{FT})
# Haven't gotten around to implement the types correctly with abstract types in FixedSizeArrays
function (::Type{M})(
        vertices::AbstractVector{Point{3, VT}}, faces::AbstractVector{FT}
    ) where {M <: HMesh, VT, FT <: Face}
    msh = PlainMesh{VT, FT}(vertices = vertices, faces = faces)
    convert(M, msh)
end
get_default(x::Union{Type, TypeVar}) = nothing
get_default(x::Type{X}) where {X <: AbstractArray} = Nothing[]

"""
generic constructor for abstract HomogenousMesh, infering the types from the keywords (which have to match the field names)
some problems with the dispatch forced me to use this method name... need to further investigate this
"""
function homogenousmesh(attribs::Dict{Symbol, Any})
    newfields = map(fieldnames(HMesh)) do name
        get(attribs, name, get_default(fieldtype(HMesh, name)))
    end
    HomogenousMesh(newfields...)
end

"""
Creates a mesh from keyword arguments, which have to match the field types of the given concrete mesh
"""
(::Type{M})(; kw_args...) where {M <: HMesh} = M(Dict{Symbol, Any}(kw_args))

"""
Creates a new mesh from a dict of `fieldname => value` and converts the types to the given meshtype
"""
function (::Type{M})(attribs::Dict{Symbol, Any}) where M <: HMesh
    newfields = map(fieldnames(HomogenousMesh)) do field
        target_type = fieldtype(M, field)
        default = fieldtype(HomogenousMesh, field) <: AbstractVector ? Nothing[] : nothing
        get(attribs, field, default)
    end
    convert(M, HomogenousMesh(newfields...))
end

"""
Creates a new mesh from an old one, with changed attributes given by the keyword arguments
"""
function (::Type{M})(mesh::AbstractMesh, attributes::Dict{Symbol, Any}) where M <: HMesh
    newfields = map(fieldnames(HomogenousMesh)) do field
        get(attributes, field, getfield(mesh, field))
    end
    HomogenousMesh(newfields...)
end
"""
Creates a new mesh from an old one, with a new constant attribute (like a color)
"""
function (::Type{HM})(mesh::AbstractMesh, constattrib::ConstAttrib) where {HM <: HMesh, ConstAttrib}
    result = Dict{Symbol, Any}()
    for (field, target_type) in zip(fieldnames(HM), HM.parameters)
        if target_type <: ConstAttrib
            result[field] = constattrib
        elseif target_type != Nothing
            result[field] = decompose(target_type, mesh)
        end
    end
    HM(result)
end
function add_attribute(m::AbstractMesh, attribute)
    attribs = attributes(m) # get all attribute values as a Dict fieldname => value
    attribs[:color] = attribute # color will probably be renamed to attribute. not sure yet...
    homogenousmesh(attribs)
end

"""
Creates a new mesh from a tuple of a geometry type and a constant attribute
"""
function (::Type{HM})(x::Tuple{P, ConstAttrib}) where {HM <: HMesh, ConstAttrib, P<:AbstractGeometry}
    any, const_attribute = x
    add_attribute(HM(any), const_attribute)
end

merge(m::AbstractVector{M}) where {M <: AbstractMesh} = merge(m...)

"""
Merges an arbitrary mesh. This function probably doesn't work for all types of meshes
parameters:
`m1` first mesh
`meshes...` other meshes
"""
function merge(m1::M, meshes::M...) where M <: AbstractMesh
    v = copy(m1.vertices)
    f = copy(m1.faces)
    attribs = deepcopy(attributes_noVF(m1))
    for mesh in meshes
        append!(f, map(f-> f .+ length(v), mesh.faces))
        append!(v, mesh.vertices)
        for (v1, v2) in zip(values(attribs), values(attributes_noVF(mesh)))
            append!(v1, v2)
        end
    end
    attribs[:vertices]  = v
    attribs[:faces]     = f
    return M(attribs)
end

"""
A mesh with one constant attribute can be merged as an attribute mesh.
Possible attributes are FSArrays
"""
function merge(
        m1::HMesh{_1, _2, _3, _4, ConstAttrib, _5, _6},
        meshes::HMesh{_1, _2, _3, _4, ConstAttrib, _5, _6}...
    ) where {_1, _2, _3, _4, ConstAttrib <: Colorant, _5, _6}
    vertices     = copy(m1.vertices)
    faces        = copy(m1.faces)
    attribs      = Dict{Symbol, Any}(filter(((k,v),) -> k != :color, attributes_noVF(m1)))
    attribs      = Dict{Symbol, Any}([(k, copy(v)) for (k,v) in attribs])
    color_attrib = RGBA{N0f8}[RGBA{N0f8}(m1.color)]
    index        = Float32[length(color_attrib)-1 for i=1:length(m1.vertices)]
    for mesh in meshes
        append!(faces, map(f-> f .+ length(vertices), mesh.faces))
        append!(vertices, mesh.vertices)
        attribsb = attributes_noVF(mesh)
        for (k,v) in attribsb
            k != :color && append!(attribs[k], copy(v))
        end
        push!(color_attrib, mesh.color)
        append!(index, Float32[length(color_attrib)-1 for i=1:length(mesh.vertices)])
    end
    delete!(attribs, :color)
    attribs[:vertices]      = vertices
    attribs[:faces]         = faces
    attribs[:attributes]    = color_attrib
    attribs[:attribute_id]  = index
    return HMesh{_1, _2, _3, _4, Nothing, typeof(color_attrib), eltype(index)}(attribs)
end

# Fast but slightly ugly way to implement mesh multiplication
# This should probably go into FixedSizeArrays.jl, Vector{FSA} * FSA
struct MeshMulFunctor{T}
    matrix::Mat{4,4,T}
end

(m::MeshMulFunctor{T})(vert) where {T} = Vec{3, T}(m.matrix*Vec{4, T}(vert..., 1))

function *(m::Mat{4,4,T}, mesh::AbstractMesh) where T
    msh = deepcopy(mesh)
    map!(MeshMulFunctor(m), msh.vertices, msh.vertices)
    msh
end

"""
The unnormalized normal of three vertices.
"""
function orthogonal_vector( v1, v2, v3 )
    a = v2 - v1
    b = v3 - v1
    return cross(a, b)
end

"""
```
normals{VT,FD,FT,FO}(vertices::Vector{Point{3, VT}},
                    faces::Vector{Face{FD,FT,FO}},
                    NT = Normal{3, VT})
```
Compute all vertex normals.
"""
function normals(vertices::AbstractVector{Point{3, VT}}, faces::AbstractVector{F},
                 NT = Normal{3, VT}) where {VT, F <: Face}
    normals_result = zeros(Point{3, VT}, length(vertices)) # initilize with same type as verts but with 0
    for face in faces
        v = vertices[face]
        # we can get away with two edges since faces are planar.
        n = orthogonal_vector(v[1], v[2], v[3])
        for i =1:length(F)
            fi = face[i]
            normals_result[fi] = normals_result[fi] + n
        end
    end
    normals_result .= NT.(normalize.(normals_result))
    return normals_result
end

"""
Slice an AbstractMesh at the specified Z axis value.
Returns a Vector of LineSegments generated from the faces at the specified
heights. Note: This will not slice in-plane faces.
"""
function slice(mesh::AbstractMesh{Point{3,VT}, Face{3, FT}}, height::Number) where {VT<:AbstractFloat,FT<:Integer}

    height_ct = length(height)
    # intialize the LineSegment array
    slice = Simplex{2,Point{2,VT}}[]

    for face in mesh.faces
        v1,v2,v3 = mesh.vertices[face]
        zmax = max(v1[3], v2[3], v3[3])
        zmin = min(v1[3], v2[3], v3[3])
        if height > zmax
            continue
        elseif zmin <= height
            if v1[3] < height && v2[3] >= height && v3[3] >= height
                p1 = v1
                p2 = v3
                p3 = v2
            elseif v1[3] > height && v2[3] < height && v3[3] < height
                p1 = v1
                p2 = v2
                p3 = v3
            elseif v2[3] < height && v1[3] >= height && v3[3] >= height
                p1 = v2
                p2 = v1
                p3 = v3
            elseif v2[3] > height && v1[3] < height && v3[3] < height
                p1 = v2
                p2 = v3
                p3 = v1
            elseif v3[3] < height && v2[3] >= height && v1[3] >= height
                p1 = v3
                p2 = v2
                p3 = v1
            elseif v3[3] > height && v2[3] < height && v1[3] < height
                p1 = v3
                p2 = v1
                p3 = v2
            else
                continue
            end

            start = Point{2,VT}(p1[1] + (p2[1] - p1[1]) * (height - p1[3]) / (p2[3] - p1[3]),
                                p1[2] + (p2[2] - p1[2]) * (height - p1[3]) / (p2[3] - p1[3]))
            finish = Point{2,VT}(p1[1] + (p3[1] - p1[1]) * (height - p1[3]) / (p3[3] - p1[3]),
                                 p1[2] + (p3[2] - p1[2]) * (height - p1[3]) / (p3[3] - p1[3]))

            push!(slice, Simplex{2,Point{2, VT}}(start, finish))
        end
    end

    return slice
end

# TODO this should be checkbounds(Bool, ...)
"""
```
checkbounds{VT,FT,FD,FO}(m::AbstractMesh{VT,Face{FD,FT,FO}})
```

Check the `Face` indices to ensure they are in the bounds of the vertex
array of the `AbstractMesh`.
"""
function Base.checkbounds(m::AbstractMesh{VT, Face{FD, FT}}) where {VT, FD, FT}
    isempty(faces(m)) && return true # nothing to worry about I guess
    flat_inds = reinterpret(FT, faces(m))
    checkbounds(Bool, vertices(m), flat_inds)
end

function unique_verts!(verts::AbstractVector{T}) where T
    table = Dict{T, Int}()
    new_idx = 0
    for (i, v) in enumerate(verts)
        if !haskey(table, v)
            new_idx += 1
            table[v] = new_idx
            # if index has moved, inplace rewrite verts
            new_idx != i && (verts[new_idx] = v)
        end
    end
    if new_idx != length(verts)
        resize!(verts, new_idx)
    end
    table
end

function reface!(point2idx, verts, uverts, faces)
    map!(faces, faces) do face
        map(face) do i
            point2idx[verts[i]]
        end
    end
end

"""
Calculate the signed volume of one tetrahedron. Be sure the orientation of your surface is right.
"""
function volume(
        vertices::AbstractVector{Point{3, VT}},
        face::Face{3,FT}
    ) where {VT,FT}
    v1, v2, v3 = vertices[face]
    sig = sign( orthogonal_vector( v1, v2, v3 ) ⋅ v1 )
    return sig * abs( v1 ⋅ ( v2 × v3 ) ) / 6
end

"""
Calculate the signed volume of all tetrahedra. Be sure the orientation of your surface is right.
"""
function volume(
        vertices::AbstractVector{Point{3, VT}},
        faces::AbstractVector{Face{3,FT}}
    ) where {VT,FT}
    return sum( x->volume( vertices, x), faces )
end

"""
    remove_overlap!(mesh::AbstractMesh)

removes non unique vertices from a mesh, and relinks the faces to point to only shared vertices.
"""
function remove_overlap!(mesh::AbstractMesh)
    verts = vertices(mesh)
    orig_verts = copy(verts)
    table = unique_verts!(verts)
    reface!(table, orig_verts, verts, faces(mesh))
    return
end

# functions related to displaying types

function show(io::IO, m::M) where M <: HMesh
    println(io, "HomogenousMesh(")
    for (key,val) in attributes(m)
        print(io, "    ", key, ": ", length(val), "x", eltype(val), ", ")
    end
    println(io, ")")
end

function (meshtype::Type{T})(c::Pyramid) where T <: AbstractMesh
    T(vertices = decompose(vertextype(T), c), faces = decompose(facetype(T), c))
end

"""
Standard way of creating a mesh from a primitive.
Just walk through all attributes of the mesh and try to decompose it.
If there are attributes missing, just hope it will get managed by the mesh constructor.
(E.g. normal calculation, which needs to have vertices and faces present)
"""
function (meshtype::Type{T})(c::GeometryPrimitive, args...) where T <: AbstractMesh
    attribs = attributes(T)
    newattribs = Dict{Symbol, Any}()
    for (fieldname, typ) in attribs
        if isdecomposable(eltype(typ), c)
            newattribs[fieldname] = decompose(eltype(typ), c, args...)
        end
    end
    return T(newattribs)
end

function to_pointn(::Type{Point{N, T}}, p::StaticVector{N2}, d = T(0)) where {T, N, N2}
    return Point{N, T}(ntuple(i-> i <= N2 ? T(p[i]) : T(d), N))
end

function (::Type{T})(c::Circle, n = 32) where T <: AbstractMesh
    newattribs = Dict{Symbol, Any}()
    VT = vertextype(T)
    verts = decompose(VT, c, n)
    N = length(verts)
    push!(verts, to_pointn(VT, origin(c))) # middle point
    middle_idx = length(verts)
    FT = facetype(T)
    faces = map(1:N) do i
        FT(i, middle_idx, i + 1)
    end
    return T(vertices = verts, faces = faces)
end

function (meshtype::Type{T})(
        c::Union{HyperCube{3,T}, HyperRectangle{3,HT}}
    ) where {T <: HMesh, HT}
    xdir = Vec{3, HT}(1, 0, 0)
    ydir = Vec{3, HT}(0, 1, 0)
    zdir = Vec{3, HT}(0, 0, 1)
    o0 = zero(Vec{3, HT})
    quads = [
        Quad(zdir, xdir, ydir), # Top
        Quad(o0,   ydir, xdir), # Bottom
        Quad(xdir, ydir, zdir), # Right
        Quad(o0,   zdir, ydir), # Left
        Quad(o0,   xdir, zdir), # Back
        Quad(ydir, zdir, xdir) #Front
    ]
    mesh = merge(map(meshtype, quads))
    v = vertices(mesh)
    w = widths(c)
    o = origin(c)
    map!(v, v) do v
        (v .* w) + o
    end
    return mesh
end

==(a::AbstractMesh, b::AbstractMesh) = false
function ==(a::M, b::M) where M <: AbstractMesh
    for ((ka, va), (kb, vb)) in zip(all_attributes(a), all_attributes(b))
        va != vb && return false
    end
    true
end