text:
using LogExpFunctions, Test

@testset "xlogx & xlogy" begin
    @test iszero(xlogx(0))
    @test xlogx(2) ≈ 2.0 * log(2.0)

    @test iszero(xlogy(0, 1))
    @test xlogy(2, 3) ≈ 2.0 * log(3.0)
    @test xlogy(2, 3.0) ≈ 2.0 * log(3.0) # promoted
end

@testset "logistic & logit" begin
    @test logistic(2)        ≈ 1.0 / (1.0 + exp(-2.0))
    @test iszero(logit(0.5))
    @test logit(logistic(2)) ≈ 2.0
end

@testset "log1psq" begin
    @test iszero(log1psq(0.0))
    @test log1psq(1.0) ≈ log1p(1.0)
    @test log1psq(2.0) ≈ log1p(4.0)
    @test log1psq(Float16(0.1)) ≈ log1p(0.01)
end

# log1pexp, log1mexp, log2mexp & logexpm1

@testset "log1pexp" begin
    @test log1pexp(2.0)    ≈ log(1.0 + exp(2.0))
    @test log1pexp(-2.0)   ≈ log(1.0 + exp(-2.0))
    @test log1pexp(10000)  ≈ 10000.0
    @test log1pexp(-10000) ≈ 0.0

    @test log1pexp(2f0)      ≈ log(1f0 + exp(2f0))
    @test log1pexp(-2f0)     ≈ log(1f0 + exp(-2f0))
    @test log1pexp(10000f0)  ≈ 10000f0
    @test log1pexp(-10000f0) ≈ 0f0
end

@testset "log1mexp" begin
    @test log1mexp(-1.0)  ≈ log1p(- exp(-1.0))
    @test log1mexp(-10.0) ≈ log1p(- exp(-10.0))
end

@testset "log2mexp" begin
    @test log2mexp(0.0)  ≈ 0.0
    @test log2mexp(-1.0) ≈ log(2.0 - exp(-1.0))
end

@testset "logexpm1" begin
    @test logexpm1(2.0)            ≈  log(exp(2.0) - 1.0)
    @test logexpm1(log1pexp(2.0))  ≈  2.0
    @test logexpm1(log1pexp(-2.0)) ≈ -2.0

    @test logexpm1(2f0)            ≈  log(exp(2f0) - 1f0)
    @test logexpm1(log1pexp(2f0))  ≈  2f0
    @test logexpm1(log1pexp(-2f0)) ≈ -2f0
end

@testset "log1pmx" begin
    @test iszero(log1pmx(0.0))
    for x in [-0.65, -0.5, -0.3, -0.1, 0.4, 1.0, 2.0]
        z = BigFloat(x)
        @test log1pmx(x) ≈ Float64(log(1 + x) - x)
    end
end

@testset "logmxp1" begin
    @test iszero(logmxp1(1.0))
    for x in [0.1, 0.35, 0.5, 2.0, 3.0]
        z = BigFloat(x)
        @test logmxp1(x) ≈ Float64(log(z) - z + 1)
    end
end

@testset "logsumexp" begin
    @test logaddexp(2.0, 3.0)     ≈ log(exp(2.0) + exp(3.0))
    @test logaddexp(10002, 10003) ≈ 10000 + logaddexp(2.0, 3.0)
    @test logaddexp(2, 3.0)     ≈ log(exp(2.0) + exp(3.0)) # promotion

    @test logsumexp([1.0, 2.0, 3.0])          ≈ 3.40760596444438
    @test logsumexp((1.0, 2.0, 3.0))          ≈ 3.40760596444438
    @test logsumexp([1.0, 2.0, 3.0] .+ 1000.) ≈ 1003.40760596444438

    let cases = [([-Inf, -Inf], -Inf),   # correct handling of all -Inf
                 ([-Inf, -Inf32], -Inf), # promotion
                 ([-Inf32, -Inf32], -Inf32), # Float32
                 ([-Inf, Inf], Inf),
                 ([-Inf, 9.0], 9.0),
                 ([Inf, 9.0], Inf),
                 ([NaN, 9.0], NaN),  # NaN propagation
                 ([NaN, Inf], NaN),  # NaN propagation
                 ([NaN, -Inf], NaN), # NaN propagation
                 ([0, 0], log(2.0))] # non-float arguments
        for (arguments, result) in cases
            @test logaddexp(arguments...) ≡ result
            @test logsumexp(arguments) ≡ result
        end
    end
end

@testset "softmax" begin
    x = [1.0, 2.0, 3.0]
    r = exp.(x) ./ sum(exp.(x))
    @test softmax([1.0, 2.0, 3.0]) ≈ r
    softmax!(x)
    @test x ≈ r
end


text:
include("radix_system_utils.jl")
include("qubo/qubo.jl")

ext: jl


text:
# queue

export ZeCommandQueue, synchronize

mutable struct ZeCommandQueue
    handle::ze_command_queue_handle_t

    context::ZeContext
    device::ZeDevice
    ordinal::Int

    function ZeCommandQueue(ctx::ZeContext, dev::ZeDevice, ordinal=1, index=1;
                            flags=0,
                            mode::ze_command_queue_mode_t=ZE_COMMAND_QUEUE_MODE_DEFAULT,
                            priority::ze_command_queue_priority_t=ZE_COMMAND_QUEUE_PRIORITY_NORMAL)
        desc_ref = Ref(ze_command_queue_desc_t(;
            ordinal=ordinal-1, index=index-1, flags, mode, priority
        ))
        handle_ref = Ref{ze_command_queue_handle_t}()
        zeCommandQueueCreate(ctx, dev, desc_ref, handle_ref)
        obj = new(handle_ref[], ctx, dev, ordinal)
        finalizer(obj) do obj
            zeCommandQueueDestroy(obj)
        end
        obj
    end
end

Base.unsafe_convert(::Type{ze_command_queue_handle_t}, queue::ZeCommandQueue) = queue.handle

synchronize(queue::ZeCommandQueue, timeout::Number=typemax(UInt32)) =
    zeCommandQueueSynchronize(queue, timeout)

## groups

export command_queue_groups, compute_groups

struct ZeCommandQueueGroups
    device::ZeDevice
end

command_queue_groups(dev::ZeDevice) = ZeCommandQueueGroups(dev)

Base.eltype(::ZeCommandQueueGroups) = ZeCommandQueueGroup

function Base.iterate(groups::ZeCommandQueueGroups, i=1)
    i >= length(groups) + 1 ? nothing : (ZeCommandQueueGroup(groups, i), i+1)
end

Base.length(groups::ZeCommandQueueGroups) = length(properties(groups))

function properties(groups::ZeCommandQueueGroups)
    count_ref = Ref{UInt32}(0)
    zeDeviceGetCommandQueueGroupProperties(groups.device, count_ref, C_NULL)

    all_props = fill(ze_command_queue_group_properties_t(), count_ref[])
    zeDeviceGetCommandQueueGroupProperties(groups.device, count_ref, all_props)

    return [(flags=props.flags,
             maxMemoryFillPatternSize=UInt(props.maxMemoryFillPatternSize),
             numQueues=Int(props.numQueues),
             ) for props in all_props[1:count_ref[]]]
end

Base.IteratorSize(::ZeCommandQueueGroups) = Base.HasLength()

struct ZeCommandQueueGroup
    groups::ZeCommandQueueGroups
    ordinal::Int
end

properties(group::ZeCommandQueueGroup) = properties(group.groups)[group.ordinal]

# short-hands
compute_groups(dev::ZeDevice) = filter(collect(command_queue_groups(dev))) do group
    properties(group).flags & oneL0.ZE_COMMAND_QUEUE_GROUP_PROPERTY_FLAG_COMPUTE != 0
end


text:
function is_equivalent_by_translation(cryst::Crystal, b1::BondRaw, b2::BondRaw)
    # Displacements between the two bonds
    D1 = b2.ri - b1.ri
    D2 = b2.rj - b1.rj
    # Round components of D1 to nearest integers
    n = round.(D1, RoundNearest)
    # If both n ≈ D1 and n ≈ D2, then the bonds are equivalent by translation
    return norm(n - D1) < cryst.symprec && norm(n - D2) < cryst.symprec
end

# Generate list of all symmetries that transform b2 into b1, along with parity
function symmetries_between_bonds(cryst::Crystal, b1::BondRaw, b2::BondRaw)
    # Fail early if two bonds describe different real-space distances
    # (dimensionless error tolerance is measured relative to the minimum lattice
    # constant ℓ)
    if b1 != b2
        ℓ = minimum(norm, eachcol(cryst.lat_vecs))
        d1 = distance(cryst, b1) / ℓ
        d2 = distance(cryst, b2) / ℓ
        if abs(d1-d2) > cryst.symprec
            return Tuple{SymOp, Bool}[]
        end
    end

    ret = Tuple{SymOp, Bool}[]
    for s in cryst.symops
        b2′ = transform(s, b2)
        if is_equivalent_by_translation(cryst, b1, b2′)
            push!(ret, (s, true))
        elseif is_equivalent_by_translation(cryst, b1, reverse(b2′))
            push!(ret, (s, false))
        end
    end
    return ret
end

"Is there a symmetry operation that transforms `b1` into either `b2` or its reverse?"
function is_related_by_symmetry(cryst::Crystal, b1::BondRaw, b2::BondRaw)
    return !isempty(symmetries_between_bonds(cryst::Crystal, b1::BondRaw, b2::BondRaw))
end

function is_related_by_symmetry(cryst::Crystal, b1::Bond, b2::Bond)
    return is_related_by_symmetry(cryst, BondRaw(cryst, b1), BondRaw(cryst, b2))
end

"Returns all bonds in `cryst` for which `bond.i == i`"
function all_bonds_for_atom(cryst::Crystal, i::Int, max_dist; min_dist=0.0)
    # be a little generous with the minimum and maximum distances
    ℓ = minimum(norm, eachcol(cryst.lat_vecs))
    max_dist += 4 * cryst.symprec * ℓ
    min_dist -= 4 * cryst.symprec * ℓ

    # columns are the reciprocal vectors
    recip_vecs = 2π * inv(cryst.lat_vecs)'

    # box_lengths[i] represents the perpendicular distance between two parallel
    # boundary planes spanned by lattice vectors a_j and a_k (where indices j
    # and k differ from i)
    box_lengths = [a⋅b/norm(b) for (a,b) = zip(eachcol(cryst.lat_vecs), eachcol(recip_vecs))]
    n_max = round.(Int, max_dist ./ box_lengths, RoundUp)

    bonds = Bond[]

    # loop over neighboring cells
    for n1 in -n_max[1]:n_max[1]
        for n2 in -n_max[2]:n_max[2]
            for n3 in -n_max[3]:n_max[3]
                n = SVector(n1, n2, n3)

                # loop over all atoms within neighboring cell
                for j in eachindex(cryst.positions)
                    b = Bond(i, j, n)
                    if min_dist <= distance(cryst, b) <= max_dist
                        push!(bonds, b)
                    end
                end
            end
        end
    end

    return bonds
end

# Calculate score for a bond. Lower would be preferred.
function _score_bond(cryst::Crystal, b::Bond)
    # Favor bonds with fewer nonzero elements in basis matrices J
    Js = basis_for_symmetry_allowed_couplings(cryst, b)
    nnz = [count(abs.(J) .> 1e-12) for J in Js]
    score = Float64(sum(nnz))

    # Favor bonds with smaller unit cell displacements. Positive
    # displacements are slightly favored over negative displacements.
    # Displacements in x are slightly favored over y, etc.
    score += norm((b.n .- 0.1) .* [0.07, 0.08, 0.09])

    # Favor smaller indices and indices where i < j
    score += 1e-2 * (b.i + b.j) + 1e-2 * (b.i < b.j ? -1 : +1)

    return score
end

"Indices of the unique elements in `a`, ordered by their first appearance."
function unique_indices(a)
    map(x->x[1], unique(x->x[2], enumerate(a)))
end

"""    reference_bonds(cryst::Crystal, max_dist)

Returns a full list of bonds, one for each symmetry equivalence class, up to
distance `max_dist`. The reference bond `b` for each equivalence class is
selected according to a scoring system that prioritizes simplification of the
elements in `basis_for_symmetry_allowed_couplings(cryst, b)`."""
reference_bonds(cryst::Crystal, max_dist) = reference_bonds(cryst, convert(Float64, max_dist))

function reference_bonds(cryst::Crystal, max_dist::Float64)
    # Bonds, one for each equivalence class
    ref_bonds = Bond[]
    for i in unique_indices(cryst.classes)
        for b in all_bonds_for_atom(cryst, i, max_dist)
            if !any(is_related_by_symmetry(cryst, b, b′) for b′ in ref_bonds)
                push!(ref_bonds, b)
            end
        end
    end

    # Sort by distance
    sort!(ref_bonds, by=b->distance(cryst, b))

    # Replace each canonical bond by the "best" equivalent bond
    return map(ref_bonds) do rb
        # Find full set of symmetry equivalent bonds
        equiv_bonds = unique([transform(cryst, s, rb) for s in cryst.symops])
        # Take the bond with lowest score
        scores = [_score_bond(cryst, b) for b in equiv_bonds]
        return equiv_bonds[findmin(scores)[2]::Int]
        # TODO: In Julia 1.7 the above two lines become
        #     return argmin(b -> _score_bond(cryst, b), equiv_bonds)
    end
end

"""
    all_symmetry_related_bonds_for_atom(cryst::Crystal, i::Int, b::Bond)

Returns a list of all bonds that start at atom `i`, and that are symmetry
equivalent to bond `b` or its reverse.
"""
function all_symmetry_related_bonds_for_atom(cryst::Crystal, i::Int, b_ref::Bond)
    bs = Bond[]
    dist = distance(cryst, b_ref)
    for b in all_bonds_for_atom(cryst, i, dist; min_dist=dist)
        if is_related_by_symmetry(cryst, b_ref, b)
            push!(bs, b)
        end
    end
    return bs
end

"""
    all_symmetry_related_bonds(cryst::Crystal, b::Bond)

Returns a list of all bonds that are symmetry-equivalent to bond `b` or its
reverse.
"""
function all_symmetry_related_bonds(cryst::Crystal, b_ref::Bond)
    bs = Bond[]
    for i in eachindex(cryst.positions)
        append!(bs, all_symmetry_related_bonds_for_atom(cryst, i, b_ref))
    end
    bs
end

"""    coordination_number(cryst::Crystal, i::Int, b::Bond)

Returns the number times that atom `i` participates in a bond equivalent to `b`.
In other words, the count of bonds that begin at atom `i` and that are
symmetry-equivalent to `b` or its reverse.

Defined as `length(all_symmetry_related_bonds_for_atom(cryst, i, b))`.
"""
function coordination_number(cryst::Crystal, i::Int, b::Bond)
    return length(all_symmetry_related_bonds_for_atom(cryst, i, b))
end


text:
"""
    TemporalPointProcess{L<:Real, M}

The common supertype of all temporal point processes (i.e. point processes on the real line) with mark type `M`.
"""
abstract type TemporalPointProcess{M} <: AbstractPointProcess{Float64,M} end

"""
    BoundedTemporalPointProcess{L<:Real, M}

Store a temporal point process `P` with pre-defined start and end times.
"""
struct BoundedTemporalPointProcess{M,P<:TemporalPointProcess{M}} <: TemporalPointProcess{M}
    pp::P
    tmin::Float64
    tmax::Float64
end

## Intensity functions

@doc raw"""
    intensity(pp, h, t, m)

Compute the conditional intensity for a temporal point process `pp` applied to history `h` and event `(t, m)`.

The conditional intensity function ``\lambda(t, m | h)`` quantifies the instantaneous risk  of an event with mark `m` occurring at time `t`[^Rasmussen_2018].

[^Rasmussen_2018]: Rasmussen, J. G. (2018), “Lecture Notes: Temporal Point Processes and the Conditional Intensity Function,” arXiv:1806.00221 [stat].
"""
function intensity(pp::TemporalPointProcess, h::TemporalHistory, t, m) end

@doc raw"""
    mark_distribution(pp, h, t)

Compute the distribution of marks for a temporal point process `pp` knowing that an event takes place at time `t` after history `h`.
"""
function mark_distribution(pp::TemporalPointProcess, h::TemporalHistory, t) end

@doc raw"""
    ground_intensity(pp, h, t)

Compute the ground intensity for a temporal point process `pp` applied to history `h` at time `t`.

The ground intensity quantifies the instantaneous risk of an event with any mark occurring at time `t`[^Rasmussen_2018]. It can be expressed as
```math
    \lambda_g(t|h) = \sum_{m \in \mathcal} \lambda(t, m|h).
```
"""
function ground_intensity(pp::TemporalPointProcess, h::TemporalHistory, t) end

@doc raw"""
    ground_intensity_bound(pp, θ, h, t)

Compute a local upper bound on the ground intensity for a temporal point process `pp` applied to history `h` at time `t`[^Rasmussen_2018].

Return a tuple of the form $(B, L)$ satisfying
```math
    \forall u \in [t, t+L), \quad \lambda_g(t|h) \leq B
```
"""
function ground_intensity_bound(pp::TemporalPointProcess, h::TemporalHistory, t) end

## Simulation with Ogata's algorithm

MeasureTheory.sampletype(::TemporalPointProcess{M}) where {M} = TemporalHistory{M}

"""
    rand(rng, pp, tmin, tmax)

Simulate a temporal point process `pp` on interval `[tmin, tmax)` using Ogata's algorithm[^Ogata_1981].

[^Ogata_1981]: Ogata, Y. (1981), “On Lewis’ simulation method for point processes,” IEEE Transactions on Information Theory, 27, 23–31. https://doi.org/10.1109/TIT.1981.1056305.
"""
function Base.rand(rng::AbstractRNG, pp::TemporalPointProcess{M}, tmin, tmax) where {M}
    h = TemporalHistory(Float64[], M[], tmin, tmax)
    t = tmin
    while t < tmax
        B, L = ground_intensity_bound(pp, h, t + eps(t))
        T = B > 0 ? rand(rng, Dists.Exponential(1 / B)) : Inf
        if T > L
            t = t + L
        elseif T <= L
            U = rand(Dists.Uniform(0, 1))
            if U < ground_intensity(pp, h, t + T) / B
                m = rand(rng, mark_distribution(pp, h, t + T))
                if t + T < tmax
                    push!(h, t + T, m)
                end
            end
            t = t + T
        end
    end
    return h
end

Base.rand(rng::AbstractRNG, tpp::BoundedTemporalPointProcess) =
    rand(rng, tpp.pp, tpp.tmin, tpp.tmax)

Base.rand(tpp::TemporalPointProcess, args...) = rand(GLOBAL_RNG, tpp, args...)

## Learning

@doc raw"""
    integrated_ground_intensity(pp, h[, t])

Compute the integrated ground intensity (or compensator) $\Lambda(t|h)$ for a temporal point process `pp` applied to history `h`:
```math
    \Lambda(h) = \int \lambda_g(t|h) \mathrm{d}t.
```

The default method uses [Quadrature.jl](https://github.com/SciML/Quadrature.jl) for numerical integration, but it should be reimplemented for specific processes if explicit integration is feasible.
"""
function integrated_ground_intensity(
    pp::TemporalPointProcess,
    h::TemporalHistory,
    t = max_time(h),
)
    par = [pp]
    f = (t, par) -> ground_intensity(par[1], h, t)
    prob = QuadratureProblem(f, min_time(h), t, par)
    sol = solve(prob, HCubatureJL())
    return sol[1]
end

@doc raw"""
    logpdf(pp, h)

Compute the log probability density function for a temporal point process `pp` applied to history `h`:
```math
    \log f(h) = \sum_i \log \lambda(t_i | h_i) - \Lambda(h)
```

The default method uses a loop over events combined with [`integrated_ground_intensity`](@ref), but it should be reimplemented for specific processes if faster computation is possible.
"""
function logpdf(pp::TemporalPointProcess, h::TemporalHistory)
    l = -integrated_ground_intensity(pp, h)
    for (t, m) in zip(event_times(h), event_marks(h))
        l += log(intensity(pp, h, t, m))
    end
    return l
end

@doc raw"""
    fit(pp0, h)

Compute the optimal parameter for a temporal point process of type `typeof(pp0)` on history `h` using maximum likelihood:
```math
    \hat{\theta} = \mathrm{argmax} \{f_{\theta}(h): \theta \in \Theta\}
```

The default method uses [GalacticOptim.jl](https://github.com/SciML/GalacticOptim.jl) for numerical optimization, but it should be reimplemented for specific processes if explicit maximization is feasible.
"""
function Dists.fit(
    pp_init::PP,
    h::TemporalHistory{M};
    adtype = SciMLBase.NoAD(),
    alg = Optim.NelderMead(),
) where {M,PP<:TemporalPointProcess}
    trans = build_transform(pp_init)
    θ_init = inverse(trans, ntfromstruct(pp_init))
    par = [nothing]
    f = (θ, par) -> -logpdf(PP(transform(trans, θ)), h)
    obj = OptimizationFunction(f, adtype)
    prob = OptimizationProblem(obj, θ_init, par)
    sol = solve(prob, alg)
    θ_opt = sol.u
    pp_opt = PP(transform(trans, θ_opt))
    return pp_opt
end

@doc raw"""
    check_residuals(pp, h)

Check whether the point process `pp` is a good fit for history `h` by applying Ogata's time rescaling method: if $(t_i)_i$ is a temporal point process with intensity $\lambda(t)$, then $(\Lambda(t_i))_i$ is a standard temporal Poisson process.
"""
function check_residuals(pp::TemporalPointProcess, h::TemporalHistory)
    Λ(t) = integrated_ground_intensity(pp, h, t)
    h_rescaled = time_change(h, Λ)
    qqplot_interevent_times(h_rescaled)
end


text:
export convexhull, convexhull!, conichull

const HAny{T} = Union{HRep{T}, HRepElement{T}}
const VAny{T} = Union{VRep{T}, VRepElement{T}}

"""
    intersect(P1::HRep, P2::HRep)

Takes the intersection of `P1` and `P2` ``\\{\\, x : x \\in P_1, x \\in P_2 \\,\\}``.
It is very efficient between two H-representations or between two polyhedron for which the H-representation has already been computed.
However, if `P1` (resp. `P2`) is a polyhedron for which the H-representation has not been computed yet, it will trigger a representation conversion which is costly.
See the [Polyhedral Computation FAQ](http://www.cs.mcgill.ca/~fukuda/soft/polyfaq/node25.html) for a discussion on this operation.

The type of the result will be chosen closer to the type of `P1`. For instance, if `P1` is a polyhedron (resp. H-representation) and `P2` is a H-representation (resp. polyhedron), `intersect(P1, P2)` will be a polyhedron (resp. H-representation).
If `P1` and `P2` are both polyhedra (resp. H-representation), the resulting polyhedron type (resp. H-representation type) will be computed according to the type of `P1`.
The coefficient type however, will be promoted as required taking both the coefficient type of `P1` and `P2` into account.
"""
function Base.intersect(p::HRep...)
    T = promote_coefficient_type(p)
    similar(p, hmap((i, x) -> convert(similar_type(typeof(x), T), x), FullDim(p[1]), T, p...)...)
end
Base.intersect(p::HRep, el::HRepElement) = p ∩ intersect(el)
Base.intersect(el::HRepElement, p::HRep) = p ∩ el

Base.intersect(hps::HyperPlane...) = hrep([hps...])
Base.intersect(hss::HalfSpace...) = hrep([hss...])
Base.intersect(h1::HyperPlane{T}, h2::HalfSpace{T}) where {T} = hrep([h1], [h2])
Base.intersect(h1::HalfSpace{T}, h2::HyperPlane{T}) where {T} = h2 ∩ h1
Base.intersect(p1::HAny{T}, p2::HAny{T}, ps::HAny{T}...) where {T} = intersect(p1 ∩ p2, ps...)
function Base.intersect(p::HAny...)
    T = promote_type(coefficient_type.(p)...)
    f(p) = convert(similar_type(typeof(p), T), p)
    intersect(f.(p)...)
end

"""
    intersect!(p::HRep, h::Union{HRepresentation, HRepElement})

Same as [`intersect`](@ref) except that `p` is modified to be equal to the intersection.
"""
function Base.intersect!(p::HRep, ::Union{HRepresentation, HRepElement})
    error("intersect! not implemented for $(typeof(p)). It probably does not support in-place modification, try `intersect` (without the `!`) instead.")
end
function Base.intersect!(p::Polyhedron, h::Union{HRepresentation, HRepElement})
    resethrep!(p, hrep(p) ∩ h)
end

"""
    convexhull(P1::VRep, P2::VRep)

Takes the convex hull of `P1` and `P2` ``\\{\\, \\lambda x + (1-\\lambda) y : x \\in P_1, y \\in P_2 \\,\\}``.
It is very efficient between two V-representations or between two polyhedron for which the V-representation has already been computed.
However, if `P1` (resp. `P2`) is a polyhedron for which the V-representation has not been computed yet, it will trigger a representation conversion which is costly.

The type of the result will be chosen closer to the type of `P1`. For instance, if `P1` is a polyhedron (resp. V-representation) and `P2` is a V-representation (resp. polyhedron), `convexhull(P1, P2)` will be a polyhedron (resp. V-representation).
If `P1` and `P2` are both polyhedra (resp. V-representation), the resulting polyhedron type (resp. V-representation type) will be computed according to the type of `P1`.
The coefficient type however, will be promoted as required taking both the coefficient type of `P1` and `P2` into account.
"""
function convexhull(p::VRep...)
    T = promote_coefficient_type(p)
    similar(p, vmap((i, x) -> convert(similar_type(typeof(x), T), x), FullDim(p[1]), T, p...)...)
end
convexhull(p::VRep, el::VRepElement) = convexhull(p, convexhull(el))
convexhull(el::VRepElement, p::VRep) = convexhull(p, el)

convexhull(ps::AbstractVector...) = vrep([ps...])
convexhull(ls::Line...) = vrep([ls...])
convexhull(rs::Ray...) = vrep([rs...])
convexhull(p::AbstractVector{T}, r::Union{Line{T}, Ray{T}}) where {T} = vrep([p], [r])
convexhull(r::Union{Line{T}, Ray{T}}, p::AbstractVector{T}) where {T} = convexhull(p, r)
convexhull(l::Line{T}, r::Ray{T}) where {T} = vrep([l], [r])
convexhull(r::Ray{T}, l::Line{T}) where {T} = convexhull(l, r)
convexhull(p1::VAny{T}, p2::VAny{T}, ps::VAny{T}...) where {T} = convexhull(convexhull(p1, p2), ps...)
function convexhull(p::VAny...)
    T = promote_type(coefficient_type.(p)...)
    f(p) = convert(similar_type(typeof(p), T), p)
    convexhull(f.(p)...)
end

"""
    convexhull!(p1::VRep, p2::VRep)

Same as [`convexhull`](@ref) except that `p1` is modified to be equal to the convex hull.
"""
function convexhull!(p::VRep, ine::VRepresentation)
    error("convexhull! not implemented for $(typeof(p)). It probably does not support in-place modification, try `convexhull` (without the `!`) instead.")
end
function convexhull!(p::Polyhedron, v::VRepresentation)
    resetvrep!(p, convexhull(vrep(p), v))
end

# conify: same than conichull except that conify(::VRepElement) returns a VRepElement and not a V-representation
conify(v::VRep) = vrep(lines(v), [collect(rays(v)); Ray.(collect(points(v)))])
conify(v::VCone) = v
conify(p::AbstractVector) = Ray(p)
conify(r::Union{Line, Ray}) = r

conichull(p...) = convexhull(conify.(p)...)

function sumpoints(::FullDim, ::Type{T}, p1, p2) where {T}
    _tout(p) = convert(similar_type(typeof(p), T), p)
    ps = [_tout(po1 + po2) for po1 in points(p1) for po2 in points(p2)]
    tuple(ps)
end
sumpoints(::FullDim, ::Type{T}, p1::Rep, p2::VCone) where {T} = change_coefficient_type.(preps(p1), T)
sumpoints(::FullDim, ::Type{T}, p1::VCone, p2::Rep) where {T} = change_coefficient_type.(preps(p2), T)

function Base.:+(p1::VRep{T1}, p2::VRep{T2}) where {T1, T2}
    T = typeof(zero(T1) + zero(T2))
    similar((p1, p2), FullDim(p1), T, sumpoints(FullDim(p1), T, p1, p2)..., change_coefficient_type.(rreps(p1, p2), T)...)
end
Base.:+(p::Rep, el::Union{Line, Ray}) = p + vrep([el])
Base.:+(el::Union{Line, Ray}, p::Rep) = p + el

# p1 has priority
function usehrep(p1::Polyhedron, p2::Polyhedron)
    hrepiscomputed(p1) && (!vrepiscomputed(p1) || hrepiscomputed(p2))
end

function hcartesianproduct(p1::HRep, p2::HRep)
    d = sum_fulldim(FullDim(p1), FullDim(p2))
    T = promote_coefficient_type((p1, p2))
    f = (i, x) -> zeropad(x, i == 1 ? FullDim(p2) : neg_fulldim(FullDim(p1)))
    similar((p1, p2), d, T, hmap(f, d, T, p1, p2)...)
end
function vcartesianproduct(p1::VRep, p2::VRep)
    d = sum_fulldim(FullDim(p1), FullDim(p2))
    T = promote_coefficient_type((p1, p2))
    # Always type of first arg
    f1 = (i, x) -> zeropad(x, FullDim(p2))
    f2 = (i, x) -> zeropad(x, neg_fulldim(FullDim(p1)))
    q1 = similar(p1, d, T, vmap(f1, d, T, p1)...)
    q2 = similar(p2, d, T, vmap(f2, d, T, p2)...)
    q1 + q2
end
cartesianproduct(p1::HRep, p2::HRep) = hcartesianproduct(p1, p2)
cartesianproduct(p1::VRep, p2::VRep) = vcartesianproduct(p1, p2)

function cartesianproduct(p1::Polyhedron, p2::Polyhedron)
    if usehrep(p1, p2)
        hcartesianproduct(p1, p2)
    else
        vcartesianproduct(p1, p2)
    end
end

"""
    *(p1::Rep, p2::Rep)

Cartesian product between the polyhedra `p1` and `p2`.
"""
Base.:(*)(p1::Rep, p2::Rep) = cartesianproduct(p1, p2)

"""
    \\(P::AbstractMatrix, p::HRep)

Transform the polyhedron represented by ``p`` into ``P^{-1} p`` by transforming each halfspace ``\\langle a, x \\rangle \\le \\beta`` into ``\\langle P^\\top a, x \\rangle \\le \\beta`` and each hyperplane ``\\langle a, x \\rangle = \\beta`` into ``\\langle P^\\top a, x \\rangle = \\beta``.
"""
Base.:(\)(P::AbstractMatrix, rep::HRep) = rep / P'

"""
    /(p::HRep, P::AbstractMatrix)

Transform the polyhedron represented by ``p`` into ``P^{-T} p`` by transforming each halfspace ``\\langle a, x \\rangle \\le \\beta`` into ``\\langle P a, x \\rangle \\le \\beta`` and each hyperplane ``\\langle a, x \\rangle = \\beta`` into ``\\langle P a, x \\rangle = \\beta``.
"""
function Base.:(/)(p::HRep{Tin}, P::AbstractMatrix) where {Tin}
    if size(P, 2) != fulldim(p)
        throw(DimensionMismatch("The number of rows of P must match the dimension of the H-representation"))
    end
    f = (i, h) -> h / P
    # FIXME For a matrix P of StaticArrays, `d` should be type stable
    d = size(P, 1)
    T = _promote_type(Tin, eltype(P))
    similar(p, d, T, hmap(f, d, T, p)...)
end

"""
    *(P::AbstractMatrix, p::VRep)

Transform the polyhedron represented by ``p`` into ``P p`` by transforming each element of the V-representation (points, symmetric points, rays and lines) `x` into ``P x``.
"""
function Base.:(*)(P::AbstractMatrix, p::VRep{Tin}) where {Tin}
    if size(P, 2) != fulldim(p)
        throw(DimensionMismatch("The number of rows of P must match the dimension of the V-representation"))
    end
    f = (i, v) -> P * v
    # For a matrix P of StaticArrays, `d` should be type stable
    d = size(P, 1)
    T = _promote_type(Tin, eltype(P))
    similar(p, d, T, vmap(f, d, T, p)...)
end


text:
using NBodyIPs, JuLIP, Test, IPFitting
using JuLIP.Potentials: evaluate_d
using IPFitting: Dat, LsqDB
using NBodyIPs: BondLengthDesc
const Lsq = IPFitting.Lsq

# generate random data
function generate_data(species, L, rmax, N, calc)
   data = Dat[]
   for n = 1:N
      at = bulk(species; cubic=true, pbc=true) * L
      rattle!(at, rand() * rmax)
      E = energy(calc, at)
      F = forces(calc, at)
      push!(data, Dat(at, "rand"; E = E, F = F))
   end
   return data
end

r0 = rnn(:Si)
calc = StillingerWeber()
data = generate_data(:Si, 2, 0.2*r0, 30, calc)

TRANSFORM = "exp( - 2 * (r/$r0 - 1) )"

rcut2 = 2 * cutoff(calc)
CUTOFF2 = "(:cos, $(rcut2-1), $rcut2)"
D2 = BondLengthDesc(TRANSFORM, CUTOFF2)

rcut3 = 2 * cutoff(calc)
CUTOFF3 = "(:cos, $(rcut3-1), $rcut3)"
D3 = BondLengthDesc(TRANSFORM, CUTOFF3)

@info("Profiling pure 3B basis assembly")
deg3 = 10
B = nbpolys(3, D3, deg3)
@show length(B)
LsqDB("", B, data[1:2]; verbose=false)
@time LsqDB("", B, data)

@info("Profiling combined 2B+3B basis assembly")
deg2, deg3 = 10, 10
B = [nbpolys(2, D2, deg2); nbpolys(3, D3, deg3)]
@show length(B)
LsqDB("", B, data[1:2], verbose=false)
@time LsqDB("", B, data)


text:
struct DriftPath{T <: SSDFloat}
    path::Vector{<:AbstractCoordinatePoint{RealQuantity}}
    timestamps::Vector{RealQuantity}
end

struct EHDriftPath{T <: SSDFloat, TT <: RealQuantity}
    e_path::Vector{<:AbstractCoordinatePoint{T}}
    h_path::Vector{<:AbstractCoordinatePoint{T}}
    timestamps_e::Vector{TT}
    timestamps_h::Vector{TT}
end

abstract type ChargeCarrier end
abstract type Electron <: ChargeCarrier end
abstract type Hole <: ChargeCarrier end

function _common_time(dp::EHDriftPath{T, TT})::TT where {T <: SSDFloat, TT<:RealQuantity}
    max(last(dp.timestamps_e), last(dp.timestamps_h))
end
_common_time(dps::Vector{<:EHDriftPath}) =
maximum(_common_time.(dps))

function _common_timestamps(dp::Union{<:EHDriftPath{T}, Vector{<:EHDriftPath{T}}}, Δt) where {T}
    range(zero(Δt), step = Δt, stop = typeof(Δt)(_common_time(dp)) + Δt)
end

@inline function get_velocity_vector(velocity_field::Interpolations.Extrapolation{<:StaticVector{3}, 3}, pt::CartesianPoint{T})::CartesianVector{T} where {T <: SSDFloat}
    return CartesianVector{T}(velocity_field(pt.x, pt.y, pt.z))
end

@inline function get_velocity_vector(velocity_field::Interpolations.Extrapolation{<:StaticVector{3}, 3}, pt::CylindricalPoint{T}) where {T <: SSDFloat}
    return CartesianVector{T}(velocity_field(pt.r, pt.φ, pt.z))
end

function _drift_charges(det::SolidStateDetector{T}, grid::Grid{T, 3}, point_types::PointTypes{T, 3},
                        starting_points::VectorOfArrays{CartesianPoint{T}}, energies::VectorOfArrays{T},
                        electric_field::Interpolations.Extrapolation{<:SVector{3}, 3},
                        Δt::RQ; max_nsteps::Int = 2000, diffusion::Bool = false, self_repulsion::Bool = false, verbose::Bool = true)::Vector{EHDriftPath{T}} where {T <: SSDFloat, RQ <: RealQuantity}

    drift_paths::Vector{EHDriftPath{T}} = Vector{EHDriftPath{T}}(undef, length(flatview(starting_points)))
    dt::T = T(to_internal_units(Δt))

    drift_path_counter::Int = 0

    for (i, start_points) in enumerate(starting_points)

        n_hits::Int = length(start_points)
        charges::Vector{T} = energies[i] ./ to_internal_units(det.semiconductor.material.E_ionisation)

        drift_path_e::Array{CartesianPoint{T}, 2} = Array{CartesianPoint{T}, 2}(undef, n_hits, max_nsteps)
        drift_path_h::Array{CartesianPoint{T}, 2} = Array{CartesianPoint{T}, 2}(undef, n_hits, max_nsteps)
        timestamps_e::Vector{T} = Vector{T}(undef, max_nsteps)
        timestamps_h::Vector{T} = Vector{T}(undef, max_nsteps)

        n_e::Int = _drift_charge!( drift_path_e, timestamps_e, det, point_types, grid, start_points, -charges, dt, electric_field, Electron, diffusion = diffusion, self_repulsion = self_repulsion, verbose = verbose )
        n_h::Int = _drift_charge!( drift_path_h, timestamps_h, det, point_types, grid, start_points,  charges, dt, electric_field, Hole, diffusion = diffusion, self_repulsion = self_repulsion, verbose = verbose )

        for i in eachindex(start_points)
            drift_paths[drift_path_counter + i] = EHDriftPath{T, T}( drift_path_e[i,1:n_e], drift_path_h[i,1:n_h], timestamps_e[1:n_e], timestamps_h[1:n_h] )
        end

        drift_path_counter += n_hits
    end

    return drift_paths
end

function modulate_surface_drift(p::CartesianVector{T})::CartesianVector{T} where {T <: SSDFloat}
    return p
end

function modulate_driftvector(sv::CartesianVector{T}, pt::CartesianPoint{T}, vdv::Vector{<:AbstractVirtualVolume{T}})::CartesianVector{T} where {T <: SSDFloat}
    for i in eachindex(vdv)
        if in(pt, vdv[i])
            return modulate_driftvector(sv, pt, vdv[i])
        end
    end
    return sv
end
modulate_driftvector(sv::CartesianVector{T}, pt::CartesianPoint{T}, vdv::Missing) where {T} = sv

@inline function _is_next_point_in_det(pt::CartesianPoint{T}, det::SolidStateDetector{T}, point_types::PointTypes{T, 3, Cylindrical})::Bool where {T <: SSDFloat}
    pt_cyl::CylindricalPoint{T} = CylindricalPoint(pt)
    pt_cyl in point_types || pt_cyl in det.semiconductor
end
@inline function _is_next_point_in_det(pt::CartesianPoint{T}, det::SolidStateDetector{T}, point_types::PointTypes{T, 3, Cartesian})::Bool where {T <: SSDFloat}
    pt in point_types || pt in det.semiconductor
end

function project_to_plane(v⃗::AbstractArray, n⃗::AbstractArray) #Vector to be projected, #normal vector of plane
    # solve (v⃗+λ*n⃗) ⋅ n⃗ = 0
    # n⃗ = n⃗ ./ norm(n⃗)
    λ = -1 * dot(v⃗, n⃗) / dot(n⃗, n⃗)
    SVector{3,eltype(v⃗)}(v⃗[1] + λ * n⃗[1], v⃗[2] + λ * n⃗[2], v⃗[3] + λ * n⃗[3])
end

function _set_to_zero_vector!(v::Vector{CartesianVector{T}})::Nothing where {T <: SSDFloat}
    for n in eachindex(v)
        v[n] = CartesianVector{T}(0,0,0)
    end
    nothing
end

function _add_fieldvector_drift!(step_vectors::Vector{CartesianVector{T}}, current_pos::Vector{CartesianPoint{T}}, done::Vector{Bool},
    electric_field::Interpolations.Extrapolation{<:StaticVector{3}, 3}, det::SolidStateDetector{T}, ::Type{S})::Nothing where {T, S}
    for n in eachindex(step_vectors)
       if !done[n]
           step_vectors[n] += get_velocity_vector(electric_field, _convert_point(current_pos[n], S))
           done[n] = (step_vectors[n] == CartesianVector{T}(0,0,0))
       end
    end
    nothing
end

function _add_fieldvector_diffusion!(step_vectors::Vector{CartesianVector{T}}, done::Vector{Bool}, length::T = T(0.5e3))::Nothing where {T <: SSDFloat}
    for n in eachindex(step_vectors)
        if done[n] continue end
        sinθ::T, cosθ::T = sincos(T(rand())*T(2π))
        sinφ::T, cosφ::T = sincos(T(rand())*T(π))
        step_vectors[n] += CartesianVector{T}( length * cosφ * sinθ, length * sinφ * sinθ, length * cosθ )
    end
    nothing
end

function _add_fieldvector_selfrepulsion!(step_vectors::Vector{CartesianVector{T}}, current_pos::Vector{CartesianPoint{T}}, done::Vector{Bool}, charges::Vector{T}, ϵ_r::T)::Nothing where {T <: SSDFloat}
    #TO DO: ignore charges that are already collected (not trapped though!)
    for n in eachindex(step_vectors)
        if done[n] continue end
        for m in eachindex(step_vectors)
            if done[m] continue end
            if m > n
                direction::CartesianVector{T} = current_pos[n] .- current_pos[m]
                if iszero(direction) # if the two charges are at the exact same position
                    continue         # don't let them self-repel each other but treat them as same change
                end                  # if diffusion is simulated, they will very likely be separated in the next step
                tmp::T = elementary_charge * inv(4 * pi * ϵ0 * ϵ_r * max(sum(direction.^2), T(1e-10))) # minimum distance is 10µm
                step_vectors[n] += charges[m] * tmp * normalize(direction)
                step_vectors[m] -= charges[n] * tmp * normalize(direction)
            end
        end
    end
    nothing
end

function _modulate_driftvectors!(step_vectors::Vector{CartesianVector{T}}, current_pos::Vector{CartesianPoint{T}}, vdv::Vector{V})::Nothing where {T <: SSDFloat, V <: AbstractVirtualVolume{T}}
    for n in eachindex(step_vectors)
        step_vectors[n] = modulate_driftvector(step_vectors[n], current_pos[n], vdv)
    end
    nothing
end
_modulate_driftvectors!(step_vectors::Vector{CartesianVector{T}}, current_pos::Vector{CartesianPoint{T}}, ::Missing) where {T <: SSDFloat} = nothing

function _get_driftvectors!(step_vectors::Vector{CartesianVector{T}}, done::Vector{Bool}, Δt::T, cdm::AbstractChargeDriftModel{T}, ::Type{Electron})::Nothing where {T <: SSDFloat}
    for n in eachindex(step_vectors)
       if !done[n]
           step_vectors[n] = getVe(SVector{3,T}(step_vectors[n]), cdm) * Δt
       end
    end
    nothing
end

function _get_driftvectors!(step_vectors::Vector{CartesianVector{T}}, done::Vector{Bool}, Δt::T, cdm::AbstractChargeDriftModel{T}, ::Type{Hole})::Nothing where {T <: SSDFloat}
    for n in eachindex(step_vectors)
       if !done[n]
           step_vectors[n] = getVh(SVector{3,T}(step_vectors[n]), cdm) * Δt
       end
    end
    nothing
end

function _check_and_update_position!(
            step_vectors::Vector{CartesianVector{T}},
            current_pos::Vector{CartesianPoint{T}},
            done::Vector{Bool},
            normal::Vector{Bool},
            drift_path::Array{CartesianPoint{T},2},
            timestamps::Vector{T},
            istep::Int,
            det::SolidStateDetector{T},
            g::Grid{T, 3, S},
            point_types::PointTypes{T, 3, S},
            startpos::AbstractVector{CartesianPoint{T}},
            Δt::T,
            verbose::Bool
        )::Nothing where {T <: SSDFloat, S}

    for n in eachindex(normal)
        done[n] = current_pos[n] == current_pos[n] + step_vectors[n]
        normal[n] = done[n] || _is_next_point_in_det(current_pos[n]+step_vectors[n], det, point_types)
    end

    if all(normal)
        #all charges are either finished or still inside the detector => drift normally
        current_pos .+= step_vectors
        drift_path[:,istep] .= current_pos
    else
        #all charges that would not be inside after the drift step
        for n in findall(.!normal)
            crossing_pos::CartesianPoint{T}, cd_point_type::UInt8, surface_normal::CartesianVector{T} =
                get_crossing_pos(det, point_types, copy(current_pos[n]), current_pos[n] + step_vectors[n])
            if cd_point_type == CD_ELECTRODE
                done[n] = true
                drift_path[n,istep] = crossing_pos
            elseif cd_point_type == CD_FLOATING_BOUNDARY
                projected_vector::CartesianVector{T} = CartesianVector{T}(project_to_plane(step_vectors[n], surface_normal))
                projected_vector = modulate_surface_drift(projected_vector)
                next_pos::CartesianPoint{T} = current_pos[n] + projected_vector
                small_projected_vector = projected_vector * T(0.001)
                i::Int = 0
                while i < 1000 && !(next_pos in det.semiconductor)
                    next_pos -= small_projected_vector
                    i += 1
                end
                if i == 1000
                    if verbose @warn("Handling of charge at floating boundary did not work as intended. Start Position (Cart): $(startpos[n])") end
                    done[n] = true
                    continue
                end
                drift_path[n,istep] = next_pos
                step_vectors *= (1 - i * T(0.001))  # scale down the step_vectors for all other charge clouds
                Δt *= (1 - i * T(0.001))            # scale down Δt for all charge clouds
                done[n] = next_pos == current_pos[n]
                current_pos[n] = next_pos
            else # if cd_point_type == CD_BULK or CD_OUTSIDE
                if verbose @warn ("Internal error for charge starting at $(startpos[n])") end
                done[n] = true
                drift_path[n,istep] = current_pos[n]
            end
        end
        #drift all other charge clouds normally according to the new Δt_min
        for n in findall(normal)
            current_pos[n] += step_vectors[n]
            drift_path[n,istep] = current_pos[n]
        end
    end
    timestamps[istep] = timestamps[istep-1] + Δt
    nothing
end

function _drift_charge!(
                            drift_path::Array{CartesianPoint{T},2},
                            timestamps::Vector{T},
                            det::SolidStateDetector{T},
                            point_types::PointTypes{T, 3, S},
                            grid::Grid{T, 3, S},
                            startpos::AbstractVector{CartesianPoint{T}},
                            charges::Vector{T},
                            Δt::T,
                            electric_field::Interpolations.Extrapolation{<:StaticVector{3}, 3},
                            ::Type{CC};
                            diffusion::Bool = false,
                            self_repulsion::Bool = false,
                            verbose::Bool = true
                        )::Int where {T <: SSDFloat, S, CC <: ChargeCarrier}

    n_hits::Int, max_nsteps::Int = size(drift_path)
    drift_path[:,1] = startpos
    timestamps[1] = zero(T)
    ϵ_r::T = T(det.semiconductor.material.ϵ_r)

    diffusion_length::T = if diffusion
        if CC == Electron
            det.semiconductor.material.diffusion_fieldvector_electrons
        else # CC == Hole
            det.semiconductor.material.diffusion_fieldvector_holes
        end
    else
        zero(T)
    end

    last_real_step_index::Int = 1
    current_pos::Vector{CartesianPoint{T}} = deepcopy(startpos)
    step_vectors::Vector{CartesianVector{T}} = Vector{CartesianVector{T}}(undef, n_hits)
    done::Vector{Bool} = broadcast(pt -> !_is_next_point_in_det(pt, det, point_types), startpos)
    normal::Vector{Bool} = deepcopy(done)

    null_step::CartesianVector{T} = CartesianVector{T}(0, 0, 0)

    @inbounds for istep in 2:max_nsteps
        last_real_step_index += 1
        _set_to_zero_vector!(step_vectors)
        _add_fieldvector_drift!(step_vectors, current_pos, done, electric_field, det, S)
        diffusion && _add_fieldvector_diffusion!(step_vectors, done, diffusion_length)
        self_repulsion && _add_fieldvector_selfrepulsion!(step_vectors, current_pos, done, charges, ϵ_r)
        _get_driftvectors!(step_vectors, done, Δt, det.semiconductor.charge_drift_model, CC)
        _modulate_driftvectors!(step_vectors, current_pos, det.virtual_drift_volumes)
        _check_and_update_position!(step_vectors, current_pos, done, normal, drift_path, timestamps, istep, det, grid, point_types, startpos, Δt, verbose)
        if all(done) break end
    end

    return last_real_step_index
end

# Point types for charge drift
const CD_ELECTRODE = 0x00
const CD_OUTSIDE = 0x01
const CD_BULK = 0x02
const CD_FLOATING_BOUNDARY = 0x04 # not 0x03, so that one could use bit operations here...

function get_crossing_pos(  det::SolidStateDetector{T}, point_types::PointTypes{T, 3, S}, pt_in::CartesianPoint{T}, pt_out::CartesianPoint{T};
                            max_n_iter::Int = 500)::Tuple{CartesianPoint{T}, UInt8, CartesianVector{T}} where {T <: SSDFloat, S}

    # check if the points are already in contacts
    if pt_in in det.contacts return (pt_in, CD_ELECTRODE, CartesianVector{T}(0,0,0)) end
    if pt_out in det.contacts return (pt_out, CD_ELECTRODE, CartesianVector{T}(0,0,0)) end

    direction::CartesianVector{T} = normalize(pt_out - pt_in)
    crossing_pos::Tuple{CartesianPoint{T}, UInt8, CartesianVector{T}} = (pt_out, CD_OUTSIDE, CartesianVector{T}(0,0,0)) # need undef version for this

    # define a Line between pt_in and pt_out
    line::ConstructiveSolidGeometry.Line{T} = ConstructiveSolidGeometry.Line{T}(pt_in, direction)

    # check if the Line intersects with a surface of a Contact
    for contact in det.contacts
        for surf in ConstructiveSolidGeometry.surfaces(contact.geometry)
            for pt in ConstructiveSolidGeometry.intersection(surf, line)
                if pt in contact &&
                    0 ≤ (pt - pt_in) ⋅ direction ≤ 1 &&              # pt within pt_in and pt_out
                    norm(pt - pt_in) < norm(crossing_pos[1] - pt_in) # pt closer to pt_in that previous crossing_pos
                    crossing_pos = (pt, CD_ELECTRODE, normalize(ConstructiveSolidGeometry.normal(surf, pt)))
                end
            end
        end
    end

    # if the Line does not intersect with a surface of a Contact, check if it does intersect with the surface of the Semiconductor
    if crossing_pos[2] & CD_OUTSIDE > 0
        tol::T = 5000 * ConstructiveSolidGeometry.csg_default_tol(T)
        # check if the Line intersects with a surface of the Semiconductor
        for surf in ConstructiveSolidGeometry.surfaces(det.semiconductor.geometry)
            for pt in ConstructiveSolidGeometry.intersection(surf, line)
                normal = normalize(ConstructiveSolidGeometry.normal(surf, pt))
                if pt + tol * normal in det.semiconductor &&         # point "before" crossing_pos should be in
                   !(pt - tol * normal in det.semiconductor) &&      # point "after" crossing_pos should be out
                    0 ≤ (pt - pt_in) ⋅ direction ≤ 1 &&              # pt within pt_in and pt_out
                    norm(pt - pt_in) < norm(crossing_pos[1] - pt_in) # pt closer to pt_in that previous crossing_pos
                    CD_POINTTYPE::UInt8 = point_types[pt] & update_bit == 0 ? CD_ELECTRODE : CD_FLOATING_BOUNDARY
                    crossing_pos = (pt, CD_POINTTYPE, normal)
                end
            end
        end
    end

    crossing_pos
end


text:
# Patrick Keast, MODERATE-DEGREE TETRAHEDRAL QUADRATURE FORMULAS
# http://mech.fsv.cvut.cz/oofem/resources/doc/oofemrefman/gaussintegrationrule_8C_source.html

function (::Type{QuadratureRule{Tetrahedron}})(quad_type::Gauss, order::Int)
    if order == 0 || order == 1
        points = [Tensors.Vec{3}([1.0/4.0, 1.0/4.0, 1.0/1.4])]
        weigths = [1.0/6.0]
    elseif order == 2
        a = ( 5. + 3. * √(5.) ) / 20.
        b = ( 5. - √(5.) ) / 20.
        w = 1. / 24.
        points = [Tensors.Vec{3}([a,b,b]),
                  Tensors.Vec{3}([b,a,b]),
                  Tensors.Vec{3}([b,b,a]),
                  Tensors.Vec{3}([b,b,b])]
        weigths = [w,w,w,w]
    elseif order == 3
        a1 = 1. / 4.
        a2 = 1. / 2.
        b2 = 1. / 6.
        w1 = -2. / 15.
        w2 = 3. / 40.
        points = [Tensors.Vec{3}([a1,a1,a1]),
                  Tensors.Vec{3}([a2,b2,b2]),
                  Tensors.Vec{3}([b2,a2,b2]),
                  Tensors.Vec{3}([b2,b2,a2]),
                  Tensors.Vec{3}([b2,b2,b2])]
        weigths = [w1,w2,w2,w2,w2]
    elseif order == 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.;
        points = [Tensors.Vec{3}([a1,a1,a1]),
                  Tensors.Vec{3}([b2,a2,a2]),
                  Tensors.Vec{3}([a2,b2,a2]),
                  Tensors.Vec{3}([a2,a2,b2]),
                  Tensors.Vec{3}([a2,a2,a2]),
                  Tensors.Vec{3}([a3,a3,b3]),
                  Tensors.Vec{3}([a3,b3,a3]),
                  Tensors.Vec{3}([a3,b3,b3]),
                  Tensors.Vec{3}([b3,a3,a3]),
                  Tensors.Vec{3}([b3,a3,b3]),
                  Tensors.Vec{3}([b3,b3,a3])]

        weigths = [w1,w2,w2,w2,w2,w3,w3,w3,w3,w3,w3]
    else
        throw(ArgumentError("Gauss rule of order $order for tetrahedron not available"))
    end
    return QuadratureRule{Tetrahedron,3,Float64}(weigths, points)
end

max_stars_repo_name: Paulms/jFEMToools

max_issues_repo_name: Paulms/jFEMToools

max_forks_repo_name: Paulms/jFEMToools


text:
include("runtests_Tmatrix_for_spheroid.jl")
include("runtests_autodiff_Tmatrix.jl")


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

function num_thetas(lshaped::AbstractLShaped, ::AbstractLShapedExecution)
    return num_thetas(num_subproblems(lshaped),
                      lshaped.aggregation,
                      scenarioproblems(lshaped.structure))
end

function timestamp(lshaped::AbstractLShaped, ::AbstractLShapedExecution)
    return lshaped.data.iterations
end

function current_decision(lshaped::AbstractLShaped, ::AbstractLShapedExecution)
    return lshaped.x
end

function incumbent_decision(::AbstractLShaped, ::Integer, regularization::AbstractRegularization, ::AbstractLShapedExecution)
    return map(regularization.ξ) do ξᵢ
        return ξᵢ.value
    end
end

function incumbent_objective(::AbstractLShaped, ::Integer, regularization::AbstractRegularization, ::AbstractLShapedExecution)
    return regularization.data.Q̃
end

function incumbent_trustregion(::AbstractLShaped, ::Integer, rd::RegularizedDecomposition, ::AbstractLShapedExecution)
    return rd.data.σ
end

function incumbent_trustregion(::AbstractLShaped, ::Integer, tr::TrustRegion, ::AbstractLShapedExecution)
    Δ = StochasticPrograms.decision(tr.decisions, tr.data.Δ)
    return Δ.value
end

function start_workers!(::AbstractLShaped, ::AbstractLShapedExecution)
    return nothing
end

function close_workers!(::AbstractLShaped, ::AbstractLShapedExecution)
    return nothing
end

function readd_cuts!(lshaped::AbstractLShaped, consolidation::Consolidation, ::AbstractLShapedExecution)
    for i in eachindex(consolidation.cuts)
        for cut in consolidation.cuts[i]
            add_cut!(lshaped, cut; consider_consolidation = false, check = false)
        end
        for cut in consolidation.feasibility_cuts[i]
            add_cut!(lshaped, cut; consider_consolidation = false, check = false)
        end
    end
    return nothing
end

function subobjectives(lshaped::AbstractLShaped, execution::AbstractLShapedExecution)
    return execution.subobjectives
end

function set_subobjectives(lshaped::AbstractLShaped, Qs::AbstractVector, execution::AbstractLShapedExecution)
    execution.subobjectives .= Qs
    return nothing
end

function model_objectives(lshaped::AbstractLShaped, execution::AbstractLShapedExecution)
    return execution.model_objectives
end

function set_model_objectives(lshaped::AbstractLShaped, θs::AbstractVector, execution::AbstractLShapedExecution)
    ids = active_model_objectives(lshaped)
    execution.model_objectives[ids] .= θs[ids]
    return nothing
end

function solve_master!(lshaped::AbstractLShaped, ::AbstractLShapedExecution)
    try
        MOI.optimize!(lshaped.master)
    catch err
        status = MOI.get(lshaped.master, MOI.TerminationStatus())
        # Master problem could not be solved for some reason.
        @unpack Q,θ = lshaped.data
        gap = abs(θ-Q)/(abs(Q)+1e-10)
        # Always print this warning
        @warn "Master problem could not be solved, solver returned status $status. The following relative tolerance was reached: $(@sprintf("%.1e",gap)). Aborting procedure."
        rethrow(err)
    end
    return MOI.get(lshaped.master, MOI.TerminationStatus())
end

function iterate!(lshaped::AbstractLShaped, ::AbstractLShapedExecution)
    # Resolve all subproblems at the current optimal solution
    Q, added = resolve_subproblems!(lshaped)
    if Q == Inf && !handle_feasibility(lshaped.feasibility)
        @warn "Stochastic program is not second-stage feasible at the current decision. Rerun procedure with feasibility_strategy = FeasibilityCuts to use feasibility cuts."
        # Early termination log
        log!(lshaped; status = MOI.INFEASIBLE)
        return MOI.INFEASIBLE
    end
    if Q == -Inf
        # Early termination log
        log!(lshaped; status = MOI.DUAL_INFEASIBLE)
        return MOI.DUAL_INFEASIBLE
    end
    lshaped.data.Q = Q
    # Update incumbent (if applicable)
    take_step!(lshaped)
    # Early gap optimality check if using level sets
    if lshaped.regularization isa LevelSet && check_optimality(lshaped, true)
        # Resolve subproblems with optimal vector
        lshaped.x .= decision(lshaped)
        resolve_subproblems!(lshaped)
        # Optimal, final log
        log!(lshaped; optimal = true)
        return MOI.OPTIMAL
    end
    # Solve master problem
    status = solve_master!(lshaped)
    if !(status ∈ AcceptableTermination)
        # Early termination log
        log!(lshaped; status = status)
        return status
    end
    # Update master solution
    update_solution!(lshaped)
    lshaped.data.θ = calculate_estimate(lshaped)
    # Handle integrality
    handle_integrality!(lshaped, lshaped.integer)
    # Log progress
    log!(lshaped)
    # Check optimality
    if check_optimality(lshaped, added)
        # Optimal, final log
        log!(lshaped; optimal = true)
        return MOI.OPTIMAL
    end
    # Calculate time spent so far and check perform time limit check
    time_spent = lshaped.progress.tlast - lshaped.progress.tinit
    if time_spent >= lshaped.parameters.time_limit
        log!(lshaped; status = MOI.TIME_LIMIT)
        return MOI.TIME_LIMIT
    end
    # Consolidate (if applicable)
    consolidate!(lshaped, lshaped.consolidation)
    # Dont return a status as procedure should continue
    return nothing
end


text:
workspace()
ENV["JuBrain_Mode"]="SingleCore"
using JuBrain
model=networkModel()

model.parameters="""
area = 20000*um^2;Cm = 1*ufarad*cm^(-2)*area;gl = 5e-5*siemens*cm^(-2)*area
El = -65*mV;EK = -90*mV;ENa = 50*mV;VT = -63*mV
g_na = 100*msiemens*cm^(-2)*area;g_kd = 30*msiemens*cm^(-2)*area
taue = 5*ms;taui = 10*ms;Ee = 0*mV;Ei = -80*mV"""
model.equations="""
dv/dt=(gl*(El-v)+ge*(Ee-v)+gi*(Ei-v)-g_na*m^3*h*(v-ENa)-g_kd*n^4*(v-EK))/Cm
dm/dt=0.32*(13-v+VT)/(exp((13-v+VT)/4)-1)*(1-m)-0.28*(v-VT-40)/(exp((v-VT-40)/5)-1)*m
dn/dt=0.032*(15-v+VT)/(exp((15-v+VT)/5)-1)*(1-n)-0.5*exp((10-v+VT)/40)*n
dh/dt=0.128*exp((17-v+VT)/18)*(1-h)-4/(1+exp((40-v+VT)/5))*h
dge/dt=-ge/taue;dgi/dt=-gi/taui"""

model.groups=Groups(groupName=["e","i"],groupSize=[3200,800])
initialize(model,groups=["e","i"],expr="v=El+5*(randn()-1)")
initialize(model,groups=["e","i"],expr="h=1.0")
initialize(model,groups=["e","i"],expr="ge=(randn()*1.5+4)*10*nS")
initialize(model,groups=["e","i"],expr="gi=(randn()*12+20)*10*nS")
spConnect(model,pre="e",post=["e","i"],expr="we=6*nS",p=0.02)
spConnect(model,pre="i",post=["e","i"],expr="wi=67*nS",p=0.02)

model.spike="v.>-20;t.>ts+3"
model.synapse="ge+=we;gi+=wi"
model.record=[1,10,100]

model.postProcs="""
@load("COBAHH_code.jld")
using MatlabPlot
figure();mplot(vRecord')"""

model.nCores=4;model.dt=0.06;model.solver="Euler4";model.synctime=0.0
buildnet(name="COBAHH_code",model=model,duration=1000)


text:
module EntropicToy

# Write your package code here.
greet() = "Hello World 3"
greet_random() = "Hello Random"

export greet,greet_random

export func, func2

"""
    func(x)

Returns double the number `x` plus `1`.
"""
func(x) = 2x + 1

"""
    func2(x,y)

Returns something else double the number `x` plus `1`.
"""
func2(x,y) = 2x + 1 + 3y

end


text:
module PaperPlots
using Plots, Measures
using Plots: PGFPlotsBackend, PyPlotBackend, GRBackend
export PGFPlotsBackend, PyPlotBackend, GRBackend
export linewidth, textwidth, scalefont, allfontsizes

export SIGPLAN

const Backend = Union{PGFPlotsBackend, PyPlotBackend, GRBackend}

include("templates.jl")

"""
Return the font size to pass to `plot` for a target font size in
the rendered pdf of

`scalefont(8, PGFPlotsBackend; width_phys = linewidth(SIGPLAN), dpi = 300)``

# Arguments
size_px : Size of the image in pixels, (i.e. size argument of Plots)
size_phy : size in
"""
function scalefont(fontsize; width_px = 800, width_phys, dpi = 100, scale)
  @show unit_width = (width_px / dpi) * inch
  @show rel_width = unit_width / width_phys
  @show fsz = fontsize * rel_width * scale
  Int(round(fsz))
end

scalefont(fontsize, ::Type{B}; width_px = 800, width_phys, dpi = dpi(B)) where {B <: Backend} =
  scalefont(fontsize; width_px = width_px, width_phys = width_phys, dpi = dpi, scale = scale(B))

""""From a small font size generate a named tuple for everything"""
function fromfontsm(fntsm; scale = 1.4)
  fntlg = Int(round(fntsm * scale))
  (titlefontsize = fntlg, guidefontsize = fntlg, tickfontsize = fntsm, legendfontsize = fntsm)
end

"""

`allfontsizes(SIGPLAN, PGFPlotsBackend; dpi = 300`
"""
function allfontsizes(venue, backend::Type; scale = 1.4, width_phys = linewidth(venue), fntsm = defaultsm(venue), kwargs...)
  fntsm_rc = scalefont(fntsm, backend; width_phys = width_phys, kwargs...)
  fromfontsm(fntsm_rc; scale = scale)
end

end


text:
# © 2016 Massachusetts Institute of Technology.  See LICENSE file for details.
type ANode{S,A}
  id::Int
  name::Compat.UTF8String
  input::Vector{ANode}
  size::Tuple{Vararg{Int}}

  arg::A

  control::Vector{ANode}

  @compat (::Type{ANode}){T}(t::Type{T}, x::T; c::Bool=false) = new{c ? :const : :load, T}(0, "", [], size(x), x, [])
  @compat (::Type{ANode}){T}(t::Type{T}, x; c::Bool=false) = new{c ? :const : :load, T}(0, "", [], size(x), T(x), [])
  @compat (::Type{ANode})(x; c::Bool=false) = ANode(typeof(x), x; c=c)

  @compat (::Type{ANode})(x::Int...) = ANode(x)
  @compat (::Type{ANode})(x::Tuple{Vararg{Int}}) = new{:input,Void}(0, "", [], x, nothing, [])

  @compat (::Type{ANode})(s::Symbol, input::Tuple{Vararg{ANode}}, size::Tuple{Vararg{Int}}, arg=nothing) =
	new{s,typeof(arg)}(0, "", collect(input), size, arg, [])
end

@compat (::Type{ANode})(name::AbstractString, x...) = name!(ANode(x...), name)

sym{S}(::ANode{S}) = S
nin(n::ANode) = length(n.input)

id(n::ANode) = n.id
name(n::ANode) = n.name
input(n::ANode) = n.input
arg(n::ANode) = n.arg
mutates(n::ANode) = false
name!(n::ANode, s::AbstractString) = (n.name = utf8(s); n)

function named_(e::Expr)
  if e.head in (:block, :const, :global)
    return Expr(e.head, map(named_, e.args)...)
  elseif e.head == :(=) && isa(e.args[1], Symbol)
    return :($(esc(e.args[1])) = name!($(esc(e.args[2])), $(string(e.args[1]))))
  elseif e.head == :line
    return e
  else
    error("unexpected expression: $e")
  end
end

macro named(arg...)
  if length(arg) == 2 && (isa(arg[1], AbstractString) || isa(arg[1], Symbol))
    name,node = arg
    return :($(esc(Symbol(name))) = name!($(esc(node)), $(string(name))))
  elseif length(arg) == 1 && isa(arg[1], Expr)
    return named_(arg[1])
  else
    error("unexpected arguments to Seep.@named")
  end
end

function order{S,A}(n::AbstractArray{ANode{S,A}})
  for i in 2:length(n)
    push!(n[i].control, n[i-1])
  end
  n[end]
end

function order(n::ANode...)
  for i in 2:length(n)
    push!(n[i].control, n[i-1])
  end
  n[end]
end

function noop(nodes::AbstractVector{ANode})
  n = ANode(:noop, (), ())
  append!(n.control, nodes)
  n
end

function noop(nodes::ANode...)
  n = ANode(:noop, (), ())
  append!(n.control, collect(ANode, nodes))
  n
end

Base.size(x::ANode) = x.size
Base.size(x::ANode, i::Integer) = i <=  length(x.size) ? x.size[i] : 1
Base.length(x::ANode) = prod(x.size)
Base.ndims(x::ANode) = length(x.size)

gradient_node(a::ANode, wrt::Union{ANode,Int}, b::ANode) = b.*gradient_node(a, wrt)

function gradient_node(a::ANode, wrt::ANode, b::ANode)
  for i in 1:length(a.input)
    if wrt == a.input[i]
      return gradient_node(a, Val{i}, b)
    end
  end
  error("Invalid gradient node")
end

connected_nodes(n::ANode...) = connected_nodes(collect(ANode, n))
connected_nodes(n::ANode, visited=Set{ANode}()) = connected_nodes(ANode[n], visited)

function remove_visited_nodes(q::Vector{ANode}, visited::Set{ANode})
  N = length(q)
  i = 1
  while i <= N && !(q[i] in visited)
    i += 1
  end

  if i <= N
    j = i
    for i in j+1:N
      if !(q[i] in visited)
        q[j] = q[i]
        j += 1
      end
    end

    resize!(q, j-1)
  end
end

function connected_nodes{A<:ANode}(n::Vector{A}, visited=Set{ANode}())
  q = ANode[n;]
  while !isempty(q)
    nn = pop!(q)
    if nn in visited
      remove_visited_nodes(q, visited)
      continue
    end
    push!(visited, nn)
    append!(q, nn.input)
    append!(q, nn.control)
  end

  return visited
end

function output_map(nodes)
  inputs = Dict{ANode,Vector{ANode}}(map(x->x=>ANode[], nodes))
  outputs = Dict{ANode,Vector{ANode}}(map(x->x=>ANode[], nodes))

  for a in nodes
    inputs[a] = unique(filter(x->haskey(inputs, x), [a.input; a.control]))
    for b in inputs[a]
      push!(outputs[b], a)
    end
  end

  for mutator in filter(mutates, nodes)
    a = mutator.input[1]
    if !haskey(outputs, a)
      # a is not a part of this subgraph.
      continue
    end
    out = outputs[a]
    @assert count(mutates, out) == 1 "Too many mutators attached to $a"
    # a has a mutator.  all of its other output must be executed before the mutator.
    for b in out
      if b == mutator continue end
      if mutator in outputs[b] continue end
      push!(inputs[mutator], b)
      push!(outputs[b], mutator)
    end
  end

  return outputs,inputs
end

toposort(n::ANode...) = toposort(collect(ANode, n))

function toposort{A<:ANode}(n::Vector{A}, limit=Union{}[])
  unsorted = connected_nodes(n, Set{ANode}(limit))
  outputs,inputs = output_map(unsorted)

  ready = Set{ANode}(keys(filter((node,input)->isempty(input), inputs)))
  union!(ready, limit)

  sorted = sizehint!(Array{ANode}(0), length(unsorted))
  while !isempty(ready)
    a = pop!(ready)

    if sym(a) != :noop && a∉limit
      push!(sorted, a)
    end

    delete!(unsorted, a)

    for b in outputs[a]
      if b in limit continue end
      @assert !(b in ready)
      if !any(x->x in unsorted, inputs[b])
        push!(ready, b)
      end
    end
  end

  @assert isempty(unsorted)

  return sorted
end

function connected_components(nodes::Vector{ANode},
  outputs::Dict{ANode,Vector{ANode}}, inputs::Dict{ANode,Vector{ANode}})
  cc = Vector{ANode}[]
  s = Set{ANode}(nodes)
  while !isempty(s)
    n = first(s)
    q = ANode[n]
    cc1 = ANode[]
    while !isempty(q)
      i = pop!(q)
      if i in s
        delete!(s, i)
        push!(cc1, i)
        append!(q, inputs[i])
        append!(q, outputs[i])
      end
    end
    push!(cc, cc1)
  end
  return cc
end

function Base.show(io::IO, n::ANode)
  if !isempty(n.name)
    print(io, "ANode{$(sym(n))}($(n.id), $(n.name))")
  else
    print(io, "ANode{$(sym(n))}($(n.id))")
  end
end

graphviz(fn::AbstractString, nn...) = open(f->graphviz(f, nn...), fn, "w")
graphviz(io::IO, nn::ANode...) = graphviz(io, collect(ANode, nn))

function graphviz{A<:ANode}(io::IO, nn::Vector{A})
  println(io, "digraph {")

  unsorted = connected_nodes(nn)

  inputs = Dict{ANode,Vector{ANode}}()
  outputs = Dict{ANode,Vector{ANode}}()
  special = Dict{ANode,Vector{ANode}}()

  id = 1
  for a in unsorted
    a.id = id
    id += 1
    inputs[a] = collect(a.input)
    append!(inputs[a], a.control)

    for b in inputs[a]
      push!(get!(outputs, b, ANode[]), a)
    end
  end

  for a in unsorted
    if mutates(a)
      @assert count(mutates, outputs[a.input[1]]) == 1 "Too many mutators attached to $(a.input[1])"
      ca = connected_nodes(a)

      for b in outputs[a.input[1]]
        if b != a && !(b in ca)
          push!(get!(special, a, ANode[]), b)
        end
      end
    end
  end

  for n in unsorted
    label = (isempty(n.name) ? "" : string(n.name, "\\n"))*"$(sym(n)) ($(n.id))"
    color = elementwise(n) ? "black" : "red"
    println(io, "  n$(n.id) [label=\"$label\",color=\"$color\"];")

    for i in 1:length(n.input)
      println(io, "    n$(n.input[i].id) -> n$(n.id);")
    end

    for i in 1:length(n.control)
      println(io, "    n$(n.control[i].id) -> n$(n.id) [style=dashed];")
    end

    if haskey(special, n)
      for i in 1:length(special[n])
        println(io, "    n$(special[n][i].id) -> n$(n.id) [style=dotted];")
      end
    end
  end

  println(io, "}")
end