text:
mutable struct EulerNet{Tv,Tt,Tw,Tdelay,Tref,TΔt,Tstrace,Tparams,Trecorder,Tplasticity,Tnoise} <: AbstractNet{Tplasticity}
    v::Tv
    n_of_neurons::Int64
    t::Tt
    posts::Array{Array{Tuple{Int64,Int64}, 1}, 1}
    pres::Array{Array{Tuple{Int64,Int64}, 1}, 1}
    weights::Tw
    spiketrace::Tstrace
    spikecandidates::Array{Int64, 1}
    Δt::TΔt
    indelay::Array{Int64, 1}
    inref::Array{Int64, 1}
    delay::Tdelay
    refractory_period::Tref
    parameters::Tparams
    recorder::Trecorder
    plasticityrule::Tplasticity
    noise::Tnoise
end
function EulerNet(net; Δt = .1f0)
    if minimum(net.parameters.delay) == 0
        error("Cannot integrate nets with zero spike transmission delay.
              delay >= Δt")
    end
    EulerNet(net.v, net.n_of_neurons,
             net.t, net.posts, net.pres, net.weights, net.spiketrace,
             Int64[],
             Δt, -ones(Int64, net.n_of_neurons), -ones(Int64, net.n_of_neurons),
             round.(Int64, net.parameters.delay/Δt),
             round.(Int64, net.parameters.refractory_period/Δt),
             net.parameters, net.recorder, net.plasticityrule, net.noise)
end
export EulerNet

"""
    deepcopynet(net;
                recorder = net.recorder,
                plasticityrule = net.plasticityrule)

Copy an EulerNet and optionally replace the `recorder` or the `plasticityrule`.
"""
function deepcopynet(net::EulerNet;
                     recorder = deepcopy(net.recorder),
                     plasticityrule = deepcopy(net.plasticityrule),
                     parameters = deepcopy(net.parameters))
    EulerNet([field == :recorder ? recorder :
         field == :plasticityrule ? plasticityrule :
         field == :parameters ? parameters :
         deepcopy(getfield(net, field)) for field in fieldnames(typeof(net))]...)
end
export deepcopynet

function addspike!(net::EulerNet, t, i)
    net.indelay[i] = getdelay(net.delay, i)
    net.inref[i] = getrefractory_period(net.refractory_period, i)
    (typeof(net.noise) == RectLinEscapeNoise) &&
    (net.noise.endofrefrperiod[i] = t + getrefractory_period(net.parameters.refractory_period, i))
    addspikepostprocess!(net, t, i)
end
@inline function eulerspike!(net::EulerNet, noise::NoNoise, i)
  if abovethreshold(net.v, net.parameters.threshold, i)
      push!(net.spikecandidates, i)
  end
end
@inline getrectlinhazardrate(threshold::Number, v, β, i) = max(0, β * (v[i] - threshold))
@inline getrectlinhazardrate(threshold::AbstractArray, v, β, i) = max(0, β * (v[i] - threshold[i]))
@inline function eulerspike!(net::EulerNet, noise::RectLinEscapeNoise, i)
  if rand() <= net.Δt*getrectlinhazardrate(net.parameters.threshold, net.v, noise.β, i)
    push!(net.spikecandidates, i)
  end
end
function updatenet!(net::EulerNet, input, tmax)
    for eventtime in net.t:net.Δt:tmax
        updatevars!(net, eventtime, input)
        net.t = eventtime
        @inbounds for i in 1:net.n_of_neurons
            eulerspike!(net, net.noise, i)
            if net.indelay[i] == 1
                broadcastspike!(net, i)
                updateweights!(net, i)
            end
            if net.indelay[i] >= 1 net.indelay[i] -= 1 end
            if net.inref[i] == 0
                net.v[i] = getreset(net.parameters.reset, i)
            end
            if net.inref[i] >= 0 net.inref[i] -= 1 end
        end
        processspikecandidates!(net)
    end
    return false
end


text:
function drop!(cur_player::AbstractPlayer, cur_piece::AbstractPiece, cur_grid::AbstractGrid, cur_direction::Integer, is_main_piece::Bool, is_step::Bool)
  can_move(cur_player) || return false

  cur_blocks = map(
    cur_block -> calc_block_coords(cur_block),
    cur_piece.blocks
  )

  did_move = false

  new_direction = 0

  for tmp_direction in sign(cur_direction):sign(cur_direction):cur_direction

    for (cur_index, cur_block) in enumerate(cur_blocks)
      cur_blocks[cur_index] = (
        first(cur_block) + sign(cur_direction),
        last(cur_block) + 0
      )
    end

    check(cur_blocks, cur_grid) || break

    did_move = true

    new_direction = tmp_direction

  end

  did_move || return false

  cur_piece.i += new_direction

  if is_main_piece && !is_step
    drop_count = abs(new_direction)

    if drop_count == 1
      score!(cur_player, "soft", -1)
    else
      score!(cur_player, "hard", drop_count)
    end
  end

  if is_main_piece
    move!(cur_player, false)

    cur_player.clock.drag =
      Nullable{UUIDs.UUID}()
  end

  true
end

function drop!(cur_player::AbstractPlayer, cur_piece::Nullable{AbstractPiece}, cur_grid::AbstractGrid, cur_direction::Integer, is_main_piece::Bool, is_step::Bool)
  false
end

ext: jl


text:
function optimize_controls!(
        m::ClimateModel;
        obj_option = "adaptive_temp", temp_goal = 1.5, budget=10., expenditure = 0.5,
        max_deployment = Dict("mitigate"=>1., "remove"=>1., "geoeng"=>1., "adapt"=>1.),
        max_slope = Dict("mitigate"=>1. /40., "remove"=>1. /40., "geoeng"=>1. /40., "adapt"=> 1. /40.),
        temp_overshoot = nothing,
        overshoot_year = 2100,
        delay_deployment = Dict(
            "mitigate"=>0,
            "remove"=>0,
            "geoeng"=>0,
            "adapt"=>0
        ),
        cost_exponent = 3.,
        print_status = false, print_statistics = false, print_raw_status = true,
    )

    # Translate JuMP printing options
    if print_status
        if print_statistics
            bool_str = "yes"
        else
            bool_str = "no"
        end
        print_int = 1
    else
        print_int = 0
        bool_str = "no"
    end

    # Shorthands
    tarr = t(m)
    Earr = E(m)
    τ = τd(m)
    dt = m.domain.dt
    t0 = tarr[1]
    tp = m.domain.present_year
    q = m.economics.baseline_emissions
    qGtCO2 = ppm_to_GtCO2(q)
    N = length(tarr)

    # Set default temperature goal for end year
    if isnothing(temp_overshoot)
        temp_overshoot = temp_goal
        odx = N
    elseif temp_overshoot >= temp_goal
        odx = argmin(abs.(tarr .- overshoot_year)) # overshoot index
    elseif temp_overshoot < temp_goal
        temp_overshoot = temp_goal
        odx = N
    end

    # Set defaults for start_deployment
    start_deployment = Dict()
    for (key, item) in delay_deployment
        start_deployment[key] = t0 + delay_deployment[key]
    end

    if typeof(cost_exponent) in [Int64, Float64]
        cost_exponent = Dict(
            "mitigate"=>cost_exponent,
            "remove"=>cost_exponent,
            "geoeng"=>cost_exponent,
            "adapt"=>cost_exponent
        )
    end

    model_optimizer = Model(optimizer_with_attributes(Ipopt.Optimizer,
        "acceptable_tol" => 1.e-8, "max_iter" => Int64(1e8),
        "acceptable_constr_viol_tol" => 1.e-3, "constr_viol_tol" => 1.e-4,
        "print_frequency_iter" => 50,  "print_timing_statistics" => bool_str,
        "print_level" => print_int,
    ))

    function fM_JuMP(α)
        if α <= 0.
            return 1000.
        else
            return α ^ cost_exponent["mitigate"]
        end
    end
    register(model_optimizer, :fM_JuMP, 1, fM_JuMP, autodiff=true)

    function fA_JuMP(α)
        if α <= 0.
            return 1000.
        else
            return α ^ cost_exponent["adapt"]
        end
    end
    register(model_optimizer, :fA_JuMP, 1, fA_JuMP, autodiff=true)

    function fR_JuMP(α)
        if α <= 0.
            return 1000.
        else
            return α ^ cost_exponent["remove"]
        end
    end
    register(model_optimizer, :fR_JuMP, 1, fR_JuMP, autodiff=true)

    function fG_JuMP(α)
        if α <= 0.
            return 1000.
        else
            return α ^ cost_exponent["geoeng"]
        end
    end
    register(model_optimizer, :fG_JuMP, 1, fG_JuMP, autodiff=true)

    scale = 5.e-3
    function Hstep(α)
        if α <= 0.
            return 0.
        elseif 0. < α <= scale/2.
            return 2*(α/scale)^2
        elseif scale/2. < α <= scale
            return 1. - 2((α - scale)/scale)^2
        elseif α > scale
            return 1.
        end
    end
    register(model_optimizer, :Hstep, 1, Hstep, autodiff=true)

    function log_JuMP(x)
        if x <= 0.
            return -1000.0
        else
            return log(x)
        end
    end
    register(model_optimizer, :log_JuMP, 1, log_JuMP, autodiff=true)

    function discounting_JuMP(t)
        if t < tp
            return 0.
        else
            return (
                (1. + m.economics.ρ) ^
                (-(t - tp))
            )
        end
    end
    register(model_optimizer, :discounting_JuMP, 1, discounting_JuMP, autodiff=true)

    # constraints on control variables
    for (key, item) in max_deployment
        if item==0.
            max_deployment[key] = item+1.e-8
        end
    end
    @variables(model_optimizer, begin
            0. <= M[1:N] <= max_deployment["mitigate"]  # emissions reductions
            0. <= R[1:N] <= max_deployment["remove"]  # negative emissions
            0. <= G[1:N] <= max_deployment["geoeng"]  # geoengineering
            0. <= A[1:N] <= max_deployment["adapt"]  # adapt
    end)

    control_vars = Dict(
        "mitigate" => M,
        "remove" => R,
        "geoeng" => G,
        "adapt" => A
    )
    controls = Dict(
        "mitigate" => m.controls.mitigate,
        "remove" => m.controls.remove,
        "geoeng" => m.controls.geoeng,
        "adapt" => m.controls.adapt
    )

    domain_idx = (tarr .> tp) # don't update past or present
    if tarr[1] == tp
        domain_idx = (tarr .>= tp) # unless present is also first timestep
    end

    M₀ = m.economics.mitigate_init
    R₀ = m.economics.remove_init
    G₀ = m.economics.geoeng_init
    A₀ = m.economics.adapt_init

    control_inits = Dict(
        "mitigate" => M₀,
        "remove" => R₀,
        "geoeng" => G₀,
        "adapt" => A₀
    )

    for (key, control) in control_vars
        if control_inits[key] !== nothing
            fix(control_vars[key][1], control_inits[key]; force = true)
            Nstart = 2
        else
            Nstart = 1
        end

        for idx in Nstart:N
            if idx <= length(tarr[.~domain_idx])
                fix(control_vars[key][idx], controls[key][idx]; force = true)
            else
                if tarr[idx] < start_deployment[key]
                    if control_inits[key] !== nothing
                        fix(control_vars[key][idx], control_inits[key]; force = true)
                    else
                        fix(control_vars[key][idx], 0.; force = true)
                    end
                end
            end
        end
    end

    # add integral function as a new variable defined by first order finite differences
    @variable(model_optimizer, cumsum_qMR[1:N]);
    for i in 1:N-1
        @constraint(
            model_optimizer, cumsum_qMR[i+1] - cumsum_qMR[i] ==
            (dt * (m.physics.r * (q[i+1] * (1. - M[i+1]) - q[1] * R[i+1])))
        )
    end
    @constraint(
        model_optimizer, cumsum_qMR[1] == (dt * (m.physics.r * (q[1] * (1. - M[1]) - q[1] * R[1])))
    );

    # add temperature kernel as new variable defined by first order finite difference
    @variable(model_optimizer, cumsum_KFdt[1:N]);
    for i in 1:N-1
        @NLconstraint(
            model_optimizer, cumsum_KFdt[i+1] - cumsum_KFdt[i] ==
            (
                dt *
                exp( (tarr[i+1] - (t0 - dt)) / τ ) * (
                    m.physics.a * log_JuMP(
                        (m.physics.c0 + cumsum_qMR[i+1]) /
                        (m.physics.c0)
                    ) - m.economics.Finf*G[i+1] )
            )
        )
    end
    @NLconstraint(
        model_optimizer, cumsum_KFdt[1] ==
        (
            dt *
            exp( dt / τ ) * (
                m.physics.a * log_JuMP(
                    (m.physics.c0 + cumsum_qMR[1]) /
                    (m.physics.c0)
                ) - m.economics.Finf*G[1] )
         )
    );

    # Add constraint of rate of changes
    present_idx = argmin(abs.(tarr .- tp))

    @variables(model_optimizer, begin
            -max_slope["mitigate"] <= dMdt[present_idx:N-1] <= max_slope["mitigate"]
            -max_slope["remove"] <= dRdt[present_idx:N-1] <= max_slope["remove"]
            -max_slope["geoeng"] <= dGdt[present_idx:N-1] <= max_slope["geoeng"]
            -max_slope["adapt"] <= dAdt[present_idx:N-1] <= max_slope["adapt"]
    end);
    for i in present_idx:N-1
        @constraint(model_optimizer, dMdt[i] == (M[i+1] - M[i]) / dt)
    end
    for i in present_idx:N-1
        @constraint(model_optimizer, dRdt[i] == (R[i+1] - R[i]) / dt)
    end
    for i in present_idx:N-1
        @constraint(model_optimizer, dGdt[i] == (G[i+1] - G[i]) / dt)
    end
    for i in present_idx:N-1
        @constraint(model_optimizer, dAdt[i] == (A[i+1] - A[i]) / dt)
    end

    ## Optimization options

    # maximize net benefits
    if obj_option == "net_benefit"
        # (in practice we solve the equivalent problem of minimizing the net cost, i.e. minus the net benefit)
        @NLobjective(model_optimizer, Min,
            sum(
                (
                    m.economics.β *
                    Earr[i] *
                    (
                        (m.physics.T0 +
                        (
                             (m.physics.a * log_JuMP(
                                        (m.physics.c0 + cumsum_qMR[i]) /
                                        (m.physics.c0)
                                    ) - m.economics.Finf*G[i]
                            ) +
                            m.physics.κ /
                            (τ * m.physics.B) *
                            exp( - (tarr[i] - (t0 - dt)) / τ ) *
                            cumsum_KFdt[i]
                        ) / (m.physics.B + m.physics.κ)
                        )^2
                    ) * ( (1 - A[i]) ) +
                    Earr[i] * m.economics.adapt_cost * fA_JuMP(A[i]) +
                    m.economics.mitigate_cost * qGtCO2[i] *
                    fM_JuMP(M[i]) +
                    m.economics.remove_cost * fR_JuMP(R[i]) +
                    Earr[i] * (
                        m.economics.geoeng_cost * fG_JuMP(G[i]) +
                        m.economics.epsilon_cost * Hstep(G[i])
                    )

                ) *
                discounting_JuMP(tarr[i]) *
                dt
            for i=1:N)
        )

    # maximize net benefits subject to a temperature goal
    elseif obj_option == "temp"
        @NLobjective(model_optimizer, Min,
            sum(
                (
                    Earr[i] * m.economics.adapt_cost * fA_JuMP(A[i]) +
                    m.economics.mitigate_cost * qGtCO2[i] *
                    fM_JuMP(M[i]) +
                    m.economics.remove_cost * fR_JuMP(R[i]) +
                    Earr[i] * (
                        m.economics.geoeng_cost * fG_JuMP(G[i]) +
                        m.economics.epsilon_cost * Hstep(G[i])
                    )
                ) *
                discounting_JuMP(tarr[i]) *
                dt
            for i=1:N)
        )

        # Subject to temperature goal (during overshoot period)
        for i in 1:odx-1
            @NLconstraint(model_optimizer,
            (m.physics.T0 +
                (
                        (m.physics.a * log_JuMP(
                                (m.physics.c0 + cumsum_qMR[i]) /
                                (m.physics.c0)
                            ) - m.economics.Finf*G[i]
                    ) +
                    m.physics.κ /
                    (τ * m.physics.B) *
                    exp( - (tarr[i] - (t0 - dt)) / τ ) *
                    cumsum_KFdt[i]
                ) / (m.physics.B + m.physics.κ)
            ) <= temp_overshoot
            )
        end
        # Subject to temperature goal (after temporary overshoot period)
        for i in odx:N
            @NLconstraint(model_optimizer,
            (m.physics.T0 +
                (
                     (m.physics.a * log_JuMP(
                                (m.physics.c0 + cumsum_qMR[i]) /
                                (m.physics.c0)
                            ) - m.economics.Finf*G[i]
                    ) +
                    m.physics.κ /
                    (τ * m.physics.B) *
                    exp( - (tarr[i] - (t0 - dt)) / τ ) *
                    cumsum_KFdt[i]
                ) / (m.physics.B + m.physics.κ)
            ) <= temp_goal
            )
        end

    # minimize control costs subject to adaptative-temperature constraint
    elseif obj_option == "adaptive_temp"

        # minimize control costs
        @NLobjective(model_optimizer, Min,
            sum(
                (
                    Earr[i] * m.economics.adapt_cost * fA_JuMP(A[i]) +
                    m.economics.mitigate_cost * qGtCO2[i] *
                    fM_JuMP(M[i]) +
                    m.economics.remove_cost * fR_JuMP(R[i]) +
                    Earr[i] * (
                        m.economics.geoeng_cost * fG_JuMP(G[i]) +
                        m.economics.epsilon_cost * Hstep(G[i])
                    )
                ) *
                discounting_JuMP(tarr[i]) *
                dt
            for i=1:N)
        )

        # Subject to adaptive temperature goal (during overshoot period)
        for i in 1:odx-1
            @NLconstraint(model_optimizer,
            ((m.physics.T0 +
                (
                     (m.physics.a * log_JuMP(
                                (m.physics.c0 + cumsum_qMR[i]) /
                                (m.physics.c0)
                            ) - m.economics.Finf*G[i]
                    ) +
                    m.physics.κ /
                    (τ * m.physics.B) *
                    exp( - (tarr[i] - (t0 - dt)) / τ ) *
                    cumsum_KFdt[i]
                ) / (m.physics.B + m.physics.κ)
            )
            * sqrt(1. - A[i]) ) <=
                temp_overshoot
            )
        end
        # Subject to adaptive temperature goal (after temporary overshoot period)
        for i in odx:N
            @NLconstraint(model_optimizer,
            ((m.physics.T0 +
                (
                     (m.physics.a * log_JuMP(
                                (m.physics.c0 + cumsum_qMR[i]) /
                                (m.physics.c0)
                            ) - m.economics.Finf*G[i]
                    ) +
                    m.physics.κ /
                    (τ * m.physics.B) *
                    exp( - (tarr[i] - (t0 - dt)) / τ ) *
                    cumsum_KFdt[i]
                ) / (m.physics.B + m.physics.κ)
            )
            * sqrt(1. - A[i]) ) <=
                temp_goal
            )
        end

    # minimize damages subject to a total discounted budget constraint
    elseif obj_option == "budget"
        @NLobjective(model_optimizer, Min,
            sum(
                m.economics.β *
                Earr[i] *
                (
                    m.physics.T0 +
                    (
                         (m.physics.a * log_JuMP(
                                    (m.physics.c0 + cumsum_qMR[i]) /
                                    (m.physics.c0)
                                ) - m.economics.Finf*G[i]
                        ) +
                        m.physics.κ /
                        (τ * m.physics.B) *
                        exp( - (tarr[i] - (t0 - dt)) / τ ) *
                        cumsum_KFdt[i]
                    ) / (m.physics.B + m.physics.κ)
                    )^2 * (1. - A[i]) *
                discounting_JuMP(t[i]) *
                dt
            for i=1:N)
        )

        @NLconstraint(model_optimizer,
            sum(
                (
                    Earr[i] * m.economics.adapt_cost * fA_JuMP(A[i]) +
                    m.economics.mitigate_cost * qGtCO2[i] *
                    fM_JuMP(M[i]) +
                    m.economics.remove_cost * fR_JuMP(R[i]) +
                    Earr[i] * (
                        m.economics.geoeng_cost * fG_JuMP(G[i]) +
                        m.economics.epsilon_cost * Hstep(G[i])
                    )
                ) *
                discounting_JuMP(tarr[i]) *
                dt
            for i=1:N) <= budget
        )

    # minimize damages subject to annual expenditure constraint (as % of GWP)
    elseif obj_option == "expenditure"
        @NLobjective(model_optimizer, Min,
            sum(
                m.economics.β *
                Earr[i] *
                (m.physics.T0 +
                    (
                         (m.physics.a * log_JuMP(
                                    (m.physics.c0 + cumsum_qMR[i]) /
                                    (m.physics.c0)
                                ) - m.economics.Finf*G[i]
                        ) +
                        m.physics.κ /
                        (τ * m.physics.B) *
                        exp( - (tarr[i] - (t0 - dt)) / τ ) *
                        cumsum_KFdt[i]
                    ) / (m.physics.B + m.physics.κ)
                )^2 * (1-A[i]) *
                discounting_JuMP(t[i]) *
                dt
            for i=1:N)
        )

        for i in 1:N
            @NLconstraint(model_optimizer,
                (
                    Earr[i] * m.economics.adapt_cost * fA_JuMP(A[i]) +
                    m.economics.mitigate_cost * qGtCO2[i] *
                    fM_JuMP(M[i]) +
                    m.economics.remove_cost * fR_JuMP(R[i]) +
                    Earr[i] * (
                        m.economics.geoeng_cost * fG_JuMP(G[i]) +
                        m.economics.epsilon_cost * Hstep(G[i])
                    )
                ) <= expenditure * Earr[i]
            )
        end
    end

    optimize!(model_optimizer)

    if print_raw_status
        print(raw_status(model_optimizer), "\n")
    end

    mitigate_values = value.(M)[domain_idx]
    mitigate_values[q[domain_idx].==0.] .= 0.
    getfield(m.controls, :mitigate)[domain_idx] = mitigate_values
    getfield(m.controls, :remove)[domain_idx] = value.(R)[domain_idx]
    getfield(m.controls, :geoeng)[domain_idx] = value.(G)[domain_idx]
    getfield(m.controls, :adapt)[domain_idx] = value.(A)[domain_idx]

    return model_optimizer
end


text:
using Coverage

cd(joinpath(@__DIR__, "..", "..")) do
    Codecov.submit(Codecov.process_folder())
    Coveralls.submit(Coveralls.process_folder())
end


text:
using .Flux.Tracker: TrackedArray, TrackedReal

if !usetyped
  unwrap(x::Union{TrackedArray,TrackedReal}) = Flux.data(x)
end


text:
struct Analysis{F} <: AbstractElement
    f::F
    kwargs::NamedTuple
end

Analysis(f; kwargs...) = Analysis(f, values(kwargs))

Spec(a::Analysis) = Spec{Any}(analyses=(a,))

(a::Analysis)(; kwargs...) = Analysis(a.f, merge(a.kwargs, values(kwargs)))

(a::Analysis)(args...; kwargs...) = a.f(args...; merge(a.kwargs, values(kwargs))...)


text:
@testset "673.number-of-longest-increasing-subsequence.jl" begin
    @test find_number_of_LIS([1, 3, 5, 4, 7]) == 2
    @test find_number_of_LIS([2, 2, 2, 2, 2]) == 5
    @test find_number_of_LIS([1, 29, 58, 28, 9, 53, 38, 42, 33, 58, 4, 9, 26, 45, 53, 5, 3,
                              55, 49, 58, 50, 13, 59, 22, 38, 20, 11, 58, 59, 37, 4, 43, 28,
                              20, 6, 4, 57, 39, 3, 1, 43, 13, 58, 25, 1, 27, 58, 27, 45,
                              58]) == 10
end


text:
using Weave

plot_bool = true

examples_dir = joinpath(dirname(dirname(@__FILE__)), "examples/HydroValleys")
docs_dir = dirname(@__FILE__)
testcases_dir = joinpath(dirname(dirname(@__FILE__)), "testcases")

weave(joinpath(examples_dir,"case3.jl"), out_path =joinpath(docs_dir,"build/examples"), doctype = "md2html", args=Dict(:testcases_dir=>testcases_dir))
weave(joinpath(examples_dir,"case3_5years.jl"), out_path =joinpath(docs_dir,"build/examples"), doctype = "md2html", args=Dict(:testcases_dir=>testcases_dir))
weave(joinpath(examples_dir,"case3_ac.jl"), out_path =joinpath(docs_dir,"build/examples"), doctype = "md2html", args=Dict(:testcases_dir=>testcases_dir))
weave(joinpath(examples_dir,"case3_soc.jl"), out_path =joinpath(docs_dir,"build/examples"), doctype = "md2html", args=Dict(:testcases_dir=>testcases_dir))


text:
immutable SharedData{T<:FloatingPoint}
	Lambda::Diagonal{T}
	sigma2::T
	diagZ::Array{T, 1}
	Z::Union(Base.LinAlg.Cholesky{T}, Array{T,2})
end
function get_shared(m, k, lik)
	Lambda = Diagonal(Float64[S(k, sqrt(lam)) for lam in m.lambdas])
	sigma2 = exp(2lik.log_sigma)
	diagZ = sigma2 ./ diag(Lambda) + diag(m.PhiTPhi)
	Z = sigma2 * inv(Lambda) + m.PhiTPhi
	try
		Z = cholfact(Z)
	catch
	end
	return SharedData(Lambda, sigma2, diagZ, Z)
end

function nll(m::LaplacianModel{OneD})
	s = get_shared(m, m.k, m.lik)
	nll(m, m.k, m.lik, s)
end

function nll(m, k, lik, s)
	n = size(m.X, 1)
	yTQy(m, k, s) / 2 + logdetQ(m, k, s) / 2 + n*log(2pi)
end

function dnlldl(m, k, lik)
	s = get_shared(m, k, lik)
	dnlldl(m, k, lik, s)
end
function dnlldnu(m, k, lik)
	s = get_shared(m, k, lik)
	dnlldnu(m, k, lik, s)
end

function dnlldsigma2(m, k, lik)
	s = get_shared(m, k, lik)
	dnlldsigma2(m, k, lik, s)
end

function dnlldl(m, k, lik, s)
	dyTQydl(m, k, s) / 2 + dlogdetQdl(m, k, s) / 2
end
function dnlldnu(m, k, lik, s)
	dyTQydnu(m, k, s) / 2 + dlogdetQdnu(m, k, s) / 2
end

function dnlldsigma2(m, k, lik, s)
	dyTQydsigma2(m, k, s) / 2 + dlogdetQdsigma2(m, k, s) / 2
end

function yTQy(m, k, s)
	(dot(m.y,m.y) - dot(m.PhiTy, s.Z \ m.PhiTy)) / s.sigma2
end

function dyTQydl(m, k, s)
	dLambdadl = Diagonal([dSdl(k, sqrt(lam)) for lam in m.lambdas])
	- dot(m.PhiTy, s.Z \ (inv(s.Lambda)^2 * dLambdadl *  (s.Z \ m.PhiTy)))
end

function dyTQydnu(m, k, s)
	dLambdadnu = Diagonal([dSdnu(k, sqrt(lam)) for lam in m.lambdas])
	- dot(m.PhiTy, s.Z \ (inv(s.Lambda)^2 * dLambdadnu *  (s.Z \ m.PhiTy)))
end

function dyTQydsigma2(m, k, s)
	2s.sigma2 * (dot(m.PhiTy, s.Z \ (inv(s.Lambda) * (s.Z \ m.PhiTy))) - yTQy(m, k, s)) / s.sigma2
end

function logdetQ(m, k, s)
	logdetZ = logdet(s.Z)
	n = size(m.X, 1)
	(n - m.jmax)*log(s.sigma2) + logdetZ + sum(log(diag(s.Lambda)))
end

function dlogdetQdl(m, k, s)
	dLambdadl = Diagonal([dSdl(k, sqrt(l)) for l in m.lambdas])
	trace(dLambdadl / s.Lambda) - s.sigma2 * sum(1 ./ (s.diagZ  .* inv(inv(s.Lambda)^2 * dLambdadl).diag))
end

function dlogdetQdnu(m, k, s)
	dLambdadnu = Diagonal([dSdnu(k, sqrt(l)) for l in m.lambdas])
	trace(dLambdadnu / s.Lambda) - s.sigma2 * sum(1 ./ (s.diagZ  .* inv(inv(s.Lambda)^2 * dLambdadnu).diag))
end

function dlogdetQdsigma2(m, k, s)
	n = size(m.X, 1)
	((n - m.jmax) / (s.sigma2) + sum(1./ (s.diagZ .* diag(s.Lambda)))) * 2s.sigma2
end


text:
using RegionTrees
using StaticArrays
using PyPlot
using Statistics
using LinearAlgebra
using SparseArrays
using VectorizedRoutines
using ProgressBars
using MAT

struct triangle2d
    pts::Array{Array{Float64,1},1}
    nodes::Union{Array{Int64,1},SVector{3,Int64}}
    bbox::HyperRectangle{2,Float64}
    A::Float64
end

struct ray2d
    origin::Union{SVector{2,Float64},Vector{Float64}}
    dir::Union{SVector{2,Float64},Vector{Float64}}
    inv_direction::Union{SVector{2,Float64},Vector{Float64}}
    sign::Union{Array{Int64,1},SVector{2,Int64}}
end

function triangle2d(
    p1::T,
    p2::T,
    p3::T,
    nlist::Union{Array{Int64,1},SArray{Tuple{3},Int64,1,3}},
) where {T}
    origin = [minimum([p1[1], p2[1], p3[1]]), minimum([p1[2], p2[2], p3[2]])]
    final = [maximum([p1[1], p2[1], p3[1]]), maximum([p1[2], p2[2], p3[2]])]
    width = final - origin
    origin = @SVector[origin[1], origin[2]]
    width = @SVector[width[1], width[2]]
    bbox = HyperRectangle(origin, width)
    A =
        0.5 * abs(
            p1[1] * p2[2] + p2[1] * p3[2] + p3[1] * p1[2] - p2[1] * p1[2] -
            p3[1] * p2[2] - p1[1] * p3[2],
        )
    return triangle2d([p1, p2, p3], nlist, bbox, A)
end

function ray2d(p1::T, p2::T) where {T}
    dir::T = p2 - p1
    invdir::T = 1.0 ./ dir
    sign = (invdir .< 0) .* 1
    return ray2d(p1, dir, invdir, sign)
end

function ray2ddir(p1::T, dir::T) where {T}
    invdir::T = 1.0 ./ dir
    sign = (invdir .< 0) .* 1
    return ray2d(p1, dir, invdir, sign)
end

@inline @inbounds function intersects(
    r::ray2d,
    b::HyperRectangle{2,Float64},
)::Float64
    parameters = [b.origin, b.origin + b.widths]
    tmin = (parameters[r.sign[1]+1][1] - r.origin[1]) * r.inv_direction[1]
    tmax = (parameters[1-r.sign[1]+1][1] - r.origin[1]) * r.inv_direction[1]
    tymin = (parameters[r.sign[2]+1][2] - r.origin[2]) * r.inv_direction[2]
    tymax = (parameters[1-r.sign[2]+1][2] - r.origin[2]) * r.inv_direction[2]
    if ((tmin > tymax) || (tymin > tmax))
        return -1.0
    end
    if (tymin > tmin)
        tmin = tymin
    end
    if (tymax < tmax)
        tmax = tymax
    end
    if (~isfinite(tmax) || ~isfinite(tmin))
        error("An ambiguous ray trajectory detected: ray touches the pixel at its boundary. ")
    end
    return norm((r.origin + r.dir * tmin) - (r.origin + r.dir * tmax))

end

@inline function intersects(
    a::HyperRectangle{N,Float64},
    b::HyperRectangle{N,Float64},
) where {N}
    maxA = a.origin + a.widths
    minA = a.origin
    maxB = b.origin + b.widths
    minB = b.origin
    return (all(maxA .>= minB) && all(maxB .> minA))

end

function inittree2d(n)
    empty::Array{Int64,1} = []
    oct = Cell(SVector(-1.0, -1), SVector(2.0, 2), empty)
    for _ = 1:n
        for leaf in allleaves(oct)
            split!(leaf)
        end
    end
    for node in allcells(oct)
        node.data = copy(node.data)
    end
    return oct
end

function inserttotree!(
    box::HyperRectangle{N,Float64},
    tree::Cell,
    index::Int64,
) where {N}
    if (isleaf(tree))
        if (intersects(box, tree.boundary))
            push!(tree.data, index)
        end
    elseif (intersects(box, tree.boundary))
        for i = 1:length(tree.children)
            inserttotree!(box, tree[i], index)
        end
    end
    return nothing
end

function possiblevoxels(r::ray2d, tree::Cell)::Vector{Int64}
    if (isleaf(tree) && intersects(r, tree.boundary) > 0.0)
        return tree.data

    elseif (~(tree.parent === nothing) && intersects(r, tree.boundary) > 0.0)
        N = length(tree.children)
        v = Vector{Int64}()
        for i = 1:N
            append!(v, possiblevoxels(r, tree[i]))
        end
        return v

    elseif ((tree.parent === nothing) && intersects(r, tree.boundary) > 0.0)
        N = length(tree.children)
        v = Vector{Int64}()
        for i = 1:N
            append!(v, possiblevoxels(r, tree[i]))
        end
        return collect(Set(v))

    else
        return []
    end

end

function elementray(r::ray2d, vv::Vector{triangle2d}, checklist::Vector{Int64})
    N = length(checklist)
    tlens = Vector{MArray{Tuple{2},MArray{Tuple{2},Float64,1,2},1,2}}()
    indices = Vector{Int64}()
    for i = 1:N
        q = intersects(r, vv[checklist[i]])
        if (~isnothing(q))
            push!(tlens, q)
            push!(indices, checklist[i])
        end
    end
    return (indices, tlens)
end

function constructmatrix(tree::Cell, vv::Vector{triangle2d}, N::Int64, theta)
    r = 1.0 # Radius of the transform.
    Na = length(theta)
    if (~isa(theta, Union{SArray,Array}))
        theta = [theta]
    end
    if (N > 1)
        span = range(r, stop = -r, length = N)
    else
        span = [0.0]
    end
    Ntotal = Na * N
    rows = Vector{Int64}()
    cols = Vector{Int64}()
    vals = Vector{Float64}()
    Nofnodes = length(Set(vcat([i.nodes for i in vv]...)))

    Nth = Threads.nthreads()
    rows = Vector{Vector{Int64}}(undef, Nth)
    cols = Vector{Vector{Int64}}(undef, Nth)
    vals = Vector{Vector{Float64}}(undef, Nth)
    for p = 1:Nth
        rows[p] = []
        cols[p] = []
        vals[p] = []
    end

    t = time()
    if (Na > 1)
        pb = ProgressBar(1:Na)
    else
        pb = 1:1
    end
    for a in pb
        dir = [cos(theta[a]), sin(theta[a])]
        aux = [-sin(theta[a]), cos(theta[a])]
        Threads.@threads for i = 1:N
            rayindex = (a - 1) * N + i
            p1 = dir + aux * span[i]
            p2 = p1 - 2 * dir * r
            or = ray2d(p1, p2)
            checklist = possiblevoxels(or, tree) #Possible elements.
            indices, tlens = elementray(or, vv, checklist)
            Nel = length(indices)
            for el = 1:Nel
                nodes = vv[indices[el]].nodes
                integrals = elementintegral(tlens[el], vv[indices[el]])
                #
                # append!(rows,rayindex*ones(Int64,3))
                # append!(cols,nodes)
                # append!(vals,integrals)

                append!(rows[Threads.threadid()], rayindex * ones(Int64, 3))
                append!(cols[Threads.threadid()], nodes)
                append!(vals[Threads.threadid()], integrals)

            end
        end
    end
    rows = vcat(rows...)
    cols = vcat(cols...)
    vals = vcat(vals...)
    M = sparse(rows, cols, vals, Ntotal, Nofnodes)
    println("Matrix constructed in ", time() - t, " seconds.")
    return M
end

function sub2ind(sizes, i::Int64, j::Int64, k::Int64)
    @assert i > 0 && i <= sizes[1]
    @assert j > 0 && j <= sizes[2]
    @assert k > 0 && k <= sizes[3]
    return (k - 1) * sizes[1] * sizes[2] + sizes[1] * (j - 1) + i
end

@inbounds function intersects(
    r::ray2d,
    t::triangle2d,
)::Union{MArray{Tuple{2},MArray{Tuple{2},Float64,1,2},1,2},Nothing}
    v = r.dir
    ips = @MVector [@MVector[NaN, NaN], @MVector[NaN, NaN]]
    cnt::Int64 = 0
    for i = 1:3
        j = mod(i, 3) + 1
        u = @SVector[t.pts[j][1] - t.pts[i][1], t.pts[j][2] - t.pts[i][2]]
        c = u[1] * v[2] - u[2] * v[1]
        a = @SVector[t.pts[i][1] - r.origin[1], t.pts[i][2] - r.origin[2]]
        b = a[1] * v[2] - a[2] * v[1]
        if (c == 0 && b == 0)
            println("Ray coincides with an edge. Ray is perturbed to avoid counting the integral twice. ")
            rnew = ray2ddir(r.origin + 1e-10 .* [-v[2], v[1]], r.dir)
            return intersects(rnew, t)
        end

        k = -b / c

        if (0 < k <= 1.0)
            cnt = cnt + 1
            ips[cnt][1] = t.pts[i][1] + u[1] * k
            ips[cnt][2] = t.pts[i][2] + u[2] * k

        end

    end
    if (cnt >= 2)
        return ips
    else
        return nothing

    end
end

@inline @inbounds function  gaussq(fvals,edge)::Float64
    val = fvals[1]*0.555555556 + fvals[2]*0.888888889 + 0.555555556*fvals[3];
    val = val*sqrt(((edge[1][1]-edge[2][1])/2 )^2 + ((edge[1][2]-edge[2][2])/2 )^2);
    return val
end

@inline @inbounds function localtoglobal(s, edge)
    Ns = length(s)
    z = zeros(MMatrix{Ns,2,Float64})
    for i = 1:Ns
        z[i, 1] =
            s[i] * (edge[2][1] - edge[1][1]) / 2 + (edge[1][1] + edge[2][1]) / 2
        z[i, 2] =
            s[i] * (edge[2][2] - edge[1][2]) / 2 + (edge[1][2] + edge[2][2]) / 2
    end
    return z
end

@inbounds function shapefun(x,n1,n2,A::Float64)
    val = 1/(2*A).*(n1[1].*n2[2]-n2[1].*n1[2] .+ (n1[2] - n2[2]).*x[:,1] + (n2[1] - n1[1]).*x[:,2] );
    return val
end

 @inbounds function shapefun!(val,x,n1,n2,A::Float64)
   val .= 1/(2*A).*(n1[1].*n2[2]-n2[1].*n1[2] .+ (n1[2] - n2[2]).*x[:,1] + (n2[1] - n1[1]).*x[:,2] );
   @assert maximum(val) <= 1.00001 && minimum(val) >= -0.00001
   # We do not want use a shapefunction out of its support.
end

function elementintegral(edge,t::triangle2d)
    v = @MVector[0.0,0.0,0.0]
    qp = localtoglobal(@SVector[-0.774596669,0,0.774596669],edge)
    fv = @MVector [0.0,0.0,0.0]
    for i = 1:3
        n1 = mod(i,3) + 1;
        n2 = mod(i+1,3) + 1;
        shapefun!(fv,qp,t.pts[n1],t.pts[n2],t.A)
        v[i] = gaussq(fv,edge)
    end

    return v
end

function nodeslookup(filename)
    file = matopen(filename)
    p=read(file,"p")
    e=read(file,"e")
    t=read(file,"t")
    t = Array{Int64}(t)
    quad = inittree2d(7) # Do 7 divisions for quadtree.
    Nelem = size(t)[2]
    trianglevector = Vector{triangle2d}(undef,Nelem)
    tt = time()
    for i =1:Nelem
        trianglevector[i] = triangle2d(p[:,t[1,i]],p[:,t[2,i]],p[:,t[3,i]],t[1:3,i])
        inserttotree!(trianglevector[i].bbox, quad, i)
    end
    println("Quadtree initialized in " ,time()-tt ," seconds.")

    return quad,trianglevector
end

function test()
    fn = "phantom.mat"
    quad, trianglevector = nodeslookup(fn)
    Npr = 400
    Na = 100
    M = constructmatrix(
        quad,
        trianglevector,
        Npr,
        Array(range(pi / 2, stop = 3 / 2 * pi, length = Na)),
    )
    file = matopen(fn)
    zn = read(file, "zn")
    zn = zn'
    y = M * zn
    rad = reshape(y, Npr, Na)
    imshow(rad, aspect = "auto")
    return nothing
end

test()

max_stars_repo_name: suurj/radon

max_issues_repo_name: suurj/radon

max_forks_repo_name: suurj/radon


text:
import Libdl: dlext, find_library
# FIXME: THis should all probably be in some kind of __init__ (see Issue #2)
# TODO: handle missing values gracefully

const ROS_VERSION = VersionNumber(ENV["ROS_VERSION"])
const ROS_DISTRO = ENV["ROS_DISTRO"]

@assert ROS_VERSION >= v"2.0" "Only ROS 2.0 is supported by this package."
#@assert ROS_DISTRO == "ardent"  # the only available version of ROS2 for now
@assert ROS_DISTRO == "galactic"  # the only available version of ROS2 for now

const AMENT_PREFIX_PATH = ENV["AMENT_PREFIX_PATH"]
const ROS_INCLUDE_PATH = joinpath(AMENT_PREFIX_PATH, "include")
const ROS_LIB_PATH = joinpath(AMENT_PREFIX_PATH, "lib")

const WRAPPED_ROS2_LIBS = ["rcutils", "rosidl_generator_c", "rmw", "rcl"]

# dynamically load the ROS2 libraries
for lib_name in WRAPPED_ROS2_LIBS
    lib_var = Symbol("lib", lib_name)
    lib_fname = "lib$lib_name.$dlext"
    @info "Loading library $lib_name"
    @eval const $lib_var = Symbol(find_library([$lib_fname]))
end


text:
using Test
using NQCDynamics
using NQCDynamics: Calculators

@testset "fill_density!" begin
    @testset "SingleState" for state ∈ (1, 2)
        electronics = SingleState(state)
        density = zeros(2,2)
        NonadiabaticDistributions.fill_density!(density, electronics)
        @test density[state, state] == 1
    end

    @testset "ElectronicPopulation" begin
        population = [0.5, 0.5]
        electronics = ElectronicPopulation(population)
        density = zeros(2, 2)
        NonadiabaticDistributions.fill_density!(density, electronics)
        @test density[1, 1] == density[2, 2] == 0.5
    end
end

@testset "initialise_density_matrix" begin
    calc = Calculators.Calculator(DoubleWell(), 1, Float64)
    electronics = SingleState(1)
    density = NonadiabaticDistributions.initialise_density_matrix(electronics, calc)
    @test density == [1 0; 0 0]
end

rs = (fill(0.0, 1, 1), fill(0.0, 1, 1, 1))
calculators = (
    Calculators.Calculator(DoubleWell(), 1, Float64),
    Calculators.Calculator(DoubleWell(), 1, 1, Float64),
    )
@testset "transform_density!" for (r, calc) ∈ zip(rs, calculators)
    density = [1.0 0; 0 0]
    NonadiabaticDistributions.transform_density!(density, calc, r, :to_adiabatic)
    @test density ≈ [0.5 0.5; 0.5 0.5]
    NonadiabaticDistributions.transform_density!(density, calc, r, :to_diabatic)
    @test density ≈ [1.0 0; 0 0]
    @test_throws ArgumentError NonadiabaticDistributions.transform_density!(density, calc, r, :blah)
end

@testset "initialise_adiabatic_density_matrix" for (r, calc) ∈ zip(rs, calculators)
    @testset "Adiabatic" begin
        electronics = SingleState(1, Adiabatic())
        density = NonadiabaticDistributions.initialise_adiabatic_density_matrix(electronics, calc, r)
        @test density[1,1] == 1
    end

    @testset "Diabatic" begin
        electronics = SingleState(1, Diabatic())
        density = NonadiabaticDistributions.initialise_adiabatic_density_matrix(electronics, calc, r)
        @test all(density .≈ 0.5)
    end
end


text:
# HEMM / SYMM
#

macro blis_linalg_hemm(typename, targetfunc, bliname)

    # Get method for the typed API backend.
    blifuncname = Symbol("bli_", cblas_typechar[ eval(typename) ], bliname)
    blifunc = getproperty(TypedBackend, blifuncname)

    return quote
        """
        BLIS-based HEMM/SYMM with approximately the same interface as `LinearAlgebra.BLAS`.
        Strided arrays are directly supported without further Julia-level dispatching.
        """
        $(esc(targetfunc))(side::AbstractChar,
                           uplo::AbstractChar,
                           alpha::$typename,
                           A::StridedMatrix{$typename},
                           B::StridedMatrix{$typename},
                           beta::$typename,
                           C::StridedMatrix{$typename}) = begin

            bli_sidea = char_to_side[side]
            bli_uploa = char_to_uplo[uplo]

            k , k_= size(A)
            m , n = size(B)
            m_, n_= size(C)

            (k == k_) || throw(DimensionMismatch("HE(SY)MM A storage is not square."))
            (m == m_ &&
             n == n_) || throw(DimensionMismatch("HE(SY)MM target buffer size mismatch."))
            (bli_sidea != BLIS_LEFT  ||
             k == m ) || throw(DimensionMismatch("HE(SY)MM contracted size mismatch."))
            (bli_sidea != BLIS_RIGHT ||
             k == n ) || throw(DimensionMismatch("HE(SY)MM contracted size mismatch."))

            $blifunc(bli_sidea,
                     bli_uploa,
                     BLIS_NO_CONJUGATE,
                     BLIS_NO_TRANSPOSE,
                     m, n,
                     [alpha],
                     A, strides(A)...,
                     B, strides(B)...,
                     [beta],
                     C, strides(C)...)
            C

        end
    end
end

@blis_linalg_hemm(Float32, hemm!, hemm!)
@blis_linalg_hemm(Float64, hemm!, hemm!)
@blis_linalg_hemm(ComplexF32, hemm!, hemm!)
@blis_linalg_hemm(ComplexF64, hemm!, hemm!)

@blis_linalg_hemm(Float32, symm!, symm!)
@blis_linalg_hemm(Float64, symm!, symm!)
@blis_linalg_hemm(ComplexF32, symm!, symm!)
@blis_linalg_hemm(ComplexF64, symm!, symm!)


text:
module BivarNormalProb

export bvn, bvnu

include("bvn.jl")

end # module


text:
using LazySets.Approximations: project
export is_intersection_empty,
       isdisjoint

# --- disjointness check for lazy sets ---

"""
    is_intersection_empty(X::LazySet{N},
                          Y::LazySet{N},
                          witness::Bool=false
                          )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}

Check whether two sets do not intersect, and otherwise optionally compute a
witness.

### Input

- `X`       -- set
- `Y`       -- another set
- `witness` -- (optional, default: `false`) compute a witness if activated

### Output

* If `witness` option is deactivated: `true` iff ``X ∩ Y = ∅``
* If `witness` option is activated:
  * `(true, [])` iff ``X ∩ Y = ∅``
  * `(false, v)` iff ``X ∩ Y ≠ ∅`` and ``v ∈ X ∩ Y``

### Algorithm

This is a fallback implementation that computes the concrete intersection,
`intersection`, of the given sets.

A witness is constructed using the `an_element` implementation of the result.
"""
function is_intersection_empty(X::LazySet{N},
                               Y::LazySet{N},
                               witness::Bool=false
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
    cap = intersection(X, Y)
    empty_intersection = isempty(cap)
    if witness
        if empty_intersection
            return (true, N[])
        else
            return (false, an_element(cap))
        end
    end
    return empty_intersection
end

# --- alias ---

"""
    isdisjoint(X, Y)

An alternative name for `is_intersection_empty(X, Y)`.
"""
const isdisjoint = is_intersection_empty

# --- AbstractHyperrectangle ---

"""
    is_intersection_empty(H1::AbstractHyperrectangle{N},
                          H2::AbstractHyperrectangle{N},
                          witness::Bool=false
                         )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}

Check whether two hyperrectangles do not intersect, and otherwise optionally
compute a witness.

### Input

- `H1` -- first hyperrectangle
- `H2` -- second hyperrectangle
- `witness` -- (optional, default: `false`) compute a witness if activated

### Output

* If `witness` option is deactivated: `true` iff ``H1 ∩ H2 = ∅``
* If `witness` option is activated:
  * `(true, [])` iff ``H1 ∩ H2 = ∅``
  * `(false, v)` iff ``H1 ∩ H2 ≠ ∅`` and ``v ∈ H1 ∩ H2``

### Algorithm

``H1 ∩ H2 ≠ ∅`` iff ``|c_2 - c_1| ≤ r_1 + r_2``, where ``≤`` is taken
component-wise.

A witness is computed by starting in one center and moving toward the other
center for as long as the minimum of the radius and the center distance.
In other words, the witness is the point in `H1` that is closest to the center
of `H2`.
"""
function is_intersection_empty(H1::AbstractHyperrectangle{N},
                               H2::AbstractHyperrectangle{N},
                               witness::Bool=false
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
    empty_intersection = false
    center_diff = center(H2) - center(H1)
    @inbounds for i in eachindex(center_diff)
        if abs(center_diff[i]) >
                radius_hyperrectangle(H1, i) + radius_hyperrectangle(H2, i)
            empty_intersection = true
            break
        end
    end

    if !witness
        return empty_intersection
    elseif empty_intersection
        return (true, N[])
    end

    # compute a witness 'v' in the intersection
    v = copy(center(H1))
    c2 = center(H2)
    @inbounds for i in eachindex(center_diff)
        if v[i] <= c2[i]
            # second center is right of first center
            v[i] += min(radius_hyperrectangle(H1, i), center_diff[i])
        else
            # second center is left of first center
            v[i] -= min(radius_hyperrectangle(H1, i), -center_diff[i])
        end
    end
    return (false, v)
end

# --- disjointness check for 1D intervals ---

"""
    is_intersection_empty(I1::Interval{N}, I2::Interval{N},
                          witness::Bool=false) where {N<:Real}

Check whether two intervals do not intersect, and otherwise optionally
compute a witness.

### Input

- `I1`      -- first interval
- `I2`      -- second interval
- `witness` -- (optional, default: `false`) compute a witness if activated

### Output

* If `witness` option is deactivated: `true` iff ``I1 ∩ I2 = ∅``
* If `witness` option is activated:
  * `(true, [])` iff ``I1 ∩ I2 = ∅``
  * `(false, v)` iff ``I1 ∩ I2 ≠ ∅`` and ``v ∈ I1 ∩ I2``

### Algorithm

``I1 ∩ I2 ≠ ∅`` iff there is a gap between the left-most point of the second
interval and the left-most point of the first interval, or vice-versa.

A witness is computed by taking the maximum over the left-most points of each
interval, which is guaranteed to belong to the intersection.
"""
function is_intersection_empty(I1::Interval{N}, I2::Interval{N},
                               witness::Bool=false) where {N<:Real}
    check = min(I2) > max(I1) || min(I1) > max(I2)
    if !witness
        return check
    else
        return (check, [max(min(I1), min(I2))])
    end
end

# --- AbstractSingleton ---

# common code for singletons
@inline function is_intersection_empty_helper_singleton(
        S::AbstractSingleton{N}, X::LazySet{N}, witness::Bool=false
       )::Union{Bool, Tuple{Bool, Vector{N}} } where {N<:Real}
    empty_intersection = element(S) ∉ X
    if witness
        return (empty_intersection, empty_intersection ? N[] : element(S))
    else
        return empty_intersection
    end
end

"""
    is_intersection_empty(X::LazySet{N},
                          S::AbstractSingleton{N},
                          witness::Bool=false
                         )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}

Check whether a convex set and a singleton do not intersect, and otherwise
optionally compute a witness.

### Input

- `X`       -- convex set
- `S`       -- singleton
- `witness` -- (optional, default: `false`) compute a witness if activated

### Output

* If `witness` option is deactivated: `true` iff ``S ∩ X = ∅``
* If `witness` option is activated:
  * `(true, [])` iff ``S ∩ X = ∅``
  * `(false, v)` iff ``S ∩ X ≠ ∅`` and
    `v` = `element(S)` ``∈ S ∩ X``

### Algorithm

``S ∩ X = ∅`` iff `element(S)` ``∉ X``.
"""
function is_intersection_empty(X::LazySet{N},
                               S::AbstractSingleton{N},
                               witness::Bool=false
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
    return is_intersection_empty_helper_singleton(S, X, witness)
end

# symmetric method
function is_intersection_empty(S::AbstractSingleton{N},
                               X::LazySet{N},
                               witness::Bool=false
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
    return is_intersection_empty_helper_singleton(S, X, witness)
end

"""
    is_intersection_empty(S1::AbstractSingleton{N},
                          S2::AbstractSingleton{N},
                          witness::Bool=false
                         )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}

Check whether two singletons do not intersect, and otherwise optionally compute
a witness.

### Input

- `S1` -- first singleton
- `S2` -- second singleton
- `witness` -- (optional, default: `false`) compute a witness if activated

### Output

* If `witness` option is deactivated: `true` iff ``S1 ∩ S2 = ∅``
* If `witness` option is activated:
  * `(true, [])` iff ``S1 ∩ S2 = ∅``
  * `(false, v)` iff ``S1 ∩ S2 ≠ ∅`` and `v` = `element(S1)` ``∈ S1 ∩ S2``

### Algorithm

``S1 ∩ S2 = ∅`` iff ``S1 ≠ S2``.
"""
function is_intersection_empty(S1::AbstractSingleton{N},
                               S2::AbstractSingleton{N},
                               witness::Bool=false
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
    empty_intersection = element(S1) != element(S2)
    if witness
        return (empty_intersection, empty_intersection ? N[] : element(S1))
    else
        return empty_intersection
    end
end

"""
    is_intersection_empty(H::AbstractHyperrectangle{N},
                          S::AbstractSingleton{N},
                          witness::Bool=false
                         )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}

Check whether a hyperrectangle and a singleton do not intersect, and otherwise
optionally compute a witness.

### Input

- `H` -- hyperrectangle
- `S` -- singleton
- `witness` -- (optional, default: `false`) compute a witness if activated

### Output

* If `witness` option is deactivated: `true` iff ``H ∩ S = ∅``
* If `witness` option is activated:
  * `(true, [])` iff ``H ∩ S = ∅``
  * `(false, v)` iff ``H ∩ S ≠ ∅`` and `v` = `element(S)` ``∈ H ∩ S``

### Algorithm

``H ∩ S = ∅`` iff `element(S)` ``∉ H``.
"""
function is_intersection_empty(H::AbstractHyperrectangle{N},
                               S::AbstractSingleton{N},
                               witness::Bool=false
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
    return is_intersection_empty_helper_singleton(S, H, witness)
end

# symmetric method
function is_intersection_empty(S::AbstractSingleton{N},
                               H::AbstractHyperrectangle{N},
                               witness::Bool=false
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
    return is_intersection_empty_helper_singleton(S, H, witness)
end

# --- Ball2 ---

"""
    is_intersection_empty(B1::Ball2{N},
                          B2::Ball2{N},
                          witness::Bool=false
                         )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:AbstractFloat}

Check whether two balls in the 2-norm do not intersect, and otherwise optionally
compute a witness.

### Input

- `B1` -- first ball in the 2-norm
- `B2` -- second ball in the 2-norm
- `witness` -- (optional, default: `false`) compute a witness if activated

### Output

* If `witness` option is deactivated: `true` iff ``B1 ∩ B2 = ∅``
* If `witness` option is activated:
  * `(true, [])` iff ``B1 ∩ B2 = ∅``
  * `(false, v)` iff ``B1 ∩ B2 ≠ ∅`` and ``v ∈ B1 ∩ B2``

### Algorithm

``B1 ∩ B2 = ∅`` iff ``‖ c_2 - c_1 ‖_2 > r_1 + r_2``.

A witness is computed depending on the smaller/bigger ball (to break ties,
choose `B1` for the smaller ball) as follows.
- If the smaller ball's center is contained in the bigger ball, we return it.
- Otherwise start in the smaller ball's center and move toward the other center
  until hitting the smaller ball's border.
  In other words, the witness is the point in the smaller ball that is closest
  to the center of the bigger ball.
"""
function is_intersection_empty(B1::Ball2{N},
                               B2::Ball2{N},
                               witness::Bool=false
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:AbstractFloat}
    center_diff_normed = norm(center(B2) - center(B1), 2)
    empty_intersection = center_diff_normed > B1.radius + B2.radius

    if !witness
        return empty_intersection
    elseif empty_intersection
        return (true, N[])
    end

    # compute a witness 'v' in the intersection
    if B1.radius <= B2.radius
        smaller = B1
        bigger = B2
    else
        smaller = B2
        bigger = B1
    end
    if center_diff_normed <= bigger.radius
        # smaller ball's center is contained in bigger ball
        v = smaller.center
    else
        # scale center difference with smaller ball's radius
        direction = (bigger.center - smaller.center)
        v = smaller.center + direction / center_diff_normed * smaller.radius
    end
    return (false, v)
end

# --- Zonotope ---

"""
    is_intersection_empty(Z::Zonotope{N}, H::Hyperplane{N}, witness::Bool=false
                         )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}

Check whether a zonotope and a hyperplane do not intersect, and otherwise
optionally compute a witness.

### Input

- `Z` -- zonotope
- `H` -- hyperplane
- `witness` -- (optional, default: `false`) compute a witness if activated

### Output

* If `witness` option is deactivated: `true` iff ``Z ∩ H = ∅``
* If `witness` option is activated:
  * `(true, [])` iff ``Z ∩ H = ∅``
  * `(false, v)` iff ``Z ∩ H ≠ ∅`` and ``v ∈ Z ∩ H``

### Algorithm

``Z ∩ H = ∅`` iff ``(b - a⋅c) ∉ \\left[ ± ∑_{i=1}^p |a⋅g_i| \\right]``,
where ``a``, ``b`` are the hyperplane coefficients, ``c`` is the zonotope's
center, and ``g_i`` are the zonotope's generators.

For witness production we fall back to a less efficient implementation for
general sets as the first argument.
"""
function is_intersection_empty(Z::Zonotope{N},
                               H::Hyperplane{N},
                               witness::Bool=false
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
    if witness
        # use less efficient implementation that supports witness production
        return invoke(is_intersection_empty,
                      Tuple{LazySet{N}, Hyperplane{N}, Bool},
                      Z, H, witness)
    end

    v = H.b - dot(H.a, Z.center)
    abs_sum = sum(abs(dot(H.a, Z.generators[:, i])) for i = 1:ngens(Z))
    return v < -abs_sum || v > abs_sum
end

# symmetric method
function is_intersection_empty(H::Hyperplane{N},
                               Z::Zonotope{N},
                               witness::Bool=false
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
    return is_intersection_empty(Z, H, witness)
end

"""
    is_intersection_empty(Z1::Zonotope{N}, Z2::Zonotope{N}, witness::Bool=false
                         ) where {N<:Real}

Check whether two zonotopes do not intersect, and otherwise optionally compute a
witness.

### Input

- `Z1`      -- zonotope
- `Z2`      -- zonotope
- `witness` -- (optional, default: `false`) compute a witness if activated

### Output

* If `witness` option is deactivated: `true` iff ``Z1 ∩ Z2 = ∅``
* If `witness` option is activated:
  * `(true, [])` iff ``Z1 ∩ Z2 = ∅``
  * `(false, v)` iff ``Z1 ∩ Z2 ≠ ∅`` and ``v ∈ Z1 ∩ Z2``

### Algorithm

``Z1 ∩ Z2 ≠ ∅`` iff ``c_1 - c_2 ∈ Z(0, (g_1, g_2))`` where ``c_i`` and ``g_i``
are the center and generators of zonotope `Zi` and ``Z(c, g)`` represents the
zonotope with center ``c`` and generators ``g``.
"""
function is_intersection_empty(Z1::Zonotope{N}, Z2::Zonotope{N},
                               witness::Bool=false) where {N<:Real}
    n = dim(Z1)
    @assert n == dim(Z2) "zonotopes need to have the same dimensions"
    Z = Zonotope(zeros(N, n), hcat(Z1.generators, Z2.generators))
    result = (center(Z1) - center(Z2)) ∈ Z
    if result
        return witness ? (true, N[]) : true
    elseif witness
        error("witness production is not supported yet")
    else
        return false
    end
end

"""
    is_intersection_empty(ls1::LineSegment{N},
                          ls2::LineSegment{N},
                          witness::Bool=false
                         )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}

Check whether two line segments do not intersect, and otherwise optionally
compute a witness.

### Input

- `ls1` -- first line segment
- `ls2` -- second line segment
- `witness` -- (optional, default: `false`) compute a witness if activated

### Output

* If `witness` option is deactivated: `true` iff ``ls1 ∩ ls2 = ∅``
* If `witness` option is activated:
  * `(true, [])` iff ``ls1 ∩ ls2 = ∅``
  * `(false, v)` iff ``ls1 ∩ ls2 ≠ ∅`` and ``v ∈ ls1 ∩ ls2``

### Algorithm

The algorithm is inspired from [here](https://stackoverflow.com/a/565282), which
again is the special 2D case of a 3D algorithm by Ronald Goldman's article on the
*Intersection of two lines in three-space* in Graphics Gems, Andrew S. (ed.), 1990.

We first check if the two line segments are parallel, and if so, if they are
collinear.
In the latter case, we check containment of any of the end points in the other
line segment.
Otherwise the lines are not parallel, so we can solve an equation of the
intersection point, if it exists.
"""
function is_intersection_empty(ls1::LineSegment{N},
                               ls2::LineSegment{N},
                               witness::Bool=false
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
    function cross(ls1::Vector{N}, ls2::Vector{N})::N where {N<:Real}
        return ls1[1] * ls2[2] - ls1[2] * ls2[1]
    end

    r = ls1.q - ls1.p
    if iszero(r)
        # first line segment is a point
        empty_intersection = ls1.q ∉ ls2
        if witness
            return (empty_intersection, empty_intersection ? N[] : ls1.q)
        else
            return empty_intersection
        end
    end

    s = ls2.q - ls2.p
    if iszero(s)
        # second line segment is a point
        empty_intersection = ls2.q ∉ ls1
        if witness
            return (empty_intersection, empty_intersection ? N[] : ls2.q)
        else
            return empty_intersection
        end
    end

    p1p2 = ls2.p - ls1.p
    u_numerator = cross(p1p2, r)
    u_denominator = cross(r, s)

    if u_denominator == 0
        # line segments are parallel
        if u_numerator == 0
            # line segments are collinear
            if ls1.p ∈ ls2
                empty_intersection = false
                if witness
                    v = ls1.p
                end
            elseif ls1.q ∈ ls2
                empty_intersection = false
                if witness
                    v = ls1.q
                end
            else
                empty_intersection = true
            end
        else
            # line segments are parallel and not collinear
            empty_intersection = true
        end
    else
        # line segments are not parallel
        u = u_numerator / u_denominator
        if u < 0 || u > 1
            empty_intersection = true
        else
            t = cross(p1p2, s) / u_denominator
            empty_intersection = t < 0 || t > 1
            if witness
                v = ls1.p + t * r
            end
        end
    end

    if witness
        return (empty_intersection, empty_intersection ? N[] : v)
    else
        return empty_intersection
    end
end

# --- Hyperplane ---

@inline function is_intersection_empty_helper_hyperplane(
        hp::Union{Hyperplane{N}, Line{N}},
        X::LazySet{N},
        witness::Bool=false
       )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
    normal_hp = hp.a
    sv_left = σ(-normal_hp, X)
    if -dot(sv_left, -normal_hp) <= hp.b
        sv_right = σ(normal_hp, X)
        empty_intersection = (hp.b > dot(sv_right, normal_hp))
    else
        empty_intersection = true
    end
    if witness
        if empty_intersection
            v = N[]
        else
            point_hp = an_element(hp)
            point_line = sv_left
            dir_line = sv_right - sv_left
            d = dot((point_hp - point_line), normal_hp) /
                dot(dir_line, normal_hp)
            v = d * dir_line + point_line
        end
        return (empty_intersection, v)
    else
        return empty_intersection
    end
end

"""
    is_intersection_empty(X::LazySet{N},
                          hp::Union{Hyperplane{N}, Line{N}},
                          [witness]::Bool=false
                         )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}

Check whether a compact set an a hyperplane do not intersect, and otherwise
optionally compute a witness.

### Input

- `X`       -- compact set
- `hp`      -- hyperplane
- `witness` -- (optional, default: `false`) compute a witness if activated

### Output

* If `witness` option is deactivated: `true` iff ``X ∩ hp = ∅``
* If `witness` option is activated:
  * `(true, [])` iff ``X ∩ hp = ∅``
  * `(false, v)` iff ``X ∩ hp ≠ ∅`` and ``v ∈ X ∩ hp``

### Notes

We assume that `X` is compact.
Otherwise, the support vector queries may fail.

### Algorithm

A compact convex set intersects with a hyperplane iff the support function in
the negative resp. positive direction of the hyperplane's normal vector ``a`` is
to the left resp. right of the hyperplane's constraint ``b``:

```math
-ρ(-a) ≤ b ≤ ρ(a)
```

For witness generation, we compute a line connecting the support vectors to the
left and right, and then take the intersection of the line with the hyperplane.
We follow
[this algorithm](https://en.wikipedia.org/wiki/Line%E2%80%93plane_intersection#Algebraic_form)
for the line-hyperplane intersection.
"""
function is_intersection_empty(X::LazySet{N},
                               hp::Union{Hyperplane{N}, Line{N}},
                               witness::Bool=false
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
    return is_intersection_empty_helper_hyperplane(hp, X, witness)
end

# symmetric function
function is_intersection_empty(hp::Union{Hyperplane{N}, Line{N}},
                               X::LazySet{N},
                               witness::Bool=false
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
    return is_intersection_empty_helper_hyperplane(hp, X, witness)
end

# disambiguation
function is_intersection_empty(hp1::Union{Hyperplane{N}, Line{N}},
                               hp2::Union{Hyperplane{N}, Line{N}},
                               witness::Bool=false
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
    return is_intersection_empty_helper_hyperplane(hp1, hp2, witness)
end

# disambiguation
function is_intersection_empty(hp::Union{Hyperplane{N}, Line{N}},
                               S::AbstractSingleton{N},
                               witness::Bool=false
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
    return is_intersection_empty_helper_singleton(S, hp, witness)
end

# symmetric method
function is_intersection_empty(S::AbstractSingleton{N},
                               hp::Union{Hyperplane{N}, Line{N}},
                               witness::Bool=false
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
    return is_intersection_empty_helper_singleton(S, hp, witness)
end

# --- HalfSpace ---

@inline function is_intersection_empty_helper_halfspace(
        hs::HalfSpace{N},
        X::LazySet{N},
        witness::Bool=false
       )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
    if !witness
        return -ρ(-hs.a, X) > hs.b
    end

    # for witness production, we compute the support vector instead
    svec = σ(-hs.a, X)
    empty_intersection = svec ∉ hs
    v = empty_intersection ? N[] : svec
    return (empty_intersection, v)
end

"""
    is_intersection_empty(X::LazySet{N},
                          hs::HalfSpace{N},
                          [witness]::Bool=false
                         )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}

Check whether a compact set an a half-space do not intersect, and otherwise
optionally compute a witness.

### Input

- `X`       -- compact set
- `hs`      -- half-space
- `witness` -- (optional, default: `false`) compute a witness if activated

### Output

* If `witness` option is deactivated: `true` iff ``X ∩ hs = ∅``
* If `witness` option is activated:
  * `(true, [])` iff ``X ∩ hs = ∅``
  * `(false, v)` iff ``X ∩ hs ≠ ∅`` and ``v ∈ X ∩ hs``

### Notes

We assume that `X` is compact.
Otherwise, the support vector queries may fail.

### Algorithm

A compact convex set intersects with a half-space iff the support vector in
the negative direction of the half-space's normal vector ``a`` is contained in
the half-space: ``σ(-a) ∈ hs``.
The support vector is thus also a witness.

Optional keyword arguments can be passed to the `ρ` function. In particular, if
`X` is a lazy intersection, options can be passed to the line search algorithm.
"""
function is_intersection_empty(X::LazySet{N},
                               hs::HalfSpace{N},
                               witness::Bool=false
                               )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
    return is_intersection_empty_helper_halfspace(hs, X, witness)
end

# symmetric function
function is_intersection_empty(hs::HalfSpace{N},
                               X::LazySet{N},
                               witness::Bool=false
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
    return is_intersection_empty_helper_halfspace(hs, X, witness)
end

"""
    is_intersection_empty(hs1::HalfSpace{N},
                          hs2::HalfSpace{N},
                          [witness]::Bool=false
                         )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}

Check whether two half-spaces do not intersect, and otherwise optionally compute
a witness.

### Input

- `hs1`     -- half-space
- `hs2`     -- half-space
- `witness` -- (optional, default: `false`) compute a witness if activated

### Output

* If `witness` option is deactivated: `true` iff ``hs1 ∩ hs2 = ∅``
* If `witness` option is activated:
  * `(true, [])` iff ``hs1 ∩ hs2 = ∅``
  * `(false, v)` iff ``hs1 ∩ hs2 ≠ ∅`` and ``v ∈ hs1 ∩ hs2``

### Algorithm

Two half-spaces do not intersect if and only if their normal vectors point in
the opposite direction and there is a gap between the two defining hyperplanes.

The latter can be checked as follows:
Let ``hs_1 : a_1⋅x = b_1`` and ``hs2 : a_2⋅x = b_2``.
Then we already know that ``a_2 = -k⋅a_1`` for some positive scaling factor
``k``.
Let ``x_1`` be a point on the defining hyperplane of ``hs_1``.
We construct a line segment from ``x_1`` to the point ``x_2`` on the defining
hyperplane of ``hs_2`` by shooting a ray from ``x_1`` with direction ``a_1``.
Thus we look for a factor ``s`` such that ``(x_1 + s⋅a_1)⋅a_2 = b_2``.
This gives us ``s = (b_2 - x_1⋅a_2) / (-k a_1⋅a_1)``.
The gap exists if and only if ``s`` is positive.

If the normal vectors do not point in opposite directions, then the defining
hyperplanes intersect and we can produce a witness as follows.
All points ``x`` in this intersection satisfy ``a_1⋅x = b_1`` and ``a_2⋅x = b_2``.
Thus we have ``(a_1 + a_2)⋅x = b_1+b_2``.
We now find a dimension where ``a_1 + a_2`` is non-zero, say, ``i``.
Then the result is a vector with one non-zero entry in dimension ``i``, defined
as ``[0, …, 0, (b_1 + b_2)/(a_1[i] + a_2[i]), 0, …, 0]``.
Such a dimension ``i`` always exists.
"""
function is_intersection_empty(hs1::HalfSpace{N},
                               hs2::HalfSpace{N},
                               witness::Bool=false
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
    a1 = hs1.a
    a2 = hs2.a
    issamedir, k = samedir(a1, -a2)
    if issamedir
        x1 = an_element(Hyperplane(a1, hs1.b))
        b2 = hs2.b
        s = (b2 - dot(x1, a2)) / (-k * dot(a1, a1))
        empty_intersection = s > 0
        if witness
            if empty_intersection
                v = N[]
            else
                # both defining hyperplanes are contained in each half-space
                v = x1
            end
        end
    else
        empty_intersection = false
        if witness
            v = zeros(N, length(a1))
            for i in 1:length(a1)
                a_sum_i = a1[i] + a2[i]
                if a_sum[i] != 0
                    v[i] = (hs1.b + hs2.b) / a_sum_i
                    break
                end
            end
        end
    end

    if !witness
        return empty_intersection
    else
        return (empty_intersection, v)
    end
end

# disambiguation
function is_intersection_empty(H::HalfSpace{N},
                               S::AbstractSingleton{N},
                               witness::Bool=false
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
    return is_intersection_empty_helper_singleton(S, H, witness)
end

# symmetric method
function is_intersection_empty(S::AbstractSingleton{N},
                               H::HalfSpace{N},
                               witness::Bool=false
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
    return is_intersection_empty_helper_singleton(S, H, witness)
end

# disambiguation
function is_intersection_empty(hp::Union{Hyperplane{N}, Line{N}},
                               hs::HalfSpace{N},
                               witness::Bool=false
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
    return is_intersection_empty_helper_halfspace(hs, hp, witness)
end

# symmetric method
function is_intersection_empty(hs::HalfSpace{N},
                               hp::Union{Hyperplane{N}, Line{N}},
                               witness::Bool=false
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
    return is_intersection_empty_helper_halfspace(hs, hp, witness)
end

# --- polyhedra ---

"""
    is_intersection_empty(P::Union{HPolyhedron{N}, AbstractPolytope{N}},
                          X::LazySet{N},
                          witness::Bool=false;
                          solver=GLPKSolverLP(method=:Simplex)
                         )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}

Check whether two polyhedra do not intersect.

### Input

- `P`         -- polyhedron
- `X`         -- another set (see the Notes section below)
- `witness`   -- (optional, default: `false`) compute a witness if activated
- `solver`    -- (optional, default: `GLPKSolverLP(method=:Simplex)`) LP solver
                 backend
- `algorithm` -- (optional, default: `"exact"`) algorithm keyword, one of:
                 * `"exact" (exact, uses a feasibility LP)
                 * `"sufficient" (sufficient, uses half-space checks)

### Output

* If `witness` option is deactivated: `true` iff ``P ∩ X = ∅``
* If `witness` option is activated:
  * `(true, [])` iff ``P ∩ X = ∅``
  * `(false, v)` iff ``P ∩ X ≠ ∅`` and ``v ∈ P ∩ X``

### Notes

For `algorithm == "exact"`, we assume that `constraints_list(X)` is defined.
For `algorithm == "sufficient"`, witness production is not supported.

### Algorithm

For `algorithm == "exact"`, see @ref(isempty(P::HPoly{N}, ::Bool)).

For `algorithm == "sufficient"`, we rely on the intersection check between the
set `X` and each constraint in `P`.
This means one support function evaluation of `X` for each constraint of `P`.
"""
function is_intersection_empty(P::Union{HPolyhedron{N}, AbstractPolytope{N}},
                               X::LazySet{N},
                               witness::Bool=false;
                               solver=GLPKSolverLP(method=:Simplex),
                               algorithm="exact"
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
    if algorithm == "sufficient"
        # sufficient check for empty intersection using half-space checks
        for Hi in constraints_list(P)
            if is_intersection_empty(X, Hi)
                if witness
                    return (true, N[])
                end
                return true
            end
        end
        if witness
            error("witness production is not supported yet")
        end
        return false
    elseif algorithm == "exact"
        # exact check for empty intersection using a feasibility LP
        return isempty(HPolyhedron([constraints_list(P); constraints_list(X)]),
                       witness; solver=solver)
    else
        error("algorithm $algorithm unknown")
    end
end

# symmetric method
function is_intersection_empty(X::LazySet{N},
                               P::Union{HPolyhedron{N}, AbstractPolytope{N}},
                               witness::Bool=false;
                               solver=GLPKSolverLP(method=:Simplex),
                               algorithm="exact"
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
    return is_intersection_empty(P, X, witness;
                                 solver=solver, algorithm=algorithm)
end

# disambiguation
function is_intersection_empty(P::Union{HPolyhedron{N}, AbstractPolytope{N}},
                               Q::Union{HPolyhedron{N}, AbstractPolytope{N}},
                               witness::Bool=false;
                               solver=GLPKSolverLP(method=:Simplex),
                               algorithm="exact"
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
    return invoke(is_intersection_empty,
                  Tuple{typeof(P), LazySet{N}, Bool},
                  P, Q, witness; solver=solver, algorithm=algorithm)
end

# disambiguation
function is_intersection_empty(P::Union{HPolyhedron{N}, AbstractPolytope{N}},
                               hs::HalfSpace{N},
                               witness::Bool=false
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
    return is_intersection_empty_helper_halfspace(hs, P, witness)
end

# symmetric method
function is_intersection_empty(hs::HalfSpace{N},
                               P::Union{HPolyhedron{N}, AbstractPolytope{N}},
                               witness::Bool=false
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
    return is_intersection_empty_helper_halfspace(hs, P, witness)
end

# disambiguation
function is_intersection_empty(P::Union{HPolyhedron{N}, AbstractPolytope{N}},
                               S::AbstractSingleton{N},
                               witness::Bool=false
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
    return is_intersection_empty_helper_singleton(S, P, witness)
end

# symmetric method
function is_intersection_empty(S::AbstractSingleton{N},
                               P::Union{HPolyhedron{N}, AbstractPolytope{N}},
                               witness::Bool=false
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
    return is_intersection_empty_helper_singleton(S, P, witness)
end

# disambiguation
function is_intersection_empty(P::Union{HPolyhedron{N}, AbstractPolytope{N}},
                               hp::Union{Hyperplane{N}, Line{N}},
                               witness::Bool=false;
                               solver=GLPKSolverLP(method=:Simplex),
                               algorithm="exact"
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
    return invoke(is_intersection_empty,
                  Tuple{typeof(P), LazySet{N}, Bool},
                  P, hp, witness, solver=solver, algorithm=algorithm)
end

# symmetric method
function is_intersection_empty(hp::Union{Hyperplane{N}, Line{N}},
                               P::Union{HPolyhedron{N}, AbstractPolytope{N}},
                               witness::Bool=false;
                               solver=GLPKSolverLP(method=:Simplex),
                               algorithm="exact"
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
    return invoke(is_intersection_empty,
                  Tuple{typeof(P), LazySet{N}, Bool},
                  P, hp, witness, solver=solver, algorithm=algorithm)
end

# --- union ---

"""
    is_intersection_empty(cup::UnionSet{N}, X::LazySet{N}, [witness]::Bool=false
                         )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}

Check whether a union of two convex sets and another set do not intersect.

### Input

- `cup` -- union of two convex sets
- `X`   -- another set

### Output

`true` iff ``\\text{cup} ∩ X = ∅``.
"""
function is_intersection_empty(cup::UnionSet{N},
                               X::LazySet{N},
                               witness::Bool=false
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where
                                  {N<:Real}
    return is_intersection_empty(UnionSetArray([cup.X, cup.Y]), X, witness)
end

# symmetric method
function is_intersection_empty(X::LazySet{N},
                               cup::UnionSet{N},
                               witness::Bool=false
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where
                                  {N<:Real}
    return is_intersection_empty(cup, X, witness)
end

# disambiguation
function is_intersection_empty(cup1::UnionSet{N},
                               cup2::UnionSet{N},
                               witness::Bool=false
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where
                                  {N<:Real}
    return is_intersection_empty(UnionSetArray([cup1.X, cup1.Y]),
                                 UnionSetArray([cup2.X, cup2.Y]), witness)
end

"""
    is_intersection_empty(cup::UnionSetArray{N},
                          X::LazySet{N},
                          [witness]::Bool=false
                         )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}

Check whether a union of a finite number of convex sets and another set do not
intersect.

### Input

- `cup` -- union of a finite number of convex sets
- `X`   -- another set

### Output

`true` iff ``\\text{cup} ∩ X = ∅``.
"""
function is_intersection_empty(cup::UnionSetArray{N},
                               X::LazySet{N},
                               witness::Bool=false
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where
                                  {N<:Real}
    result = true
    w = N[]
    for Y in array(cup)
        if witness
            result, w = is_intersection_empty(Y, X, witness)
        else
            result = is_intersection_empty(Y, X, witness)
        end
        if !result
            break
        end
    end
    return witness ? (result, w) : result
end

# symmetric method
function is_intersection_empty(X::LazySet{N},
                               cup::UnionSetArray{N},
                               witness::Bool=false
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where
                                  {N<:Real}
    return is_intersection_empty(cup, X, witness)
end

# disambiguation
function is_intersection_empty(cup1::UnionSet{N},
                               cup2::UnionSetArray{N},
                               witness::Bool=false
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where
                                  {N<:Real}
    return is_intersection_empty(UnionSetArray([cup1.X, cup1.Y]), cup2, witness)
end

# disambiguation
function is_intersection_empty(cup1::UnionSetArray{N},
                               cup2::UnionSet{N},
                               witness::Bool=false
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where
                                  {N<:Real}
    return is_intersection_empty(cup1, UnionSetArray([cup2.X, cup2.Y]), witness)
end

# disambiguation
function is_intersection_empty(cup1::UnionSetArray{N},
                               cup2::UnionSetArray{N},
                               witness::Bool=false
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where
                                  {N<:Real}
    result = true
    w = N[]
    for X in array(cup1)
        for Y in array(cup2)
            if witness
                result, w = is_intersection_empty(Y, X, witness)
            else
                result = is_intersection_empty(Y, X, witness)
            end
            if !result
                break
            end
        end
    end
    return witness ? (result, w) : result
end

# --- Universe ---

"""
    is_intersection_empty(U::Universe{N},
                          X::LazySet{N},
                          [witness]::Bool=false
                         )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}

Check whether a universe and another set do not intersect.

### Input

- `U` -- universe
- `X` -- another set

### Output

`true` iff ``X ≠ ∅``.
"""
function is_intersection_empty(U::Universe{N},
                               X::LazySet{N},
                               witness::Bool=false
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where
                                  {N<:Real}
    result = isempty(X)
    if result
        return witness ? (result, N[]) : result
    else
        return witness ? (result, an_element(X)) : result
    end
end

# symmetric method
function is_intersection_empty(X::LazySet{N},
                               U::Universe{N},
                               witness::Bool=false
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where
                                  {N<:Real}
    return is_intersection_empty(U, X, witness)
end

# disambiguation
function is_intersection_empty(U::Universe{N},
                               ::Universe{N},
                               witness::Bool=false
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where
                                  {N<:Real}
    return witness ? (false, an_element(U)) : false
end
function is_intersection_empty(P::AbstractPolyhedron{N},
                               U::Universe{N},
                               witness::Bool=false
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where
                                  {N<:Real}
    return invoke(is_intersection_empty,
                  Tuple{Universe{N}, LazySet{N}, Bool},
                  U, P, witness)
end
function is_intersection_empty(U::Universe{N},
                               P::Union{AbstractPolytope{N}, HPolyhedron{N}},
                               witness::Bool=false
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where
                                  {N<:Real}
    return invoke(is_intersection_empty,
                  Tuple{Universe{N}, LazySet{N}, Bool},
                  U, P, witness)
end
function is_intersection_empty(P::Union{AbstractPolytope{N}, HPolyhedron{N}},
                               U::Universe{N},
                               witness::Bool=false
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where
                                  {N<:Real}
    return invoke(is_intersection_empty,
                  Tuple{Universe{N}, LazySet{N}, Bool},
                  U, P, witness)
end
function is_intersection_empty(U::Universe{N},
                               P::AbstractPolytope{N},
                               witness::Bool=false
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where
                                  {N<:Real}
    return invoke(is_intersection_empty,
                  Tuple{Universe{N}, LazySet{N}, Bool},
                  U, P, witness)
end
function is_intersection_empty(S::AbstractSingleton{N},
                               U::Universe{N},
                               witness::Bool=false
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where
                                  {N<:Real}
    return invoke(is_intersection_empty,
                  Tuple{Universe{N}, LazySet{N}, Bool},
                  U, S, witness)
end
function is_intersection_empty(U::Universe{N},
                               S::AbstractSingleton{N},
                               witness::Bool=false
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where
                                  {N<:Real}
    return invoke(is_intersection_empty,
                  Tuple{Universe{N}, LazySet{N}, Bool},
                  U, S, witness)
end
function is_intersection_empty(hs::HalfSpace{N},
                               U::Universe{N},
                               witness::Bool=false
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where
                                  {N<:Real}
    return invoke(is_intersection_empty,
                  Tuple{Universe{N}, LazySet{N}, Bool},
                  U, hs, witness)
end
function is_intersection_empty(U::Universe{N},
                               hp::Union{Hyperplane{N}, Line{N}},
                               witness::Bool=false
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where
                                  {N<:Real}
    return invoke(is_intersection_empty,
                  Tuple{Universe{N}, LazySet{N}, Bool},
                  U, hs, witness)
end
function is_intersection_empty(hp::Union{Hyperplane{N}, Line{N}},
                               U::Universe{N},
                               witness::Bool=false
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where
                                  {N<:Real}
    return invoke(is_intersection_empty,
                  Tuple{Universe{N}, LazySet{N}, Bool},
                  U, hs, witness)
end
function is_intersection_empty(U::Universe{N},
                               hs::HalfSpace{N},
                               witness::Bool=false
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where
                                  {N<:Real}
    return invoke(is_intersection_empty,
                  Tuple{Universe{N}, LazySet{N}, Bool},
                  U, hs, witness)
end

# --- Complement ---

"""
    is_intersection_empty(C::Complement{N},
                          X::LazySet{N},
                          [witness]::Bool=false
                         )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}

Check whether the complement of a convex set and another set do not intersect.

### Input

- `C` -- complement of a convex set
- `X` -- convex set

### Output

* If `witness` option is deactivated: `true` iff ``X ∩ C = ∅``
* If `witness` option is activated:
  * `(true, [])` iff ``X ∩ C = ∅``
  * `(false, v)` iff ``X ∩ C ≠ ∅`` and ``v ∈ X ∩ C``

### Algorithm

We fall back to `X ⊆ C.X`, which can be justified as follows:

```math
    X ∩ Y^C = ∅ ⟺ X ⊆ Y
```
"""
function is_intersection_empty(C::Complement{N},
                               X::LazySet{N},
                               witness::Bool=false
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where
                                  {N<:Real}
    return ⊆(X, C.X, witness)
end

# symmetric method
function is_intersection_empty(X::LazySet{N},
                               C::Complement{N},
                               witness::Bool=false
                              )::Union{Bool, Tuple{Bool, Vector{N}}} where
                                  {N<:Real}
    return is_intersection_empty(C, X, witness)
end

"""
    is_intersection_empty(X::CartesianProductArray{N, S},
                          P::HPolyhedron{N}) where {N<:Real, S<:LazySet{N}}

Check whether a Cartesian product array intersects with a H-polyhedron.

### Input

- `X` -- Cartesian product array of convex sets
- `P` -- H-polyhedron

### Output

`true` iff ``X ∩ Y = ∅ ``.
"""
function is_intersection_empty(X::CartesianProductArray{N},
                               P::HPolyhedron{N}) where {N<:Real}
    cpa_low_dim, vars, block_structure = get_constrained_lowdimset(X, P)

    return isdisjoint(cpa_low_dim, project(P, vars))
end

# symmetric method
function is_intersection_empty(P::HPolyhedron{N},
                               X::CartesianProductArray{N, S}) where {N<:Real, S<:LazySet{N}}
        is_intersection_empty(X, P)
end

"""
    is_intersection_empty(X::CartesianProductArray{N}, Y::CartesianProductArray{N}) where {N}

Check whether two Cartesian products of a finite number of convex sets do not intersect.

### Input

- `X` -- Cartesian product array of convex sets
- `Y` -- Cartesian product array of convex sets

### Output

`true` iff ``X ∩ Y = ∅ ``.
"""
function is_intersection_empty(X::CartesianProductArray{N}, Y::CartesianProductArray{N}) where {N}
    @assert same_block_structure(array(X), array(Y)) "block structure has to be the same"

    for i in 1:length(X.array)
        if isdisjoint(X.array[i], Y.array[i])
            return true
        end
    end

    return false
end


text:
# Work-inefficient inclusive scan
# - uses shared memory to reduce
#
# Based on http://http.developer.nvidia.com/GPUGems3/gpugems3_ch39.html

using CUDAdrv, CUDAnative, CuArrays

function cpu_accumulate!(op::Function, data::Matrix{T}) where {T}
    cols = size(data,2)
    for col in 1:cols
        accum = zero(T)
        rows = size(data,1)
        for row in 1:size(data,1)
            accum = op(accum, data[row,col])
            data[row,col] = accum
        end
    end
end

function gpu_accumulate!(op::Function, data::CuDeviceMatrix{T}) where {T}
    col = blockIdx().x
    cols = gridDim().x

    row = threadIdx().x
    rows = blockDim().x

    if col <= cols && row <= rows
        shmem = @cuDynamicSharedMem(T, 2*rows)
        shmem[row] = data[row,col]
        sync_threads()

        # parallel reduction
        pin, pout = 1, 0
        offset = 1
        while offset < rows
            pout = 1 - pout
            pin = 1 - pin
            if row > offset
                shmem[pout * rows + row] =
                    op(shmem[pin * rows + row],
                       shmem[pin * rows + row - offset])
            else
                 shmem[pout * rows + row] =
                    shmem[pin * rows + row]
            end
            sync_threads()
            offset *= UInt32(2)
        end
        shmem[pin * rows + row] = shmem[pout * rows + row]
        sync_threads()

        # write back results
        data[row,col] = shmem[row]
    end

    return
end

rows = 5
cols = 4

a = rand(Int, rows, cols)

cpu_a = copy(a)
cpu_accumulate!(+, cpu_a)

gpu_a = CuArray(a)
@cuda blocks=cols threads=rows shmem=2*rows*sizeof(eltype(a)) gpu_accumulate!(+, gpu_a)

using Test

@test cpu_a ≈ Array(gpu_a)

# FURTHER IMPROVEMENTS:
# - work efficiency
# - avoid memory bank conflcits
# - large array support


text:
module InMemoryDatasets

using TableTraits,IteratorInterfaceExtensions
using Reexport
using Compat
using Printf
using PrettyTables, REPL
using Markdown
using PooledArrays
@reexport using Missings, InvertedIndices
@reexport using Statistics
@reexport using Dates
import DataAPI,
       DataAPI.All,
       DataAPI.Between,
       DataAPI.Cols,
       DataAPI.describe,
       DataAPI.innerjoin,
       DataAPI.outerjoin,
       DataAPI.rightjoin,
       DataAPI.leftjoin,
       DataAPI.semijoin,
       DataAPI.antijoin,
       DataAPI.nrow,
       DataAPI.ncol,
       # DataAPI.crossjoin,
       Tables,
       Tables.columnindex

const IMD = InMemoryDatasets
export
      # types
      IMD,
      @c_str,
      Characters,
      AbstractDataset,
      DatasetColumns,
      DatasetColumn,
      SubDataset,
      SubDatasetColumn,
      Dataset,
      GatherBy,
      GroupBy,
      Between,
      splitter,
      # functions
      nrow,
      ncol,
      rename!,
      rename,
      duplicates,
      getformat,
      setformat!,
      removeformat!,
      content,
      completecases,
      dropmissing,
      dropmissing!,
      flatten,
      flatten!,
      repeat!,
      select,
      select!,
      mapcols,
      insertcols!,
      mask,
      compare,
      groupby!,
      groupby,
      gatherby,
      describe,
      issorted!,
      unsort!,
      ungroup!,
      modify,
      modify!,
      combine,
      setinfo!,
      getinfo,
      # allowmissing!,
      # from byrow operations
      byrow,
      nunique,
      # from stat
      lag,
      lag!,
      lead,
      lead!,
      stdze,
      stdze!,
      rescale,
      topk,
      cummax,
      cummax!,
      cummin,
      cummin!,
      ffill!,
      ffill,
      bfill!,
      bfill,
      # from join
      innerjoin,
      outerjoin,
      leftjoin,
      leftjoin!,
      # rightjoin,
      antijoin,
      semijoin,
      antijoin!,
      semijoin!,
      closejoin,
      closejoin!,
      update,
      update!

include("other/index.jl")
include("characters/characters.jl")
include("other/utils.jl")
include("stat/non_hp_stat.jl")
include("stat/hp_stat.jl")
include("stat/stat.jl")
include("abstractdataset/abstractdataset.jl")
include("abstractdataset/dscol.jl")
# create dataset
include("dataset/constructor.jl")
# get elements
include("dataset/getindex.jl")
# set elements
include("dataset/setindex.jl")
# delete and append observations
include("dataset/del_and_append.jl")
# concatenate
include("dataset/cat.jl")

# byrow operations
include("byrow/util.jl")
include("byrow/row_functions.jl")
include("byrow/hp_row_functions.jl")
include("byrow/byrow.jl")
include("byrow/doc.jl")

# other functions
include("dataset/other.jl")
include("subdataset/subdataset.jl")
include("datasetrow/datasetrow.jl")
include("other/broadcasting.jl")

# modifying dataset
include("dataset/modify.jl")
include("dataset/combine.jl")
include("abstractdataset/selection.jl")
# sorting
include("sort/util.jl")
include("sort/qsort.jl")
include("sort/int.jl")
include("sort/int_any.jl")
include("sort/sortperm.jl")
include("sort/sort.jl")
include("sort/gatherby.jl")
include("sort/groupby.jl")

# transpose
include("dataset/transpose.jl")

# joins
include("join/join.jl")
include("join/join_dict.jl")
include("join/closejoin.jl")
include("join/update.jl")
include("join/main.jl")

include("abstractdataset/iteration.jl")
include("abstractdataset/prettytables.jl")
include("abstractdataset/show.jl")
include("datasetrow/show.jl")

include("abstractdataset/io.jl")

include("other/tables.jl")

# taking care of missings in other packages
include("missings/missings.jl")
# ds stat
include("stat/ds_stat.jl")
# precompile
# include("precompile/precompile.jl")
include("precompile/warmup.jl")
include("precompile/create_sysimage.jl")
# _precompile_()

function __init__()
   if Threads.nthreads() == 1
         if get(ENV, "IMD_WARN_THREADS", "1") == "1"
               @warn "Julia started with single thread, to enable multithreaded functionalities in InMemoryDatasets.jl start Julia with multiple threads."
         end
   end
end
end


text:
using StaticTrafficAssignment
using DataFrames

nnodes, ftnode, linkdf, trips, geometry, bestsolution = readtntpdata("examples/data/SiouxFalls/")

nzones = size(trips)[1]
zonedf = DataFrame(:id => 1:nzones, :issource => ones(Bool, nzones), :issink => ones(Bool, nzones), :thruallowed => (1:nzones) .>= ftnode)

network = SimpleNetwork(nnodes, linkdf, zonedf; idkey=:id, issrckey=:issource, issnkkey=:issink, thrukey=:thruallowed,  upnkey=:init_node, dwnkey=:term_node)

uecostfn = CostFunctionUE(network, BPR)
socostfn = CostFunctionSO(network, BPR)

@time aonsoln = allornothing(network, trips, uecostfn);

@time uesoln_msa = msa(network, trips, uecostfn);
@time sosoln_msa = msa(network, trips, socostfn);

@time uesoln_fw = frankwolfe(network, trips, uecostfn);
@time sosoln_fw = frankwolfe(network, trips, socostfn);

@time uesoln_cfw = conjugatefrankwolfe(network, trips, uecostfn);
@time sosoln_cfw = conjugatefrankwolfe(network, trips, socostfn);

@time uesoln_B = algorithmB(network, trips, uecostfn; λ=0.1);
@time sosoln_B = algorithmB(network, trips, socostfn; λ=0.05);


text:
# "smart" summary and REPL printing

function Base.summary(sys::System)
    return "System (base power $(get_base_power(sys)))"
end

function Base.show(io::IO, sys::System)
    println(io, "base_power=$(get_base_power(sys))")
    show(io, sys.data)
end

function Base.show(io::IO, ::MIME"text/plain", sys::System)
    println(io, "System")
    println(io, "======")
    println(io, "System Units Base: $(get_units_base(sys))")
    println(io, "Base Power: $(get_base_power(sys))")
    println(io, "Base Frequency: $(get_frequency(sys))\n")
    show(io, MIME"text/plain"(), sys.data)
end

function Base.show(io::IO, ::MIME"text/html", sys::System)
    println(io, "
System
")
    println(io, "
Base Power: $(get_base_power(sys))
")
    show(io, MIME"text/html"(), sys.data)
end

function Base.summary(tech::DeviceParameter)
    return "$(typeof(tech))"
end

function Base.show(io::IO, ::MIME"text/plain", data::PowerSystemTableData)
    println(io, "$(typeof(data)):")
    println(io, "  directory:  $(data.directory)")
    if !isnothing(data.timeseries_metadata_file)
        println(io, "  timeseries_metadata_file:  $(data.timeseries_metadata_file)")
    end
    println(io, "  base_power:  $(data.base_power)")
    for (field, df) in data.category_to_df
        print(io, "  $field:  ")
        println(io, "$(summary(df))")
    end
end