text:
using Test
using Caching

# Make file
FILE, FILEid = mktemp()
close(FILEid)
_tmpdir = tempdir()

# Test Cache, @cache
@testset "Cache" begin
    # Function arguments
    _a_float = 0.1234
    _a_float_2 = 0.5678
    _an_int = 1234

    # Define functions
    global foo(x) = x
    bar(x; y=1) = x + y

    # Test caching of simple functions
    foo_c1 = Cache(foo)
    @test typeof(foo_c1) <: AbstractCache
    @test typeof(foo_c1) <: Cache
    @test typeof(foo_c1.cache) <: Dict
    @test typeof(foo_c1.offsets) <: Dict
    @test typeof(foo_c1.filename) <: AbstractString
    @test typeof(foo_c1.name) <: AbstractString
    @test typeof(foo_c1.func) <: Function
    @test foo_c1.func == foo
    @test foo(_a_float) == foo_c1(_a_float)
    @test length(foo_c1.cache) == 1

    # Test functions with keyword arguments
    bar_c1 = Cache(bar)
    @test typeof(bar_c1) <: AbstractCache
    @test typeof(bar_c1) <: Cache
    @test typeof(bar_c1.cache) <: Dict
    @test typeof(bar_c1.offsets) <: Dict
    @test typeof(bar_c1.filename) <: AbstractString
    @test typeof(bar_c1.name) <: AbstractString
    @test typeof(bar_c1.func) <: Function
    @test bar_c1.func == bar
    @test bar(_a_float) == bar_c1(_a_float)
    @test bar(_a_float; y=_a_float_2) == bar_c1(_a_float; y=_a_float_2)
    @test length(bar_c1.cache) == 2

	# Test macro support
    foo_m1 = @eval @cache foo $FILE  # no type annotations
    foo_m2 = @eval @cache foo::Int $FILE  # standard type annotation
	for _foo in (foo_m1, foo_m2)
    	@test typeof(_foo) <: AbstractCache
    	@test typeof(_foo) <: Cache
        @test typeof(_foo.cache) <: Dict
        @test typeof(_foo.offsets) <: Dict
        @test typeof(_foo.filename) <: AbstractString
        @test typeof(_foo.name) <: AbstractString
        @test typeof(_foo.func) <: Function
        @test _foo.func == foo
		@test _foo(_an_int) == foo(_an_int)	 # all should work with Int's
		if !(_foo === foo_m1)
		# `foo_m2` fails with arguments other than `Int`
		    @test_throws AssertionError _foo(_a_float)
	        @test _foo.cache isa Dict{UInt64, Int}
		else
			@test _foo(_a_float) == foo(_a_float)
			@test _foo.cache isa Dict{UInt64, Any}
		end
	end

    dc = @eval @cache foo $(joinpath(_tmpdir,"somefile.bin"))
    for i in 1:3 dc(i); end
    @persist! dc
    @empty! dc
    @test length(dc.offsets) == 3
    @test isempty(dc.cache)
    for i in 4:6 dc(i); end
    for i in 1:6
        @test dc(i) == foo(i) == i
    end
    @empty! dc true
end

# syncache!, @syncache!
@testset "syncache!, @syncache!" begin
    # Define functions
    foo(x) = x+1  # define function

    n1 = 5
    fc = @eval @cache foo $FILE  # make a cache object
    [fc(i) for i in 1:n1]  # populate the memorycache
    @persist! fc  # write to disk the cache
    @empty! fc  # delete the memory cache
    @test isempty(fc.cache)
    @syncache! fc "disk"  # load cache from disk
    @test isfile(fc.filename)
    @test length(fc.cache) == n1
    @empty! fc  # 5 entries on disk, 0 in memory

    n2 = 3
    [fc(-i) for i in 1:n2]  # populate the memory cache
    @syncache! fc "memory"  # write memory cache to disk
    @test length(fc.cache) == n2
    @empty! fc
    @test isempty(fc.cache)
    @syncache! fc "disk"    # load cache from disk
    @test length(fc.cache) == n1 + n2
    @empty! fc true  # remove everything

    [fc(i) for i in 1:n1]  # populate the memory cache
    @syncache! fc "memory"  # write to disk
    @empty! fc
    [fc(-i) for i in 1:n1]  # populate the memory cache
    @syncache! fc "both"     # sync both memory and disk
    @test length(fc.offsets) == 2*n1
    @test length(fc.cache) == 2*n1
    @empty! fc true
    @test !isfile(fc.filename)
end

# empty!, @empty!
@testset "empty!, @empty!" begin
    # Define functions
    foo(x) = x
    bar(x; y=1) = x + y

    dc =@eval @cache foo $FILE
    dc(1)
    @test length(dc.cache) == 1
    @persist! dc "somefile.bin"
    @test length(dc.offsets) == 1
    @test isfile("somefile.bin")
    @empty! dc true  # remove offsets and file
    @test isempty(dc.cache)
    @test isempty(dc.offsets)
    @test !isfile("somefile.bin")
end

# persist!, @persist!
@testset "persist!, @persist!" begin
    # Define functions
    foo(x) = x
    global dc
    dc = @eval @cache foo $FILE
    @test dc.filename == abspath(FILE)
    @test isempty(dc.offsets)
    N = 3
    for i in 1:N dc(i); end  # add N entries

    @persist! dc
    @test length(dc.offsets) == N
    @test isfile(FILE)
    @test dc.filename == abspath(FILE)
    buf = open(read, FILE)
    rm(FILE)
    FILE2 = joinpath(_tmpdir, "some_other_file.bin")
    @eval @persist! dc $FILE2
    @test length(dc.offsets) == N
    @test isfile(FILE2)
    @test dc.filename == abspath(FILE2)
    buf2 = open(read, FILE2)
    rm(FILE2)

    @test buf == buf2  # sanity check
end

# Compression
@testset "Compression" begin
    files = joinpath.(_tmpdir, ["somefile.bin",
                                "somefile.bin.gz",
                                "somefile.bin.gzip",
                                "somefile.bin.bz2",
                                "somefile.bin.bzip2"])
    for _file in files
        global foo(x) = x
        dc = @eval @cache foo $_file
        data = [1,2,3]
        dc(data)
        @persist! dc
        @empty! dc
        try
            @test dc(data) == foo(data)
            @test true
        catch
            @test false
        end
        @empty! dc true
    end
end

# Size constraints
@testset "Size constraints: number of objects" begin
    file = joinpath(_tmpdir, "somefile.bin")
    foo(x) = x
    dc = @eval @cache foo $file 3  # 3 objects max
    dc(1)
    dc(2)
    dc(3)
    #--> Cache is full
    dc(4)  #--> 1 is removed
    @test !(1 in values(dc.cache)) &&
        all(i in values(dc.cache) for i in 2:4)
    @persist! dc
    @empty! dc
    for i in 5:6 dc(i) end
    #--> 2,3,4 on disk, 5,6 in memory
    @syncache! dc #--> brings 4 in memory and 5,6 on disk
    @test all(i in values(dc.cache) for i in 4:6)
    @test length(dc) == 3 # sanity check
    @empty! dc  #--> nothing in memory
    @syncache! dc "disk"  #--> load from disk max entries i.e. 3
    @test all(i in values(dc.cache) for i in 4:6)
    @empty! dc true  # remove everything
end

@testset "Size constraints: bytes of memory" begin
    file = joinpath(_tmpdir, "somefile.bin")
    foo(x) = x
    dc = @eval @cache foo $file 1.0  # 1024 bytes=1.0 KiB max
    for i in 1:128 dc(i) end
    #--> Cache is full (128 x 8bytes/number=1024 bytes)
    dc(129)  #--> 1 is removed
    @test !(1 in values(dc.cache)) &&
        all(i in values(dc.cache) for i in 2:129)
    @persist! dc
    @empty! dc  #--> 2,...,129 on disk, nothing in memory
    for i in 130:130+126 dc(i) end  # write 127 entries
    #--> 130,..,256 in memory, 2,...,129 on disk
    @syncache! dc #--> brings 129 in memory and 130,...,256 on disk
    @test all(i in values(dc.cache) for i in 129:130+126)
    @test length(dc) == 128 # sanity check
    @empty! dc  #--> nothing in memory
    @syncache! dc "disk"  #--> load from disk max entries i.e. 128
    @test all(i in values(dc.cache) for i in 129:130+126)
    @empty! dc true  # remove everything
end

@testset "@cache <func. def.>" begin
    @test !@isdefined foo
    @cache foo=x->x
    @test @isdefined foo
    @test typeof(foo.cache) == Dict{UInt,Any}
    @test foo(1) == 1

    T = Int
    @test !@isdefined bar
    @eval @cache bar=x::$T->x
    @test @isdefined bar
    @eval @test typeof(bar.cache) == Dict{UInt,$T}
    @test bar(1) == 1

    @test !@isdefined baz
    @eval @cache function baz(x)::$T x end
    @test @isdefined baz
    @eval @test typeof(baz.cache) == Dict{UInt,$T}
    @test baz(1) == 1
end

# Hashing function
struct custom_type{T}
    x::T
end
@testset "Hashing function" begin
    num = 1
    str = "a string"
    # test hashing function for various objects
    for object in (Int(num), UInt(num), Float64(num),
                   [num num], (num, num),
                   str, [str, str], (str,str),
                   custom_type(num), custom_type(str),
                   [custom_type(num)], (custom_type(num),),
                   custom_type(custom_type(num))
                  )
        object2 = deepcopy(object)
        @test arghash(object) == arghash(object2)
    end
end

# show methods
@testset "Show methods" begin
    buf = IOBuffer()
    foo(x) = x
    dc = @eval @cache foo $FILE
    try
        show(buf, dc)
        @test true
    catch
        @test false
    end
end


text:
include("testUncertain.jl")
include("testQuantities.jl")

ext: jl


text:
outcome = "bart"

include(abspath("src", "04_full-extent.jl"))

## SCBD
# Prep values
inds_obs = _indsobs(Y)
Yobs = _Yobs(Y, inds_obs)
Ytransf = _Ytransf(Yobs)
BDstats = betadiv(Ytransf)

# Get SCBD values
scbd = vec(BDstats.SCBDj)

# Visualize
histogram(scbd; bins=20)
scatter(scbd)

## Species counts
# Count sites with occurrences per species
specounts = vec(sum(Yobs; dims=1))

# Visualize
histogram(specounts; bins=20)
occupancy_plot = scatter(
    specounts;
    formatter=:plain,
    legend=:none,
    xlabel="Species ID number",
    ylabel="Species occupancy",
)

## Combine stuff
scbd_plot = scatter(
    specounts,
    scbd;
    xlabel="Species occupancy",
    ylabel="SCBD value",
    formatter=:plain,
    label=:none,
)
scatter!([NaN]; label="r = $(round(cor(specounts, scbd), digits = 3))")

## Save figures
savefig(occupancy_plot, joinpath(".", "fig", "bart", "x_sp-occupancy.png"))
savefig(scbd_plot, joinpath(".", "fig", "bart", "x_sp-scbd.png"))

max_stars_repo_name: gabrieldansereau/betadiversity-hotspots

max_issues_repo_name: gabrieldansereau/betadiversity-hotspots

max_forks_repo_name: gabrieldansereau/betadiversity-hotspots


text:
using Base.Test
using MIPVerify
using MIPVerify: UnrestrictedPerturbationFamily, BlurringPerturbationFamily
using MIPVerify: LInfNormBoundedPerturbationFamily
isdefined(:TestHelpers) || include("../../../TestHelpers.jl")
using TestHelpers: batch_test_adversarial_example

@testset "Conv + FC + Softmax" begin

### Parameters for neural net
batch = 1
in1_height = 8
in1_width = 8
stride1_height = 2
stride1_width = 2
strides1 = (1, stride1_height, stride1_width, 1)
pooled1_height = round(Int, in1_height/stride1_height, RoundUp)
pooled1_width = round(Int, in1_width/stride1_width, RoundUp)
in1_channels = 1
filter1_height = 2
filter1_width = 2
out1_channels = 2

A_height = 5
A_width = pooled1_height*pooled1_width*out1_channels

B_height = 3
B_width = A_height

srand(5)
input_size = (batch, in1_height, in1_width, in1_channels)
x0 = rand(input_size)

kernelc1 = rand(filter1_height, filter1_width, in1_channels, out1_channels)*2-1
biasc1 = rand(out1_channels)*2-1

kernelf1 = rand(-10:10, A_width, A_height)
biasf1 = rand(-10:10, A_height)

kernelf2 = rand(B_width, B_height)*2-1
biasf2 = rand(B_height)*2-1

nn = Sequential(
    [
        Conv2d(kernelc1, biasc1), MaxPool(strides1), ReLU(),
        Flatten(4),
        Linear(kernelf1, biasf1), ReLU(),
        Linear(kernelf2, biasf2)
    ],
    "tests.integration.generated_weights.conv+fc+softmax"
)

pp_blur = BlurringPerturbationFamily((5, 5))
pp_unrestricted = UnrestrictedPerturbationFamily()

expected_objective_values = Dict(
    (1, pp_blur, 1, 0) => 0,
    (1, pp_blur, Inf, 0) => 0,
    (2, pp_unrestricted, 1, 0) => 2.98266,
    (2, pp_unrestricted, 1, 0.1) => 3.03465,
    (2, pp_unrestricted, Inf, 0) => 0.235631,
    (2, LInfNormBoundedPerturbationFamily(0.23), Inf, 0) => NaN,
    (2, LInfNormBoundedPerturbationFamily(0.236), Inf, 0) => 0.235631,
    (2, pp_unrestricted, Inf, 0.1) => 0.240124,
    (2, pp_blur, 1, 0) => NaN,
    (2, pp_blur, Inf, 0) => NaN,
    (3, pp_unrestricted, Inf, 0) => 0.00288325,
    (3, pp_unrestricted, Inf, 1) => 0.0110296,
    (3, pp_blur, 1, 0) => 0.261483,
    (3, pp_blur, 1, 10) => NaN,
    (3, pp_blur, Inf, 0) => 0.0105534,
    (3, pp_blur, Inf, 10) => NaN,
    ([2, 3], pp_unrestricted, Inf, 0) => 0.00288325,
    ([2, 3], pp_unrestricted, Inf, 1) => 0.0110296,
    ([2, 3], pp_blur, 1, 0) => 0.261483,
    ([2, 3], pp_blur, Inf, 0) => 0.0105534
)

if Pkg.installed("Gurobi") == nothing
    # Skip some tests if Gurobi is not installed.
    selected_test_keys = [
        (2, pp_unrestricted, Inf, 0),
        (2, LInfNormBoundedPerturbationFamily(0.23), Inf, 0),
        (2, LInfNormBoundedPerturbationFamily(0.236), Inf, 0),
        (3, pp_blur, 1, 0),
        ([2, 3], pp_blur, 1, 0)
    ]
    expected_objective_values = filter((k, v) -> k in selected_test_keys, expected_objective_values)
end

batch_test_adversarial_example(nn, x0, expected_objective_values)

end


text:
function APMParser()
    return APMParser(())
end

function APMParser(arg0::DataContext)
    return APMParser((DataContext,), arg0)
end

function parse(obj::APMParser, arg0::InputStream)
    return jcall(obj, "parse", APMFile, (InputStream,), arg0)
end

function parse(obj::APMParser, arg0::InputStream, arg1::JString)
    return jcall(obj, "parse", APMFile, (InputStream, JString), arg0, arg1)
end

function parse(obj::APMParser, arg0::JString)
    return jcall(obj, "parse", APMFile, (JString,), arg0)
end

function with_conventions(obj::APMParser, arg0::IERSConventions)
    return jcall(obj, "withConventions", APMParser, (IERSConventions,), arg0)
end

function with_data_context(obj::APMParser, arg0::DataContext)
    return jcall(obj, "withDataContext", APMParser, (DataContext,), arg0)
end

function with_international_designator(obj::APMParser, arg0::jint, arg1::jint, arg2::JString)
    return jcall(obj, "withInternationalDesignator", APMParser, (jint, jint, JString), arg0, arg1, arg2)
end

function with_mission_reference_date(obj::APMParser, arg0::AbsoluteDate)
    return jcall(obj, "withMissionReferenceDate", APMParser, (AbsoluteDate,), arg0)
end

function with_mu(obj::APMParser, arg0::jdouble)
    return jcall(obj, "withMu", APMParser, (jdouble,), arg0)
end

function with_simple_eop(obj::APMParser, arg0::jboolean)
    return jcall(obj, "withSimpleEOP", APMParser, (jboolean,), arg0)
end


text:
for (cfunction, jfunction, ntask) in (
        (CreateDIChan, :digital_input, :DITask),
        (CreateDOChan, :digital_output, :DOTask))
    @eval function $jfunction(channel::ASCIIString)
      t = $ntask()
      $jfunction(t, channel)
      t
    end
    @eval function $jfunction(t::$ntask, channel::ASCIIString)
        catch_error( $cfunction(t.th,
                pointer(channel),
                pointer(""),
                Val_ChanPerLine) )
        nothing
    end
end

@doc """
`digital_output(channel) -> task`

`digital_output(task, channel)`

create a digital output channel, and a new NIDAQ task if one is not specified
""" digital_output

@doc """
`digital_input(channel) -> task`

`digital_input(task, channel)`

create a digital input channel, and a new task if one is not specified
""" digital_input

read_digital_cfunctions = Dict{Type,Function}(
    UInt8 => ReadDigitalU8,
    UInt16 => ReadDigitalU16,
    UInt32 => ReadDigitalU32 )

function read(t::DITask, precision::DataType, num_samples_per_chan::Integer = -1)
    num_channels = getproperties(t)["NumChans"][1]
    num_samples_per_chan_read = Int32[0]
    if num_samples_per_chan==-1;  num_samples_per_chan=1024;  end
    data = Array(precision, num_samples_per_chan*num_channels)
    catch_error( read_digital_cfunctions[precision](t.th,
        convert(Int32,num_samples_per_chan),
        1.0,
        reinterpret(Bool32,Val_GroupByChannel),
        pointer(data),
        convert(UInt32,num_samples_per_chan*num_channels),
        pointer(num_samples_per_chan_read),
        reinterpret(Ptr{Bool32},C_NULL)) )
    data = data[1:num_samples_per_chan_read[1]*num_channels]
    num_channels==1 ? data : reshape(data, (num_samples_per_chan, convert(Int64,num_channels)))
end

for (cfunction, types) in (
        (WriteDigitalU8,  UInt8),
        (WriteDigitalU16, UInt16),
        (WriteDigitalU32, UInt32))
    @eval function write(t::DOTask, data::Matrix{$types})
        num_samples_per_chan::Int32 = size(data, 1)
        data = reshape(data, length(data))
        num_samples_per_chan_written = Int32[0]
        catch_error( $cfunction(t.th,
            num_samples_per_chan,
            reinterpret(Bool32,UInt32(false)),
            1.0,
            reinterpret(Bool32,Val_GroupByChannel),
            pointer(data),
            pointer(num_samples_per_chan_written),
            reinterpret(Ptr{Bool32},C_NULL)) )
        num_samples_per_chan_written[1]
    end
    @eval write(t::DOTask, data::Vector{$types}) =
        write(t, reshape(data, (length(data),1)))
end


text:
import Base: show, push!, start, next, done

## Store

abstract Store

function push!(store::Store, k::Key, v::Value)
    error("Method `push!` is not implemented for store of type $(typeof(store))")
end

function prepare!(store::Store)
    error("Method `prepare` is not implemented for store " *
          "of type $(typeof(store))")
end

function start(store::Store)
    error("Iteration is not implemented for store " *
          "of type $(typeof(store))")
end

# MemStore

abstract MemStore <: Store

type ArrayMemStore <: MemStore
    data::Vector{PokeRecord}
    next_idx::Int
    ArrayMemStore(capacity=1_000_000) = new(Array(PokeRecord, capacity), 1)
end

function show(io::IO, store::ArrayMemStore)
    print(io, "ArrrayMemStore($(store.next_idx-1)/$(length(store.data)))")
end

function push!(store::ArrayMemStore, k::Key, v::Value)
    store.data[store.next_idx] = PokeRecord(k, v)
    store.next_idx += 1
end

function prepare!(store::ArrayMemStore)
    sort!(store.data)
end

function start(store::ArrayMemStore)
    it = uniquesorted(store.data)
    s = start(it)
    return (it, s)  # note: putting both - iterator and actual state into
                    # ArrayMemStore's "state"
end

function next(store::ArrayMemStore, state)
    it, s = state
    x, new_s = next(it, s)
    return x, (it, new_s)
end

function done(store::ArrayMemStore, state)
    it, s = state
    return done(it, s)
end

# FileStore

type FileStore <: Store
    path::AbstractString
end

function push!(store::FileStore)
    error("FileStore doesn't support direct pushing, " *
          "use `createdump` or `mergedump` instead")
end

function prepare!(store::FileStore)
    # do nothing
end

function start(store::FileStore)
    # TODO: we open file, but don't close in iterator
    # need to handle it somehow
    io = open(store.path)
    it = PokeIterator(io)
    s = start(it)
    return (it, s)
end

function next(store::FileStore, iter_state::Tuple)
    it, s = iter_state
    x, new_s = next(it, s)
    return x, (it, new_s)
end

function done(store::FileStore, iter_state::Tuple)
    it, s = iter_state
    return done(it, s)
end

function createdump(dumpf::IOStream, memstore::MemStore)
    for rec in memstore
        writeobj(dumpf, rec)
    end
end

function mergedump(outf::IOStream, filestore::FileStore, memstore::MemStore)
    # important: memstore first since it has more recent records
    merged = mergesorted(memstore, filestore)
    for rec in merged
        writeobj(outf, rec)
    end
end

# TODO:
# 1) block read for PokeIterator
# 2) mergesorted for multiple iterators (min(map(head, iterators)); move that iter; if done, remove from list)

function perf_test()
    N = 1_000_000
    println("$(N * 110 / 1024 / 1024) Mb")
    K = [rand(UInt8, 10) for i=1:N]
    V = [rand(UInt8, 100) for i=1:N]
    KV = collect(zip(K, V))
    A = Array(Tuple{Vector{UInt8},Vector{UInt8}}, N)
    SD = create_memstore()
    @time @inbounds for i=1:N A[i] = KV[i] end
    @time sort(A, lt=(x,y) -> (bytestring(x[1]) < (bytestring(y[1]))))
    @time sort(A)
    @time @inbounds for i=1:N SD[K[i]] = V[i] end
    @time @inbounds for i=1:N SD[K[i]] = V[i] end
    @time collect(SD)
end

function main()
    memstore = ArrayMemStore()
    @time for i=1:1_000_000
        push!(memstore, rand(UInt8, 10), rand(UInt8, 100))
    end
    @time prepare!(memstore)
    @time open("/tmp/dump", "w") do dumpf
        createdump(dumpf, memstore)
    end

    fstore = FileStore("/tmp/new_dump")
    @time open("/tmp/dump", "w") do new_dumpf
        @profile mergedump(new_dumpf, fstore, memstore)
    end
end


text:
include("../utils.jl")

@testset "Multinomial doubly stochastic Matrices" begin
    M = MultinomialDoubleStochastic(3)
    @test manifold_dimension(M) == 4
    @test repr(M) == "MultinomialDoubleStochastic(3)"
    p = ones(3, 3) ./ 3
    X = zeros(3, 3)
    @test is_point(M, p)
    @test is_vector(M, p, X)
    pf1 = [0.1 0.9 0.1; 0.1 0.9 0.1; 0.1 0.1 0.9] #not sum 1
    @test_throws CompositeManifoldError is_point(M, pf1, true)
    pf2r = [0.1 0.9 0.1; 0.8 0.05 0.15; 0.1 0.05 0.75]
    @test_throws DomainError is_point(M, pf2r, true)
    @test_throws CompositeManifoldError is_point(M, pf2r', true)
    pf3 = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0] # contains nonpositive entries
    @test_throws CompositeManifoldError is_point(M, pf3, true)
    Xf2c = [-0.1 0.0 0.1; -0.2 0.1 0.1; 0.2 -0.1 -0.1] #nonzero columns
    @test_throws CompositeManifoldError is_vector(M, p, Xf2c, true)
    @test_throws DomainError is_vector(M, p, Xf2c', true)
    @test representation_size(M) == (3, 3)
    pE = similar(p)
    embed!(M, pE, p)
    pE2 = embed(M, p)
    @test pE == p
    @test pE2 == p
    @test_throws DomainError project(M, -ones(3, 3))
    @test project(M, p) == p
    p2 = [0.1 0.2 0.7; 0.2 0.7 0.1; 0.7 0.1 0.2]
    p3 = [0.1 0.4 0.5; 0.4 0.5 0.1; 0.5 0.1 0.4]

    X2 = [-0.1 0.0 0.1; 0.0 0.2 -0.2; 0.1 -0.2 0.1]
    @test is_vector(
        M,
        p2,
        vector_transport_to(M, p, X2, p2, ProjectionTransport());
        atol=10^-15,
    )
    X3 = [1.0 1.0 1.0; 0.0 0.0 0.0; 0.0 0.0 0.0]
    @test inner(M, p, X3, X3) == 9.0

    types = [Matrix{Float64}]
    TEST_FLOAT32 && push!(types, Matrix{Float32})
    TEST_STATIC_SIZED && push!(types, MMatrix{4,2,Float64,8})

    for T in types
        pts = [convert(T, p), convert(T, p2), convert(T, p3)]
        @testset "Type $T" begin
            test_manifold(
                M,
                pts,
                test_injectivity_radius=false,
                test_reverse_diff=false,
                test_forward_diff=false,
                test_project_tangent=true,
                test_exp_log=false,
                test_default_vector_transport=true,
                vector_transport_methods=[ProjectionTransport()],
                default_inverse_retraction_method=nothing,
                default_retraction_method=ProjectionRetraction(),
                is_tangent_atol_multiplier=20,
                is_point_atol_multiplier=20,
                test_inplace=true,
            )
        end
    end
end


text:
################################################################################
#                               Type Conversions                               #
################################################################################
const RangeLike = Union{AbstractRange, AbstractVector, ClosedInterval}

# if no plot type based conversion is defined, we try using a trait
function convert_arguments(T::PlotFunc, args...; kw...)
    ct = conversion_trait(T)
    try
        convert_arguments(ct, args...; kw...)
    catch e
        if e isa MethodError
            try
                convert_arguments_individually(T, args...)
            catch ee
                if ee isa MethodError
                    error(
                        """
                        `Makie.convert_arguments` for the plot type $T and its conversion trait $ct was unsuccessful.

                        The signature that could not be converted was:
                        $(join("::" .* string.(typeof.(args)), ", "))

                        Makie needs to convert all plot input arguments to types that can be consumed by the backends (typically Arrays with Float32 elements).
                        You can define a method for `Makie.convert_arguments` (a type recipe) for these types or their supertypes to make this set of arguments convertible (See http://makie.juliaplots.org/stable/documentation/recipes/index.html).

                        Alternatively, you can define `Makie.convert_single_argument` for single arguments which have types that are unknown to Makie but which can be converted to known types and fed back to the conversion pipeline.
                        """
                    )
                else
                    rethrow(ee)
                end
            end
        else
            rethrow(e)
        end
    end
end

# in case no trait matches we try to convert each individual argument
# and reconvert the whole tuple in order to handle missings centrally, e.g.
function convert_arguments_individually(T::PlotFunc, args...)
    # convert each single argument until it doesn't change type anymore
    single_converted = recursively_convert_argument.(args)
    # if the type of args hasn't changed this function call didn't help and we error
    if typeof(single_converted) == typeof(args)
        throw(MethodError(convert_arguments, (T, args...)))
    end
    # otherwise we try converting our newly single-converted args again because
    # now a normal conversion method might work again
    convert_arguments(T, single_converted...)
end

function recursively_convert_argument(x)
    newx = convert_single_argument(x)
    if typeof(newx) == typeof(x)
        x
    else
        recursively_convert_argument(newx)
    end
end

################################################################################
#                          Single Argument Conversion                          #
################################################################################

# if no specific conversion is defined, we don't convert
convert_single_argument(x) = x

# replace missings with NaNs
function convert_single_argument(a::AbstractArray{<:Union{Missing, <:Real}})
    [ismissing(x) ? NaN32 : convert(Float32, x) for x in a]
end

# same for points
function convert_single_argument(a::AbstractArray{<:Union{Missing, <:Point{N}}}) where N
    [ismissing(x) ? Point{N, Float32}(NaN32) : Point{N, Float32}(x) for x in a]
end

################################################################################
#                                  PointBased                                  #
################################################################################

"""
Wrap a single point or equivalent object in a single-element array.
"""
function convert_arguments(::PointBased, position::VecTypes{N, <: Number}) where N
    ([convert(Point{N, Float32}, position)],)
end

function convert_arguments(::PointBased, positions::AbstractVector{<: VecTypes{N, <: Number}}) where N
    (elconvert(Point{N, Float32}, positions),)
end

function convert_arguments(::PointBased, positions::SubArray{<: VecTypes, 1})
    # TODO figure out a good subarray solution
    (positions,)
end

"""
Enables to use scatter like a surface plot with x::Vector, y::Vector, z::Matrix
spanning z over the grid spanned by x y
"""
function convert_arguments(::PointBased, x::AbstractVector, y::AbstractVector, z::AbstractMatrix)
    (vec(Point3f.(x, y', z)),)
end
"""
    convert_arguments(P, x, y, z)::(Vector)

Takes vectors `x`, `y`, and `z` and turns it into a vector of 3D points of the values
from `x`, `y`, and `z`.
`P` is the plot Type (it is optional).
"""
convert_arguments(::PointBased, x::RealVector, y::RealVector, z::RealVector) = (Point3f.(x, y, z),)

"""
    convert_arguments(P, x)::(Vector)

Takes an input GeometryPrimitive `x` and decomposes it to points.
`P` is the plot Type (it is optional).
"""
convert_arguments(p::PointBased, x::GeometryPrimitive) = convert_arguments(p, decompose(Point, x))

function convert_arguments(::PointBased, pos::AbstractMatrix{<: Number})
    (to_vertices(pos),)
end

convert_arguments(P::PointBased, x::AbstractVector{<:Real}, y::AbstractVector{<:Real}) = (Point2f.(x, y),)

convert_arguments(P::PointBased, x::AbstractVector{<:Real}, y::AbstractVector{<:Real}, z::AbstractVector{<:Real}) = (Point3f.(x, y, z),)

"""
    convert_arguments(P, y)::Vector
Takes vector `y` and generates a range from 1 to the length of `y`, for plotting on
an arbitrary `x` axis.

`P` is the plot Type (it is optional).
"""
convert_arguments(P::PointBased, y::RealVector) = convert_arguments(P, keys(y), y)

"""
    convert_arguments(P, x, y)::(Vector)

Takes vectors `x` and `y` and turns it into a vector of 2D points of the values
from `x` and `y`.

`P` is the plot Type (it is optional).
"""
#convert_arguments(::PointBased, x::RealVector, y::RealVector) = (Point2f.(x, y),)
convert_arguments(P::PointBased, x::ClosedInterval, y::RealVector) = convert_arguments(P, LinRange(extrema(x)..., length(y)), y)
convert_arguments(P::PointBased, x::RealVector, y::ClosedInterval) = convert_arguments(P, x, LinRange(extrema(y)..., length(x)))

"""
    convert_arguments(P, x)::(Vector)

Takes an input `Rect` `x` and decomposes it to points.

`P` is the plot Type (it is optional).
"""
function convert_arguments(P::PointBased, x::Rect2)
    # TODO fix the order of decompose
    return convert_arguments(P, decompose(Point2f, x)[[1, 2, 4, 3]])
end

function convert_arguments(P::PointBased, mesh::AbstractMesh)
    return convert_arguments(P, decompose(Point3f, mesh))
end

function convert_arguments(PB::PointBased, linesegments::FaceView{<:Line, P}) where {P<:AbstractPoint}
    # TODO FaceView should be natively supported by backends!
    return convert_arguments(PB, collect(reinterpret(P, linesegments)))
end

function convert_arguments(P::PointBased, rect::Rect3)
    return (decompose(Point3f, rect),)
end

function convert_arguments(P::Type{<: LineSegments}, rect::Rect3)
    f = decompose(LineFace{Int}, rect)
    p = connect(decompose(Point3f, rect), f)
    return convert_arguments(P, p)
end

function convert_arguments(::Type{<: Lines}, rect::Rect3)
    points = unique(decompose(Point3f, rect))
    push!(points, Point3f(NaN)) # use to seperate linesegments
    return (points[[1, 2, 3, 4, 1, 5, 6, 2, 9, 6, 8, 3, 9, 5, 7, 4, 9, 7, 8]],)
end
"""

    convert_arguments(PB, LineString)

Takes an input `LineString` and decomposes it to points.
"""
function convert_arguments(PB::PointBased, linestring::LineString)
    return convert_arguments(PB, decompose(Point, linestring))
end

"""
    convert_arguments(PB, Union{Array{<:LineString}, MultiLineString})

Takes an input `Array{LineString}` or a `MultiLineString` and decomposes it to points.
"""
function convert_arguments(PB::PointBased, linestring::Union{Array{<:LineString}, MultiLineString})
    arr = copy(convert_arguments(PB, linestring[1])[1])
    for ls in 2:length(linestring)
        push!(arr, Point2f(NaN))
        append!(arr, convert_arguments(PB, linestring[ls])[1])
    end
    return (arr,)
end

"""

    convert_arguments(PB, Polygon)

Takes an input `Polygon` and decomposes it to points.
"""
function convert_arguments(PB::PointBased, pol::Polygon)
    arr = copy(convert_arguments(PB, pol.exterior)[1])
    push!(arr, arr[1]) # close exterior
    if !isempty(pol.interiors)
        push!(arr, Point2f(NaN))
        for interior in pol.interiors
            inter = convert_arguments(PB, interior)[1]
            append!(arr, inter)
            # close interior + separate!
            push!(arr, inter[1], Point2f(NaN))
        end
    end
    return (arr,)
end

"""

    convert_arguments(PB, Union{Array{<:Polygon}, MultiPolygon})

Takes an input `Array{Polygon}` or a `MultiPolygon` and decomposes it to points.
"""
function convert_arguments(PB::PointBased, mp::Union{Array{<:Polygon}, MultiPolygon})
    arr = copy(convert_arguments(PB, mp[1])[1])
    for p in 2:length(mp)
        push!(arr, Point2f(NaN))
        append!(arr, convert_arguments(PB, mp[p])[1])
    end
    return (arr,)
end

################################################################################
#                                 SurfaceLike                                  #
################################################################################

function edges(v::AbstractVector)
    l = length(v)
    if l == 1
        return [v[1] - 0.5, v[1] + 0.5]
    else
        # Equivalent to
        # mids = 0.5 .* (v[1:end-1] .+ v[2:end])
        # borders = [2v[1] - mids[1]; mids; 2v[end] - mids[end]]
        borders = [0.5 * (v[max(1, i)] + v[min(end, i+1)]) for i in 0:length(v)]
        borders[1] = 2borders[1] - borders[2]
        borders[end] = 2borders[end] - borders[end-1]
        return borders
    end
end

function adjust_axes(::DiscreteSurface, x::AbstractVector{<:Number}, y::AbstractVector{<:Number}, z::AbstractMatrix)
    x̂, ŷ = map((x, y), size(z)) do v, sz
        return length(v) == sz ? edges(v) : v
    end
    return x̂, ŷ, z
end

adjust_axes(::SurfaceLike, x, y, z) = x, y, z

"""
    convert_arguments(SL::SurfaceLike, x::VecOrMat, y::VecOrMat, z::Matrix)

If `SL` is `Heatmap` and `x` and `y` are vectors, infer from length of `x` and `y`
whether they represent edges or centers of the heatmap bins.
If they are centers, convert to edges. Convert eltypes to `Float32` and return
outputs as a `Tuple`.
"""
function convert_arguments(SL::SurfaceLike, x::AbstractVecOrMat{<: Number}, y::AbstractVecOrMat{<: Number}, z::AbstractMatrix{<: Union{Number, Colorant}})
    return map(el32convert, adjust_axes(SL, x, y, z))
end
function convert_arguments(SL::SurfaceLike, x::AbstractVecOrMat{<: Number}, y::AbstractVecOrMat{<: Number}, z::AbstractMatrix{<:Number})
    return map(el32convert, adjust_axes(SL, x, y, z))
end

convert_arguments(::SurfaceLike, x::AbstractMatrix, y::AbstractMatrix) = (x, y, zeros(size(y)))

"""
    convert_arguments(P, x, y, z)::Tuple{ClosedInterval, ClosedInterval, Matrix}

Takes 2 ClosedIntervals's `x`, `y`, and an AbstractMatrix `z`, and converts the closed range to
linspaces with size(z, 1/2)
`P` is the plot Type (it is optional).
"""
function convert_arguments(P::SurfaceLike, x::ClosedInterval, y::ClosedInterval, z::AbstractMatrix)
    convert_arguments(P, to_linspace(x, size(z, 1)), to_linspace(y, size(z, 2)), z)
end

"""
    convert_arguments(P, Matrix)::Tuple{ClosedInterval, ClosedInterval, Matrix}

Takes an `AbstractMatrix`, converts the dimesions `n` and `m` into `ClosedInterval`,
and stores the `ClosedInterval` to `n` and `m`, plus the original matrix in a Tuple.

`P` is the plot Type (it is optional).
"""
function convert_arguments(sl::SurfaceLike, data::AbstractMatrix)
    n, m = Float32.(size(data))
    convert_arguments(sl, 0f0 .. n, 0f0 .. m, el32convert(data))
end

function convert_arguments(ds::DiscreteSurface, data::AbstractMatrix)
    n, m = Float32.(size(data))
    convert_arguments(ds, edges(1:n), edges(1:m), el32convert(data))
end

function convert_arguments(SL::SurfaceLike, x::AbstractVector{<:Number}, y::AbstractVector{<:Number}, z::AbstractVector{<:Number})
    if !(length(x) == length(y) == length(z))
        error("x, y and z need to have the same length. Lengths are $(length.((x, y, z)))")
    end

    xys = tuple.(x, y)
    if length(unique(xys)) != length(x)
        c = StatsBase.countmap(xys)
        cdup = filter(x -> x[2] > 1, c)
        error("Found duplicate x/y coordinates: $cdup")
    end

    x_centers = sort(unique(x))
    any(isnan, x_centers) && error("x must not have NaN values.")
    y_centers = sort(unique(y))
    any(isnan, y_centers) && error("x must not have NaN values.")
    zs = fill(NaN32, length(x_centers), length(y_centers))
    foreach(zip(x, y, z)) do (xi, yi, zi)
        i = searchsortedfirst(x_centers, xi)
        j = searchsortedfirst(y_centers, yi)
        @inbounds zs[i, j] = zi
    end
    convert_arguments(SL, x_centers, y_centers, zs)
end

"""
    convert_arguments(P, x, y, f)::(Vector, Vector, Matrix)

Takes vectors `x` and `y` and the function `f`, and applies `f` on the grid that `x` and `y` span.
This is equivalent to `f.(x, y')`.
`P` is the plot Type (it is optional).
"""
function convert_arguments(sl::SurfaceLike, x::AbstractVector{T1}, y::AbstractVector{T2}, f::Function) where {T1, T2}
    if !applicable(f, x[1], y[1])
        error("You need to pass a function with signature f(x::$T1, y::$T2). Found: $f")
    end
    T = typeof(f(x[1], y[1]))
    z = similar(x, T, (length(x), length(y)))
    z .= f.(x, y')
    return convert_arguments(sl, x, y, z)
end

################################################################################
#                                  VolumeLike                                  #
################################################################################

"""
    convert_arguments(P, Matrix)::Tuple{ClosedInterval, ClosedInterval, ClosedInterval, Matrix}

Takes an array of `{T, 3} where T`, converts the dimesions `n`, `m` and `k` into `ClosedInterval`,
and stores the `ClosedInterval` to `n`, `m` and `k`, plus the original array in a Tuple.

`P` is the plot Type (it is optional).
"""
function convert_arguments(::VolumeLike, data::AbstractArray{T, 3}) where T
    n, m, k = Float32.(size(data))
    return (0f0 .. n, 0f0 .. m, 0f0 .. k, el32convert(data))
end

function convert_arguments(::VolumeLike, x::RangeLike, y::RangeLike, z::RangeLike, data::AbstractArray{T, 3}) where T
    return (x, y, z, el32convert(data))
end
"""
    convert_arguments(P, x, y, z, i)::(Vector, Vector, Vector, Matrix)

Takes 3 `AbstractVector` `x`, `y`, and `z` and the `AbstractMatrix` `i`, and puts everything in a Tuple.

`P` is the plot Type (it is optional).
"""
function convert_arguments(::VolumeLike, x::AbstractVector, y::AbstractVector, z::AbstractVector, i::AbstractArray{T, 3}) where T
    (x, y, z, el32convert(i))
end

"""
    convert_arguments(P, x, y, z, f)::(Vector, Vector, Vector, Matrix)

Takes `AbstractVector` `x`, `y`, and `z` and the function `f`, evaluates `f` on the volume
spanned by `x`, `y` and `z`, and puts `x`, `y`, `z` and `f(x,y,z)` in a Tuple.

`P` is the plot Type (it is optional).
"""
function convert_arguments(::VolumeLike, x::AbstractVector, y::AbstractVector, z::AbstractVector, f::Function)
    if !applicable(f, x[1], y[1], z[1])
        error("You need to pass a function with signature f(x, y, z). Found: $f")
    end
    _x, _y, _z = ntuple(Val(3)) do i
        A = (x, y, z)[i]
        reshape(A, ntuple(j-> j != i ? 1 : length(A), Val(3)))
    end
    return (x, y, z, el32convert.(f.(_x, _y, _z)))
end

################################################################################
#                                <:Lines                                       #
################################################################################

function convert_arguments(::Type{<: Lines}, x::Rect2)
    # TODO fix the order of decompose
    points = decompose(Point2f, x)
    return (points[[1, 2, 4, 3, 1]],)
end

################################################################################
#                                <:LineSegments                                #
################################################################################

"""
Accepts a Vector of Pair of Points (e.g. `[Point(0, 0) => Point(1, 1), ...]`)
to encode e.g. linesegments or directions.
"""
function convert_arguments(::Type{<: LineSegments}, positions::AbstractVector{E}) where E <: Union{Pair{A, A}, Tuple{A, A}} where A <: VecTypes{N, T} where {N, T}
    (elconvert(Point{N, Float32}, reinterpret(Point{N, T}, positions)),)
end

function convert_arguments(::Type{<: LineSegments}, x::Rect2)
    # TODO fix the order of decompose
    points = decompose(Point2f, x)
    return (points[[1, 2, 2, 4, 4, 3, 3, 1]],)
end

################################################################################
#                                    <:Text                                    #
################################################################################

"""
    convert_arguments(x)::(String)

Takes an input `AbstractString` `x` and converts it to a string.
"""
# convert_arguments(::Type{<: Text}, x::AbstractString) = (String(x),)

################################################################################
#                                    <:Mesh                                    #
################################################################################

"""
    convert_arguments(Mesh, x, y, z)::GLNormalMesh

Takes real vectors x, y, z and constructs a mesh out of those, under the assumption that
every 3 points form a triangle.
"""
function convert_arguments(
        T::Type{<:Mesh},
        x::RealVector, y::RealVector, z::RealVector
    )
    convert_arguments(T, Point3f.(x, y, z))
end
"""
    convert_arguments(Mesh, xyz::AbstractVector)::GLNormalMesh

Takes an input mesh and a vector `xyz` representing the vertices of the mesh, and
creates indices under the assumption, that each triplet in `xyz` forms a triangle.
"""
function convert_arguments(
        MT::Type{<:Mesh},
        xyz::AbstractVector
    )
    faces = connect(UInt32.(0:length(xyz)-1), GLTriangleFace)
    # TODO support faceview natively
    return convert_arguments(MT, xyz, collect(faces))
end

function convert_arguments(::Type{<:Mesh}, mesh::GeometryBasics.Mesh{N}) where {N}
    # Make sure we have normals!
    if !hasproperty(mesh, :normals)
        n = normals(mesh)
        # Normals can be nothing, when it's impossible to calculate the normals (e.g. 2d mesh)
        if n !== nothing
            mesh = GeometryBasics.pointmeta(mesh, decompose(Vec3f, n))
        end
    end
    return (GeometryBasics.mesh(mesh, pointtype=Point{N, Float32}, facetype=GLTriangleFace),)
end

function convert_arguments(
        MT::Type{<:Mesh},
        meshes::AbstractVector{<: Union{AbstractMesh, AbstractPolygon}}
    )
    return (meshes,)
end

function convert_arguments(
        MT::Type{<:Mesh},
        xyz::Union{AbstractPolygon, AbstractVector{<: AbstractPoint{2}}}
    )
    return convert_arguments(MT, triangle_mesh(xyz))
end

function convert_arguments(MT::Type{<:Mesh}, geom::GeometryPrimitive)
    # we convert to UV mesh as default, because otherwise the uv informations get lost
    # - we can still drop them, but we can't add them later on
    return (GeometryBasics.uv_normal_mesh(geom),)
end

"""
    convert_arguments(Mesh, x, y, z, indices)::GLNormalMesh

Takes real vectors x, y, z and constructs a triangle mesh out of those, using the
faces in `indices`, which can be integers (every 3 -> one triangle), or GeometryBasics.NgonFace{N, <: Integer}.
"""
function convert_arguments(
        T::Type{<: Mesh},
        x::RealVector, y::RealVector, z::RealVector,
        indices::AbstractVector
    )
    return convert_arguments(T, Point3f.(x, y, z), indices)
end

"""
    convert_arguments(Mesh, vertices, indices)::GLNormalMesh

Takes `vertices` and `indices`, and creates a triangle mesh out of those.
See [`to_vertices`](@ref) and [`to_triangles`](@ref) for more information about
accepted types.
"""
function convert_arguments(
        ::Type{<:Mesh},
        vertices::AbstractArray,
        indices::AbstractArray
    )
    m = normal_mesh(to_vertices(vertices), to_triangles(indices))
    (m,)
end

################################################################################
#                             Function Conversions                             #
################################################################################

function convert_arguments(P::PlotFunc, r::AbstractVector, f::Function)
    ptype = plottype(P, Lines)
    to_plotspec(ptype, convert_arguments(ptype, r, f.(r)))
end

function convert_arguments(P::PlotFunc, i::AbstractInterval, f::Function)
    x, y = PlotUtils.adapted_grid(f, endpoints(i))
    ptype = plottype(P, Lines)
    to_plotspec(ptype, convert_arguments(ptype, x, y))
end

# The following `tryrange` code was copied from Plots.jl
# https://github.com/JuliaPlots/Plots.jl/blob/15dc61feb57cba1df524ce5d69f68c2c4ea5b942/src/series.jl#L399-L416

# try some intervals over which the function may be defined
function tryrange(F::AbstractArray, vec)
    rets = [tryrange(f, vec) for f in F] # get the preferred for each
    maxind = maximum(indexin(rets, vec)) # get the last attempt that succeeded (most likely to fit all)
    rets .= [tryrange(f, vec[maxind:maxind]) for f in F] # ensure that all functions compute there
    rets[1]
end

function tryrange(F, vec)
    for v in vec
        try
            tmp = F(v)
            return v
        catch
        end
    end
    error("$F is not a Function, or is not defined at any of the values $vec")
end

# OffsetArrays conversions
function convert_arguments(sl::SurfaceLike, wm::OffsetArray)
  x1, y1 = wm.offsets .+ 1
  nx, ny = size(wm)
  x = range(x1, length = nx)
  y = range(y1, length = ny)
  v = parent(wm)
  return convert_arguments(sl, x, y, v)
end

################################################################################
#                               Helper Functions                               #
################################################################################

to_linspace(interval, N) = range(minimum(interval), stop = maximum(interval), length = N)

"""
Converts the elemen array type to `T1` without making a copy if the element type matches
"""
elconvert(::Type{T1}, x::AbstractArray{T2, N}) where {T1, T2, N} = convert(AbstractArray{T1, N}, x)
float32type(x::Type) = Float32
float32type(::Type{<: RGB}) = RGB{Float32}
float32type(::Type{<: RGBA}) = RGBA{Float32}
float32type(::Type{<: Colorant}) = RGBA{Float32}
float32type(x::AbstractArray{T}) where T = float32type(T)
float32type(x::T) where T = float32type(T)
el32convert(x::AbstractArray) = elconvert(float32type(x), x)
el32convert(x::AbstractArray{Float32}) = x
el32convert(x::Observable) = lift(el32convert, x)
el32convert(x) = convert(float32type(x), x)

function el32convert(x::AbstractArray{T, N}) where {T<:Union{Missing, <: Number}, N}
    return map(x) do elem
        return (ismissing(elem) ? NaN32 : convert(Float32, elem))::Float32
    end::Array{Float32, N}
end
"""
    to_triangles(indices)

Convert a representation of triangle point indices `indices` to its canonical representation as a `Vector{Makie.GLTriangleFace}`. `indices` can be any of the following:

- An `AbstractVector{Int}`, containing groups of 3 1-based indices,
- An `AbstractVector{UIn32}`, containing groups of 3 0-based indices,
- An `AbstractVector` of `TriangleFace` objects,
- An `AbstractMatrix` of `Integer`s, where each row is a triangle.
"""
function to_triangles(x::AbstractVector{Int})
    idx0 = UInt32.(x .- 1)
    return to_triangles(idx0)
end

function to_triangles(idx0::AbstractVector{UInt32})
    reinterpret(GLTriangleFace, idx0)
end

function to_triangles(faces::AbstractVector{TriangleFace{T}}) where T
    elconvert(GLTriangleFace, faces)
end

function to_triangles(faces::AbstractMatrix{T}) where T <: Integer
    let N = Val(size(faces, 2)), lfaces = faces
        broadcast(1:size(faces, 1), N) do fidx, n
            to_ndim(GLTriangleFace, ntuple(i-> lfaces[fidx, i], n), 0.0)
        end
    end
end

"""
    to_vertices(v)

Converts a representation of vertices `v` to its canonical representation as a
`Vector{Point3f}`. `v` can be:

- An `AbstractVector` of 3-element `Tuple`s or `StaticVector`s,

- An `AbstractVector` of `Tuple`s or `StaticVector`s, in which case exta dimensions will
  be either truncated or padded with zeros as required,

- An `AbstractMatrix`"
  - if `v` has 2 or 3 rows, it will treat each column as a vertex,
  - otherwise if `v` has 2 or 3 columns, it will treat each row as a vertex.
"""
function to_vertices(verts::AbstractVector{<: VecTypes{3, T}}) where T
    vert3f0 = T != Float32 ? Point3f.(verts) : verts
    return reinterpret(Point3f, vert3f0)
end

function to_vertices(verts::AbstractVector{<: VecTypes})
    to_vertices(to_ndim.(Point3f, verts, 0.0))
end

function to_vertices(verts::AbstractMatrix{<: Number})
    if size(verts, 1) in (2, 3)
        to_vertices(verts, Val(1))
    elseif size(verts, 2) in (2, 3)
        to_vertices(verts, Val(2))
    else
        error("You are using a matrix for vertices which uses neither dimension to encode the dimension of the space. Please have either size(verts, 1/2) in the range of 2-3. Found: $(size(verts))")
    end
end

function to_vertices(verts::AbstractMatrix{T}, ::Val{1}) where T <: Number
    N = size(verts, 1)
    if T == Float32 && N == 3
        reinterpret(Point{N, T}, elconvert(T, vec(verts)))
    else
        let N = Val(N), lverts = verts
            broadcast(1:size(verts, 2), N) do vidx, n
                to_ndim(Point3f, ntuple(i-> lverts[i, vidx], n), 0.0)
            end
        end
    end
end

function to_vertices(verts::AbstractMatrix{T}, ::Val{2}) where T <: Number
    let N = Val(size(verts, 2)), lverts = verts
        broadcast(1:size(verts, 1), N) do vidx, n
            to_ndim(Point3f, ntuple(i-> lverts[vidx, i], n), 0.0)
        end
    end
end

################################################################################
#                            Attribute conversions                             #
################################################################################

"""
    to_color(color)
Converts a `color` symbol (e.g. `:blue`) to a color RGBA.
"""
to_color(color) = convert_attribute(color, key"color"())

"""
    to_colormap(cm[, N = 20])

Converts a colormap `cm` symbol (e.g. `:Spectral`) to a colormap RGB array, where `N` specifies the number of color points.
"""
to_colormap(colormap) = convert_attribute(colormap, key"colormap"())
to_rotation(rotation) = convert_attribute(rotation, key"rotation"())
to_font(font) = convert_attribute(font, key"font"())
to_align(align) = convert_attribute(align, key"align"())
to_textsize(textsize) = convert_attribute(textsize, key"textsize"())

convert_attribute(x, key::Key, ::Key) = convert_attribute(x, key)
convert_attribute(s::SceneLike, x, key::Key, ::Key) = convert_attribute(s, x, key)
convert_attribute(s::SceneLike, x, key::Key) = convert_attribute(x, key)
convert_attribute(x, key::Key) = x

convert_attribute(p, ::key"highclip") = to_color(p)
convert_attribute(p::Nothing, ::key"highclip") = p
convert_attribute(p, ::key"lowclip") = to_color(p)
convert_attribute(p::Nothing, ::key"lowclip") = p
convert_attribute(p, ::key"nan_color") = to_color(p)

struct Palette{N}
   colors::SArray{Tuple{N},RGBA{Float32},1,N}
   i::Ref{UInt8}
   Palette(colors) = new{length(colors)}(SVector{length(colors)}(to_color.(colors)), zero(UInt8))
end
Palette(name::Union{String, Symbol}, n = 8) = Palette(to_colormap(name, n))

function convert_attribute(p::Palette{N}, ::key"color") where {N}
    p.i[] = p.i[] == N ? one(UInt8) : p.i[] + one(UInt8)
    p.colors[p.i[]]
end

convert_attribute(c::Colorant, ::key"color") = convert(RGBA{Float32}, c)
convert_attribute(c::Symbol, k::key"color") = convert_attribute(string(c), k)
function convert_attribute(c::String, ::key"color")
    return parse(RGBA{Float32}, c)
end

# Do we really need all colors to be RGBAf?!
convert_attribute(c::AbstractArray{<: Colorant}, k::key"color") = el32convert(c)
convert_attribute(c::AbstractArray, k::key"color") = to_color.(c)

convert_attribute(c::AbstractArray, ::key"color", ::key"heatmap") = el32convert(c)

convert_attribute(c::Tuple, k::key"color") = convert_attribute.(c, k)
convert_attribute(p::AbstractPattern, k::key"color") = p

function convert_attribute(c::Tuple{T, F}, k::key"color") where {T, F <: Number}
    RGBAf(Colors.color(to_color(c[1])), c[2])
end

convert_attribute(b::Billboard{Float32}, ::key"rotations") = to_rotation(b.rotation)
convert_attribute(b::Billboard{Vector{Float32}}, ::key"rotations") = to_rotation.(b.rotation)
convert_attribute(r::AbstractArray, ::key"rotations") = to_rotation.(r)
convert_attribute(r::StaticVector, ::key"rotations") = to_rotation(r)
convert_attribute(r, ::key"rotations") = to_rotation(r)

convert_attribute(c, ::key"markersize", ::key"scatter") = to_2d_scale(c)
convert_attribute(c, k1::key"markersize", k2::key"meshscatter") = to_3d_scale(c)

convert_attribute(x, ::key"uv_offset_width") = Vec4f(x)
convert_attribute(x::AbstractVector{Vec4f}, ::key"uv_offset_width") = x

to_2d_scale(x::Number) = Vec2f(x)
to_2d_scale(x::VecTypes) = to_ndim(Vec2f, x, 1)
to_2d_scale(x::Tuple{<:Number, <:Number}) = to_ndim(Vec2f, x, 1)
to_2d_scale(x::AbstractVector) = to_2d_scale.(x)

to_3d_scale(x::Number) = Vec3f(x)
to_3d_scale(x::VecTypes) = to_ndim(Vec3f, x, 1)
to_3d_scale(x::AbstractVector) = to_3d_scale.(x)

convert_attribute(c::Number, ::key"glowwidth") = Float32(c)
convert_attribute(c, ::key"glowcolor") = to_color(c)
convert_attribute(c, ::key"strokecolor") = to_color(c)
convert_attribute(c::Number, ::key"strokewidth") = Float32(c)

convert_attribute(x::Nothing, ::key"linestyle") = x

"""
    `AbstractVector{<:AbstractFloat}` for denoting sequences of fill/nofill. e.g.

[0.5, 0.8, 1.2] will result in 0.5 filled, 0.3 unfilled, 0.4 filled. 1.0 unit is one linewidth!
"""
convert_attribute(A::AbstractVector, ::key"linestyle") = A

"""
    A `Symbol` equal to `:dash`, `:dot`, `:dashdot`, `:dashdotdot`
"""
convert_attribute(ls::Union{Symbol,AbstractString}, ::key"linestyle") = line_pattern(ls, :normal)

function convert_attribute(ls::Tuple{<:Union{Symbol,AbstractString},<:Any}, ::key"linestyle")
    line_pattern(ls[1], ls[2])
end

function line_pattern(linestyle, gaps)
    pattern = line_diff_pattern(linestyle, gaps)
    isnothing(pattern) ? pattern : float.([0.0; cumsum(pattern)])
end

"The linestyle patterns are inspired by the LaTeX package tikZ as seen here https://tex.stackexchange.com/questions/45275/tikz-get-values-for-predefined-dash-patterns."

function line_diff_pattern(ls::Symbol, gaps = :normal)
    if ls == :solid
        nothing
    elseif ls == :dash
        line_diff_pattern("-", gaps)
    elseif ls == :dot
        line_diff_pattern(".", gaps)
    elseif ls == :dashdot
        line_diff_pattern("-.", gaps)
    elseif ls == :dashdotdot
        line_diff_pattern("-..", gaps)
    else
        error(
            """
            Unkown line style: $ls. Available linestyles are:
            :solid, :dash, :dot, :dashdot, :dashdotdot
            or a sequence of numbers enumerating the next transparent/opaque region.
            This sequence of numbers must be cumulative; 1 unit corresponds to 1 line width.
            """
        )
    end
end

function line_diff_pattern(ls_str::AbstractString, gaps = :normal)
    dot = 1
    dash = 3
    check_line_pattern(ls_str)

    dot_gap, dash_gap = convert_gaps(gaps)

    pattern = Float64[]
    for i in 1:length(ls_str)
        curr_char = ls_str[i]
        next_char = i == lastindex(ls_str) ? ls_str[firstindex(ls_str)] : ls_str[i+1]
        # push dash or dot
        if curr_char == '-'
            push!(pattern, dash)
        else
            push!(pattern, dot)
        end
        # push the gap (use dot_gap only between two dots)
        if (curr_char == '.') && (next_char == '.')
            push!(pattern, dot_gap)
        else
            push!(pattern, dash_gap)
        end
    end
    pattern
end

"Checks if the linestyle format provided as a string contains only dashes and dots"
function check_line_pattern(ls_str)
    isnothing(match(r"^[.-]+$", ls_str)) &&
        throw(ArgumentError("If you provide a string as linestyle, it must only consist of dashes (-) and dots (.)"))

    nothing
end

function convert_gaps(gaps)
  error_msg = "You provided the gaps modifier $gaps when specifying the linestyle. The modifier must be `∈ ([:normal, :dense, :loose])`, a real number or a collection of two real numbers."
  if gaps isa Symbol
      gaps in [:normal, :dense, :loose] || throw(ArgumentError(error_msg))
      dot_gaps  = (normal = 2, dense = 1, loose = 4)
      dash_gaps = (normal = 3, dense = 2, loose = 6)

      dot_gap  = getproperty(dot_gaps, gaps)
      dash_gap = getproperty(dash_gaps, gaps)
  elseif gaps isa Real
      dot_gap = gaps
      dash_gap = gaps
  elseif length(gaps) == 2 && eltype(gaps) <: Real
      dot_gap, dash_gap = gaps
  else
      throw(ArgumentError(error_msg))
  end
  (dot_gap = dot_gap, dash_gap = dash_gap)
end

function convert_attribute(f::Symbol, ::key"frames")
    f == :box && return ((true, true), (true, true))
    f == :semi && return ((true, false), (true, false))
    f == :none && return ((false, false), (false, false))
    throw(MethodError("$(string(f)) is not a valid framestyle. Options are `:box`, `:semi` and `:none`"))
end
convert_attribute(f::Tuple{Tuple{Bool,Bool},Tuple{Bool,Bool}}, ::key"frames") = f

convert_attribute(c::Tuple{<: Number, <: Number}, ::key"position") = Point2f(c[1], c[2])
convert_attribute(c::Tuple{<: Number, <: Number, <: Number}, ::key"position") = Point3f(c)
convert_attribute(c::VecTypes{N}, ::key"position") where N = Point{N, Float32}(c)

"""
    Text align, e.g.:
"""
convert_attribute(x::Tuple{Symbol, Symbol}, ::key"align") = Vec2f(alignment2num.(x))
convert_attribute(x::Vec2f, ::key"align") = x

const _font_cache = Dict{String, NativeFont}()

"""
    font conversion

a string naming a font, e.g. helvetica
"""
function convert_attribute(x::Union{Symbol, String}, k::key"font")
    str = string(x)
    get!(_font_cache, str) do
        str == "default" && return to_font("Dejavu Sans")

        # check if the string points to a font file and load that
        if isfile(str)
            font = FreeTypeAbstraction.try_load(str)
            if isnothing(font)
                error("Could not load font file $str")
            else
                return font
            end
        end

        fontpath = assetpath("fonts")
        font = FreeTypeAbstraction.findfont(str; additional_fonts=fontpath)
        if font === nothing
            @warn("Could not find font $str, using Dejavu Sans")
            if "dejavu sans" == lowercase(str)
                # since we fall back to dejavu sans, we need to check for recursion
                error("Recursion encountered; DejaVu Sans cannot be located in the font path $fontpath")
            end
            return to_font("dejavu sans")
        end
        return font
    end
end
convert_attribute(x::Vector{String}, k::key"font") = convert_attribute.(x, k)
convert_attribute(x::NativeFont, k::key"font") = x

"""
    rotation accepts:
    to_rotation(b, quaternion)
    to_rotation(b, tuple_float)
    to_rotation(b, vec4)
"""
convert_attribute(s::Quaternionf, ::key"rotation") = s
convert_attribute(s::Quaternion, ::key"rotation") = Quaternionf(s.data...)
function convert_attribute(s::VecTypes{N}, ::key"rotation") where N
    if N == 4
        Quaternionf(s...)
    elseif N == 3
        rotation_between(Vec3f(0, 0, 1), to_ndim(Vec3f, s, 0.0))
    elseif N == 2
        rotation_between(Vec3f(0, 1, 0), to_ndim(Vec3f, s, 0.0))
    else
        error("The $N dimensional vector $s can't be converted to a rotation.")
    end
end

function convert_attribute(s::Tuple{VecTypes, AbstractFloat}, ::key"rotation")
    qrotation(to_ndim(Vec3f, s[1], 0.0), s[2])
end
convert_attribute(angle::AbstractFloat, ::key"rotation") = qrotation(Vec3f(0, 0, 1), Float32(angle))
convert_attribute(r::AbstractVector, k::key"rotation") = to_rotation.(r)
convert_attribute(r::AbstractVector{<: Quaternionf}, k::key"rotation") = r

convert_attribute(x, k::key"colorrange") = x==nothing ? nothing : Vec2f(x)

convert_attribute(x, k::key"textsize") = Float32(x)
convert_attribute(x::AbstractVector, k::key"textsize") = convert_attribute.(x, k)
convert_attribute(x::AbstractVector{T}, k::key"textsize") where T <: Number = el32convert(x)
convert_attribute(x::AbstractVector{T}, k::key"textsize") where T <: VecTypes = elconvert(Vec2f, x)
convert_attribute(x, k::key"linewidth") = Float32(x)
convert_attribute(x::AbstractVector, k::key"linewidth") = el32convert(x)

# ColorBrewer colormaps that support only 8 colors require special handling on the backend, so we show them here.
const colorbrewer_8color_names = String.([
    :Accent,
    :Dark2,
    :Pastel2,
    :Set2
])

const plotutils_names = String.(union(
    keys(PlotUtils.ColorSchemes.colorschemes),
    keys(PlotUtils.COLORSCHEME_ALIASES),
    keys(PlotUtils.MISC_COLORSCHEMES)
))

const all_gradient_names = Set(vcat(plotutils_names, colorbrewer_8color_names))

"""
    available_gradients()

Prints all available gradient names.
"""
function available_gradients()
    println("Gradient Symbol/Strings:")
    for name in sort(collect(all_gradient_names))
        println("    ", name)
    end
end

"""
Reverses the attribute T upon conversion
"""
struct Reverse{T}
    data::T
end

function convert_attribute(r::Reverse, ::key"colormap", n::Integer=20)
    reverse(to_colormap(r.data, n))
end

function convert_attribute(cs::ColorScheme, ::key"colormap", n::Integer=20)
    return to_colormap(cs.colors, n)
end

"""
    to_colormap(b, x)

An `AbstractVector{T}` with any object that [`to_color`](@ref) accepts.
"""
convert_attribute(cm::AbstractVector, ::key"colormap", n::Int=length(cm)) = to_colormap(to_color.(cm), n)

function convert_attribute(cm::AbstractVector{<: Colorant}, ::key"colormap", n::Int=length(cm))
    colormap = length(cm) == n ? cm : resample(cm, n)
    return el32convert(colormap)
end

"""
Tuple(A, B) or Pair{A, B} with any object that [`to_color`](@ref) accepts
"""
function convert_attribute(cs::Union{Tuple, Pair}, ::key"colormap", n::Int=2)
    return to_colormap([to_color.(cs)...], n)
end

function convert_attribute(cs::Tuple{<: Union{Symbol, AbstractString}, Real}, ::key"colormap", n::Int=30)
    return RGBAf.(to_colormap(cs[1]), cs[2]) # We need to rework this to conform to the backend interface.
end

function convert_attribute(cs::NamedTuple{(:colormap, :alpha, :n), Tuple{Union{Symbol, AbstractString}, Real, Int}}, ::key"colormap")
    return RGBAf.(to_colormap(cs.colormap, cs.n), cs.alpha)
end

to_colormap(x, n::Integer) = convert_attribute(x, key"colormap"(), n)

"""
A Symbol/String naming the gradient. For more on what names are available please see: `available_gradients()`.
For now, we support gradients from `PlotUtils` natively.
"""
function convert_attribute(cs::Union{String, Symbol}, ::key"colormap", n::Integer=40)
    cs_string = string(cs)
    if cs_string in all_gradient_names
        if cs_string in colorbrewer_8color_names # special handling for 8 color only
            return to_colormap(ColorBrewer.palette(cs_string, 8), n)
        else                                    # cs_string must be in plotutils_names
            return to_colormap(PlotUtils.get_colorscheme(Symbol(cs_string)).colors, n)
        end
    else
        error(
            """
            There is no color gradient named $cs.
            See `Makie.available_gradients()` for the list of available gradients,
            or look at http://makie.juliaplots.org/dev/generated/colors#Colormap-reference.
            """
        )
    end
end

function Makie.convert_attribute(cg::PlotUtils.ContinuousColorGradient, ::key"colormap", n::Integer=length(cg.values))
    # PlotUtils does not always give [0, 1] range, so we adapt to what it has
    return getindex.(Ref(cg), LinRange(first(cg.values), last(cg.values), n))
end

function Makie.convert_attribute(cg::PlotUtils.CategoricalColorGradient, ::key"colormap", n::Integer = length(cg.colors) * 20)
    # PlotUtils does not always give [0, 1] range, so we adapt to what it has
    return vcat(fill.(cg.colors.colors, Ref(n ÷ length(cg.colors)))...)
end

"""
    to_volume_algorithm(b, x)

Enum values: `IsoValue` `Absorption` `MaximumIntensityProjection` `AbsorptionRGBA` `AdditiveRGBA` `IndexedAbsorptionRGBA`
"""
function convert_attribute(value, ::key"algorithm")
    if isa(value, RaymarchAlgorithm)
        return Int32(value)
    elseif isa(value, Int32) && value in 0:5
        return value
    elseif value == 7
        return value # makie internal contour implementation
    else
        error("$value is not a valid volume algorithm. Please have a look at the docstring of `to_volume_algorithm` (in the REPL, `?to_volume_algorithm`).")
    end
end

"""
Symbol/String: iso, absorption, mip, absorptionrgba, indexedabsorption
"""
function convert_attribute(value::Union{Symbol, String}, k::key"algorithm")
    vals = Dict(
        :iso => IsoValue,
        :absorption => Absorption,
        :mip => MaximumIntensityProjection,
        :absorptionrgba => AbsorptionRGBA,
        :indexedabsorption => IndexedAbsorptionRGBA,
        :additive => AdditiveRGBA,
    )
    convert_attribute(get(vals, Symbol(value)) do
        error("$value is not a valid volume algorithm. It must be one of $(keys(vals))")
    end, k)
end

const _marker_map = Dict(
    :rect => '■',
    :star5 => '★',
    :diamond => '◆',
    :hexagon => '⬢',
    :cross => '✚',
    :xcross => '❌',
    :utriangle => '▲',
    :dtriangle => '▼',
    :ltriangle => '◀',
    :rtriangle => '▶',
    :pentagon => '⬟',
    :octagon => '⯄',
    :star4 => '✦',
    :star6 => '🟋',
    :star8 => '✷',
    :vline => '┃',
    :hline => '━',
    :+ => '+',
    :x => 'x',
    :circle => '●'
)

"""
    available_marker_symbols()

Displays all available marker symbols.
"""
function available_marker_symbols()
    println("Marker Symbols:")
    for (k, v) in _marker_map
        println("    ", k, " => ", v)
    end
end

"""
    FastPixel()

Use

```julia
scatter(..., marker=FastPixel())
```

For significant faster plotting times for large amount of points.
Note, that this will draw markers always as 1 pixel.
"""
struct FastPixel end

to_spritemarker(x::FastPixel) = x
to_spritemarker(x::Circle) = x
to_spritemarker(::Type{<: Circle}) = Circle(Point2f(0), 1f0)
to_spritemarker(::Type{<: Rect}) = Rect(Vec2f(0), Vec2f(1))
to_spritemarker(x::Rect) = x

"""
    to_spritemarker(b, marker::Char)

Any `Char`, including unicode
"""
to_spritemarker(marker::Char) = marker

"""
Matrix of AbstractFloat will be interpreted as a distancefield (negative numbers outside shape, positive inside)
"""
to_spritemarker(marker::Matrix{<: AbstractFloat}) = el32convert(marker)

"""
Any AbstractMatrix{<: Colorant} or other image type
"""
to_spritemarker(marker::AbstractMatrix{<: Colorant}) = marker

"""
A `Symbol` - Available options can be printed with `available_marker_symbols()`
"""
function to_spritemarker(marker::Symbol)
    if haskey(_marker_map, marker)
        return to_spritemarker(_marker_map[marker])
    else
        @warn("Unsupported marker: $marker, using ● instead")
        return '●'
    end
end

to_spritemarker(marker::String) = marker
to_spritemarker(marker::AbstractVector{Char}) = String(marker)

"""
Vector of anything that is accepted as a single marker will give each point it's own marker.
Note that it needs to be a uniform vector with the same element type!
"""
function to_spritemarker(marker::AbstractVector)
    marker = to_spritemarker.(marker)
    if isa(marker, AbstractVector{Char})
        String(marker)
    else
        marker
    end
end

convert_attribute(value, ::key"marker", ::key"scatter") = to_spritemarker(value)
convert_attribute(value, ::key"isovalue", ::key"volume") = Float32(value)
convert_attribute(value, ::key"isorange", ::key"volume") = Float32(value)

function convert_attribute(value::Symbol, ::key"marker", ::key"meshscatter")
    if value == :Sphere
        return normal_mesh(Sphere(Point3f(0), 1f0))
    else
        error("Unsupported marker: $(value)")
    end
end

function convert_attribute(value::AbstractGeometry, ::key"marker", ::key"meshscatter")
    return normal_mesh(value)
end


text:
#= An assembly of functions enabling the generation of moment
   equations for any chemical reaction network with any type of
   (infinitely differentiable) propensities up to arbitrary order =#

function fact(i)

    #= Calculate a multi-variate factorial of moment vector i,
       i.e., if i=(a, b, c), then fact(i) = a!b!c! =#

    fact = 1
    for j in i
        fact *= factorial(j)
    end
    return fact

end

"""
    generate_central_moment_eqs(rn::Union{ReactionSystem, ReactionSystemMod}, m_order::Int, q_order=nothing;
                                combinatoric_ratelaw=true, smap=speciesmap(rn))

Given a [`ReactionSystem`](https://catalyst.sciml.ai/stable/api/catalyst_api/#ModelingToolkit.ReactionSystem)
or [`ReactionSystemMod`](@ref), return the [`CentralMomentEquations`](@ref) of the system generated up to `m_order`.

Notes:
- if `q_order` is not specified by the user, it is assumed that the reaction network
  contains *only* polynomial propensity functions and hence `q_order` is determined
  automatically as in [`generate_raw_moment_eqs`](@ref). However, `q_order` must
  be specified if non-polynomial propensities are included. Note that the expansion
  order ``q`` denotes the highest order of central moments which will be included
  in the ODEs [(due to the Taylor expansion of propensity functions)](@ref central_moment_eqs).
- `combinatoric_ratelaw=true` uses binomials in calculating the propensity functions
  of a `ReactionSystem`, see the notes for [`ModelingToolkit.jumpratelaw`]
  (https://mtk.sciml.ai/stable/systems/ReactionSystem/#ModelingToolkit.jumpratelaw).
  *Note* that this field is irrelevant using `ReactionSystemMod` as then the
  propensities are defined directly by the user.
- `smap` sets the variable ordering in the moment equations (which index corresponds to which species
  in the reaction network). By default, this is consistent with the internal system ordering
  accessible with [`speciesmap`](@ref).
"""
function generate_central_moment_eqs(rn::Union{ReactionSystem, ReactionSystemMod},
                                     m_order::Int, q_order::Int=0;
                                     combinatoric_ratelaw=true, smap=speciesmap(rn))

    N = numspecies(rn) # no. of molecular species in the network
    R = numreactions(rn) # no. of reactions in the network

    # propensity functions of all reactions in the network
    a = propensities(rn; combinatoric_ratelaw)

    # quite messy way to check whether all propensity functions are polynomials
    # and extract the moment expansion order automatically (if not set by the user)
    if iszero(q_order)
        try
            _, _, poly_order = polynomial_propensities(a, rn)
            q_order = poly_order + m_order - 1
        catch e
            error("non-polynomial rates (OR A BUG): please specify q_order.\n" * string(e))
        end
    end

    # net stoichiometric matrix
    S = netstoichmat(rn; smap)

    # iterator over all moments from lowest to highest moment order
    iter_all = construct_iter_all(N, q_order)
    # iterator over central moments up to order m
    iter_m = filter(x -> 1 < sum(x) <= m_order, iter_all)
    # iterator over central moments of order greater than m up to q_order
    iter_q = filter(x -> m_order < sum(x) <= q_order, iter_all)
    # iterator over the first order central moments
    iter_1 = filter(x -> sum(x) == 1, iter_all)

    #= Define the first order raw moments μ and
       central moments Mᵢ as symbolic variables
       using the functionality of ModelingToolkit.jl =#

    @parameters t # need to define time as a parameter

    μ = define_μ(N, 1, iter_1)
    M = define_M(N, q_order)

    #= Obtain all derivatives of the propensity functions with respect
    to all molecular species up to order defined by q_order.
    For example, Da[2][(1, 2)] is equivalent to (∂³/∂n₁∂n₂²) a₂(μ) =#

    Da = []
    # define the number of molecules of each species as a variable (required for the differentiation)
    n = species(rn)
    dict_n_to_μ = Dict(zip(n, values(μ)))

    # save derivatives in a dictionary
    derivs = Dict()
    for r in 1:R
        for k in iter_all
            derivs[k] = a[r]
            for i in 1:N
                for d in 1:k[i]
                    derivs[k] = expand_derivatives(Differential(n[i])(derivs[k]), true) # simplify = true/false
                end
            end
            derivs[k] = substitute(derivs[k], dict_n_to_μ)
        end
        push!(Da, copy(derivs))
    end

    # generate the equations for the first raw moments (means)
    du = []
    for r in 1:R
        suma = 0
        for j in iter_all
            suma = suma + Da[r][j]*M[j]*1//fact(j)
            #suma = simplify(suma)
            # here // gives fractions (otherwise it's a float number)
            # only issue is that it does not latexify properly
        end

        for i in 1:N
            if r == 1
                push!(du, S[i,r]*suma)
            else
                du[i] = S[i, r]*suma + du[i]
            end
            du[i] = expand(du[i])
        end
    end

    # generate the equations for the central moments up to order m

    dM = Dict()
    for i in iter_m
        dM[i] = 0
        for r in 1:R
            iter_j = filter(x -> all(x .<= i) && sum(x) <= sum(i)-1, iter_all)
            for j in iter_j
                factor_j = 1.0
                for k in 1:N
                    factor_j *= expected_coeff(S[k, r], i[k]-j[k]) * binomial(i[k], j[k])
                end
                # m+1 -> MA_order
                iter_k = filter(x -> sum(x) <= q_order-sum(j), iter_all)
                suma = 0.0
                for k in iter_k
                    suma += Da[r][k]*M[j.+k]*1//fact(k)
                end
                dM[i] += factor_j*suma
            end
        end
        for j in 1:N
            if i[j] > 0
                dM[i] -= i[j]*du[j]*M[i.-iter_1[j]]
            end
        end
        dM[i] = expand(dM[i])
    end

    iv = get_iv(rn)
    D = Differential(iv)
    eqs = Equation[]
    for i in 1:N
        push!(eqs, D(μ[iter_1[i]]) ~ du[i])
    end
    for i in iter_m
        push!(eqs, D(M[i]) ~ dM[i])
    end

    vars = extract_variables(eqs, N, q_order)
    odename = Symbol(nameof(rn), "_central_moment_eqs_m", m_order, "_q", q_order)
    odes = ODESystem(eqs, iv, vars, get_ps(rn); name=odename)

    CentralMomentEquations(
        odes,
        μ,
        M,
        N,
        m_order,
        q_order,
        iter_all,
        iter_m,
        iter_q,
        iter_1
    )
end


text:
@time using MianalyzerBackend: TimeSeriesController, CorrelationsController, BasicStatisticsController, GaussianMixturesController, KMeansController, GMeansController, HierachicalClusteringsController

function handler(f, jsonpayload)
  return f(jsonpayload)
end

function routes(path, payload)
  if path == "/time_series"
    return handler(TimeSeriesController.call, payload)
  elseif path == "/correlations"
    return handler(CorrelationsController.call, payload)
  elseif path == "/basic_statistics"
    return handler(BasicStatisticsController.call, payload)
  elseif path == "/gaussian_mixtures"
    return handler(GaussianMixturesController.call, payload)
  elseif path == "/kmeans"
    return handler(KMeansController.call, payload)
  elseif path == "/gmeans"
    return handler(GMeansController.call, payload)
  elseif path == "/hierachical_clustering"
    return handler(HierachicalClusteringsController.call, payload)
  else
    throw(ErrorException("No route matches."))
  end
end


text:
"""
Dynamically import the top-level module given by `sym`, and return it if
possible.
"""
function dynamic_import_toplevel_module(sym)::Union{Module, Nothing}
    info("dynamic import: $sym")
    try
        eval(Main, :(import $sym))
        getfield(Main, sym)
    catch
        nothing
    end
end


text:
import Test

Test.@testset "PredictMDExtra.jl" begin
    import InteractiveUtils # stdlib
    import Pkg # stdlib
    import Test # stdlib

    @debug(string("Julia depot paths: "), Base.DEPOT_PATH)
    @debug(string("Julia load paths: "), Base.LOAD_PATH)

    @debug(string("Attempting to import PredictMDExtra...",))
    import PredictMDExtra
    @debug(string("Successfully imported PredictMDExtra.",))
    @debug(string("PredictMDExtra version: "),PredictMDExtra.version(),)
    @debug(string("PredictMDExtra package directory: "),PredictMDExtra.package_directory(),)

    @debug(string("Julia depot paths: "), Base.DEPOT_PATH)
    @debug(string("Julia load paths: "), Base.LOAD_PATH)

    Test.@testset "Unit tests              " begin
        testmodulea_filename::String = joinpath("PredictMDExtraTestModuleA", "PredictMDExtraTestModuleA.jl")
        testmoduleb_filename::String  = joinpath(
            "PredictMDExtraTestModuleB", "directory1", "directory2", "directory3",
            "directory4", "directory5", "PredictMDExtraTestModuleB.jl",
            )
        testmodulec_filename::String  = joinpath(mktempdir(), "PredictMDExtraTestModuleC.jl")
        rm(testmodulec_filename; force = true, recursive = true)
        open(testmodulec_filename, "w") do io
            write(io, "module PredictMDExtraTestModuleC end")
        end
        include(testmodulea_filename)
        include(testmoduleb_filename)
        include(testmodulec_filename)
        include(joinpath("test_package_directory.jl"))
        include(joinpath("test_package_list.jl"))
        include(joinpath("test_registry_url_list.jl"))
        include(joinpath("test_version.jl"))
    end

    Test.@testset "Import required packages" begin
        include(joinpath("test_import_required_packages.jl"))
    end

    Test.@testset "Test import_all()" begin
        include(joinpath("test_import_all.jl"))
    end
end


text:
import Base: print, show, copy, deepcopy, empty!, isempty, print, show

# Annotations
# ===========

"""
The `Annotations` type is basically a container for `Dict`s with the annotations of a multiple sequence alignment.
`Annotations` was designed for storage of annotations of the **Stockholm format**.

MIToS also uses MSA annotations to keep track of:

- **Modifications** of the MSA (`MIToS_...`) as deletion of sequences or columns.
- Positions numbers in the original MSA file (**column mapping:** `ColMap`)
- Position of the residues in the sequence (**sequence mapping:** `SeqMap`)
"""
type Annotations
    file::OrderedDict{ASCIIString, ByteString}
    sequences::Dict{Tuple{ASCIIString,ASCIIString},ASCIIString}
    columns::Dict{ASCIIString,ASCIIString}
    residues::Dict{Tuple{ASCIIString,ASCIIString},ASCIIString}
end

call(::Type{Annotations}) = Annotations(OrderedDict{ASCIIString, ByteString}(),
                                        Dict{Tuple{ASCIIString,ASCIIString},ASCIIString}(),
                                        Dict{ASCIIString, ASCIIString}(),
                                        Dict{Tuple{ASCIIString,ASCIIString},ASCIIString}())

"""
Creates empty MSA `Annotations` of length 0 using `sizehint!`
"""
empty(::Type{Annotations}) = Annotations( sizehint!(OrderedDict{ASCIIString, ByteString}(), 0), # There is not sizehint! for OrderedDict right now
                                         sizehint!(Dict{Tuple{ASCIIString,ASCIIString},ASCIIString}(), 0),
                                         sizehint!(Dict{ASCIIString, ASCIIString}(), 0),
                                         sizehint!(Dict{Tuple{ASCIIString,ASCIIString},ASCIIString}(), 0))

# Filters
# -------

# This function is useful because of the Julia issue 12495
_filter(str::ASCIIString, mask) = ASCIIString( str.data[mask] )

"""
For filter column and sequence mapping of the format: ",,,,10,11,,12"
"""
_filter_mapping(str_map::ASCIIString, mask) = join(split(str_map, ',')[mask], ',')

"""
`filtersequences!(data::Annotations, ids::IndexedArray, mask::AbstractArray{Bool,1})` is useful for deleting sequence annotations.
`ids` should be an `IndexedArray` with the `seqname`s of the annotated sequences and `mask` should be a logical vector.
"""
function filtersequences!(data::Annotations, ids::IndexedArray, mask::AbstractArray{Bool,1})
    if length(data.sequences) > 0 || length(data.residues) > 0
        del = ids[ !mask ]
        for key in keys(data.residues)
            if key[1] in del
                delete!(data.residues, key)
            end
        end
        for key in keys(data.sequences)
            if key[1] in del
                delete!(data.sequences, key)
            end
        end
        data.sequences = sizehint!(data.sequences, length(data.sequences))
        data.residues = sizehint!(data.residues, length(data.residues))
    end
    data
end

"""
`filtercolumns!(data::Annotations, mask)` is useful for deleting column annotations (creating a subset in place).
"""
function filtercolumns!(data::Annotations, mask)
    if length(data.residues) > 0
        for (key,value) in data.residues
            data.residues[key] = _filter(value, mask)
        end
    end
    if length(data.columns) > 0
        for (key,value) in data.columns
            data.columns[key] = _filter(value, mask)
        end
    end
    if length(data.sequences) > 0
        for (key,value) in data.sequences
            if key[2] == "SeqMap"
                data.sequences[key] = _filter_mapping(value, mask)
            end
        end
    end
    filecolumnmapping = get(data.file, "ColMap", "")
    if filecolumnmapping != ""
        data.file["ColMap"] = _filter_mapping(filecolumnmapping, mask)
    end
    data
end

# Copy, deepcopy and empty!
# -------------------------

for fun in [:copy, :deepcopy]
    @eval $(fun)(ann::Annotations) = Annotations( $(fun)( ann.file ), $(fun)( ann.sequences ), $(fun)( ann.columns ), $(fun)( ann.residues ))
end

function empty!(ann::Annotations)
    empty!(ann.file)
    empty!(ann.sequences)
    empty!(ann.columns)
    empty!(ann.residues)
    ann
end

isempty(ann::Annotations) = isempty(ann.file) && isempty(ann.sequences) && isempty(ann.columns) && isempty(ann.residues)

# ncolumns
# --------

"""
`ncolumns(ann::Annotations)` returns the number of columns/residues with annotations.
This function returns `-1` if there is not annotations per column/residue.
"""
function ncolumns(ann::Annotations)
    for value in values(ann.columns)
        return(length(value))
    end
    for value in values(ann.residues)
        return(length(value))
    end
    -1
end

# Getters & Setters
# -----------------

for (fun, field) in [ (:getannotfile,   :(ann.file)),
                     (:getannotcolumn, :(ann.columns))]
    @eval begin
        $(fun)(ann::Annotations) = $(field)
        $(fun)(ann::Annotations, feature::ASCIIString) = getindex($(field), feature)
        $(fun)(ann::Annotations, feature::ASCIIString, default::ASCIIString) = get($(field), feature, default)
    end
end

for (fun, field) in [ (:getannotsequence, :(ann.sequences)),
                     (:getannotresidue,  :(ann.residues))]
    @eval begin
        $(fun)(ann::Annotations) = $(field)
        $(fun)(ann::Annotations, seqname::ASCIIString, feature::ASCIIString) = getindex($(field), (seqname, feature))
        $(fun)(ann::Annotations, seqname::ASCIIString, feature::ASCIIString, default::ASCIIString) = get($(field), (seqname, feature), default)
    end
end

function _test_feature_name(feature::ASCIIString)
    length(feature) <= 50 || throw(ErrorException("Feature name has a limit of 50 characters."))
    ismatch(r"\s", feature) && throw(ErrorException("Feature name must not have spaces."))
end

function setannotfile!(ann::Annotations, feature::ASCIIString, annotation::ByteString)
    _test_feature_name(feature)
    previous = get(ann.file, feature, "")
    ann.file[feature] = previous != "" ? string(previous, '\n', annotation) : annotation
end

function setannotsequence!(ann::Annotations, seqname::ASCIIString, feature::ASCIIString, annotation::ASCIIString)
    _test_feature_name(feature)
    previous = get(ann.sequences, (seqname, feature), "")
    ann.sequences[(seqname, feature)] = previous != "" ? string(previous, '\n', annotation) : annotation
end

function setannotcolumn!(ann::Annotations, feature::ASCIIString, annotation::ASCIIString)
    _test_feature_name(feature)
    len = ncolumns(ann)
    if (len == -1) || (len == length(annotation))
        setindex!(ann.columns, annotation, feature)
    else
        throw(DimensionMismatch("You should have exactly 1 char per column ($len columns/residues)"))
    end
end

function setannotresidue!(ann::Annotations, seqname::ASCIIString, feature::ASCIIString, annotation::ASCIIString)
    _test_feature_name(feature)
    len = ncolumns(ann)
    if (len == -1) || (len == length(annotation))
        setindex!(ann.residues, annotation, (seqname, feature))
    else
        throw(DimensionMismatch("You should have exactly 1 char per residue ($len columns/residues)"))
    end
end

@doc "`getannotfile(ann[, feature[,default]])` returns per file annotation for `feature`" getannotfile
@doc "`getannotcolumn(ann[, feature[,default]])` returns per column annotation for `feature`" getannotcolumn
@doc "`getannotsequence(ann[, seqname, feature[,default]])` returns per sequence annotation for `(seqname, feature)`" getannotsequence
@doc "`getannotresidue(ann[, seqname, feature[,default]])` returns per residue annotation for `(seqname, feature)`" getannotresidue
@doc "`setannotfile!(ann, feature, annotation)` stores per file `annotation` for `feature`" setannotfile!
@doc "`setannotcolumn!(ann, feature, annotation)` stores per column `annotation` (1 char per column) for `feature`" setannotcolumn!
@doc "`setannotsequence!(ann, seqname, feature, annotation)` stores per sequence `annotation` for `(seqname, feature)`" setannotsequence!
@doc "`setannotresidue!(ann, seqname, feature, annotation)` stores per residue `annotation` (1 char per residue) for `(seqname, feature)`" setannotresidue!

# MIToS modification annotations
# ===============================

"""
Annotates on file annotations the modifications realized by MIToS on the MSA
"""
function annotate_modification!(ann::Annotations, modification::ASCIIString)
    setannotfile!(ann, string("MIToS_", Dates.now()), modification)
    true # generally used on bool context: annotate && annotate_modification!(...
end

"""
Deletes all the MIToS annotated modifications
"""
function delete_annotated_modifications!(ann::Annotations)
    for key in keys(ann.file)
        if length(key) > 6 && key[1:6] == "MIToS_"
            delete!(ann.file, key)
        end
    end
end

using Base.Markdown

"""
Prints MIToS annotated modifications
"""
function printmodifications(ann::Annotations)
    for (key,value) in ann.file
        if length(key) > 6 && key[1:6] == "MIToS_"
            list_k = split(key, '_')
            println("-------------------")
            println(list_k[2])
            println('\n', value)
        end
    end
end

# Show & Print Annotations
# ========================

function _printfileannotations(io::IO, ann::Annotations)
    if !isempty(ann.file)
        for (key, value) in ann.file
            for val in split(value, '\n')
                println(io, string("#=GF ", key, "   ", val))
            end
        end
    end
end

function _printcolumnsannotations(io::IO, ann::Annotations)
    if !isempty(ann.columns)
        for (key, value) in ann.columns
            println(io, string("#=GC ", key, "\t\t\t", value))
        end
    end
end

function _printsequencesannotations(io::IO, ann::Annotations)
    if !isempty(ann.sequences)
        for (key, value) in ann.sequences
            for val in split(value, '\n')
                println(io, string("#=GS ", key[1], '\t', key[2], ' ', val))
            end
        end
    end
end

function _printresiduesannotations(io::IO, ann::Annotations)
    if !isempty(ann.residues)
        for (key, value) in ann.residues
            println(io, string("#=GR ", key[1], '\t', key[2], '\t', value))
        end
    end
end

function print(io::IO, ann::Annotations)
    _printfileannotations(io, ann)
    _printsequencesannotations(io, ann)
    _printresiduesannotations(io, ann)
    _printcolumnsannotations(io, ann)
end

show(io::IO, ann::Annotations) = print(io, ann)


text:
using AFL
using Documenter

makedocs(;
    modules=[AFL],
    authors="Joost <sevenfourtwo@protonmail.com> and contributors",
    repo="https://github.com/sevenfourtwo/AFL.jl/blob/{commit}{path}#L{line}",
    sitename="AFL.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://sevenfourtwo.github.io/AFL.jl",
        assets=String[],
    ),
    pages=[
        "index.md",
        "example.md",
        "internals.md",
        "api.md",
    ],
)

deploydocs(;
    repo="github.com/sevenfourtwo/AFL.jl",
)


text:
"""
    mlmnet_perms(data, lambdas, alphas;
                  method = "ista", isNaive =false,
                  permFun = shuffle_rows,
                  hasXIntercept=true, hasZIntercept=true,
                  toXReg=trues(data.p),
                  toZReg=trues(data.q),
                  toXInterceptReg::Bool=false, toZInterceptReg::Bool=false,
                  toNormalize=true, isVerbose=true,
                  stepsize=0.01, setStepsize=true,
                  funArgs...)

Permutes response matrix Y in RawData object and then calls the mlmnet core
function.

# Arguments

- data = RawData object
- lambdas = 1d array of floats consisting of the total penalties in descending
  order. If they are not in descending order, they will be sorted.
- alphas = 1d array of floats consisting of the penalty ratio that
  determines the mix of penalties between L1 and L2

# Keyword arguments

- methods = function name that applies the Elastic-net penalty estimate method;
  default is `ista`, and the other methods are `fista`, `fista_bt`, `admm` and `cd`
- isNaive = boolean flag indicating whether to solve the Naive or non-Naive
  Elastic-net problem
- permFun = function used to permute `Y`. Defaults to `shuffle_rows`
  (shuffles rows of `Y`).
- hasXIntercept = boolean flag indicating whether or not to include an `X`
  intercept (row main effects). Defaults to `true`.
- hasZIntercept = boolean flag indicating whether or not to include a `Z`
  intercept (column main effects). Defaults to `true`.
- toXReg = 1d array of bit flags indicating whether or not to regularize each
  of the `X` (row) effects. Defaults to 2d array of `true`s with length
  equal to the number of `X` effects (equivalent to `data.p`).
- toZReg = 1d array of bit flags indicating whether or not to regularize each
  of the `Z` (column) effects. Defaults to 2d array of `true`s with length
  equal to the number of `Z` effects (equivalent to `data.q`).
- toXInterceptReg = boolean flag indicating whether or not to regularize the
  `X` intercept Defaults to `false`.
- toZInterceptReg = boolean flag indicating whether or not to regularize the
  `Z` intercept. Defaults to `false`.
- toNormalize = boolean flag indicating if the columns of `X` and `Z`
  should be standardized (to mean 0, standard deviation 1). Defaults to `true`.
- isVerbose = boolean flag indicating whether or not to print messages.
  Defaults to `true`.
- setStepsize = boolean flag indicating whether the fixed step size should be
  calculated (for `ista!` and `fista!`). Defaults to `true`.
- stepsize = float; step size for updates (irrelevant for coordinate
  descent and when `setStepsize` is set to `true` for `ista!` and `fista!`).
  Defaults to `0.01`.
- funArgs = variable keyword arguments to be passed into `fun`

# Value

An MLMnet object

# Some notes

The default method for choosing the fixed step size for `fista!` or `ista!`
is to use the reciprocal of the product of the maximum eigenvalues of
`X*transpose(X)` and `Z*transpose(Z)`. This is computed when `fista!` or
`ista!` is passed into the `fun` argument and `setStepsize` is set to `true`.
If `setStepsize` is set to `false`, the value of the `stepsize` argument will
be used as the fixed step size. Note that obtaining the eigenvalues when `X`
and/or `Z` are very large may exceed computational limitations.

Specifying a good starting step size (`stepsize`) and multiplying factor
(`gamma`) when `fista_bt!` is passed into the `fun` argument can be difficult.
Shrinking the step size too gradually can result in slow convergence. Doing so
too quickly can cause the criterion to diverge. We have found that setting
`stepsize` to 0.01 often works well in practice; choice of `gamma` appears to
be less consequential.

"""
function mlmnet_perms(data::RawData,
                      lambdas::AbstractArray{Float64,1}, alphas::AbstractArray{Float64,1};
                      method::String = "ista", isNaive::Bool=false,
                      permFun::Function=shuffle_rows,
                      hasXIntercept::Bool=true, hasZIntercept::Bool=true,
                      toXReg::BitArray{1}=trues(data.p),
                      toZReg::BitArray{1}=trues(data.q),
                      toXInterceptReg::Bool=false,
                      toZInterceptReg::Bool=false,
                      toNormalize::Bool=true, isVerbose::Bool=true,
                      stepsize::Float64=0.01, setStepsize=true, funArgs...)

    # Create RawData object with permuted Y
    dataPerm = RawData(Response(permFun(get_Y(data))), data.predictors)

    # Run penalty on the permuted data
    return mlmnet(dataPerm, lambdas, alphas;
                    method = method, isNaive=isNaive,
                    hasXIntercept=hasXIntercept, hasZIntercept=hasZIntercept,
                    toXReg=toXReg, toZReg=toZReg,
                    toXInterceptReg=toXInterceptReg, toZInterceptReg=toZInterceptReg,
                    toNormalize=toNormalize, isVerbose=isVerbose,
                    stepsize=stepsize, setStepsize=setStepsize,
                    funArgs...)

end

"""
mlmnet_perms(data, lambdas;
             method = "ista", isNaive =false,
             permFun = shuffle_rows,
             hasXIntercept=true, hasZIntercept=true,
             toXReg=trues(data.p),
             toZReg=trues(data.q),
             toXInterceptReg::Bool=false, toZInterceptReg::Bool=false,
             toNormalize=true, isVerbose=true,
             stepsize=0.01, setStepsize=true,
             funArgs...)

"""

function mlmnet_perms(data::RawData,
                      lambdas::AbstractArray{Float64,1};
                      method::String = "ista", isNaive::Bool=false,
                      permFun::Function=shuffle_rows,
                      hasXIntercept::Bool=true, hasZIntercept::Bool=true,
                      toXReg::BitArray{1}=trues(data.p),
                      toZReg::BitArray{1}=trues(data.q),
                      toXInterceptReg::Bool=false,
                      toZInterceptReg::Bool=false,
                      toNormalize::Bool=true, isVerbose::Bool=true,
                      stepsize::Float64=0.01, setStepsize=true, funArgs...)

    # Create RawData object with permuted Y
    dataPerm = RawData(Response(permFun(get_Y(data))), data.predictors)

    alphas = [1.0] # default LASSO, 𝛼 = 1

    # Run L1-L2 penalties on the permuted data
  rslts =  mlmnet(dataPerm, lambdas, alphas;
                            method = method, isNaive=isNaive,
                            hasXIntercept=hasXIntercept, hasZIntercept=hasZIntercept,
                            toXReg=toXReg, toZReg=toZReg,
                            toXInterceptReg=toXInterceptReg, toZInterceptReg=toZInterceptReg,
                            toNormalize=toNormalize, isVerbose=isVerbose,
                            stepsize=stepsize, setStepsize=setStepsize,
                            funArgs...)

return rslts

end


text:
# Patrick Keast, MODERATE-DEGREE TETRAHEDRAL QUADRATURE FORMULAS
# http://mech.fsv.cvut.cz/oofem/resources/doc/oofemrefman/gaussintegrationrule_8C_source.html
function _get_gauss_tetdata(n::Int)
    if n == 1
        a = 1. / 4.
        w = 1. / 6.
        xw = [a a a w]
    elseif n == 2
        a = ( 5. + 3. * √(5.) ) / 20.
        b = ( 5. - √(5.) ) / 20.
        w = 1. / 24.
        xw = [a b b w
              b a b w
              b b a w
              b b b w]
    elseif n == 3
        a1 = 1. / 4.
        a2 = 1. / 2.
        b2 = 1. / 6.
        w1 = -2. / 15.
        w2 = 3. / 40.
        xw = [a1 a1 a1 w1
              a2 b2 b2 w2
              b2 a2 b2 w2
              b2 b2 a2 w2
              b2 b2 b2 w2]
    elseif n == 4
        a1 = 1. / 4.;
        w1 = -74. / 5625.;

        a2 = 5. / 70.;
        b2 = 11. / 14.;
        w2 = 343. / 45000.;

        a3 = ( 1. + √(5. / 14.) ) / 4.;
        b3 = ( 1. - √(5. / 14.) ) / 4.;
        w3 = 28. / 1125.;

        xw = [a1 a1 a1 w1
              b2 a2 a2 w2
              a2 b2 a2 w2
              a2 a2 b2 w2
              a2 a2 a2 w2
              a3 a3 b3 w3
              a3 b3 a3 w3
              a3 b3 b3 w3
              b3 a3 a3 w3
              b3 a3 b3 w3
              b3 b3 a3 w3]
    else
        throw(ArgumentError("unsupported order for tetraheder gauss-legendre integration"))
    end
    return xw
end


text:
using BlockHankelMatrices
using Test

@testset "BlockHankelMatrices.jl" begin
    # Write your tests here.
end