text:
using BinaryBuilder

name = "czmq"
version = v"4.2.1"

# Collection of sources required to build ZMQ
sources = [
    GitSource("https://github.com/zeromq/czmq.git", "4a50c2153586cf510d6cc3fcfbb9f5ea2e02c419"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/czmq
sh autogen.sh

# Hint to find libstc++, required to link against C++ libs when using C compiler
if [[ "${target}" == *-linux-* ]]; then
    if [[ "${nbits}" == 32 ]]; then
        export CFLAGS="-Wl,-rpath-link,/opt/${target}/${target}/lib"
    else
        export CFLAGS="-Wl,-rpath-link,/opt/${target}/${target}/lib64"
    fi
fi

atomic_patch -p1 ../patches/ntohll.patch

./configure --prefix=$prefix \
    --build=${MACHTYPE} \
    --host=${target} \
    --enable-drafts \
    --disable-static \
    CXXFLAGS="-O2 -fms-extensions"
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libczmq", :libczmq),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
    Dependency("ZeroMQ_jll", v"4.3.2")
]

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
module KroneckerProducts
using LinearAlgebra
using LinearAlgebra: checksquare
import LinearAlgebra: \, *, /, rdiv!, ldiv!
include("lu.jl")
using LazyInverses

const AbstractMatOrFac{T} = Union{AbstractMatrix{T}, Factorization{T}}
LinearAlgebra.factorize(F::Factorization) = F

export kronecker, ⊗, KroneckerProduct, FactorizedKroneckerProduct, kronecker_power

# has field for temporary storage needed for non-allocating multiplies and solves
# currently made struct mutable and automatically store latest temporary in structure
# IDEA: could allocate two potentially large arrays for MatMatMul,
# and create temporary views for MatVecMul
mutable struct KroneckerProduct{T, A<:Tuple, V} <: Factorization{T}
    factors::A
    temporaries::V
end
const FactorizedKroneckerProduct{T} = KroneckerProduct{T, <:Tuple{Vararg{Factorization}}}

Base.eltype(K::KroneckerProduct{T}) where {T} = T
Base.IndexStyle(::Type{<:KroneckerProduct}) = IndexCartesian()

function KroneckerProduct(A::Tuple{Vararg{Union{Number, AbstractMatOrFac}}})
    T = promote_type(eltype.(A)...)
    temporaries = nothing
    T_temp = Union{Nothing, NTuple{length(A), Vector}}
    KroneckerProduct{T, typeof(A), T_temp}(A, temporaries)
end

# num_cols is number of columns to be pre-allocated for mul!, if only matrix-vector
# multiplies are required, num_cols = 1 is sufficient
# eltype_temp is the element type of the temporaries to be allocate
function KroneckerProduct(A::Tuple{Vararg{Union{Number, AbstractMatOrFac}}},
                          num_cols::Int, eltype_temp = promote_type(eltype.(A)...))
    temporaries = initialize_temporaries(A, eltype_temp, num_cols)
    T = promote_type(eltype.(A)...)
    KroneckerProduct{T, typeof(A), typeof(temporaries)}(A, temporaries)
end

function kronecker end # smart constructor
const ⊗ = kronecker
kronecker(A::Union{Number, AbstractMatOrFac}) = A
kronecker(A::Union{Number, AbstractMatOrFac}...) = KroneckerProduct(A)
kronecker(A::Tuple) = KroneckerProduct(A)
kronecker(K::KroneckerProduct) = K

function Base.copy(K::KroneckerProduct)
    factors = tuple(copy.(K.factors)...)
    temporaries = isnothing(K.temporaries) ? nothing : tuple(copy.(K.temporaries)...)
    typeof(K)(factors, temporaries)
end

# IDEA: this makes sure we have a flat Kronecker hierarchy? could have separate "flatten" function
# only makes sense to break if we have a subset, which is a power
kronecker(A::Union{Number, AbstractMatOrFac}, K::KroneckerProduct) = kronecker(A, K.factors...)
kronecker(K::KroneckerProduct, A::Union{Number, AbstractMatOrFac}) = kronecker(K.factors..., A)
kronecker(K::KroneckerProduct, L::KroneckerProduct) = kronecker(K.factors..., L.factors...)

# create new product by applying function to each component matrix
kronecker(f, K::KroneckerProduct) = kronecker((f(A) for A in K.factors)...) #kronecker(f.(K.factors)...)

kronecker_power(A, pow::Integer) = kronecker((A for _ in 1:pow)...)

Base.size(K::KroneckerProduct) = (size(K, 1), size(K, 2))
Base.size(K::KroneckerProduct, dim::Integer) = prod(A->size(A, dim), K.factors)

function collect(K::KroneckerProduct)
    f(A) = A isa Factorization ? AbstractMatrix(A) : A
    order(K) == 1 ? f(K.factors[1]) : kron(f.(K.factors)...)
end
LinearAlgebra.Matrix(K::KroneckerProduct) = collect(K)

issquare(A::AbstractMatOrFac) = size(A, 1) == size(A, 2)
# IDEA: meta programming
allsquare(K::KroneckerProduct) = all(issquare, K.factors)
LinearAlgebra.issymmetric(K::KroneckerProduct) = all(issymmetric, K.factors)
LinearAlgebra.isposdef(K::KroneckerProduct) = all(isposdef, K.factors)
# KroneckerProduct with factorizations
LinearAlgebra.issuccess(K::KroneckerProduct) = all(issuccess, K.factors)

# isfactorized(W::Woodbury) = all(isfactorized, (W.A, W.C))
# checks if all component matrices of a KroneckerProduct have the same size
samesize(K::KroneckerProduct) = samesize(K.factors)
samesize(A...) = samesize(A)
samesize(A::Tuple) = all(==(size(A[1])), A)

# checks if all component matrices are the same
ispower(K::KroneckerProduct) = isequiv(K) || all(==(K.factors[1]), K.factors)
# same as ispower but only does pointer check (constant for each factor)
isequiv(K::KroneckerProduct) = all(F -> F === K.factors[1], K.factors)
order(M::AbstractMatrix) = 1
order(K::KroneckerProduct) = length(K.factors)

function LinearAlgebra.det(K::KroneckerProduct)
    logabsdet_K, sign_K = logabsdet(K)
    exp(logabsdet_K) * sign_K
end
function LinearAlgebra.logdet(K::KroneckerProduct)
    logabsdet_K, sign_K = logabsdet(K)
    if sign_K ≈ 1
        logabsdet_K
    else
        throw(DomaingError("determinant non-positive: sign = $sign_K"))
    end
end
function LinearAlgebra.logabsdet(K::KroneckerProduct)
    n = checksquare(K)
    if !allsquare(K) # implies that there is a rank deficient factor
        real(eltype(K))(-Inf), zero(eltype(K))
    elseif isequiv(K) # if all factors are equivalent (have same pointer), use more efficient formula
        A = K.factors[1]
        logabsdet_A, sign_A = logabsdet(A)
        p = length(K.factors)
        k = p * (n ÷ size(A, 1))
        k * logabsdet_A, sign_A^k
    else
        logabsdet_K = zero(eltype(K))
        sign_K = one(eltype(K))
        for A in K.factors
            logabsdet_A, sign_A = logabsdet(A)
            k = n ÷ size(A, 1)
            logabsdet_K += k * logabsdet_A
            sign_K *= sign_A^k
        end
        return logabsdet_K, sign_K
    end
end

function LinearAlgebra.tr(K::KroneckerProduct)
    n = checksquare(K)
    allsquare(K) ? prod(tr, K.factors) : sum(K[i] for i in diagind(K))
end

# TODO: delete this?
function Base.inv(K::KroneckerProduct)
    checksquare(K)
    allsquare(K) ? ⊗(inv, K) : throw(SingularException(1))
end

function LazyInverses.inverse(K::KroneckerProduct)
    checksquare(K)
    allsquare(K) ? ⊗(inverse, K) : throw(SingularException(1))
end
LinearAlgebra.pinv(K::KroneckerProduct) = ⊗(pinv, K)
LazyInverses.pseudoinverse(K::KroneckerProduct) = ⊗(pseudoinverse, K)

LinearAlgebra.factorize(K::KroneckerProduct) = ⊗(factorize, K)
LinearAlgebra.adjoint(K::KroneckerProduct) = ⊗(adjoint, K)
LinearAlgebra.transpose(K::KroneckerProduct) = ⊗(transpose, K)
LinearAlgebra.conj(K::KroneckerProduct) = ⊗(conj, K)

# useful for least-squares problems
function LinearAlgebra.qr(K::KroneckerProduct, piv = NoPivot())
    kronecker(A->qr(A, piv), K)
end
function LinearAlgebra.cholesky(K::KroneckerProduct, piv = Val(false); check = true)
    kronecker(A->cholesky(A, piv, check = check), K)
end

############################ multiplication ####################################
# mixed product property
function *(K1::KroneckerProduct, K2::KroneckerProduct)
    size(K1, 2) == size(K2, 1) || throw(DimensionMismatch("size(K1, 2) = $(size(K1, 2)) ≠ $(size(K2, 1)) = size(K2, 1)"))
    if all((AB) -> size(AB[1]) == size(AB[2]), zip(K1.factors, K2.factors)) # if the components have the same sizes
        kronecker(K1.factors .* K2.factors)
    else # fallback
        Matrix(K1) * Matrix(K2)
    end
end

*(a::Number, K::KroneckerProduct) = ⊗(a * K.factors[1], K.factors[2:end]...)
*(K::KroneckerProduct, a::Number) = ⊗(K.factors[1:end-1]..., K.factors[end] * a)

function LinearAlgebra.lmul!(a::Number, K::KroneckerProduct)
    lmul!(a, K.factors[1])
    return K
end

function LinearAlgebra.rmul!(K::KroneckerProduct, a::Number)
    rmul!(K.factors[end], a)
    return K
end

################################ vec trick #####################################
function *(K::KroneckerProduct, x::AbstractVector)
    T = promote_type(eltype(K), eltype(x))
    y = zeros(T, size(K, 1))
    mul!(y, K, x)
end
function *(K::KroneckerProduct, X::AbstractMatrix)
    T = promote_type(eltype(K), eltype(X))
    Y = zeros(T, size(K, 1), size(X, 2))
    mul!(Y, K, X)
end
*(x::Adjoint{<:Number, <:AbstractVector}, K::KroneckerProduct) = adjoint(K'*x')
*(x::AbstractMatrix, K::KroneckerProduct) = adjoint(K'*x')
\(K::KroneckerProduct, x::AbstractVecOrMat) = pseudoinverse(factorize(K)) * x
/(x::AbstractVecOrMat, K::KroneckerProduct) = x * pseudoinverse(factorize(K))

function LinearAlgebra.dot(x::AbstractVecOrMat, K::KroneckerProduct, y::AbstractVecOrMat)
    dot(x, K*y) # fallback for now
end

function *(x::AbstractMatrix, K::KroneckerProduct, y::AbstractVecOrMat)
    x*(K*y)
end

function LinearAlgebra.mul!(y::AbstractVecOrMat, K::KroneckerProduct, x::AbstractVecOrMat,
                            α::Real = 1, β::Real = 0)
    _mul_check_args(y, K, x)
    X = x
    temporaries = set_temporaries!(K, x, y)
    for (A, z) in Iterators.reverse(zip(K.factors, temporaries)) # this loop takes up the majority of the runtime
        X = reshape(X, size(A, 2), :)
        Z = reshape(z, :, size(A, 1))
        mul!(Z, transpose(X), transpose(A)) # same as Z = transpose(A*X)
        X = Z # pointer change
    end
    Kx = (x isa AbstractVector) ? vec(X) : transpose(reshape(X, (size(x, 2), :)))
    @. y = α * Kx + β * y
end

function LinearAlgebra.ldiv!(y::AbstractVecOrMat, K::KroneckerProduct, x::AbstractVecOrMat)
    mul!(y, pseudoinverse(K), x)
end

# we can do this because temporaries are used to carry out the intermediate multiplications
function LinearAlgebra.ldiv!(K::KroneckerProduct, x::AbstractVecOrMat)
    mul!(x, pseudoinverse(K), x)
end

function _mul_check_args(y::AbstractVecOrMat, K::KroneckerProduct, x::AbstractVecOrMat)
    size(K, 2) ≠ size(x, 1) && throw(DimensionMismatch("size(K, 2) = $(size(K, 2)) ≠ $(size(x, 1)) = size(x, 1)"))
    size(K, 1) ≠ size(y, 1) && throw(DimensionMismatch("size(K, 1) = $(size(K, 1)) ≠ $(size(y, 1)) = size(y, 1)"))
    size(x, 2) ≠ size(y, 2) && throw(DimensionMismatch("size(x, 2) = $(size(x, 2)) ≠ $(size(y, 2)) = size(y, 2)"))
    return true
end

# returns temporaries necessary for efficient multiplication y = K*x with KroneckerProduct
function set_temporaries!(K::KroneckerProduct, x::AbstractVecOrMat, y::AbstractVecOrMat)
    temporaries = K.temporaries
    allocated = !isnothing(temporaries)
    if allocated
        correct_type = eltype(temporaries[1]) == eltype(y) # is a little more restrictive than necessary, i.e. real could be cast into complex
        correct_length = length(y) == length(temporaries[1])
        allocated = correct_type && correct_length
    end
    # need to allocate temporaries if size is not compatible, or element type differs
    if !allocated
        temporaries = initialize_temporaries(K, y)
    end
    return K.temporaries = temporaries
end

# allocates temporaries for multiplication and solve,
# if external vectors have different element types or different sizes,
# need to allocate new temporaries
# does not allocate new temporaries if existing ones satisfy the type and size
# requirement for the current multiplication
function initialize_temporaries(K::KroneckerProduct, y::AbstractVecOrMat)
    initialize_temporaries(K.factors, eltype(y), size(y, 2))
end
function initialize_temporaries(factors::Tuple, T = promote_type(eltype.(A)...), num_cols::Int = 1)
    # first, calculate maximum length of temporary array
    n = prod(A->size(A, 2), factors) * num_cols
    max_n = n
    for A in Iterators.reverse(factors)
        n = size(A, 1) * (n ÷ size(A, 2))
        max_n = max(n, max_n)
    end
    # second, allocate two arrays of this size and create views into them alternatingly
    t1, t2 = zeros(T, max_n), zeros(T, max_n) # only need to allocate two sufficiently large arrays
    temporaries = []
    n = prod(A->size(A, 2), factors) * num_cols
    for (i, A) in enumerate(Iterators.reverse(factors))
        n = size(A, 1) * (n ÷ size(A, 2))
        t = @view t1[1:n]
        # @time temporaries[i] = t
        push!(temporaries, t)
        t1, t2 = t2, t1
    end
    return tuple(Iterators.reverse(temporaries)...)
end

############################### indexing #######################################
# test if components support getindex ...
function Base.getindex(K::KroneckerProduct, i::Integer, j::Integer)
    @boundscheck checkbounds(K, i, j)
    val = one(eltype(K))
    n, m = size(K)
    @inbounds for A in K.factors
        n ÷= size(A, 1); m ÷= size(A, 2);
        val *= A[cld(i, n), cld(j, m)]  # can this be replaced with fld1, mod1?
        i, j = mod1(i, n), mod1(j, m) # or some kind of shifted remainder?
    end
    return val
end

# we can delete this if we subtype AbstractMatrix
function Base.checkbounds(K::KroneckerProduct, i::Integer, j::Integer)
    (1 ≤ i ≤ size(K, 1) && 1 ≤ j ≤ size(K, 2)) || throw(BoundsError(K, [i,j]))
end

end # module

ext: jl


text:
using BinaryProvider # requires BinaryProvider 0.3.0 or later

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
products = [
    LibraryProduct(prefix, ["libscsindir"], :indirect),
    LibraryProduct(prefix, ["libscsdir"], :direct),
]

# Download binaries from hosted location
bin_prefix = "https://github.com/JuliaOpt/SCSBuilder/releases/download/v2.0.2-NOGOMP"

# Listing of files generated by BinaryBuilder:
download_info = Dict(
    Linux(:aarch64, :glibc) => ("$bin_prefix/SCSBuilder.v2.0.2.aarch64-linux-gnu.tar.gz", "9ac63df8629d95588d84f0c8219d7c15cb229eab8fb78d528a4e75fbd2f5c136"),
    Linux(:armv7l, :glibc, :eabihf) => ("$bin_prefix/SCSBuilder.v2.0.2.arm-linux-gnueabihf.tar.gz", "1c1241faa174ff47c2c290793b6fec1a9f0db290fda57c7826f9b8942b41943b"),
    Linux(:i686, :glibc) => ("$bin_prefix/SCSBuilder.v2.0.2.i686-linux-gnu.tar.gz", "7731d8a2e644bf1e070e013eb100d66e49f17796b8903c3350a806946f7ba752"),
    Windows(:i686) => ("$bin_prefix/SCSBuilder.v2.0.2.i686-w64-mingw32.tar.gz", "354cce228b102e9a09a727d44c1f437ad49e97077ba4abff96b22d6d129c751f"),
    Linux(:powerpc64le, :glibc) => ("$bin_prefix/SCSBuilder.v2.0.2.powerpc64le-linux-gnu.tar.gz", "7a2847ad3f8621a6927e4281d8ee7f2704a3339cfed41eb1b227a997a59cadf5"),
    MacOS(:x86_64) => ("$bin_prefix/SCSBuilder.v2.0.2.x86_64-apple-darwin14.tar.gz", "dc7f0d7bb3f2ee0072cb81653729ac2f5123abfa0735898cd44bc39a2952d248"),
    Linux(:x86_64, :glibc) => ("$bin_prefix/SCSBuilder.v2.0.2.x86_64-linux-gnu.tar.gz", "d5503f18e7d504e4a496ae4ea79ef326fab78983cf42ebd2a691e0a682790d96"),
    FreeBSD(:x86_64) => ("$bin_prefix/SCSBuilder.v2.0.2.x86_64-unknown-freebsd11.1.tar.gz", "3ef5005d53ca7dee21842856ca5628e0efc244083e125e939304a14092d54858"),
    Windows(:x86_64) => ("$bin_prefix/SCSBuilder.v2.0.2.x86_64-w64-mingw32.tar.gz", "289504d5df960c3434e5877f16cc42949e32280eb0ec684020bf6f9336c9316f"),
)

custom_library = false
if haskey(ENV,"JULIA_SCS_LIBRARY_PATH")
    custom_products = [LibraryProduct(ENV["JULIA_SCS_LIBRARY_PATH"],product.libnames,product.variable_name) for product in products]
    if all(satisfied(p; verbose=verbose) for p in custom_products)
        products = custom_products
        custom_library = true
    else
        error("Could not install custom libraries from $(ENV["JULIA_SCS_LIBRARY_PATH"]).\nTo fall back to BinaryProvider call delete!(ENV,\"JULIA_SCS_LIBRARY_PATH\") and run build again.")
    end
end

# Install unsatisfied or updated dependencies:
unsatisfied = any(!satisfied(p; verbose=verbose) for p in products)
if !custom_library && haskey(download_info, platform_key())
    url, tarball_hash = download_info[platform_key()]
    if unsatisfied || !isinstalled(url, tarball_hash; prefix=prefix)
        # Download and install binaries
        install(url, tarball_hash; prefix=prefix, force=true, verbose=verbose)
    end
elseif unsatisfied
    # If we don't have a BinaryProvider-compatible .tar.gz to download, complain.
    # Alternatively, you could attempt to install from a separate provider,
    # build from source or something even more ambitious here.
    error("Your platform $(triplet(platform_key())) is not supported by this package!")
end

# Write out a deps.jl file that will contain mappings for our products
write_deps_file(joinpath(@__DIR__, "deps.jl"), products, verbose=verbose)


text:
include("JToricBuilder.jl")

InstallDevelopmentVersionNConvex()

max_stars_repo_name: fingolfin/JToric

max_issues_repo_name: fingolfin/JToric

max_forks_repo_name: fingolfin/JToric


text:
#=
Created on Thursday 2 April 2020
Last update:

@author: Michiel Stock
michielfmstock@gmail.com

Experiments on the honeybee spillover dataset.
=#

using DataFrames, CSV
using BSON
using Plots, Statistics

interactions = CSV.read("data/Magrach2017_network.csv")
first(interactions, 5)

names(interactions)

sites = unique(interactions.site_id)

pollinators = unique(interactions.pollinator_species);
plants = unique(interactions.plant_species);

make_inverse_index = l -> Dict(n => i for (i, n) in enumerate(l))

pollinators_inv_index = make_inverse_index(pollinators);
plants_inv_index = make_inverse_index(plants);

println("Plants:")
for p in plants
    println("- $p")
end

println("Pollinators:")
for p in pollinators
    println("- $p")
end

function make_incidence(interactions, pl_inv_index, pol_inv_index)
    Y = zeros(length(pl_inv_index), length(pol_inv_index))
    for (pl, pol, int) in zip(interactions[:,:plant_species], interactions[:,:pollinator_species], interactions[:,:transect])
        i, j = pl_inv_index[pl], pol_inv_index[pol]
        Y[i,j] += int
    end
    return Y
end

Y = make_incidence(interactions, plants_inv_index, pollinators_inv_index)

networks_after = Dict{String, Array{Float64,2}}()
networks_during = Dict{String, Array{Float64,2}}()

for site in interactions.site_id
    interactions_site = interactions[interactions.site_id.==site,:]
    Y_after = make_incidence(interactions_site[interactions_site.period.=="after",:],
                            plants_inv_index, pollinators_inv_index)
    Y_during = make_incidence(interactions_site[interactions_site.period.=="during",:],
                            plants_inv_index, pollinators_inv_index)

    networks_after[site] = Y_after
    networks_during[site] = Y_during
end

for (time, networks) in zip(["after", "during"], [networks_after, networks_during])
    println("$time")
    for (site, Y) in networks
        println("\t- $site: $(sum(Y)) ($(sum(Y.>0)))")
    end
end

using OTSIN

Ps = [(Y .+ 1e-4) |> OTSIN.normalize for Y in values(networks_after)];

γ = 0.010

@time M_ab = fitM(Ps, fix_a=true, fix_b=true, maxitter=20, γ=γ);

@time M_a = fitM(Ps, fix_a=true, fix_b=false, γ=γ);

@time M_b = fitM(Ps, fix_a=false, fix_b=true, γ=γ);

@time M = fitM(Ps, fix_a=false, fix_b=false, γ=γ);

plot(heatmap(M, title="M"), heatmap(M_a, title="M (a fixed)"),
        heatmap(M_b, title="M (b fixed)"), heatmap(M_ab, title="M (a, b fixed)"))
savefig("figures/M_honeybee_spillover.svg")

Ps_during = [Y |> OTSIN.normalize for Y in values(networks_during)];

kl_ot_ab = Float64[]
kl_ot_a = Float64[]
kl_ot_b = Float64[]
kl_ot_softmax = Float64[]
kl_ot_neutral = Float64[]

for P in Ps_during
    a, b = marginals(P)
    P = P[a.>0, b.>0]
    M = M_ab[a.>0, b.>0]
    a, b = marginals(P)
    # optimal transport both
    Q = optimaltransport(M, a, b)
    push!(kl_ot_ab, KL(P, Q))
    Q .= optimaltransport(M, a, nothing)
    push!(kl_ot_a, KL(P, Q))
    Q .= optimaltransport(M, nothing, b)
    push!(kl_ot_b, KL(P, Q))
    Q .= optimaltransport(M, nothing, nothing)
    push!(kl_ot_softmax, KL(P, Q))
    Q .= optimaltransport(a, b)
    push!(kl_ot_neutral, KL(P, Q))
end

Ps_after = Ps

BSON.@save "results/spillover.bson" Ps_during Ps_after M M_a M_b M_ab kl_ot_ab  kl_ot_a kl_ot_b kl_ot_softmax kl_ot_neutral γ

max_stars_repo_name: PoisotLab/OTSIN

max_issues_repo_name: PoisotLab/OTSIN

max_forks_repo_name: PoisotLab/OTSIN


text:
module PSSFSSLen
using Reexport

@reexport using Unitful
export mm, cm, inch, mil, PSSFSSLength

const mm = unit(1u"mm")
const cm = unit(1u"cm")
const inch = unit(1u"inch")
const mil = unit(1u"mil")
const PSSFSSLength = Union{typeof(mm),typeof(cm),typeof(inch),typeof(mil)}

end


text:
using MeshCat
using MeshCatMechanisms
rbd = MeshCatMechanisms.rbd

vis =  Visualizer()
render(vis)

urdf = "models/bob.urdf"
robot = MeshCatMechanisms.parse_urdf(urdf)
mvis = MechanismVisualizer(robot, URDFVisuals(urdf), vis)

neck = rbd.findjoint(robot, "neck")
robotstate = rbd.MechanismState(robot)
rbd.set_configuration!(mvis, neck, -0.0)

left_gripper = rbd.findjoint(robot, "rack_to_left_fork")
rbd.set_configuration!(mvis, left_gripper, 0.05)

right_gripper = rbd.findjoint(robot, "rack_to_right_fork")
rbd.set_configuration!(mvis, right_gripper, -0.05)

rbd.set_configuration!(robotstate, neck, -1.0)
rbd.configuration(robotstate)


text:
import PolyJuMP

@testset "Creating polynomial with empty MonomialVector" begin
    @polyvar x
    X = emptymonovec(typeof(x^2))
    model = SOSModel()
    v = PolyJuMP.Variable(SOSPoly(X), false, false)
    #@inferred JuMP.add_variable(model, v) # FIXME broken
    @test JuMP.add_variable(model, v) == 0
end
@testset "GramMatrix JuMP.value" begin
    model = JuMP.Model()
    @variable model α
    @variable model β
    @polyvar x y
    q = GramMatrix([α β; β α], [x, y])
    JuMP.fix(α, 2)
    JuMP.fix(β, 3)
    @test_broken JuMP.value(q) == 2x^2 + 2y^2 + 6x*y
end
@testset "Container of GramMatrix" begin
    model = JuMP.Model()
    using DynamicPolynomials
    @polyvar x y
    X = monomials([x, y], 0:2)
    for cone in [SOSPoly(X), SDSOSPoly(X), DSOSPoly(X)]
        p = @variable(model, [1:2], cone)
        @test p[1].basis.monomials == X
        @test p[2].basis.monomials == X
        @test eltype(p) <: GramMatrix{JuMP.VariableRef}
    end
end


text:
include("maxflow.jl")

using MatrixNetworks
using SparseArrays
using LinearAlgebra
using DataStructures
using MAT

# a list of edges which forms a matching
# size of matrix edges is (matching_size, 2)
mutable struct matching
    edges::Matrix{Int64}
end

# a list of matchings whose union formes the expander
# each matching is a matrix of size [div(n, 2), 2]
mutable struct expander
    matchings::Vector{matching} end

"""
Generate an unit vector in R^n which is orthogonal to all ones vector
n: the dimension of the vector

Return r, which is a vector satisfying r' * ones(n) = zeros(n)
"""
function gen_rand_unitv(n::Int)
    # Generate one vector which is element-wise Gaussian
    r = randn((n))

    # remove the component which is on the direction of ones(n)
    vones = ones(n)
    r = r - r' * vones * vones / n

    # return the normalized vector
    return r / norm(r)
end

"""
Function for cut player to choose a bisection based on previous matchings
returned by matching player.

H = collections of matchings returned by matching player in previous rounds
r = a random unit vector which is orthogonal to all-ones vector

Return:
S = a set of vertices which induce the bisection
"""
function find_bisection(H::expander, r::Vector{Float64})
    n = length(r)
    newr = zeros(n)

    # Here we Compute M_i(M_{i-1}(...(M_2(M_1(r)))))
    # where M(r) means projecting distribution r onto matching M, i.e.
    # do a one-step lazy random walk on M starting from distribution r
    for i = 1:length(H.matchings)
        M = H.matchings[i]
        for j = 1:size(M.edges, 1)
            u = M.edges[j, 1]
            v = M.edges[j, 2]

            # for a one-step lazy random walk, w.p. 1/2 it will stay
            #                                  w.p. 1/2 it will traverse to the neighbor
            newr[u] = (r[u] + r[v]) / 2
            newr[v] = newr[v]
        end
        r = newr
    end

    # sort the vertices according to r
    # intuitively, if in \cup_i M_i two vertices are well-connected(a lot of paths between),
    # then their scores will be mixed well during the projections.
    # Therefore by splitting them according to score, it's more possible to find
    # some bisection which induces a sparse cut
    id = sortperm(r)
    S = id[1:div(n, 2)]
    return S
end

"""
Given a set S which is a subset of [N], return [N] - S
"""
function set_complement(N::Int, S::Vector{Int})
    V = collect(1:N)
    S_bar = setdiff(V, S)
    return S_bar
end

"""
Given the adjacency matrix A of a graph G and a seed set S,
return the expansion of S w.r.t. G.

Here Expansion(S, G) := E(S, V(G) - S) / min(|S|, |V(G) - S|)
"""
function compute_expansion(A::SparseMatrixCSC, S::Vector{Int64})
    N = size(A, 1)
    S_bar = set_complement(N, S)

    # find the non-zeros entries of Sparse Adjacency Matrix A indexed by [S, S_bar]
    # which is the number of edges of the cut (S, S_bar)
    _, _, E = findnz(A[S, S_bar])
    return sum(E) / min(length(S), length(S_bar))
end

"""
Given the bisection (S, S_bar) of G, we build a graph G' where we assign capacity
Cap to each edge in G and add one super source s, one super sink t. Further we
add one undirected edge from s to each vertex v in S with capacity 1, one undirected
edge from each vertex v in S_bar to t with capacity 1 as well.

A = N by N adjacnecy matrix of undirected graph G
S = a seed set of vertices, which is a subset of vertices of G
Cap = the capacity we assign to each edge of G when building G'

Return:
C = the capacity matrix of G'
s = super sourse's index
t = super sink's index

"""
function build_graph(A::SparseMatrixCSC, S::Vector{Int64}, Cap::Int64)
    N = size(A, 1)
    S_bar = set_complement(N, S)

    # Here for simplicity we let s = 1, t = N + 2 and shift the index of each vertex
    # of G by 1.
    s = 1
    t = N + 2
    svec = zeros(N)
    svec[S] .= 1
    tvec = zeros(N)
    tvec[S_bar] .= 1
    C = [spzeros(1, 1) sparse(svec') spzeros(1, 1);
        sparse(svec) Cap * A sparse(tvec);
        spzeros(1, 1) sparse(tvec') spzeros(1, 1)]
    return C, s, t
end

"""
Given a max s-t flow and a bisection (S, S_bar) of original graph G, we decompose
the s-t flow into a set of s-t paths. Further by picking the first vertex belonging
to S and the last vertex belonging to S_bar, we map each s-t path to a path between
one vertex from S and one vertex from S_bar. Thus we get a matching between
S and S_bar which is routable in G.

F = Flow matrix of the max s-t flow.
    F[u, v] > 0 means there exists F[u, v] unit flow flowing from u to v
    F[u, v] < 0 means there exists -F[u, v] unit flow flowing from v to u

s = source's index
t = sink's index
S = a set of vertices which induces the bisection

Return:
M =  a matching between S and S_bar which is routable in G

Notes: For convenience, we let s = 1 and t = N + 2 when building auxiliary graphs
for computing matching. So here it's always true that s = 1, t = N + 2.
"""
function flow_decomposition(F::SparseMatrixCSC, s::Int64, t::Int64, S::Vector{Int64})
    N = size(F, 1) - 2
    S_bar = set_complement(N, S)

    # target matching size = N // 2
    matching_size = div(N, 2)

    # store the result matching
    M = zeros(Int, matching_size, 2)

    # indicator variables for S and S_bar
    in_S = zeros(N + 2)
    in_S_bar = zeros(N + 2)
    in_S[S .+ 1] .= 1
    in_S_bar[S_bar .+ 1] .= 1

    # Compute Neighbors and degree for each vertex based on the flow matrix
    # Notice here neighbors include both in and out neighbors
    Neighs, d = ConstructAdj(F, N + 2)

    # pointers for Neighbor lists used for BFS
    Neighs_ptrs = ones(Int, N + 2)

    # queue used for BFS, each item is a pair, which is
    # (starting vertex, ending vertex) of a path
    q = Queue{Pair{Int, Int}}()

    # initialize each path with (s, s)
    for _ = 1:length(S)
        enqueue!(q, s => s)
    end

    edges = 0
    while !isempty(q)
        x = dequeue!(q)
        u = x.first
        from = x.second
        uNeighs = Neighs[u]
        @assert Neighs_ptrs[u] <= d[u]

        # find the first outgoing neighbor with positive flow remaining
        v = uNeighs[Neighs_ptrs[u]]
        while F[u, v] <= 0
            Neighs_ptrs[u] += 1
            @assert Neighs_ptrs[u] <= d[u]
            v = uNeighs[Neighs_ptrs[u]]
        end

        # since finally we want to match vertices in S with vertices in S_bar
        # we change the starting vertex at the first time we meet any vertex in S
        if from == s
            @assert in_S[v] == 1
            from = v
        end

        # use one unit flow
        F[u, v] -= 1

        # if we reach the sink, it means we have already find one s-t path,
        # store it in M
        if v == t
            edges += 1
            M[edges, 1] = from
            M[edges, 2] = u
        else
        # if not, then we push the temporary path into the queue again
            enqueue!(q, v => from)
        end
    end
    @assert edges == matching_size println(string(edges) * " " * string(matching_size))
    return M
end

"""
Given the bisection (S, S_bar), binary search capacity c such that
the matching between (S, S_bar) can be routed in G when we assign each
edge with capacity c and can't be routed when the capacity is c/2.

The max s-t flow we get when the capacity of each edge is assigned to
c can give us one G-routable matching between (S, S_bar).

The min s-t cut according to the max s-t flow we get when the capacity
of each edge is assigned to c/2 can give us one sparse cut with expansion
at most 2/c.

A = adjacency matrix of G
S = a set of vertices which induces the bisection

Return:
cut_S = a set of vertices which induces sparse cut with expansion at most 2/c
M = a matching between (S, S_bar) which is routable in G when capacity is c
"""
function flow_match(A::SparseMatrixCSC, S::Vector{Int})
    N = size(A, 1)
    lb = 0
    ub = Int(ceil(log2(N)))

    # binary search the smallest capacity (within factor 2) c such that
    # a matching between (S, S_bar) can be routed in G
    while lb < ub
        mid = div(lb + ub, 2)
        B, s, t = build_graph(A, S, 2^mid)
        F = maxflow(B, s, t)
        if F.flowvalue >= div(N, 2)
            ub = mid
        else
            lb = mid + 1
        end
    end
    # compute the matching
    # Here we need to compute the flow when the capacity of each edge is c
    B, s, t = build_graph(A, S, 2^ub)
    F = maxflow(B, s, t)
    M = flow_decomposition(F.F, s, t, S) .- 1

    # compute the sparse cut
    # Here we need to compute the flow when the capacity of each edge is c/2
    power = max(0, ub - 1)
    B, s, t = build_graph(A, S, 2^power)
    F = maxflow(B, s, t)

    # Let S be the set of vertices reachable from sourse in the residual graph
    # which induces the desired min cut
    cut_S = source_nodes(F) .- 1
    cut_S = setdiff(cut_S, 0)
    return cut_S, M
end

"""
Cut-matching game:
In round i,
Cut Player: Based on the matching M_1, M_2, ..., M_{i-1} returned by matching player
in the first (i-1) rounds. Find a new bisection (S, S_bar) for this round.

Matching Player: Given the bisection (S, S_bar) chosen by the cut player, find
a matching between S and S_bar, which is routable in G with congestion c_i.
Meanwhile output a sparse cut whose expansion is at most c_i/2.

After T = Omega(log2(N)^2) rounds, we find one expander which is routable
in G with congestion (T * max_i c_i), and one sparse cut with expansion at most
(2 / max_i c_i). This certifies that phi(G) (expansion of G) is at most 1/T *
phi(S).

We can show that the union of M_1, M_2, ..., M_T is an expander with high probability
for some T = O(log2(N)^2)

A = adjacency matrix of G

"""
function cut_matching(A::SparseMatrixCSC, T::Int)
    # store the set of matchings which forms an expander
    Expender = expander(Vector{matching}())
    N = size(A, 1)

    # store the minimum expansion
    min_expansion = N^2

    # store the set inducing the sparsest cut
    min_S = []

    # initialize the min_S, min_expansion with S = {v}, i.e. single vertex
    for i = 1:N
        S = [i]
        expansion = compute_expansion(A, S)
        if expansion < min_expansion
            min_expansion = expansion
            min_S = S
        end
    end

    # run cut-matching game for T rounds
    for _ = 1:T
        # Cut Player: give the bisection
        r = gen_rand_unitv(N)
        bisec_S = find_bisection(Expender, r)

        # Matching Player: give the matching and sparse cut
        cut_S, M = flow_match(A, bisec_S)

        # update the min_expansion and min_S with new sparse cut
        if length(cut_S) > 0 && length(cut_S) < N
            expansion = compute_expansion(A, cut_S)
            if expansion < min_expansion
                min_expansion = expansion
                min_S = cut_S
            end
        end

        #store the new matching into expander
        push!(Expender.matchings, matching(M))
    end
    return min_S, min_expansion
end

max_stars_repo_name: luotuoqingshan/Sparsest-Cut

max_issues_repo_name: luotuoqingshan/Sparsest-Cut

max_forks_repo_name: luotuoqingshan/Sparsest-Cut


text:
"""
A struct representing an Itoprocess.
### Members
* `t0` - The initial time of the `ItoProcess`
* `value` - The initial value of the `ItoProcess`
* `drift` - The drift of the `ItoProcess`.
* `stochastic` - The stochastic process.
"""
mutable struct ItoProcess{R<:Real}
  t0::R
  value::R
  drift::UnivariateFunction
  stochastic::ItoIntegral
end

"""
    evolve!(itoprocess::ItoProcess, stochastic::Real, new_time::Real)
Evolve the `ItoProcess` forward. This changes the time (`t0`) as well as the `value`.
### Inputs
* `itoprocess` - The `ItoProcess` to be evolved.
* `stochastic` - A draw from the `ItoIntegral`.
* `new_time` - The new time.
### Return
Nothing. It changes the input `ItoProcess`
"""
function evolve!(itoprocess::ItoProcess, stochastic::Real, new_time::Real)
  old_const = itoprocess.value
  realised_drift = evaluate_integral(itoprocess.drift,itoprocess.t0, new_time)
  itoprocess.value = old_const + realised_drift + stochastic
  itoprocess.t0 = new_time
end

"""
    evolve!(itoprocesses::Union{Dict{Symbol,ItoProcess{R}},Dict{Symbol,ItoProcess}}, stochastics::Union{Dict{Symbol,R},Dict{Symbol,Real}}, new_time::Real) where R<:Real
Evolve the `ItoProcess` forward. This changes the time (`t0`) as well as the `value`.
### Inputs
* `itoprocesses` - A `Dict` of `ItoProcess`es
* `stochastic` - A `Dict` of draws from the `ItoIntegral`.
* `new_time` - The new time.
### Return
Nothing. It changes the input `ItoProcess`
"""
function evolve!(itoprocesses::Union{Dict{Symbol,ItoProcess{R}},Dict{Symbol,ItoProcess}}, stochastics::Union{Dict{Symbol,R},Dict{Symbol,Real}}, new_time::Real) where R<:Real
  for k in keys(itoprocesses)
    evolve!(itoprocesses[k], stochastics[k], new_time)
  end
end

"""
    evolve_covar_and_ito_processes!(itoprocesses::Union{Dict{Symbol,ItoProcess{R}},Dict{Symbol,ItoProcess}}, covar::ForwardCovariance, new_time::Real; number_generator::NumberGenerator) where R<:Real
Evolve the `ItoProcess`es forward.
### Inputs
* `itoprocesses` - A `Dict` of `ItoProcess`es
* `covar` - The covariance matrix.
* `new_time` - The new time.
* `number_generator` - The number generator
### Return
* The updated `ItoProcess`es
* The Covariance matrix
"""
function evolve_covar_and_ito_processes!(itoprocesses::Union{Dict{Symbol,ItoProcess{R}},Dict{Symbol,ItoProcess}}, covar::ForwardCovariance, new_time::Real; number_generator::NumberGenerator) where R<:Real
  if new_time - covar.to_ < 1000*eps()
    stochastic_draws  = Dict{Symbol,Real}(covar.covariance_labels_ .=> Float64.(zeros(length(covar.covariance_labels_))))
    evolve!(itoprocesses, stochastic_draws, new_time)
    return itoprocesses, covar
  end
  new_covar = ForwardCovariance(covar, covar.to_, new_time)
  stochastic_draws = get_draws(new_covar, 1; number_generator = number_generator)[1]
  evolve!(itoprocesses, stochastic_draws, new_time)
  return itoprocesses, new_covar
end

"""
    evolve_covar_and_ito_processes!(itoprocesses::Union{Dict{Symbol,ItoProcess{R}},Dict{Symbol,ItoProcess}}, covar::SimpleCovariance, new_time::Real; number_generator::NumberGenerator) where R<:Real
Evolve the `ItoProcess`es forward.
### Inputs
* `itoprocesses` - A `Dict` of `ItoProcess`es
* `covar` - The covariance matrix.
* `new_time` - The new time.
* `number_generator` - The number generator
### Return
* The updated `ItoProcess`es
* The Covariance matrix
"""
function evolve_covar_and_ito_processes!(itoprocesses::Union{Dict{Symbol,ItoProcess{R}},Dict{Symbol,ItoProcess}}, covar::SimpleCovariance, new_time::Real; number_generator::NumberGenerator) where R<:Real
  if new_time - covar.to_ < 1000*eps()
    update!(covar, covar.from_, new_time) # Note that as newtime and to are about the same I am just trying to adjust the from time here.
    stochastic_draws  = Dict{Symbol,Real}(covar.covariance_labels_ .=> Float64.(zeros(length(covar.covariance_labels_))))
    evolve!(itoprocesses, stochastic_draws, new_time)
    return itoprocesses, covar
  end
  update!(covar, covar.to_, new_time)
  stochastic_draws = get_draws(covar, 1; number_generator = number_generator)[1]
  evolve!(itoprocesses, stochastic_draws, new_time)
  return itoprocesses, covar
end

"""
    make_ito_process_syncronous_time_series(ito_processes::Union{Dict{Symbol,ItoProcess{R}},Dict{Symbol,ItoProcess}},
                                                     covar::Union{ForwardCovariance,SimpleCovariance}, timegap::Real, total_number_of_ticks::Integer; ito_twister = MersenneTwister(2)) where R<:Real
Evolve the `ItoProcess`es forward.
### Inputs
* `itoprocesses` - A `Dict` of `ItoProcess`es
* `covar` - The covariance matrix.
* `timegap` - The time gap between ticks.
* `total_number_of_ticks` - The total number of ticks.
* `ito_twister` The `MersenneTwister` used for the RNG.
### Return
* A DataFrame with the ticks.
"""
function make_ito_process_syncronous_time_series(ito_processes::Union{Dict{Symbol,ItoProcess{T}},Dict{Symbol,ItoProcess}},
                                                     covar::Union{ForwardCovariance,SimpleCovariance}, timegap::R, total_number_of_ticks::Integer; ito_twister = MersenneTwister(2)) where R<:Real where T<:Real
    assets = collect(keys(ito_processes))
    d1 = DataFrame(Time = Array{R}([]),Name = Array{Symbol}([]),Value = Array{R}([]))
    at_time = covar.to_
    number_generator = Mersenne(ito_twister, length(assets))
    for i in 1:total_number_of_ticks
      next_tick = at_time + timegap
      ito_processes, covar = evolve_covar_and_ito_processes!(ito_processes, covar, next_tick; number_generator = number_generator)
      for stock in assets
        d2 = Dict([:Time, :Name, :Value] .=> [next_tick, stock, ito_processes[stock].value])
        d1 = push!(d1,d2)
      end
      at_time = next_tick
    end
    return d1
end

"""
    make_ito_process_non_syncronous_time_series(ito_processes::Union{Dict{Symbol,ItoProcess{R}},Dict{Symbol,ItoProcess}},
                                                covar::Union{ForwardCovariance,SimpleCovariance}, update_rates::Union{OrderedDict{Symbol,D},Dict{Symbol,D},OrderedDict{Symbol,Distribution},Dict{Symbol,Distribution}},
                                                total_number_of_ticks::Integer; timing_twister::MersenneTwister = MersenneTwister(1), ito_twister = MersenneTwister(2)) where R<:Real where D<:Distribution
Evolve the `ItoProcess`es forward.
### Inputs
* `itoprocesses` - A `Dict` of `ItoProcess`es
* `covar` - The covariance matrix.
* `update_rates` - The update rates of the exponential waiting times between ticks for each asset.
* `total_number_of_ticks` - The total number of ticks.
* `ito_twister` The `MersenneTwister` used for the RNG.
### Return
* A DataFrame with the ticks.
"""
function make_ito_process_non_syncronous_time_series(ito_processes::Union{Dict{Symbol,ItoProcess{R}},Dict{Symbol,ItoProcess}},
                                                     covar::Union{ForwardCovariance,SimpleCovariance}, update_rates::Union{OrderedDict{Symbol,D},Dict{Symbol,D},OrderedDict{Symbol,Distribution},Dict{Symbol,Distribution}},
                                                     total_number_of_ticks::Integer; timing_twister::MersenneTwister = MersenneTwister(1), ito_twister = MersenneTwister(2)) where R<:Real where D<:Distribution
    update_rates = OrderedDict(update_rates)
    assets = collect(keys(ito_processes))
    d1 = DataFrame(Time = Array{R}([]),Name = Array{Symbol}([]),Value = Array{R}([]))
    at_time = covar.to_
    number_generator = Mersenne(ito_twister, length(assets))
    for i in 1:total_number_of_ticks
      starts = vcat(rand.(Ref(timing_twister), values(update_rates), 1)...)
      next_tick = at_time + minimum(starts)

      ito_processes, covar = evolve_covar_and_ito_processes!(ito_processes, covar, next_tick; number_generator = number_generator)

      # Moving everything up to the tick.
      what_stock = collect(keys(update_rates))[findall(abs.(starts .- minimum(starts)) .< 1e-15)[1]]
      d2 = Dict([:Time, :Name, :Value] .=> [next_tick, what_stock, ito_processes[what_stock].value])
      d1 = push!(d1,d2)
      at_time = next_tick
    end
    return d1
end


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

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

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

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

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

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

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

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

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

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

"""
    meta(x::MetaObject)

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

metafree(x::T) where T = x

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


text:
include("../src/infrastructure/infrastructure.jl")

using Documenter

makedocs(
    sitename = "Invoicing",
    format = Documenter.HTML(),
    pages = [
        "Invoicing" => "index.md",
        "1 - Domain" => "chapter1.md",
        "2 - API" => "chapter2.md",
        "3 - Infrastructure" => "chapter3.md",
        "4 - Examples" => "chapter4.md"
    ]
)


text:
export generate_tensor_form,
       multivariate_tensor

using Base.Cartesian

"""
    generate_tensor_form(B::ImplicitForm{N}; algorithm::String="loop") where {N}

Compute the Bernstein expansion in tensorial form given its implicit form.

### Input

- `B`         -- Bernstein expansion in implicit form
- `algorithm` -- (optional, default: `"loop"`) the implementation used for this
                 calculation, valid options are:

    * `"loop"`
    * `"kron"`
"""
function generate_tensor_form(B::ImplicitForm{N}; algorithm::String="loop") where {N}
    # unpack fields
    v = VectorOfArray(B.array) # FIXME : move to type?
    generate_tensor_form(v, algorithm=algorithm)
end

function generate_tensor_form(v::VectorOfArray{N, D, T}; algorithm::String="loop") where {N, D, T}
    if algorithm == "loop"
        # preallocate output array
        A = Array{N, length(v)}(Tuple(length(vi) for vi in v))
        _loop_tensor!(A, v)
    elseif algorithm == "kron"
        # preallocate output array
        A = Array{N, 1}() # FIXME: array size matters?
        _kron_tensor!(A, v)
    else
        error("the tensorial expansion algorithm $algorithm is unknown")
    end
    return A
end

"""
```
    _kron_tensor!(B::ImplicitForm{N}) where {N}
```

Compute the Bernstein expansion in tensorial form given its implicit form, using
the Kronecker power.

### Input

- `B` -- Bernstein expansion in implicit form

### Output

A one-dimensional array of the same numeric type as the given implicit form.

### Algorithm

This implementation uses Julia's Kronecker product `kron` function in a loop.
"""
function _kron_tensor!(A::AbstractArray{T, N}, v::VectorOfArray) where {T, N}
    A = v[1] # FIXME: this is not using the given A
    @inbounds for i in 2:length(v)
        A = kron(A, v[i])
    end
end

"""
```
    _loop_tensor!(A::AbstractArray{T, N}, v::VectorOfArray) where {T, N}
```

Compute the Bernstein expansion in tensorial form given its implicit form, using
a nested loop.

### Input

- `B` -- Bernstein expansion in implicit form

### Output

A multi-dimensional array containing the coefficients.

### Algorithm

This implementation uses Julia's `Base.Cartesian` to generate a set of nested
loops to compute the element ``A[i_1, …, i_n]``.
"""
@generated function _loop_tensor!(A::AbstractArray{T, N}, v::VectorOfArray) where {T, N}
    quote
    cache = Vector{Float64}($N-2)
    @inbounds for $(Symbol(:i_, N)) in eachindex(v[$N])
                  # FIXME: split me in several lines!
                  @nloops $(N-1) i i -> eachindex(v[i]) d -> d == 1 ? nothing : d == $N-1 ? cache[d-1] = v[i_{d}, d] * v[i_{d+1}, d+1] : cache[d-1] = cache[d] * v[i_{d}, d] begin (@nref $N A i) = v[i_1, 1] * cache[1] end
              end
    end
end

"""
    multivariate_tensor(k::Vector{Int64}, l::Vector{Int64},
                          low::Vector{T}, high::Vector{T};
                          algorithm::String="loop")

Compute the Bernstein coefficients of a multivariate monomial in the tensor form.

### Input

- `k`    -- vector of degrees for each monomial
- `l`    -- vector of degrees of the Bernstein polynomial for each monomial
- `low`  -- the lower bounds of the interval where the Bernstein coefficients are computed
- `high` -- the upper bounds of the interval the Bernstein coefficients are computed
- `algorithm` -- (optional, default: `"loop"`) the implementation used for this
                 calculation, valid options are:

    * `"loop"`
    * `"kron"`
"""
function multivariate_tensor(k::Vector{Int64}, l::Vector{Int64},
                             low::Vector{T}, high::Vector{T};
                             algorithm::String="loop") where {T}

    n = length(low)
    berncoeffs = VectorOfArray([univariate(k[i], l[i], low[i], high[i]) for i in 1:n])
    generate_tensor_form(berncoeffs, algorithm=algorithm)
end


text:
# Julia wrapper for header: /usr/local/ffmpeg/2.4/include/libavdevice/avdevice.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    avdevice_version,
    avdevice_configuration,
    avdevice_license,
    avdevice_register_all,
    av_input_audio_device_next,
    av_input_video_device_next,
    av_output_audio_device_next,
    av_output_video_device_next,
    avdevice_app_to_dev_control_message,
    avdevice_dev_to_app_control_message,
    avdevice_capabilities_create,
    avdevice_capabilities_free,
    avdevice_list_devices,
    avdevice_free_list_devices

function avdevice_version()
    ccall((:avdevice_version,libavdevice),Uint32,())
end

function avdevice_configuration()
    ccall((:avdevice_configuration,libavdevice),Ptr{Uint8},())
end

function avdevice_license()
    ccall((:avdevice_license,libavdevice),Ptr{Uint8},())
end

function avdevice_register_all()
    ccall((:avdevice_register_all,libavdevice),Void,())
end

function av_input_audio_device_next(d)
    ccall((:av_input_audio_device_next,libavdevice),Ptr{AVInputFormat},(Ptr{AVInputFormat},),d)
end

function av_input_video_device_next(d)
    ccall((:av_input_video_device_next,libavdevice),Ptr{AVInputFormat},(Ptr{AVInputFormat},),d)
end

function av_output_audio_device_next(d)
    ccall((:av_output_audio_device_next,libavdevice),Ptr{AVOutputFormat},(Ptr{AVOutputFormat},),d)
end

function av_output_video_device_next(d)
    ccall((:av_output_video_device_next,libavdevice),Ptr{AVOutputFormat},(Ptr{AVOutputFormat},),d)
end

function avdevice_app_to_dev_control_message(s,_type::AVAppToDevMessageType,data,data_size::Csize_t)
    ccall((:avdevice_app_to_dev_control_message,libavdevice),Cint,(Ptr{AVFormatContext},AVAppToDevMessageType,Ptr{Void},Csize_t),s,_type,data,data_size)
end

function avdevice_dev_to_app_control_message(s,_type::AVDevToAppMessageType,data,data_size::Csize_t)
    ccall((:avdevice_dev_to_app_control_message,libavdevice),Cint,(Ptr{AVFormatContext},AVDevToAppMessageType,Ptr{Void},Csize_t),s,_type,data,data_size)
end

function avdevice_capabilities_create(caps,s,device_options)
    ccall((:avdevice_capabilities_create,libavdevice),Cint,(Ptr{Ptr{AVDeviceCapabilitiesQuery}},Ptr{AVFormatContext},Ptr{Ptr{AVDictionary}}),caps,s,device_options)
end

function avdevice_capabilities_free(caps,s)
    ccall((:avdevice_capabilities_free,libavdevice),Void,(Ptr{Ptr{AVDeviceCapabilitiesQuery}},Ptr{AVFormatContext}),caps,s)
end

function avdevice_list_devices(s,device_list)
    ccall((:avdevice_list_devices,libavdevice),Cint,(Ptr{AVFormatContext},Ptr{Ptr{AVDeviceInfoList}}),s,device_list)
end

function avdevice_free_list_devices(device_list)
    ccall((:avdevice_free_list_devices,libavdevice),Void,(Ptr{Ptr{AVDeviceInfoList}},),device_list)
end


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   Modia3D.ContactDetectionMPR (Modia3D/contactDetection/ContactDetectionMPR/_module.jl)
#

# Collision detection algorithm based on the MPR algorithm

mutable struct SupportPoint
  p::SVector{3,Float64}    # support point
  n::SVector{3,Float64}    # support normal unit vector
  a::SVector{3,Float64}    # point on shapeA
  b::SVector{3,Float64}    # point on shapeB
  function SupportPoint(p::SVector{3,Float64},n::SVector{3,Float64},a::SVector{3,Float64},b::SVector{3,Float64})
    new(p,n,a,b)
  end
end

function getSupportPoint(shapeA::Modia3D.Composition.Object3D, shapeB::Composition.Object3D, n::SVector{3,Float64}; scale::Float64=1.0)
  a = Modia3D.supportPoint(shapeA.data.geo, shapeA.r_abs, shapeA.R_abs, n)
  b = Modia3D.supportPoint(shapeB.data.geo, shapeB.r_abs, shapeB.R_abs, -n)
  return SupportPoint((a-b).*scale,n,a,b)
end

function barycentric(r1::SupportPoint,r2::SupportPoint,r3::SupportPoint,r4::SupportPoint)
  r21 = r2.p - r1.p
  r31 = r3.p - r1.p
  r41 = r4.p - r1.p

  # Determine normal to triangle (and test that r4 is on the triangle plane)
  n = cross(r21,r31)
  e = normalize(n)
  #@assert(abs(dot(e,r41)) <= 1e-10)

  # Determine normals to r21 and e, as well as r31 and e
  n21 = cross(r21,e)
  n31 = cross(r31,e)

  # Determine barycentric coordinates b21, b31:
  #   n21*(r41 = b21*r21 + b31*r31)    ->  n21*r21 = 0
  #   n31*(r41 = b21*r21 + b31*r31)    ->  n31*r31 = 0
  b21 = dot(n31,r41) / dot(n31,r21)
  b31 = dot(n21,r41) / dot(n21,r31)

  r4.a = r1.a + b21*(r2.a - r1.a) + b31*(r3.a - r1.a)
  r4.b = r1.b + b21*(r2.b - r1.b) + b31*(r3.b - r1.b)

  ### only works if r4.n goes through zero
  #=
  b1 = dot(cross(r2.p,r3.p),r4.n)
  b2 = dot(cross(r3.p,r1.p),r4.n)
  b3 = dot(cross(r1.p,r2.p),r4.n)
  b_ges = b1 + b2 + b3

  r4.a .= (b1*r1.a + b2*r2.a + b3*r3.a) / b_ges
  r4.b .= (b1*r1.b + b2*r2.b + b3*r3.b) / b_ges
  =#
end

###########      Phase 1, Minkowski Portal Refinement      ###################

getCentroid(obj::Composition.Object3D, centroidObj::SVector{3,Float64}) = (obj.r_abs + obj.R_abs'*centroidObj)

# checks if centers of shapeA and shapeB are overlapping
# belongs to construction of r0
function checkCentersOfShapesOverlapp(r0::SupportPoint, neps::Float64, shapeA::Composition.Object3D, shapeB::Composition.Object3D)
  if norm(r0.p) <= neps
    error("MPR: Too large penetration (prerequisite of MPR violated). Centers are overlapping. Look at shapeA = ", shapeA, " shapeB = ", shapeB)
  end
end

function checkIfShapesArePlanar(r0::SupportPoint,r1::SupportPoint,r2::SupportPoint,n2::SVector{3,Float64}, neps::Float64,
                                shapeA::Composition.Object3D,shapeB::Composition.Object3D)
  # r3 is in the direction of plane normal that contains triangle r0-r1-r2
  n3 = cross(r1.p-r0.p, r2.p-r0.p)
  # the triangle r0-r1-r2 has degenerated into a line segment
  if norm(n3) <= neps
    # change search direction for r2
    # because we are still interested in distances if shapes are not intersecting
    n2 = -n2
    r2 = getSupportPoint(shapeA, shapeB, n2)
    if abs(dot((r2.p-r1.p),n2)) <= neps
      # Shape is purely planar. Computing the shortest distance for a planar shape
      # requires an MPR 2D algorithm (using lines instead of triangles as portals).
      # However, this is not implemented and therefore the shortest distance cannot be computed
      error("MPR: Shapes are planar and MPR2D is not supported. abs(dot((r2.p-r1.p),n2)). Look at shapeA = ", shapeA, " shapeB = ", shapeB)
    end
    # new normal to the triangle plane (r0-r1-r2_new)
    n3 = cross(r1.p-r0.p, r2.p-r0.p)   # |n3| > 0 guaranteed, due to construction
  end

  if dot(n3,r0.p) >= neps
    n3 = -n3
  end

  # check if portal triangle r1-r2-r3 has degenerated into a line segment <--> points r1,r2,r3 are on the same line
  r3 = getSupportPoint(shapeA, shapeB, Basics.normalizeVector(n3))
  n3b = cross(r2.p-r1.p, r3.p-r1.p)
  if norm(n3b) <= neps
    # change search direction for r3
    r3 = getSupportPoint(shapeA, shapeB, -r3.n)
    if abs(dot((r3.p-r1.p),r3.n)) <= neps
      # Shape is purely planar. Computing the shortest distance for a planar shape
      # requires an MPR 2D algorithm (using lines instead of triangles as portals).
      # However, this is not implemented and therefore the shortest distance cannot be computed
      error("MPR: Shapes are planar and MPR2D is not supported. r1, r2, r3 are on the same ray. abs(dot((r3.p-r1.p),r3.n)) <= neps. Look at shapeA = ", shapeA, " shapeB = ", shapeB)
    end
  end

  return (r2, r3, n2, n3)
end

###########      Phase 2, Minkowski Portal Refinement      ###################
# loop around to "ensure" the tetrahedron r0,r1,r2 and r3 encloses the origin
# stimmt so nicht wirklich, muss ich nochmal nachlesen!!!
# Der Ursprung muss nicht enthalten sein!!!
function tetrahedronEncloseOrigin(r0::SupportPoint,r1::SupportPoint,r2::SupportPoint,r3::SupportPoint,
                                  neps::Float64, niter_max::Int64,
                                  shapeA::Composition.Object3D,shapeB::Composition.Object3D, scale::Float64)
  aux = SVector(0.0, 0.0, 0.0)
  success = false
  for i in 1:niter_max
    aux = cross(r1.p-r0.p,r3.p-r0.p)
    if dot(aux,r0.n) < -neps
       r2 = r3
       r3 = getSupportPoint(shapeA,shapeB,Basics.normalizeVector(aux), scale=scale)
       continue
    end
    aux = cross(r3.p-r0.p,r2.p-r0.p)
    if dot(aux,r0.n) < -neps
       r1 = r3
       r3 = getSupportPoint(shapeA,shapeB,Basics.normalizeVector(aux), scale=scale)
       continue
    end
    success = true
    break
  end
  if success != true
    error("MPR: Max. number of iterations is reached in phase2. Look at shapeA = ", shapeA, " shapeB = ", shapeB)
  end
  # pointInTetrahedron(r0.p,r1.p,r2.p,r3.p,SVector(0.0, 0.0, 0.0))
  return (r1, r2, r3)
end

###########      Phase 3, Minkowski Portal Refinement      ###################
# construction of r4
function constructR4(r0::SupportPoint,r1::SupportPoint,r2::SupportPoint,r3::SupportPoint,
                     neps::Float64, shapeA::Composition.Object3D,shapeB::Composition.Object3D, scale::Float64)
  n4 = cross(r2.p-r1.p, r3.p-r1.p)
  if norm(n4) <= neps
    r3 = getSupportPoint(shapeA, shapeB, -r3.n, scale=scale) # change search direction
    if abs(dot((r3.p-r1.p),r3.n)) <= neps
      # Shape is purely planar. Computing the shortest distance for a planar shape
      # requires an MPR 2D algorithm (using lines instead of triangles as portals).
      # However, this is not implemented and therefore the shortest distance cannot be computed
      error("MPR: Shapes are planar and MPR2D is not supported. abs(dot((r3.p-r1.p),r3.n)). Look at shapeA = ", shapeA, " shapeB = ", shapeB)
    end
    n4 = cross(r2.p-r1.p, r3.p-r1.p)   # |n4| > 0 guaranteed, due to construction
  end
  if dot(n4,r0.p) >= neps
    n4 = -n4
  end
  r4 = getSupportPoint(shapeA, shapeB, Basics.normalizeVector(n4), scale=scale)

  return (r3, r4, n4)
end

isNextPortal(r0,r1,r2,r4) = dot( cross( (r2.p-r0.p), (r4.p-r0.p) ), r0.n ) <= 0.0 &&
                            dot( cross( (r4.p-r0.p), (r1.p-r0.p) ), r0.n ) <= 0.0

# Construction of three baby tetrahedrons and decide which one will replace
# the current tetrahedron r1r2 - r2r3 - r3r1 - current tetrahedron edges
# leads to 3 new tetrahedrons:
#   r1r2r4
#   r2r3r4
#   r3r1r4
function createBabyTetrahedrons(r0::SupportPoint,r1::SupportPoint,r2::SupportPoint,r3::SupportPoint,r4::SupportPoint,
                                shapeA::Composition.Object3D,shapeB::Composition.Object3D)
  if isNextPortal(r0,r1,r2,r4)
     r3 = r4
  elseif isNextPortal(r0,r2,r3,r4)
     r1 = r4
  elseif isNextPortal(r0,r3,r1,r4)
     r2 = r4
  else
     # Compute the closest distance to the portal and return it:
     (r4.p, distance) = signedDistanceToPortal(r0.p,r1.p,r2.p,r3.p)
     barycentric(r1,r2,r3,r4)
     error("MPR: Numerical issues with distance computation between shapeA = ", shapeA, " and shapeB = ", shapeB,". Used distance = ",distance)
     return (distance,r4.a,r4.b,r4.n,r1.a,r1.b,r2.a,r2.b,r3.a,r3.b)
  end
  return (r1,r2,r3)
end

function scaleVector(scale, ri)
  return ri.*scale
end

function skalarization(r0,r1,r2,r3)
  x = maximum([abs.(r0.p) abs.(r1.p) abs.(r2.p) abs.(r3.p)])
  scale = 1/x
  r0.p = scaleVector(scale, r0.p)
  r1.p = scaleVector(scale, r1.p)
  r2.p = scaleVector(scale, r2.p)
  r3.p = scaleVector(scale, r3.p)
  return (r0,r1,r2,r3,scale)
end

# MPR - Minkowski Portal Refinement algorithm
# construction of points r0 is in the interior of Minkowski Difference and points r1,r2,r3,r4 are on the boundary of Minkowski Difference
# Phase 1
#   Construct a tetrahedron such one point is in the interior r0 and three points r1,r2,r3 are on the boundary of the Minkowski Difference
#   Construction of r0 and of initial portal triangle points r1, r2, r3
#   Phase 1.1: construction of r0
#   Phase 1.2: construction of initial r1
#     Termination Condition 1
#   Phase 1.3: construction of initial r2
#   Phase 1.4: construction of initial r3
# Phase 2
#   make sure the origin ray passes through portal triangle r1-r2-r3, otherwise it is not sure that the portal's normal points in the right direction
#   if the origin ray passes through the portal, --> the normal of all other triangles of the tetrahedron point away from the origin ray
#   replace one of the points in the portal with a new one until the ray from r0 passes through the portal
#   loop around to "ensure" the tetrahedron r0,r1,r2 and r3 encloses the origin
# Phase 3
#   search iteratively in the direction of the portal's normal
#   construct a new portal, such the origin passes still the portal, and it is closer to the surface of the Minkowski Difference
#   than the old portal
#   Phase 3.1: construct r4,
#   Phase 3.2: check if r4 is close to the origin
#     Termination Condition 2
#     Termination Condition 3
#   Phase 3.3: construct baby tetrahedrons with r1,r2,r3,r4 and create a new portal
function mpr(ch::Composition.ContactDetectionMPR_handler, shapeA::Composition.Object3D, shapeB::Modia3D.Composition.Object3D, geoA::Modia3D.AbstractSolidGeometry, geoB::Modia3D.AbstractSolidGeometry)
  tol_rel = ch.tol_rel
  tol_rel = 10e-6
  niter_max = ch.niter_max
  neps = ch.neps

  ###########      Phase 1, Minkowski Portal Refinement      ###################
  # Construction of r0 and initial portal triangle points r1, r2, r3

  ### Phase 1.1: construction of r0 ###
  # r0 is a point inside the Minkowski Difference (centroidA-centroidB)
  # the direction of the origin ray r0 is -r0.p
  centroidA = getCentroid(shapeA, Modia3D.centroid(shapeA.data.geo))
  centroidB = getCentroid(shapeB, Modia3D.centroid(shapeB.data.geo))
  r0 = SupportPoint(centroidA-centroidB, -(centroidA-centroidB), SVector{3,Float64}(0.0,0.0,0.0), SVector{3,Float64}(0.0,0.0,0.0))
  # check if centers of shapes are overlapping
  checkCentersOfShapesOverlapp(r0, neps, shapeA, shapeB)

  ### Phase 1.2: construction of initial r1 ###
  # r1 is the farthest point in the direction to the origin
  # first portal point should point in the direction of the origin ray (-r0.p)
  # therefore choose search direction -r0.p
  r1 = getSupportPoint(shapeA, shapeB, Basics.normalizeVector(r0.n))

  ### Phase 1.3: construction of initial r2 ###
  n2 = cross(r0.n, r1.p)
  n2abs = norm(n2)
  ## TERMINATION CONDITION 1 ##
  if n2abs <= neps
    # r0 || r1, r0 ist parallel zu r1
    # centers of shapes (r0.p and r1.p) are on the same line segment
    # the origin of Minkowski Difference is on the line from r0 with direction r1-r0
    # e.g. any collision/or distance between two spheres
    #println("TC 1")
    distance = dot(r1.p,normalize(r0.p))
    return (distance,r1.a,r1.b,r1.n, nothing, nothing, nothing, nothing, nothing, nothing)
  else
    # normalize n2
    n2 = n2/n2abs
  end
  r2 = getSupportPoint(shapeA, shapeB, n2)

  ### Phase 1.4: construction of initial r3 ###
  # r3 is in the direction of plane normal that contains triangle r0-r1-r2
  (r2, r3, n2, n3) = checkIfShapesArePlanar(r0, r1, r2, n2, neps, shapeA, shapeB)

  (r0,r1,r2,r3,scale) = skalarization(r0,r1,r2,r3)

  ###########      Phase 2, Minkowski Portal Refinement      ###################
  # loop around to "ensure" the tetrahedron r0,r1,r2 and r3 encloses the origin
  (r1,r2,r3) = tetrahedronEncloseOrigin(r0,r1,r2,r3,neps,niter_max,shapeA,shapeB, scale)
  # doesRayIntersectPortal(r1.p,r2.p,r3.p, r0.p,neps) # Portal.A, Portal.B, Portal.C, point, neps

  ###########      Phase 3, Minkowski Portal Refinement      ###################
  for i in 1:niter_max
    ### Phase 3.1: construct r4 ###
    # Find support point using the tetrahedron face
    (r3,r4,n4) = constructR4(r0,r1,r2,r3,neps,shapeA,shapeB, scale)

    ### Phase 3.2: check if r4 is close to the origin ###
    # check if the new point r4 is already on the plane of the triangle r1,r2,r3,
    # we're already as close as we can get to the origin

    ## TERMINATION CONDITION 2 ##
    if norm(cross(r4.p,r4.n)) < tol_rel
      #println("TC 2")
      #if !analyzeFinalPortal(r1.p, r2.p, r3.p, r4.p, neps)
        # error("shapeA = ", shapeA, " shapeB = ", shapeB)
      #end
      distance = -dot(r4.n, r4.p)
      barycentric(r1,r2,r3,r4)
      # println("r4.a = ", r4.a, ", r4.b = ", r4.b, " ,r1.a = ", r1.a, " ,r1.b = ", r1.b," ,r2.a = ", r2.a, " ,r2.b = ",r2.b," ,r3.a = ",r3.a," ,r3.b = ",r3.b)
      (sphereWithOthers,contactPointA, contactPointB) = handleSphereWithOtherShapes(shapeA, shapeB, r4.n, distance)
      if sphereWithOthers
        return (distance,contactPointA,contactPointB,r4.n,r1.a,r1.b,r2.a,r2.b,r3.a,r3.b)
      else
        return (distance,r4.a,r4.b,r4.n,r1.a,r1.b,r2.a,r2.b,r3.a,r3.b)
      end

    ## TERMINATION CONDITION 3 ##
    elseif abs(dot(r4.p-r1.p, r4.n)) < tol_rel
      #println("TC 3")

      #doesRayIntersectPortal(r1.p,r2.p,r3.p, r4.p,neps)
      #println("r1.p = ", r1.p , " r2.p = ", r2.p ," r3.p = ", r3.p)
      #println("r4.p = ", r4.p)
      #println(" ")
      #if !analyzeFinalPortal(r1.p, r2.p, r3.p, r4.p, neps)
        # error("shapeA = ", shapeA, " shapeB = ", shapeB)
      #end
      (r4.p, distance) = signedDistanceToPortal(r0.p,r1.p,r2.p,r3.p)
      barycentric(r1,r2,r3,r4)
      # println("r4.n = ", r4.n)

      (sphereWithOthers,contactPointA, contactPointB) = handleSphereWithOtherShapes(shapeA, shapeB, r4.n, distance)
      if sphereWithOthers
        #rContact = (contactPointA + contactPointB)/2
        #Printf.@printf("TC3: contact normal = [%.3g,%.3g,%.3g], contact position = [%.3g,%.3g,%.3g]\n",
        #             r4.n[1],r4.n[2],r4.n[3],rContact[1],rContact[2],rContact[3])

        return (distance,contactPointA,contactPointB,r4.n,r1.a,r1.b,r2.a,r2.b,r3.a,r3.b)
      else
        return (distance,r4.a,r4.b,r4.n,r1.a,r1.b,r2.a,r2.b,r3.a,r3.b)
      end
      #println("r4.a = ", r4.a, ", r4.b = ", r4.b, " ,r1.a = ", r1.a, " ,r1.b = ", r1.b," ,r2.a = ", r2.a, " ,r2.b = ",r2.b," ,r3.a = ",r3.a," ,r3.b = ",r3.b)

    end

    #### Phase 3.3: construct baby tetrahedrons with r1,r2,r3,r4 and create a new portal ###
    # Construction of three baby tetrahedrons
    (r1,r2,r3) = createBabyTetrahedrons(r0,r1,r2,r3,r4,shapeA,shapeB)

#=
    r0RayOnPortal = isr0RayOnPortal(r0.p,r1.p,r2.p,r3.p)
    if r0RayOnPortal[1] != true
       println("r0RayOnPortal = ", r0RayOnPortal, ": Wrong baby-tetrahedron selected. " , i)
    end
    =#

  end
  error("MPR: Should never get here! Computation failed. Look at shapeA = ", ModiaMath.instanceName(shapeA),
       " shapeB = ",  ModiaMath.instanceName(shapeB))
end

function mpr(ch::Composition.ContactDetectionMPR_handler, shapeA::Composition.Object3D,
              shapeB::Modia3D.Composition.Object3D,
              geoA::Modia3D.Solids.SolidSphere, geoB::Modia3D.Solids.SolidSphere)
  neps = ch.neps
  radiusA = shapeA.data.geo.Dx*0.5
  radiusB = shapeB.data.geo.Dx*0.5
  centroidSphereA = getCentroid(shapeA, Modia3D.centroid(shapeA.data.geo))
  centroidSphereB = getCentroid(shapeB, Modia3D.centroid(shapeB.data.geo))
  n = centroidSphereB - centroidSphereA
  distanceCentroids = norm(n)
  if distanceCentroids <= neps
    error("Centers of two spheres are overlapping. SphereA = ", shapeA, " sphereB = ", shapeB)
  else
    normal = n/distanceCentroids
  end
  distance = distanceCentroids - radiusA - radiusB
  contactPointShapeA = centroidSphereA + normal*radiusA
  contactPointShapeB = centroidSphereB - normal*radiusB
  return (distance,contactPointShapeA,contactPointShapeB,normal, nothing, nothing, nothing, nothing, nothing, nothing)
end

function handleSphereWithOtherShapes(shapeA::Composition.Object3D, shapeB::Composition.Object3D, normal::SVector{3,Float64}, distance::Float64)
  # normal is pointing from shapeA to shapeB
  if typeof(shapeA.data.geo) != Modia3D.Solids.SolidSphere && typeof(shapeB.data.geo) != Modia3D.Solids.SolidSphere
    return (false, nothing, nothing)
  elseif typeof(shapeA.data.geo) == Modia3D.Solids.SolidSphere
    centroidSphere = getCentroid(shapeA, Modia3D.centroid(shapeA.data.geo))
    radius = shapeA.data.geo.Dx*0.5
    contactPointSphere = centroidSphere + radius*normal
    contactPointOtherShape = contactPointSphere + distance*normal
    return (true, contactPointSphere, contactPointOtherShape)
  elseif typeof(shapeB.data.geo) == Modia3D.Solids.SolidSphere
    normal = -normal
    centroidSphere = getCentroid(shapeB, Modia3D.centroid(shapeB.data.geo))
    radius = shapeB.data.geo.Dx*0.5
    contactPointSphere = centroidSphere + radius*normal
    contactPointOtherShape = contactPointSphere + distance*normal # distance is negative (otherwise direction of normal must be changed)
    return (true, contactPointOtherShape, contactPointSphere)
  end
end


text:
function tempIntegrate!(x, tspan, ps; avoidErr = false)
    h0 = 0.5 / ps.crp.TU
    maxh = 10.0 / ps.crp.TU

    # Set utype
    cSc = ps.sp.isp*9.81*ps.crp.TU / (ps.crp.LU*1000.0)
    λv = norm(view(x,11:13))
    S = computeS(x, λv, cSc)
    if S > ps.ϵ
        ps.utype = 0
    elseif S < -ps.ϵ
        ps.utype = 2
    else
        ps.utype = 1
    end

    # Initialize integrator
    STM = true
    if length(x) == 210
        integ = ATO_Integraion.initODE((dz,z,p,t)->cr3bpEomIndirectWithSTM!(dz,z,p,t),
                                       x,
                                       tspan,
                                       ATO_Integraion.Vern9ME();
                                       ps = ps,
                                       RelTol = 1e-14,
                                       AbsTol = 1e-14,
                                       maxh = maxh)
        x0sm = Vector{Float64}(undef, 14)
        x0sm .= view(x,1:14)
        integ_sw = ATO_Integraion.initODE((dz,z,p,t)->cr3bpEomIndirect!(dz,z,p,t),
                                          x0sm,
                                          tspan,
                                          ATO_Integraion.Vern9ME();
                                          ps = ps,
                                          RelTol = 1e-14,
                                          AbsTol = 1e-14,
                                          maxh = maxh)
    elseif length(x) == 14
        STM = false
        integ = ATO_Integraion.initOde((dz,z,p,t)->cr3bpEomIndirect!(dz,z,p,t),
                                       x,
                                       tspan,
                                       ATO_Integraion.Vern9ME();
                                       ps = ps,
                                       RelTol = 1e-14,
                                       AbsTol = 1e-14,
                                       maxh = maxh)
    end

    # Begin integration loop
    done = false
    switching = false
    reducing = false
    sw_cnt = 0
    while !done

        # Incriment Switching Counter
        sw_cnt += 1

        # Take Step
        ATO_Integration.step!(integ)

        # Check for Switching
        switching, reducing = Check_4_Switching(integ.ukp1, integ.ps)

        # Detect Switching Time
        if switching && sw_cnt > ss_iters # Switching
            if STM
                Update_IntegSW!(integ_sw, integ)
                tsw = Detect_t_Switch!(integ_sw, mn_iters, mr_iters, s_tol, avoid_err)
            else
                tsw = Detect_t_Switch!(integ, mn_iters, mr_iters, s_tol, avoid_err)
            end
            if tsw < 0 || tsw == integ.tk
                switching = false
                reducing = false

                # Make sure utype is corect, reduce step size, and restep
                _Correct_Utype!(integ)
                integ.h = h0
                ATO_Integration.step!(integ)
            end
        end

        # Handle Switching or Reducing
        if switching && sw_cnt > ss_iters # Switching

            # Reset Switching Counter
            sw_cnt = 0

            # Detect Switching Time
            #if STM
            #    Update_IntegSW!(integ_sw, integ)
            #    tsw = Detect_t_Switch!(integ_sw, mn_iters, mr_iters, s_tol, avoid_err)
            #else
            #    tsw = Detect_t_Switch!(integ, mn_iters, mr_iters, s_tol, avoid_err)
            #end

            # Step to tsw
            ATO_Integration.step!(integ, tsw - integ.tk)

            # Update Utype
            _Update_Utype!(integ)

            # Accept Step
            ATO_Integration.accept_step!(integ)

            # Update Step Size
            integ.h = integ.h * reduce_fact
            ATO_Integration.init_h!(integ)

            # Propagate STM over Discontinuity
            if params.ϵ == 0.0 && STM
                _Prop_STM_Discont!(integ)
            end

        elseif reducing # Reducing
            integ.h = integ.h * reduce_fact

        else # Continue to Next Step W/O Doing Anything
            ATO_Integration.accept_step!(integ)
        end

        # Handle Closeness to tf
        if integ.tk >= integ.tf
            done = true
        elseif integ.tk + integ.h >= integ.tf
            integ.h = integ.tf - integ.tk
        end
    end

    # Update x before returning
    @inbounds for i in 1:length(x)
        x[i] = integ.uk[i]
    end
end

function Check_4_Switching(x::AbstractArray, params)

    # Set Bools
    switching = false
    reducing = false

    # Compute Skp1
    Skp1 = compute_S(x, params)

    # Check for Switching
    if params.utype == 2 # Thrust is On
        if Skp1 < -params.ϵ # No Switching
            switching = false
        elseif Skp1 > -params.ϵ && Skp1 < params.ϵ # Switching to Medium
            switching = true
        else # Reducing Step Size or Switching
            if params.ϵ == 0.0
                switching = true
            else
                reducing = true
            end
        end
    elseif params.utype == 0 # Thrust is Off

        if Skp1 > params.ϵ # No Switching
            switching = false
        elseif Skp1 > -params.ϵ && Skp1 < params.ϵ # Switching to Medium
            switching = true
        else # Reducing Step Size or Switching
            if params.ϵ == 0.0
                switching = true
            else
                reducing = true
            end
        end
    else # Medium Thrust
        if !(Skp1 > -params.ϵ && Skp1 < params.ϵ) # Switching On or Off
            switching = true
        end
    end

    return switching, reducing
end

function Update_IntegSW!(integ_sw::T,integ::U) where {T,U<:ATO_Integration.ODEInt}
    integ_sw.tk = integ.tk
    integ_sw.tkp1 = integ.tkp1
    integ_sw.h = integ.h
    integ_sw.qold = integ_sw.qold
    integ_sw.ps.utype = integ.ps.utype
    @inbounds for i in 1:14
        integ_sw.uk[i] = integ.uk[i]
        integ_sw.ukp1[i] = integ.ukp1[i]
    end
end

function Detect_t_Switch!(integ::ATO_Integration.ODEInt, mn_iters, mr_iters, s_tol, avoid_err)

    # Initialize tsw
    tkp1 = integ.tkp1
    tsw = integ.tk

    # Initialize Switching Function
    Sjp1 = 0.0
    if integ.ps.utype == 1
        Sjp1 = compute_S(integ.ukp1, integ.ps)
    end
    _fsw(tsw) = _Switch_Function!(tsw, Sjp1, integ)

    # Begin Newton's Itterations
    stop = false
    succ = false
    niters = 0
    while !stop
        # Incriment Iteration Counter
        niters += 1

        # Compute Switching Function
        fsw = _fsw(tsw)

        # Compute 1/Sdot(tsw)
        if tsw > integ.tk
            dSinv = _Compute_dSinv(integ.ukp1, integ.ps)
        else
            dSinv = _Compute_dSinv(integ.uk, integ.ps)
        end

        # Update tsw
        tsw -= fsw*dSinv

        # Check for proper convergence
        if tsw < integ.tk || tsw > integ.tkp1 || niters > mn_iters
            stop = true
        elseif abs(fsw) < s_tol
            stop = true
            succ = true
        end
    end

    # Begin Regula-Falsi is Newtons Failed
    if !succ
        # Initialize
        stop = false
        niters = 0
        a = integ.tk
        b = tkp1
        fa = _fsw(a)
        fb = _fsw(b)

        # Check for root bracketing
        if fa * fb > 0
            if avoid_err
                #if abs(fa) < abs(fb)
                #    tsw = a
                #else
                #    tsw = b
                #end
                tsw = -1.0
                stop = true
            else
                throw(ErrorException("Root not bracketed for Regula-Falsi switching detection!"))
            end
        end

        # Regula Falsi Loop
        while !stop

            # Incriment Iteration Counter
            niters += 1

            # Compute c
            c = b - (fb*(b - a))/(fb - fa)
            fc = _fsw(c)

            # Update
            if fa*fc < 0
                b = c
                fb = fc
            else
                a = c
                fa = fc
            end

            # Update tsw and check for proper convergence
            if abs(fa) < abs(fb)
                fsw = fa
                tsw = a
            else
                fsw = fb
                tsw = b
            end
            if abs(fsw) < s_tol
                stop = true
                succ = true
            elseif niters > mr_iters
                stop = true
            end
        end
    end

    return tsw
end

function _Switch_Function!(tsw, Sjp1, integ::ATO_Integration.ODEInt)

    # Determine Switching Direction
    ϵ = -integ.ps.ϵ
    if integ.ps.utype == 2 || (integ.ps.utype == 1 && Sjp1 < -integ.ps.ϵ)
        ϵ *= -1.0
    end

    # Take Step
    h = tsw - integ.tk
    if h > 0
        ATO_Integration.step!(integ, h)
        Stsw = compute_S(integ.ukp1, integ.ps)
    else
        Stsw = compute_S(integ.uk, integ.ps)
    end

    # Compute Switching Function
    return Stsw + ϵ
end

function _Compute_dSinv(x::AbstractVector, params)
    λ_v = sqrt(x[11]^2 + x[12]^2 + x[13]^2)
    c_nsc = params.sp.Isp * 9.81 * params.eom.TU / (params.eom.LU * 1000)

    # Compute 1/Sdot(tsw)
    dSinv = (x[8] - 2*x[12])*x[11] + (x[9] + 2*x[11])*x[12] + x[10]*x[13]
    dSinv = λ_v*x[7]/(c_nsc*dSinv)

    return dSinv
end

function _Update_Utype!(integ::ATO_Integration.ODEInt)
    if (integ.ps.utype == 2 || integ.ps.utype == 0) && integ.ps.ϵ != 0.0
        integ.ps.utype = 1
    elseif integ.ps.ϵ != 0.0
        Stsw = compute_S(integ.ukp1, integ.ps)
        if abs(Stsw + integ.ps.ϵ) < abs(Stsw - integ.ps.ϵ)
            integ.ps.utype = 2
        else
            integ.ps.utype = 0
        end
    else
        if integ.ps.utype == 2
            integ.ps.utype = 0
        else
            integ.ps.utype = 2
        end
    end
end

function _Correct_Utype!(integ::ATO_Integration.ODEInt)
    S = compute_S(integ.uk, integ.ps)
    if S < -integ.ps.ϵ
        integ.ps.utype = 2
    elseif S > integ.ps.ϵ
        integ.ps.utype = 0
    else
        integ.ps.utype = 1
    end
end

function _Prop_STM_Discont!(integ::ATO_Integration.ODEInt)

    # Revert Utype of Parameters Copy to Previous Value
    integ.ps.utype == 2 ? integ.ps.utype = 0 : integ.ps.utype = 2

    # Initialize Requirements
    temp = MVector{14, Float64}(undef)
    dy⁻ = MVector{14,Float64}(undef)
    dy⁺ = MVector{14,Float64}(undef)

    # Compute dy⁻
    @inbounds for i in 1:14; temp[i] = integ.uk[i]; end
    CR3BPIMF_EOM!(dy⁻, temp, integ.ps, 0.0)

    # Compute dy⁺
    integ.ps.utype == 2 ? integ.ps.utype = 0 : integ.ps.utype = 2
    @inbounds for i in 1:14; temp[i] = integ.uk[i]; end
    CR3BPIMF_EOM!(dy⁺, temp, integ.ps, 0.0)

    # Compute ∂S and dSinv
    λ_v = sqrt(integ.uk[11]^2 + integ.uk[12]^2 + integ.uk[13]^2)
    c_nsc = integ.ps.sp.Isp * 9.81 * integ.ps.eom.TU / (integ.ps.eom.LU * 1000)

    @views @. temp .= [0.0;0.0;0.0;0.0;0.0;0.0;
                       λ_v*c_nsc / (integ.uk[7]^2);
                       0.0;0.0;0.0;
                       integ.uk[11:13] * (-c_nsc / (λ_v*integ.uk[7]));
                       -1.0]
    dSinv = (integ.uk[8] - 2*integ.uk[12])*integ.uk[11] +
            (integ.uk[9] + 2*integ.uk[11])*integ.uk[12] +
            (integ.uk[10])*integ.uk[13]
    dSinv = λ_v*integ.uk[7] / (c_nsc*dSinv)

    # Compute Φ and Ψ
    uSTM = reshape(view(integ.uk, 15:210), 14, 14)
    Φ = SizedMatrix{14,14,Float64}(undef)
    Ψ = SizedMatrix{14,14,Float64}(undef)
    Threads.@threads for row in 1:14
        for col in 1:14
            Φ[row, col] = uSTM[row, col]
            Ψ[row, col] = I[row,col] + (dy⁺[row] - dy⁻[row]) * temp[col]*dSinv
        end
    end

    # Multiplying Ψ and Φ
    Threads.@threads for row in 1:14
        for col in 1:14
            sum = 0.0
            for k in 1:14
                sum += Ψ[row, k] * Φ[k, col]
            end
            uSTM[row, col] = sum
        end
    end
end

function compute_S(x::AbstractVector, params)
    λ_v = sqrt(x[11]^2 + x[12]^2 + x[13]^2)
    c_nsc   = params.sp.Isp * 9.81 * params.eom.TU / (params.eom.LU * 1000)

    return (-λ_v * c_nsc / x[7] - x[14]  + 1.0)
end