text:
module DynamicMovementPrimitives
using LinearAlgebra, RecipesBase
import OrdinaryDiffEq

export DMP, DMPopts, centraldiff,fit, solve, force, acceleration, solve_canonical, kernel_vector, plotdmp, plotdmpphase, euler

export DMP2dofopts, DMP2dof

"""
`DMPopts(Nbasis,αx,αz) = DMPopts(Nbasis, αx, αz, βz = αz/4)`\n
Holds parameters for fitting a DMP
# Fields
`Nbasis,αx,αz,βz,sched_sig,fitmethod`\n
`sched_sig` can be chosen as `:canonical` (default), `:time` or `position`\n
`fitmethod` can be chosen as `:lwr` or `:leastsquares`

See example file or the paper by Ijspeert et al. 2013
"""
struct DMPopts
    Nbasis::Int
    αx::Float64
    αz::Float64
    βz::Float64
    sched_sig::Symbol
    fitmethod::Symbol
end

DMPopts(Nbasis,αx,αz,sched_sig::Symbol=:canonical,fitmethod::Symbol = :lwr) = DMPopts(Nbasis,αx,αz,αz/4,sched_sig,fitmethod)

abstract type AbstractDMP end

"""
The result of fitting a DMP
#Fields
`opts,g,y,ẏ,ÿ,w,τ,c,σ2`\n
See example file or the paper by Ijspeert et al. 2013
"""
mutable struct DMP <: AbstractDMP
    opts::DMPopts
    g::Vector{Float64}
    y::Matrix{Float64}
    ẏ::Matrix{Float64}
    ÿ::Matrix{Float64}
    t::AbstractVector
    w::Matrix
    τ::Float64
    c::VecOrMat{Float64}
    σ2::VecOrMat{Float64}
end

include("utilities.jl")

"""
`solve_canonical(αx,τ,T::AbstractVector)`\n
`solve_canonical(αx,τ,T::Real)`\n
`solve_canonical(dmp::DMP [,t])`
"""
solve_canonical(αx,τ,T::AbstractVector) = exp.(-αx/τ.*T)
solve_canonical(αx,τ,T::Real)           = solve_canonical(αx,τ,(0:T-1))
solve_canonical(dmp::DMP,t)             = solve_canonical(dmp.opts.αx, dmp.τ, t)
solve_canonical(dmp::DMP)               = solve_canonical(dmp.opts.αx, dmp.τ, dmp.t)

function get_sched_sig(s,αx,τ,t,y,g)
    if s == :canonical
        return solve_canonical(αx,τ,t)
    elseif s == :position
        return g'.-y
    elseif s == :time
        return t
    end
    warn("Scheduling signal $s unknown")
    return solve_canonical(αx,τ,t)
end
get_sched_sig(dmp::DMP) = get_sched_sig(dmp.opts.sched_sig,dmp.opts.αx,dmp.τ,dmp.t,dmp.g)

"""
`fit(y, ẏ, ÿ, opts, g=y[end])`

Fits a DMP to data\n
`y, ẏ, ÿ` are position, velocity and acceleration respectively, `T×n` matrices where `T` is the number of time steps and `n` is the number of degrees of freedom.

see also `solve`, `plotdmp`
"""
function fit(y,ẏ,ÿ,t,opts,g=y[end,:][:])
    if opts.sched_sig != :canonical
        error("Scheduling signal $(opts.sched_sig) currently not supported")
    end
    T       = t[end]
    N       = size(y,1)
    n       = isa(y,Matrix) ? size(y,2) : 1
    τ       = T/3 # After three time constants we have reached 1-exp(-3) ≈ 0.95
    Nbasis  = opts.Nbasis
    βz      = opts.βz
    αz      = opts.αz
    αx      = opts.αx
    y0      = _1(y)
    x       = get_sched_sig(opts.sched_sig,αx,τ,t,y,g)
    ft      = τ^2*ÿ .- αz*(βz*(g'.-y).-τ*ẏ)
    ξ       = x.*(g-y0)'
    c,σ2    = get_centers_linear(Nbasis,x)
    if opts.sched_sig != :position
        Ψ = kernel_matrix(x,c,σ2)
    end
    w = zeros(Nbasis,n)
    for i = 1:n #joints
        if opts.sched_sig == :position
            Ψ = kernel_matrix(x[:,i],c[:,i],σ2[:,i])
        end
        if opts.fitmethod == :leastsquares
            w = Ψ\(ft./ξ)
        else # LWR
            for j = 1:Nbasis
                sTΓ = ξ[:,i].*Ψ[:,j]
                w[j,i] = dot(sTΓ,ft[:,i])/dot(sTΓ,ξ[:,i])
            end
        end
    end
    return DMP(opts, g, y, ẏ, ÿ,t, w, τ,c,σ2)
end

function force_single(d::AbstractDMP,x::Number, i)
    # ODE Point case
    y0 = _1(d)
    Ψ  = kernel_vector(x, d.c, d.σ2)
    f = dot(Ψ,d.w[:,i]) * x*(d.g[i]-y0[i])
end

function force_single(d::AbstractDMP,x::Number)
    # Point case
    y0 = _1(d)
    Ψ  = kernel_matrix([x], d.c, d.σ2)
    f  = (Ψ*d.w)[:] .* (x*(d.g-y0))
end

function force_single(d::AbstractDMP,x::AbstractVector)
    # Trajectory case
    y0 = _1(d)
    Ψ    = kernel_matrix(x, d.c, d.σ2)
    f = Ψ*d.w .* (x.*(d.g-y0)')
end

function force_multiple(d::AbstractDMP,x::Number, i)
    # ODE Point case
    y0 = _1(d)
    Ψ  = kernel_vector(x, d.c[:,i], d.σ2[:,i])
    f  = dot(Ψ,d.w[:,i]) * x*(d.g[i]-y0[i])
end

function force_multiple(d::AbstractDMP,x::AbstractVector)
    # Point case
    n = size(d.c,2) # Number of DOF
    y0 = _1(d)
    f = zeros(n)
    for i = 1:n
        Ψ    = kernel_matrix([x[i]], d.c[:,i], d.σ2[:,i])
        f[i] = dot(Ψ,d.w[:,i]) * x[i].*(d.g[i]-y0[i])
    end
    f
end

function force_multiple(d::AbstractDMP,x::AbstractMatrix)
    # Trajectory case
    T,n = size(x)
    y0 = _1(d)
    f = zeros(T,n)
    for i = 1:n
        Ψ    = kernel_matrix(x[:,i], d.c[:,i], d.σ2[:,i])
        f[:,i] = Ψ*d.w[:,i] .* (x[:,i].*(d.g[i]-y0[i])')
    end
    f
end

"""
`force(dmp,x)`

Calculate the forcing term for `dmp` when the phase variable is `x`\n
The return value will be an `n` Vector or `T×n` Matrix depending on `typeof(x)`
"""
function force(d::AbstractDMP, x)
    if d.opts.sched_sig == :position
        force_multiple(d,x)
    else
        force_single(d,x)
    end
end

function force(d::AbstractDMP, x,i)
    if d.opts.sched_sig == :position
        force_multiple(d,x,i)
    else
        force_single(d,x,i)
    end
end

function acceleration(d::DMP, y::Number,ẏ::Number,x::Number,g::Number,i=1)
    f = force(d,x,i)
    (d.opts.αz*(d.opts.βz*(g-y)-d.τ*ẏ)+f)/d.τ^2
end

function acceleration(d::DMP, y::AbstractVector,ẏ::AbstractVector,x ,g::AbstractVector = d.g)
    f = force(d,x)
    (d.opts.αz*(d.opts.βz*(g-y)-d.τ*ẏ)+f)/d.τ^2
end

function acceleration(d::DMP, y::AbstractMatrix,ẏ::AbstractMatrix,x::AbstractVecOrMat,g::AbstractVector)
    f = force(d,x)
    (d.opts.αz*(d.opts.βz*(g'.-y)-d.τ*ẏ)+f)/d.τ^2
end

function solve_canonical(dmp::DMP, t, y0, g, solver; kwargs...)
    T,n = size(dmp.y)
    αx  = dmp.opts.αx
    τ   = dmp.τ
    z   = zeros(T,n)
    y   = zeros(T,n)
    x   = zeros(T)
    for i = 1:n
        function time_derivative(dstate, state, p, t)
            local z   = state[1]
            local y   = state[2]
            local x   = state[3]
            dstate[1] = acceleration(dmp, y, z, x,g[i],i)
            dstate[2] = z
            dstate[3] = -αx/τ * x
        end
        state0 = [dmp.ẏ[1,i]; y0[i]; 1.]
        prob = OrdinaryDiffEq.ODEProblem(time_derivative,state0,(t[1],t[end]))
        sol = OrdinaryDiffEq.solve(prob,solver;saveat=t,dt=t[2], kwargs...)
        z[:,i] = sol[1,:]
        y[:,i] = sol[2,:]
        x      = sol[3,:] # TODO: se till att denna är samma för alla DOF
    end
    t,y,z,x
end

function solve_position(dmp, t, y0, g, solver; kwargs...)
    T,n = size(dmp.y)
    αx  = dmp.opts.αx
    τ   = dmp.τ
    z   = zeros(T,n)
    y   = zeros(T,n)
    for i = 1:n
        function time_derivative(dstate, state, p, t)
            local z   = state[1]
            local y   = state[2]
            dstate[1] = acceleration(dmp, y, z, g[i]-y,g[i],i)
            dstate[2] = z
        end
        state0  = [0; y0[i]]
        prob = OrdinaryDiffEq.ODEProblem(time_derivative,state0,(t[1],t[end]))
        sol = OrdinaryDiffEq.solve(prob,solver;saveat=t,dt=t[2], kwargs...)
        z[:,i] = sol[1,:]
        y[:,i] = sol[2,:]
    end
    t,y,z,y
end

function solve_time(dmp, t, y0, g, solver; kwargs...)
    T,n     = size(dmp.y)
    αx      = dmp.opts.αx
    τ       = dmp.τ
    z       = zeros(T,n)
    y       = zeros(T,n)
    for i = 1:n
        function time_derivative(dstate, state, p, t)
            local z   = state[1]
            local y   = state[2]
            dstate[1] = acceleration(dmp, y, z, t ,g[i],i)
            dstate[2] = z
        end
        state0  = [0; y0[i]]
        prob = OrdinaryDiffEq.ODEProblem(time_derivative,state0,(t[1],t[end]))
        sol = OrdinaryDiffEq.solve(prob,solver;saveat=t,dt=t[2], kwargs...)
        z[:,i] = sol[1,:]
        y[:,i] = sol[2,:]
    end
    t,y,z,t
end

"""
`t,y,z,x = solve(dmp, t = 0:length(dmp.t)-1; y0 = _1(dmp), g = dmp.g, solver=ode45)`

`t` time vector

## Keyword arguments: \n
`y0` start position, defaults to the initial point in training data from `dmp`
`g` goal, defaults to goal from `dmp`\n
`solver` the ode solver to use, see https://github.com/JuliaLang/ODE.jl \n
The default solver is `solver=ode54`, a faster alternative is `solver=euler`

see also `plotdmp`
"""
function solve(dmp::AbstractDMP, t = dmp.t; y0 = _1(dmp), g = dmp.g, solver=OrdinaryDiffEq.Tsit5(), kwargs...)
    if dmp.opts.sched_sig == :position
        return solve_position(dmp, t, y0, g, solver; kwargs...)
    elseif dmp.opts.sched_sig == :time
        return solve_time(dmp, t, y0, g, solver; kwargs...)
    end
    solve_canonical(dmp, t, y0, g, solver; kwargs...)
end

"""
    A = linearize(dmp::AbstractDMP, y, ẏ, x, y0)

Linearize the dynamics of the DMP around the point (y,ẏ,x). The resulting linear system is on the form
`q̈ = A*[y; ẏ; x]`
"""
function linearize(dmp::AbstractDMP, y, ẏ, x, y0)
    error("g not taken care of")
    n = size(y,1)
    Ø = 0I
    ϕ = kernel_vector(dmp,x)
    ∇ϕ = - ϕ * (x-dmp.c)/dmp.σ2
    ∇ϕwx = ∇ϕ*dmp.w*x + ϕ*dmp.w
    A = -dmp.αz/dmp.τ * [Ø I Ø; dmp.βz/dmp.τ*I I Ø; Ø Ø dmp.αx*I] # Linear part
    A[n+1:2n, 2n+1:3n] = ∇ϕwx * (dmp.g-y0) / dmp.τ^2# Linearized force term
end

include("DMP2dof.jl")
end # module


text:
module MyMCMCSampler
# export simpleRandomGraph,randomGraph, singlePotential, multiPotential, makeDataDict, sampleSingleChain, sampleAllChains, ThirdOrderExpFamily

using Reexport
@reexport using Mamba
using Distributions
using DataFrames
using ForwardDiff
using Distributed
using Parameters
using Combinatorics: combinations

###################### Types #############################
struct ThirdOrderExpFamily <: ContinuousMultivariateDistribution
    edge_index::AbstractArray # an array of node index in each edge
    edge_attr::AbstractArray # edge potential
    edge_type::AbstractArray # array of edge type
    node_attr::AbstractArray # singleton potential
    nnode::Integer # number of variable node
    nfactor::Integer # number of factor node
    alpha::Float64  #  coefficient before the |x|^4_2 base measure
end

ThirdOrderExpFamily(edge_index, edge_attr, edge_type, node_attr, nnode::Integer, nfactor::Integer) = ThirdOrderExpFamily(edge_index, edge_attr, edge_type, node_attr, nnode, nfactor, 0.2)
ThirdOrderExpFamily(edge_index, edge_attr, edge_type, node_attr) = ThirdOrderExpFamily(edge_index, edge_attr, edge_type, node_attr, length(node_attr), length(edge_type))
ThirdOrderExpFamily(edge_index, edge_attr, edge_type, node_attr, alpha::Real) = ThirdOrderExpFamily(edge_index, edge_attr, edge_type, node_attr, length(node_attr), length(edge_type), alpha)

struct SingleThirdOrderExpFamily <: ContinuousMultivariateDistribution
    node_attr::AbstractArray # singleton potential
    alpha::Float64
end

SingleThirdOrderExpFamily(node_attr) = SingleThirdOrderExpFamily(node_attr, 0.2)

################### Spin Glass Distributions ######################3
struct DiscreteSpinGlass <: DiscreteMultivariateDistribution
    edge_index::AbstractArray # an array of node index in each edge
    edge_attr::AbstractArray # edge potential
    edge_type::AbstractArray # edge type
    node_attr::AbstractArray # singleton potential
    nnode::Integer # number of variable node
    nfactor::Integer # number of factor node
    beta::Real # inverse temperature
end

# DiscreteSpinGlass(edge_index,edge_attr,node_attr,nnode::Integer) = DiscreteSpinGlass(edge_index,edge_attr,node_attr,nnode,nfactor)
DiscreteSpinGlass(edge_index, edge_attr, edge_type, node_attr, beta::Real) = DiscreteSpinGlass(edge_index, edge_attr,edge_type, node_attr, length(node_attr), length(edge_attr), beta)
DiscreteSpinGlass(edge_index, edge_attr, edge_type, node_attr) = DiscreteSpinGlass(edge_index, edge_attr, edge_type, node_attr, 1.0)
# DiscreteSpinGlass(;kwargs) = DiscreteSpinGlass(;kwargs...)
# DiscreteSpinGlass(edge_index,edge_attr,node_attr)=DiscreteSpinGlass(edge_index,edge_attr,node_attr,length(node_attr),length(edge_attr))

function add_singleton_factor(d::DiscreteSpinGlass)
    @unpack edge_index, edge_attr, edge_type, node_attr, beta,nnode,nfactor = d
    nfactor += nnode
    for i in 1:nnode
        push!(edge_index, [i])
        push!(edge_type,1)
        push!(edge_attr, node_attr[i])
    end
    return DiscreteSpinGlass(edge_index, edge_attr, edge_type, node_attr, nnode, nfactor,beta)
end

function add_singleton_factor(d::ThirdOrderExpFamily)
    @unpack edge_index, edge_attr, edge_type, node_attr, alpha,nnode,nfactor = d
    nfactor += nnode
    for i in 1:nnode
        push!(edge_index, [i])
        push!(edge_type,1)
        push!(edge_attr, node_attr[i])
    end
    return ThirdOrderExpFamily(edge_index, edge_attr, edge_type, node_attr, nnode, nfactor,alpha)
end

function has_singleton_factor(d::Union{ThirdOrderExpFamily, DiscreteSpinGlass})
    return (1 in unique(d.edge_type))
end

###################### Include #############################
include("thirdexpfam.jl")
include("mcmcutils.jl")
include("random.jl")
include("bp.jl")
# include("gaussian.jl")

###################### Exports #########################
export
    ThirdOrderExpFamily,
    SingleThirdOrderExpFamily,
    DiscreteSpinGlass,
    sampleSingleChain,
    sampleMultiChain,
    sampleMultiChain2,
    potential,
    dist2dict,
    getstats,
    isvalid,
    rand_regular_graph,
    add_singleton_factor,
    random_ws_factor_graph,
    random_ws_factor_structure,
    incident_matrix,
    computemarginal,
    get_neighbours,
    bp_3spin,
    bp

end # module MyMCMCSampler

ext: jl


text:
function Div(A, B )
	C = A ./ B;
	if typeof(B) <: Array
		for i=1:length(C)
			if (A[i] == 0) & (B[i] == 0)
				C[i] = 0
			end
		end
	elseif typeof(B) <: Number
		for i=1:length(C)
			if (A[i] == 0) & (B == 0)
				C[i] = 0
			end
		end
	end
	return C
end

max_stars_repo_name: barzan/dbseer

max_issues_repo_name: barzan/dbseer

max_forks_repo_name: barzan/dbseer


text:
@inline function ExaTronKernels.dscal(n::Int,da::Float64,dx::oneDeviceArray{Float64,1},incx::Int)
    tx = get_local_id()

    # Ignore incx for now.
    if tx <= n
        @inbounds dx[tx] = da*dx[tx]
    end
    barrier()

    return
end


text:
@with_kw mutable struct MPCOptions{T}
    "Maximum horizon of the MPC"
    h::T=3.0

    "Maximum number of MPC steps"
    M::Int=300
end

mutable struct MPCStatistics{T}
    iter::Int
    t::T
    dt::Vector{T}
    traj::Vector{Algames.KnotPoint}
end

function MPCStatistics()
    iter = 0
    t = 0.0
    dt = zeros(0)
    traj = Vector{Algames.KnotPoint}()
    return MPCStatistics(iter,t,dt,traj)
end

function Algames.reset!(stats::MPCStatistics)
    stats.iter = 0
    stats.t = 0.0
    stats.dt = zeros(0)
    stats.traj = Vector{Algames.KnotPoint}()
    return nothing
end

function Algames.record!(stats::MPCStatistics, z::Algames.KnotPoint)
    stats.iter += 1
    stats.t += z.dt
    push!(stats.dt, z.dt)
    push!(stats.traj, z)
    return nothing
end

mutable struct TrajStatistics{T}
    iter::Int
    probsize::ProblemSize
    t::T
    dt::Vector{T}
    traj::Vector{Vector{Algames.Traj}}
end

function TrajStatistics(probsize)
    iter = 0
    t = 0.0
    dt = zeros(0)
    traj = Vector{Vector{Algames.Traj}}()
    return TrajStatistics(iter, probsize,t,dt,traj)
end

function Algames.reset!(stats::TrajStatistics)
    stats.iter = 0
    stats.t = 0.0
    stats.dt = zeros(0)
    stats.traj = Vector{Vector{Algames.Traj}}()
    return nothing
end

function Algames.record!(stats::TrajStatistics, dt::T,
    trajs::Vector{<:Algames.Traj}) where {T}
    @assert length(trajs) == stats.probsize.p
    stats.iter += 1
    stats.t += dt
    push!(stats.dt, dt)
    push!(stats.traj, trajs)
    return nothing
end


text:
using JuMP, GLPK

function SM(state::DynamicIndividualState, static_state::StaticIndividualState,
    depth::Int64; allow_nothing::Bool = false,
    allow_overfarming::Bool = false, sim_to_end::Bool = false)
    A, B = get_possible_decisions(state, static_state,
        allow_nothing = allow_nothing, allow_overfarming = allow_overfarming)

    (Base.ctpop_int(A) == 0x00 || Base.ctpop_int(B) == 0x00 || depth == 0) &&
        return sim_to_end ?
            NashResult(sum(battle_scores(state, static_state, 100)
                / 100) - 0.5, no_strat, no_strat) :
            NashResult(battle_score(state, static_state) - 0.5,
                no_strat, no_strat)

    payoffs = zeros(Float64, Base.ctpop_int(A), Base.ctpop_int(B))

    state_1, state_2 = state, state
    odds = 1.0
    chance = get_chance(state)
    if chance == UInt32(5)
        state_1 = DynamicIndividualState(state[0x01], state[0x02], state.data - UInt32(9702))
        state_2 = DynamicIndividualState(state[0x01], state[0x02], state.data - UInt32(9261))
        odds = 0.5
    else
        agent = chance < UInt32(3) ? 0x01 : 0x02
        move = isodd(chance) ? static_state[agent].charged_move_1 :
            static_state[agent].charged_move_2
        data = apply_buff(state.data, move, agent)
        state_1 = DynamicIndividualState(state[0x01], state[0x02],
            data - chance * UInt32(2205))
        state_2 = DynamicIndividualState(state[0x01], state[0x02],
            state.data - chance * UInt32(2205))
        odds = move.buffChance / 100
    end
    for i = 0x01:Base.ctpop_int(A), j = 0x01:Base.ctpop_int(B)
        if odds > 0.99
            @inbounds payoffs[i, j] = SM(play_turn(state_1, static_state,
                get_decisions(A, B, i, j)), static_state, depth - 1,
                allow_nothing = allow_nothing,
                allow_overfarming = allow_overfarming,
                sim_to_end = sim_to_end).payoff
        else
            @inbounds payoffs[i, j] = odds * SM(play_turn(state_1, static_state,
                get_decisions(A, B, i, j)), static_state, depth - 1,
                allow_nothing = allow_nothing,
                allow_overfarming = allow_overfarming,
                sim_to_end = sim_to_end).payoff +
                (1 - odds) * SM(play_turn(state_2, static_state,
                get_decisions(A, B, i, j)), static_state, depth - 1,
                allow_nothing = allow_nothing,
                allow_overfarming = allow_overfarming,
                sim_to_end = sim_to_end).payoff
        end
    end
    return nash(payoffs)
end

function solve_battle(s::DynamicIndividualState, static_s::StaticIndividualState, depth::Int64;
    allow_nothing::Bool = false, allow_overfarming::Bool = false, sim_to_end = false)
    value = 0.0
    decision = 0, 0
    strat = IndividualStrategy([], [], [])
    while true
        s = resolve_chance(s, static_s)
        A, B = get_possible_decisions(s, static_s,
            allow_nothing = allow_nothing, allow_overfarming = allow_overfarming)

        (Base.ctpop_int(A) == 0x00 || Base.ctpop_int(B) == 0x00) &&
            return value, strat
        if Base.ctpop_int(A) == 0x01 && Base.ctpop_int(B) == 0x01
            decision = get_decisions(A, B, 0x01, 0x01)
        else
            nash_result = SM(s, static_s, depth, sim_to_end = sim_to_end)
            value = nash_result.payoff
            d1, d2 = rand(rb_rng), rand(rb_rng)
            decision1, decision2 = UInt8(length(nash_result.row_strategy)),
                UInt8(length(nash_result.column_strategy))
            j = 0.0
            for i = 0x01:(UInt8(length(nash_result.row_strategy)) - 0x01)
                @inbounds j += nash_result.row_strategy[i]
                if d1 < j
                    decision1 = i
                    break
                end
            end
            j = 0.0
            for i = 0x01:(UInt8(length(nash_result.column_strategy)) - 0x01)
                @inbounds j += nash_result.column_strategy[i]
                if d2 < j
                    decision2 = i
                    break
                end
            end
            decision = get_decisions(A, B, decision1, decision2)
        end
        s = play_turn(s, static_s, decision)
        push!(strat.decisions, decision)
        push!(strat.scores, value + 0.5)
        push!(strat.hps, (get_hp(s[0x01]), get_hp(s[0x02])))
    end
end


text:
using FoldsCUDABenchmarks
SUITE = FoldsCUDABenchmarks.setup()


text:
using StructuralCausalModels, Test

ProjDir = @__DIR__

d_str = "dag{k1 -> k2;k2 -> y;v -> x2;w -> k1;x1 -> v;x1 -> w;x2 -> y;x3 -> w;x3 -> y}"

dag = DAG("test_desc", d_str)

dagitty_str = to_dagitty(dag)

fname = joinpath(ProjDir, "test_descendents_03.dot")
to_graphviz(dag, fname)
#Sys.isapple() && run(`open -a GraphViz.app $(fname)`)

bs = basis_set(dag)
println()

bs[6] |> display
@show d_separation(dag, :v, :x3; c=[:x1])
@show d_separation(dag, :v, :x3)
@show d_separation(dag, :x1, :x3; c=[:x1]) # <====
println()

bs[7] |> display
@show d_separation(dag, :v, :w; c=[:x1, :x3])
@show d_separation(dag, :v, :w; c=[:x1])
@show d_separation(dag, :v, :w; c=[:x3])
println()

bs[8] |> display
@show d_separation(dag, :v, :k1; c=[:x1, :w])
@show d_separation(dag, :v, :k1; c=[:x1])
@show d_separation(dag, :v, :k1; c=[:w])
println()

bs[9] |> display
@show d_separation(dag, :v, :y; c=[:x1, :x2, :x3, :k2])
@show d_separation(dag, :v, :y; c=[:k2, :x2, :x3])
@show d_separation(dag, :v, :y; c=[:k1, :x2, :x3])
@show d_separation(dag, :v, :y; c=[:w, :x2, :x3])
@show d_separation(dag, :v, :y; c=[:x1, :x2])
println()

bs[10] |> display
@show d_separation(dag, :x2, :x3; c=[:v])
@show d_separation(dag, :x2, :x3)
println()


text:
# # Langmuir turbulence example
#
# This example implements a Langmuir turbulence simulation reported in section
# 4 of
#
# > [Wagner et al., "Near-inertial waves and turbulence driven by the growth of swell", Journal of Physical Oceanography (2021)](https://journals.ametsoc.org/view/journals/phoc/51/5/JPO-D-20-0178.1.xml)
#
# This example demonstrates
#
#   * How to run large eddy simulations with surface wave effects via the
#     Craik-Leibovich approximation.
#
#   * How to specify time- and horizontally-averaged output.

# ## Install dependencies
#
# First let's make sure we have all required packages installed.

# ```julia
# using Pkg
# pkg"add Oceananigans, Plots"
# ```

using Oceananigans
using Oceananigans.Units: minute, minutes, hours

# ## Model set-up
#
# To build the model, we specify the grid, Stokes drift, boundary conditions, and
# Coriolis parameter.
#
# ### Domain and numerical grid specification
#
# We use a modest resolution and the same total extent as Wagner et al. 2021,

grid = RectilinearGrid(size=(32, 32, 32), extent=(128, 128, 64))

# ### The Stokes Drift profile
#
# The surface wave Stokes drift profile prescribed in Wagner et al. 2021,
# corresponds to a 'monochromatic' (that is, single-frequency) wave field.
#
# A monochromatic wave field is characterized by its wavelength and amplitude
# (half the distance from wave crest to wave trough), which determine the wave
# frequency and the vertical scale of the Stokes drift profile.

using Oceananigans.BuoyancyModels: g_Earth

 amplitude = 0.8 # m
wavelength = 60 # m
wavenumber = 2π / wavelength # m⁻¹
 frequency = sqrt(g_Earth * wavenumber) # s⁻¹

## The vertical scale over which the Stokes drift of a monochromatic surface wave
## decays away from the surface is `1/2wavenumber`, or
const vertical_scale = wavelength / 4π

## Stokes drift velocity at the surface
const Uˢ = amplitude^2 * wavenumber * frequency # m s⁻¹

# The `const` declarations ensure that Stokes drift functions compile on the GPU.
# To run this example on the GPU, write `architecture = GPU()` in the constructor
# for `NonhydrostaticModel` below.
#
# The Stokes drift profile is

uˢ(z) = Uˢ * exp(z / vertical_scale)

# which we'll need for the initial condition.
#
# !!! info "The Craik-Leibovich equations in Oceananigans"
#     Oceananigans implements the Craik-Leibovich approximation for surface wave effects
#     using the _Lagrangian-mean_ velocity field as its prognostic momentum variable.
#     In other words, `model.velocities.u` is the Lagrangian-mean ``x``-velocity beneath surface
#     waves. This differs from models that use the _Eulerian-mean_ velocity field
#     as a prognostic variable, but has the advantage that ``u`` accounts for the total advection
#     of tracers and momentum, and that ``u = v = w = 0`` is a steady solution even when Coriolis
#     forces are present. See the
#     [physics documentation](https://clima.github.io/OceananigansDocumentation/stable/physics/surface_gravity_waves/)
#     for more information.
#
# The vertical derivative of the Stokes drift is

∂z_uˢ(z, t) = 1 / vertical_scale * Uˢ * exp(z / vertical_scale)

# Finally, we note that the time-derivative of the Stokes drift must be provided
# if the Stokes drift and surface wave field undergoes _forced_ changes in time.
# In this example, the Stokes drift is constant
# and thus the time-derivative of the Stokes drift is 0.

# ### Boundary conditions
#
# At the surface at ``z=0``, Wagner et al. 2021 impose

Qᵘ = -3.72e-5 # m² s⁻², surface kinematic momentum flux

u_boundary_conditions = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵘ))

# Wagner et al. 2021 impose a linear buoyancy gradient `N²` at the bottom
# along with a weak, destabilizing flux of buoyancy at the surface to faciliate
# spin-up from rest.

Qᵇ = 2.307e-9 # m³ s⁻², surface buoyancy flux
N² = 1.936e-5 # s⁻², initial and bottom buoyancy gradient

b_boundary_conditions = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵇ),
                                                bottom = GradientBoundaryCondition(N²))

# !!! info "The flux convention in Oceananigans"
#     Note that Oceananigans uses "positive upward" conventions for all fluxes. In consequence,
#     a negative flux at the surface drives positive velocities, and a positive flux of
#     buoyancy drives cooling.

# ### Coriolis parameter
#
# Wagner et al. (2021) use

coriolis = FPlane(f=1e-4) # s⁻¹

# which is typical for mid-latitudes on Earth.

# ## Model instantiation
#
# We are ready to build the model. We use a fifth-order Weighted Essentially
# Non-Oscillatory (WENO) advection scheme and the `AnisotropicMinimumDissipation`
# model for large eddy simulation. Because our Stokes drift does not vary in ``x, y``,
# we use `UniformStokesDrift`, which expects Stokes drift functions of ``z, t`` only.

model = NonhydrostaticModel(; grid,
                            advection = WENO5(),
                            timestepper = :RungeKutta3,
                            tracers = :b,
                            buoyancy = BuoyancyTracer(),
                            coriolis = coriolis,
                            closure = AnisotropicMinimumDissipation(),
                            stokes_drift = UniformStokesDrift(∂z_uˢ=∂z_uˢ),
                            boundary_conditions = (u=u_boundary_conditions, b=b_boundary_conditions))

# ## Initial conditions
#
# We make use of random noise concentrated in the upper 4 meters
# for buoyancy and velocity initial conditions,

Ξ(z) = randn() * exp(z / 4)
nothing # hide

# Our initial condition for buoyancy consists of a surface mixed layer 33 m deep,
# a deep linear stratification, plus noise,

initial_mixed_layer_depth = 33 # m
stratification(z) = z < - initial_mixed_layer_depth ? N² * z : N² * (-initial_mixed_layer_depth)

bᵢ(x, y, z) = stratification(z) + 1e-1 * Ξ(z) * N² * model.grid.Lz

# The simulation we reproduce from Wagner et al. (2021) is zero Lagrangian-mean velocity.
# This initial condition is consistent with a wavy, quiescent ocean suddenly impacted
# by winds. To this quiescent state we add noise scaled by the friction velocity to ``u`` and ``w``.

u★ = sqrt(abs(Qᵘ))
uᵢ(x, y, z) = u★ * 1e-1 * Ξ(z)
wᵢ(x, y, z) = u★ * 1e-1 * Ξ(z)

set!(model, u=uᵢ, w=wᵢ, b=bᵢ)

# ## Setting up the simulation

simulation = Simulation(model, Δt=45.0, stop_time=4hours)

# We use the `TimeStepWizard` for adaptive time-stepping
# with a Courant-Freidrichs-Lewy (CFL) number of 1.0,

wizard = TimeStepWizard(cfl=1.0, max_change=1.1, max_Δt=1minute)

simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(10))

# ### Nice progress messaging
#
# We define a function that prints a helpful message with
# maximum absolute value of ``u, v, w`` and the current wall clock time.

using Printf

function print_progress(simulation)
    u, v, w = simulation.model.velocities

    ## Print a progress message
    msg = @sprintf("i: %04d, t: %s, Δt: %s, umax = (%.1e, %.1e, %.1e) ms⁻¹, wall time: %s\n",
                   iteration(simulation),
                   prettytime(time(simulation)),
                   prettytime(simulation.Δt),
                   maximum(abs, u), maximum(abs, v), maximum(abs, w),
                   prettytime(simulation.run_wall_time))

    @info msg

    return nothing
end

simulation.callbacks[:progress] = Callback(print_progress, IterationInterval(10))

# ## Output
#
# ### A field writer
#
# We set up an output writer for the simulation that saves all velocity fields,
# tracer fields, and the subgrid turbulent diffusivity.

output_interval = 5minutes

fields_to_output = merge(model.velocities, model.tracers, (; νₑ=model.diffusivity_fields.νₑ))

simulation.output_writers[:fields] =
    JLD2OutputWriter(model, fields_to_output,
                     schedule = TimeInterval(output_interval),
                     filename = "langmuir_turbulence_fields.jld2",
                     overwrite_existing = true)

# ### An "averages" writer
#
# We also set up output of time- and horizontally-averaged velocity field and
# momentum fluxes,

u, v, w = model.velocities
b = model.tracers.b

 U = Average(u, dims=(1, 2))
 V = Average(v, dims=(1, 2))
 B = Average(b, dims=(1, 2))
wu = Average(w * u, dims=(1, 2))
wv = Average(w * v, dims=(1, 2))

simulation.output_writers[:averages] =
    JLD2OutputWriter(model, (; U, V, B, wu, wv),
                     schedule = AveragedTimeInterval(output_interval, window=2minutes),
                     filename = "langmuir_turbulence_averages.jld2",
                     overwrite_existing = true)

# ## Running the simulation
#
# This part is easy,

run!(simulation)

# # Making a neat movie
#
# We look at the results by loading data from file with FieldTimeSeries,
# and plotting vertical slices of ``u`` and ``w``, and a horizontal
# slice of ``w`` to look for Langmuir cells.

using Plots

time_series = (;
     w = FieldTimeSeries("langmuir_turbulence_fields.jld2", "w"),
     u = FieldTimeSeries("langmuir_turbulence_fields.jld2", "u"),
     B = FieldTimeSeries("langmuir_turbulence_averages.jld2", "B"),
     U = FieldTimeSeries("langmuir_turbulence_averages.jld2", "U"),
     V = FieldTimeSeries("langmuir_turbulence_averages.jld2", "V"),
    wu = FieldTimeSeries("langmuir_turbulence_averages.jld2", "wu"),
    wv = FieldTimeSeries("langmuir_turbulence_averages.jld2", "wv"))

times = time_series.w.times
xw, yw, zw = nodes(time_series.w)
xu, yu, zu = nodes(time_series.u)
nothing # hide

# This utility is handy for calculating nice contour intervals:

function nice_divergent_levels(c, clim; nlevels=20)
    levels = range(-clim, stop=clim, length=nlevels)
    cmax = maximum(abs, c)
    clim < cmax && (levels = vcat([-cmax], levels, [cmax]))
    return (-clim, clim), levels
end
nothing # hide

# Finally, we're ready to animate.

@info "Making an animation from the saved data..."

anim = @animate for (i, t) in enumerate(times)

    @info "Drawing frame $i of $(length(times))..."

    ## Get fields at ith save point
     wᵢ = time_series.w[i]
     uᵢ = time_series.u[i]
     Bᵢ = time_series.B[i][1, 1, :]
     Uᵢ = time_series.U[i][1, 1, :]
     Vᵢ = time_series.V[i][1, 1, :]
    wuᵢ = time_series.wu[i][1, 1, :]
    wvᵢ = time_series.wv[i][1, 1, :]

    ## Extract slices
    k = searchsortedfirst(grid.zᵃᵃᶠ[:], -8)
    wxy = wᵢ[:, :, k]
    wxz = wᵢ[:, 1, :]
    uxz = uᵢ[:, 1, :]

    wlims, wlevels = nice_divergent_levels(interior(w), 0.03)
    ulims, ulevels = nice_divergent_levels(interior(w), 0.05)

    B_plot = plot(Bᵢ, zu,
                  label = nothing,
                  legend = :bottom,
                  xlabel = "Buoyancy (m s⁻²)",
                  ylabel = "z (m)")

    U_plot = plot([Uᵢ Vᵢ], zu,
                  label = ["\$ \\bar u \$" "\$ \\bar v \$"],
                  legend = :bottom,
                  xlabel = "Velocities (m s⁻¹)",
                  ylabel = "z (m)")

    wu_label = "\$ \\overline{wu} \$"
    wv_label = "\$ \\overline{wv} \$"

    fluxes_plot = plot([wuᵢ, wvᵢ], zw,
                       label = [wu_label wv_label],
                       legend = :bottom,
                       xlabel = "Momentum fluxes (m² s⁻²)",
                       ylabel = "z (m)")

    wxy_plot = contourf(xw, yw, wxy';
                        color = :balance,
                        linewidth = 0,
                        aspectratio = :equal,
                        clims = wlims,
                        levels = wlevels,
                        xlims = (0, grid.Lx),
                        ylims = (0, grid.Ly),
                        xlabel = "x (m)",
                        ylabel = "y (m)")

    wxz_plot = contourf(xw, zw, wxz';
                              color = :balance,
                          linewidth = 0,
                        aspectratio = :equal,
                              clims = wlims,
                             levels = wlevels,
                              xlims = (0, grid.Lx),
                              ylims = (-grid.Lz, 0),
                             xlabel = "x (m)",
                             ylabel = "z (m)")

    uxz_plot = contourf(xu, zu, uxz';
                              color = :balance,
                          linewidth = 0,
                        aspectratio = :equal,
                              clims = ulims,
                             levels = ulevels,
                              xlims = (0, grid.Lx),
                              ylims = (-grid.Lz, 0),
                             xlabel = "x (m)",
                             ylabel = "z (m)")

    wxy_title = @sprintf("w(x, y, t) (m s⁻¹) at z=-8 m and t = %s ", prettytime(t))
    wxz_title = @sprintf("w(x, z, t) (m s⁻¹) at y=0 m and t = %s", prettytime(t))
    uxz_title = @sprintf("u(x, z, t) (m s⁻¹) at y=0 m and t = %s", prettytime(t))

    plot(wxy_plot, B_plot, wxz_plot, U_plot, uxz_plot, fluxes_plot,
         layout = Plots.grid(3, 2, widths=(0.7, 0.3)), size = (900.5, 1000.5),
         title = [wxy_title "" wxz_title "" uxz_title ""])
end

mp4(anim, "langmuir_turbulence.mp4", fps = 8) # hide


text:
## Compute negative vertex of a polynomial where,
## a vertex of a polynomial is a term corresponding to a vertex of its Newton polytope,'
## and the sign of a vertex (a term) is the sign of its coefficient.
##
## Here, Polymake computes the correspondence between vertices and exponent_vectors.
## see NegativeVertex.old.jl for an alternative approach, where is julia and C++ libflint who computes/gives
## such a correspondence.

const ver_sufix = "_vertices.cvs"
const map_sufix = "_ver_pt_map.cvs"

"""
$(TYPEDEF)

A polynomial `p`, its Newton polytope `Newtonpolyt` (as Polymake big object, with the vertices computed) and the map vertices -> points `vertex_point_map` (vector storing the position of vertex in polynomial `p`).

Methods to create a `PolyPolyt`:
```julia
PolyPolyt(p::MPolyElem) # Requires computing the vertices of the Newton polytope.
PolyPolyt(p::MPolyElem, vertices, ver_pt_map) # Loads vertices and ver-to-pt map.
PolyPolyt(p::MPolyElem, name) # Loads them from files generated by save_polypolyt(name::String, pp)
```
"""
struct PolyPolyt{T}
    p::MPolyElem{T}
    Newtonpolyt::Polymake.BigObject
    vertex_point_map::Vector{Int}
end

"""
$(SIGNATURES)

Returns the Newton polytope of a `PolyPolyt`.
"""
Newtonpolytope(pp::PolyPolyt) = pp.Newtonpolyt

"""
$(SIGNATURES)

Returns the ver-to-pt map of a `PolyPolyt`.
"""
vertex_point_map(pp::PolyPolyt) = pp.vertex_point_map

"""
$(SIGNATURES)

Given a vertex `i` of a `PolyPolyt`, returns the position of `i` in ver-to-pt map.
That is, `vertex_point_map(pp)[vertex_index(pp, i)] == i`.
Or more interestingly,
```julia
verticesof(Newtonpolytope(pp))[:, vertex_index(pp, i)] == collect(exponent_vectors(pp.p))[:,i]
```
"""
vertex_index(pp::PolyPolyt, i::Int) = searchsortedfirst(vertex_point_map(pp), i - 1)

function Base.convert(::Type{PolyPolyt}, p)
    pt_config = Polymake.polytope.PointConfiguration(POINTS=polymake_homog(permutedims(exponents_matrix(p))))
    @computing ver_pt_map = Vector(pt_config.VERTEX_POINT_MAP) "vertices of Newton polytope"
    return PolyPolyt(p,
                     Polymake.polytope.Polytope(VERTICES=pt_config.CONVEX_HULL.VERTICES),
                     ver_pt_map)
end

PolyPolyt(p::MPolyElem) = convert(PolyPolyt, p)

PolyPolyt(p::MPolyElem, vertices, ver_pt_map) =
    PolyPolyt(p,
              Polymake.polytope.Polytope(VERTICES=vertices),
              vec(ver_pt_map))

PolyPolyt(p::MPolyElem, name::AbstractString) =
    PolyPolyt(p,
              readdlm(name * ver_sufix, Int),
              readdlm(name * map_sufix, Int))

# TODO save and load the polynomial as well
# PolyPolyt(p::AbstractString, vertices::AbstractString, ver_pt_map::AbstractString) where T =
#     PolyPolyt(...,
#               Polymake.polytope.Polytope(VERTICES=DF.readdlm(vertices, Int)),
#               readdlm(ver_pt_map, Int))

"""
$(SIGNATURES)

Save vertices of `pp::PolyPolyt` in file `name * ver_sufix` and ver-to-pt map in `name * map_sufix`
"""
function save_polypolyt(name::AbstractString, pp::PolyPolyt)
    open(name * ver_sufix, "w") do io
        writedlm(io, Newtonpolytope(pp).VERTICES)
    end
    open(name * map_sufix, "w") do io
        writedlm(io, vertex_point_map(pp))
    end
end

negvertices(pp::PolyPolyt, isneg_fun=_isnegative) =
    filter(i -> isneg_fun(coeff(pp.p, i)), vertex_point_map(pp) .+ 1)
negvertices(p, name::AbstractString, isneg_fun=_isnegative) = negvertices(PolyPolyt(p, name), isneg_fun=isneg_fun)
negvertices(p, isneg_fun=_isnegative) = negvertices(PolyPolyt(p), isneg_fun=isneg_fun)

# posvertices(pp::PolyPolyt, ispos_fun=_ispositive) =
#     filter(i -> ispos_fun(coeff(pp.p, i)), vertex_point_map(pp) .+ 1)
# posvertices(p, name::AbstractString, isneg_fun=_isnegative) = posvertices(PolyPolyt(p, name),isneg_fun=isneg_fun)
# posvertices(p, isneg_fun=_isnegative) = posvertices(PolyPolyt(p),isneg_fun=isneg_fun)


text:
import MyterialColors: orange_light, teal, purple_light, pink

"""
Definition of several type of columns for progress bars.
Used in progress.jl.
"""

# ---------------------------------------------------------------------------- #
#                                    columns                                   #
# ---------------------------------------------------------------------------- #
abstract type AbstractColumn <: AbstractRenderable end

setwidth!(col::AbstractColumn, width::Int) = nothing

# ---------------------------------------------------------------------------- #
#                                 TEXT COLUMNS                                 #
# ---------------------------------------------------------------------------- #
# ---------------------------- description column ---------------------------- #
mutable struct DescriptionColumn <: AbstractColumn
    job::ProgressJob
    segments::Vector{Segment}
    measure::Measure
    text::String

    function DescriptionColumn(job::ProgressJob; style::String = orange_light)
        seg = Segment(job.description, style)
        return new(job, [seg], seg.measure, seg.text)
    end
end

update!(col::DescriptionColumn, args...)::String = col.text

# ----------------------------- separator column ----------------------------- #
struct SeparatorColumn <: AbstractColumn
    job::ProgressJob
    segments::Vector{Segment}
    measure::Measure
    text::String

    function SeparatorColumn(job::ProgressJob)
        seg = Segment("●", pink)
        return new(job, [seg], Measure(1, 1), seg.text)
    end
end

update!(col::SeparatorColumn, args...)::String = col.text

struct SpaceColumn <: AbstractColumn
    job::ProgressJob
    segments::Vector{Segment}
    measure::Measure
    text::String

    function SpaceColumn(job::ProgressJob; width = 1)
        seg = Segment(" "^width, pink)
        return new(job, [seg], seg.measure, seg.text)
    end
end

update!(col::SpaceColumn, args...)::String = col.text

# ----------------------------- completed column ----------------------------- #
struct CompletedColumn <: AbstractColumn
    job::ProgressJob
    segments::Vector{Segment}
    measure::Measure
    text::String
    padwidth::Int
    style::String

    function CompletedColumn(job::ProgressJob; style::String = "white bold")
        if isnothing(job.N)
            seg = Segment(" "^10)
            return new(job, [seg], seg.measure, "", 0, style)
        else
            width = length(string(job.N)) * 2 + 1
            seg = Segment(" "^width)
            text = apply_style(
                "{white bold}/{/white bold}{(.1, .8, .4) underline}$(job.N){/(.1, .8, .4) underline}",
            )
            return new(job, [seg], seg.measure, text, length(digits(job.N)), style)
        end
    end
end

function update!(col::CompletedColumn, color::String, args...)::String
    isnothing(col.job.N) && return apply_style(string(col.job.i), col.style)
    return apply_style(lpad(string(col.job.i), col.padwidth) * col.text, color * " bold")
end

# ----------------------------- percentage column ---------------------------- #
struct PercentageColumn <: AbstractColumn
    job::ProgressJob
    segments::Vector{Segment}
    measure::Measure

    function PercentageColumn(job::ProgressJob)
        seg = Segment(" "^4) # "xxx %
        return new(job, [seg], seg.measure)
    end
end

function update!(col::PercentageColumn, args...)::String
    isnothing(col.job.N) && return ""
    frac = int(col.job.i / col.job.N * 100)
    p = string(frac)
    return "\e[2m" * (frac == 100 ? p : lpad(p, 3)) * "%\e[0m"
end

# ----------------------------- downloaded column ---------------------------- #
struct DownloadedColumn <: AbstractColumn
    job::ProgressJob
    segments::Vector{Segment}
    measure::Measure
    tot_size::String

    function DownloadedColumn(job::ProgressJob)
        tot_size = get_file_format(job.N)
        seg = Segment(" "^(length(tot_size) * 2 + 1))
        return new(job, [seg], seg.measure, tot_size)
    end
end

function update!(col::DownloadedColumn, args...)::String
    isnothing(col.job.N) && return ""
    completed = get_file_format(col.job.i)

    return lpad(completed * "/" * col.tot_size, col.measure.w)
end

# ---------------------------------------------------------------------------- #
#                                TIMING COLUMNS                                #
# ---------------------------------------------------------------------------- #
# ------------------------------ elapsed column ------------------------------ #

struct ElapsedColumn <: AbstractColumn
    job::ProgressJob
    segments::Vector{Segment}
    measure::Measure
    style::String
    padwidth::Int

    ElapsedColumn(job::ProgressJob; style = purple_light) =
        new(job, [], Measure(6 + 9, 1), style, 6)
end

function update!(col::ElapsedColumn, args...)::String
    isnothing(col.job.startime) && return " "^(col.measure.w)
    elapsedtime = (now() - col.job.startime).value  # in ms

    # format elapsed message
    msg = if elapsedtime < 1000
        string(elapsedtime, "ms")
    elseif elapsedtime < (60 * 1000)
        # under a minute
        string(round(elapsedtime / 1000; digits = 2), "s")
    else
        # show minutes
        string(round(elapsedtime / (60 * 1000); digits = 2), "min")
    end

    msg = lpad(truncate(msg, col.padwidth), col.padwidth)
    return apply_style("elapsed: $(msg)", col.style)
end

# -------------------------------- ETA column -------------------------------- #

struct ETAColumn <: AbstractColumn
    job::ProgressJob
    segments::Vector{Segment}
    measure::Measure
    style::String
    padwidth::Int

    ETAColumn(job::ProgressJob; style = teal) = new(job, [], Measure(7 + 11, 1), style, 7)
end

function update!(col::ETAColumn, args...)::String
    isnothing(col.job.startime) && return " "^(col.measure.w)
    isnothing(col.job.N) && return " "^(col.measure.w)

    # get remaining time in ms
    elapsed = (now() - col.job.startime).value  # in ms
    perc = col.job.i / col.job.N
    remaining = elapsed * (1 - perc) / perc

    # format elapsed message
    msg = if remaining < 1000
        string(round(remaining; digits = 0), "ms")
    elseif remaining < (60 * 1000)
        # under a minute
        string(round(remaining / 1000; digits = 2), "s")
    else
        # show minutes
        string(round(remaining / (60 * 1000); digits = 2), "min")
    end

    msg = lpad(truncate(msg, col.padwidth), col.padwidth)

    return apply_style("remaining: $(msg)", col.style)
end

# ---------------------------------------------------------------------------- #
#                               PROGRESS COLUMNS                               #
# ---------------------------------------------------------------------------- #

# ------------------------------ progress column ----------------------------- #
mutable struct ProgressColumn <: AbstractColumn
    job::ProgressJob
    segments::Vector{Segment}
    measure::Measure
    nsegs::Int

    ProgressColumn(job::ProgressJob) = new(job, Vector{Segment}(), Measure(0, 0), 0)
end

function setwidth!(col::ProgressColumn, width::Int)
    col.measure = Measure(width, 1)
    return col.nsegs = width
end

function update!(col::ProgressColumn, color::String, args...)::String
    completed = int(col.nsegs * col.job.i / col.job.N)
    remaining = col.nsegs - completed
    return apply_style(
        "{" *
        color *
        " bold}" *
        '━'^(completed) *
        "{/" *
        color *
        " bold}" *
        " "^(remaining),
    )
end

# ------------------------------ spinner columns ----------------------------- #

SPINNERS = Dict(
    :dot => Dict(
        "period" => 100,
        "frames" => [
            "( ●    )",
            "(  ●   )",
            "(   ●  )",
            "(    ● )",
            "(     ●)",
            "(    ● )",
            "(   ●  )",
            "(  ●   )",
            "( ●    )",
            "(●     )",
        ],
    ),
    :circle => Dict("period" => 125, "frames" => ["◐", "◓", "◑", "◒"]),
    :toggle => Dict("period" => 250, "frames" => ["⦾⦿", "⦿⦾"]),
    :toggle2 => Dict("period" => 250, "frames" => ["◯", "⬤"]),
    :bar => Dict(
        "period" => 100,
        "frames" => [
            "(=   )",
            "(==  )",
            "(=== )",
            "( ===)",
            "(  ==)",
            "(   =)",
            "(   =)",
            "(  ==)",
            "( ===)",
            "(====)",
            "(=== )",
            "(==  )",
            "(=   )",
        ],
    ),
    :greek => Dict("period" => 350, "frames" => ["ϴ", "Ω", "Φ", "Ο"]),
)

mutable struct SpinnerColumn <: AbstractColumn
    job::ProgressJob
    segments::Vector{Segment}
    measure::Measure
    frames::Vector{String}
    Δt::Float64                 # how frequently to update display, in milliseconds
    frameidx::Int
    nframes::Int
    lastupdated::Int
    lasttext::String

    function SpinnerColumn(
        job::ProgressJob;
        spinnertype::Symbol = :dot,
        style = "bold blue",
    )
        spinnerdata = copy(SPINNERS[spinnertype])

        spinnerdata["frames"] =
            map(frame -> apply_style(frame, style), spinnerdata["frames"])

        seg = Segment(spinnerdata["frames"][1])

        return new(
            job,
            [seg],
            seg.measure,
            spinnerdata["frames"],
            spinnerdata["period"],
            1,
            length(spinnerdata["frames"]),
            0,
            spinnerdata["frames"][1],
        )
    end
end

function update!(col::SpinnerColumn, args...)::String
    col.job.started || return " "^(col.measure.w)
    col.job.finished && return "{green bold}✔{/green bold}"

    t = (now() - col.job.startime).value

    if t - col.lastupdated ≥ col.Δt
        col.lastupdated = t
        col.frameidx = col.frameidx == col.nframes ? 1 : col.frameidx + 1
        col.lasttext = col.frames[col.frameidx]
    end

    return col.lasttext
end

# ---------------------------------------------------------------------------- #
#                                COLUMNS PRESETS                               #
# ---------------------------------------------------------------------------- #

function get_columns(columnsset::Symbol)::Vector{DataType}
    return if columnsset == :minimal
        [DescriptionColumn, ProgressColumn]
    elseif columnsset == :default
        [
            DescriptionColumn,
            SeparatorColumn,
            ProgressColumn,
            SeparatorColumn,
            CompletedColumn,
            PercentageColumn,
        ]
    elseif columnsset == :spinner
        [DescriptionColumn, SpaceColumn, SpinnerColumn, SpaceColumn, CompletedColumn]
    elseif columnsset == :detailed
        # extensive
        [
            DescriptionColumn,
            SeparatorColumn,
            ProgressColumn,
            SeparatorColumn,
            CompletedColumn,
            PercentageColumn,
            SeparatorColumn,
            ElapsedColumn,
            ETAColumn,
        ]
    else
        @warn "Columns name not recognized: $columnsset"
        get_columns(:minimal)
    end
end


text:
import HTTP
import JSON

const _httpjson_github_owner = "bcbi"
const _httpjson_github_repo = "ModelSanitizer.jl"

function _httpjson__get_github_pull_request_title_without_number_authenticated(travis_pull_request::AbstractString)::String
    _travis_pull_request::String = lowercase(strip(travis_pull_request))
    if length(_travis_pull_request) == 0 || _travis_pull_request == "false"
        return ""
    else
        _github_token::String = strip(ENV["GITHUB_TOKEN"])
        pull_request_json = HTTP.request(
            "GET",
            "https://api.github.com/repos/$(github_owner)/$(github_repo)/pulls/$(_travis_pull_request)",
            [
                "User-Agent" => "bcbi-bot",
                "Authorization" => "token $(_github_token)",
                ],
            )
        pull_request_dict = JSON.parse(String(deepcopy(pull_request_json.body)))
        pull_request_title::String = strip(pull_request_dict["title"])
        return pull_request_title
    end
end

function _httpjson_get_github_pull_request_title_without_number_authenticated(d::AbstractDict)
    travis_pull_request::String = lowercase(strip(get(d, "TRAVIS_PULL_REQUEST", "false")))
    result = _httpjson__get_github_pull_request_title_without_number_authenticated(travis_pull_request)
    return result
end

function _httpjson_get_github_pull_request_title_without_number_authenticated()::String
    result::String = _httpjson_get_github_pull_request_title_without_number_authenticated(ENV)
    return result
end

function _httpjson_get_github_pull_request_title_authenticated(d::AbstractDict)::String
    travis_pull_request::String = lowercase(strip(get(d, "TRAVIS_PULL_REQUEST", "false")))
    if length(travis_pull_request) == 0 || travis_pull_request == "false"
        return ""
    else
        result_without_number::String = _httpjson_get_github_pull_request_title_without_number_authenticated(d)
        result_with_number::String = string(travis_pull_request, ": ", result_without_number)
    end
end

function _httpjson_get_github_pull_request_title_authenticated()::String
    result::String = _httpjson_get_github_pull_request_title_authenticated(ENV)
    return result
end


text:
using Statistics: median
using RobustModels: mean_and_sem

funcs = (;
          mean=mean,
          std=std,
          var=var,
          sem=sem,
          mean_and_std=mean_and_std,
          mean_and_var=mean_and_var,
          mean_and_sem=mean_and_sem
        )

xorig = [
  0.529,
 -0.921,
 -0.247,
  0.689,
 -0.254,
 -1.087,
  0.067,
  0.679,
  0.089,
 -0.278,
  1.327,
 -0.571,
 -0.559,
 -0.491,
  0.237,
  1.196,
 -0.807,
  0.35 ,
  0.259,
  1.006,
  0.257,
]
#x = randn(100)

x = copy(xorig)
x[1:2] .= 50

@testset "robust univariate statistics: MEstimator{$(lossname)}" for lossname in ("L2Loss", "L1Loss", "HuberLoss", "TukeyLoss")
    typeloss = getproperty(RobustModels, Symbol(lossname))
    est = MEstimator{typeloss}()

    @testset "method: $(name)" for (name, func) in pairs(funcs)
        if lossname == "L1Loss"
            @test_throws ArgumentError func(est, x)
            continue
        end

        # Compute robust statistics
        s = func(est, x)
        if lossname == "L2Loss"
            @test all(isapprox.(s, func(x); rtol=1e-7))
        elseif lossname == "TukeyLoss"
            @test all(s .<= func(x))
            # the estimate is better than when the outliers are removed...
            @test_skip all(func(xorig) .<= s)
        else
            @test all(func(xorig) .<= s .<= func(x))
        end
    end
end

@testset "robust univariate statistics: Bounded estimator $(typeest)" for typeest in (SEstimator, MMEstimator, TauEstimator)
    est = typeest{TukeyLoss}()

    @testset "method: $(name)" for (name, func) in pairs(funcs)
        s = func(est, x)
#        println("statistics $(name): $(round.(s; digits=4)) ≈ $(round.(func(xorig); digits=4)) (with outliers removed)")
        @test all(isapprox.(s, func(xorig); rtol=2))
    end
end

# define the Estimator for the next tests
est = MEstimator{L2Loss}()

x1 = reshape(xorig, (21, 1))
x2 = reshape(xorig, (1, 21))
x3 = reshape(xorig, (1, 3, 7))
x4 = reshape(xorig, (7, 3, 1))

@testset "robust univariate statistics: Array size: $(size(a))" for a in (x, x1, x2, x3, x4)
    @testset "method: $(name)" for (name, func) in pairs(funcs)
        @testset "dims=$(dims)" for dims in (1, 2, (1,), (1,2), (3,1), 4, (:))
            s = func(est, a; dims=dims)
            if name in (:mean, :std, :var)
                res = func(a; dims=dims)
                @test all(isapprox.(s, res; rtol=1e-7, nans=true))
            elseif name in (:mean_and_std, :mean_and_var)
                # tuples are not allowed for the `dims` arg in StatsBase.mean_and_var
                if isa(dims, Tuple) && length(dims) > 1
                    continue
                end
                res = if dims===(:); func(a) else func(a, first(dims)) end
                @test all(isapprox.(s, res; rtol=1e-7, nans=true))
            end
        end
    end
end

d = Base.values(Dict(:a=>1, :b=>2, :c=>3))
g = (i for i in 1:3)

@testset "robust univariate statistics: iterable $(typeof(a).name.name)" for a in (d, g)
    @test all(isapprox.(mean(est, a), mean(a); rtol=1e-7))
end


text:
function _create_model(rd::RoutingData)
    m = Model(rd.solver)
    set_silent(m) # TODO: add a parameter for this
    return m
end


text:
data_list = ["Cernet" , "Globalcenter" , "Goodnet" , "Gridnet" , "gscale" , "Janetbackbone" ,
"Rediris" , "Renater2010" , "Sinet" , "Sprint" , "SprintINET" , "Surfnet" ,
"SwitchL3" , "TataNld" , "Tinet" , "Uninett2011" , "Uunet" , "VtlWavenet2011" , "Xeex" , "Xspedius"]

function gen_gravity(data)
    hosts = readlines("data/$data.hosts")
    open("data/$data.demands", "w") do fout
        for k in 1:60
            #α = 2^(23 + k)
            α = 2^(24 + floor(k/10))
            weigths = map(x->rand(), hosts)
            list = [ α * weigths[i] * weigths[j]
                     for i in 1:length(hosts)
                     for j in 1:length(hosts) ]
            println(fout, join(list, ' '))
        end
    end
end
for i in data_list
    gen_gravity(i)
end

max_stars_repo_name: ylxdzsw/TEexp

max_issues_repo_name: ylxdzsw/TEexp

max_forks_repo_name: ylxdzsw/TEexp


text:
import Base: run, show
export update_build_image, DockerRunner, run, runshell

# The docker image we use
const BUILD_IMAGE = "staticfloat/julia_crossbuild:x64"
BUILD_IMAGE_UPDATED = false

"""
    should_update_build_image(force::Bool)

Helper function to determine whether we shuold or should not attempt to update
the build image containing all cross-compiler toolchains. To override automatic
updating, set the environment variable `BINBUILD_IMAGE_AUTOUPDATE` to `"true"`
or `"false"`.  To force updating always, set `force` to `true`.
"""
function should_update_build_image(force::Bool)
    global BUILD_IMAGE_UPDATED

    # If we're forcing the issue, then return true
    if force
        return true
    end

    # If we've already updated, then don't bother to
    if BUILD_IMAGE_UPDATED
        return false
    end

    # If the user has explicitly preferred us not to, then don't
    const no_synonyms = ["n", "no", "false"]
    if lowercase(get(ENV, "BINBUILD_IMAGE_AUTOUPDATE", "y")) in no_synonyms
        return false
    end

    # Otherwise, DEWIT!
    return true
end

"""
    update_build_image(; verbose::Bool = false, force::Bool = false)

Updates the build image containing all cross-compiler toolchains.  Checks for
updates upon first attempt to run a build by default.  Set `force` to `true`
to force an update, and to disable automatic updates (for instance if you are
sitting in a coffee shop desperately trying to test your build scripts and you
don't particularly feel like waiting for download the new 4GB+ image that
`@staticfloat` just pushed up) set the environment variable
`BINBUILD_IMAGE_AUTOUPDATE` to `"false"` before importing `BinaryBuilder`.
"""
function update_build_image(; verbose::Bool = false, force::Bool = false)
    global BUILD_IMAGE_UPDATED
    if should_update_build_image(force)
        msg = """
        Updating build image $BUILD_IMAGE, this may take a few minutes. To
        disable automatic image updating, set the `BINBUILD_IMAGE_AUTOUPDATE`
        environment variable to `"false"` in your shell environment or in your
        `~/.juliarc.jl` file.
        """
        info(replace(strip(msg), "\n", " "))

        oc = OutputCollector(`docker pull $BUILD_IMAGE`; verbose=verbose)
        did_succeed = wait(oc)
        if !did_succeed
            error("Could not update build image $BUILD_IMAGE")
        end

        BUILD_IMAGE_UPDATED = true
    end
end

"""
    DockerRunner

A `DockerRunner` represents an "execution context", an object that bundles all
necessary information to run commands within the docker container that contains
our crossbuild environment.  Use `run()` to actually run commands within the
`DockerRunner`, and `runshell()` as a quick way to get an interactive shell
within the crossbuild environment.
"""
type DockerRunner
    cmd_prefix::Cmd
    platform::Symbol
end

function show(io::IO, x::DockerRunner)
    write(io, "$(x.platform) DockerRunner")
end

"""
    getuid()

Wrapper function around libc's `getuid()` function
"""
function getuid()
    return ccall((:getuid, :libc), UInt32, ())
end

"""
    getgid()

Wrapper function around libc's `getgid()` function
"""
function getgid()
    return ccall((:getgid, :libc), UInt32, ())
end

"""
    target_envs(target::String)

Given a `target` (this term is used interchangeably with `triplet`), generate a
`Dict` mapping representing all the environment variables to be set within the
build environment to force compiles toward the defined target architecture.

Examples of things set are `PATH`, `CC`, `RANLIB`, as well as nonstandard
things like `target`.
"""
function target_envs(target::AbstractString)
    target_tool = tool -> "/opt/$(target)/bin/$(target)-$(tool)"
    mapping = Dict(
        "PATH" => "/opt/$(target)/bin:/usr/local/bin:/usr/local/sbin:/usr/bin:/usr/sbin:/bin:/sbin",
        "AR" => target_tool("ar"),
        "CC" => target_tool("gcc"),
        "CXX" => target_tool("g++"),
        "FC" => target_tool("gfortran"),
        "LD" => target_tool("ld"),
        "NM" => target_tool("nm"),
        "AS" => target_tool("as"),
        "RANLIB" => target_tool("ranlib"),
        "STRIP" => target_tool("strip"),
        "INSTALL_NAME_TOOL" => target_tool("install_name_tool"),
        "LIBTOOL" => target_tool("libtool"),
        "LIPO" => target_tool("lipo"),
        "OTOOL" => target_tool("otool"),
        "target" => target,

        # Autotools really appreciates being able to build stuff for the
        # host system, so we set this to ease its woes
        "CC_FOR_BUILD" => "/opt/x86_64-linux-gnu/bin/gcc",
    )

    return mapping
end

"""
    DockerRunner(; prefix::Prefix = global_prefix,
                   platform::Symbol = platform_key(),
                   volume_mapping::Vector = [])

Creates a `DockerRunner` object to run commands within the environment defined
by the given `prefix`, `platform` and `volume_mapping`s.  The `prefix` given
will be mounted into the Docker image at the same path as it exists within the
host file system, and any extra volumes that should be mapped in can be done
so with `volume_mapping`, which expects tuples of paths in a similar spirit to
`docker`'s `-v` option.
"""
function DockerRunner(;prefix::Prefix = BinaryProvider.global_prefix,
                       platform::Symbol = platform_key(),
                       volume_mapping::Vector = [])
    # We are using `docker run` to provide isolation
    cmd_prefix = `docker run --rm -i`

    # The volumes we'll always map into the docker container
    push!(volume_mapping, (prefix.path, prefix.path))
    for v in volume_mapping
        cmd_prefix = `$cmd_prefix -v $(v[1]):$(v[2])`
    end

    # The environment variables we'll set
    env_mapping = target_envs(platform_triplet(platform))
    for v in env_mapping
        cmd_prefix = `$cmd_prefix -e $(v[1])=$(v[2])`
    end

    # Set our user id and group id to match the outside world
    cmd_prefix = `$cmd_prefix --user=$(getuid()):$(getgid())`

    # Don't run with a confined seccomp profile
    cmd_prefix = `$cmd_prefix --security-opt seccomp=unconfined`

    # Manually set DESTDIR environment variable
    cmd_prefix = `$cmd_prefix -e DESTDIR=$(prefix.path)`

    # Actually update the build image, if we need to
    update_build_image()

    return DockerRunner(cmd_prefix, platform)
end

"""
    run(dr::DockerRunner, cmd::Cmd, logpath::String; verbose::Bool = false)

Given a `DockerRunner`, runs `cmd` within the docker environment, storing any
output logs into `logpath` and returning `false` if the command did not
complete successfully.  This command will also mount in the current directory
into the docker environment.
"""
function run(dr::DockerRunner, cmd::Cmd, logpath::AbstractString; verbose::Bool = false)
    # Create the directory where we'll store logs, if we need to
    mkpath(dirname(logpath))

    # Run the command
    d = pwd()
    user_cmd = `$(dr.cmd_prefix) -w $(d) -v $(d):$(d) $BUILD_IMAGE $cmd`

    oc = OutputCollector(user_cmd; verbose=verbose)
    did_succeed = wait(oc)

    # Write out the logfile, regardless of whether it was successful or not
    mkpath(dirname(logpath))
    open(logpath, "w") do f
        # First write out the actual command, then the command output
        println(f, cmd)
        print(f, merge(oc))
    end

    # Return whether we succeeded or not
    return did_succeed
end

"""
    runshell(dr::DockerRunner)

Open an interactive session inside a Docker environment defined by a
`DockerRunner` object with the current directory mapped in as well.
"""
function runshell(dr::DockerRunner)
    d = pwd()
    user_cmd = `$(dr.cmd_prefix) -w $(d) -v $(d):$(d) -t $BUILD_IMAGE bash`
    run(user_cmd)
end

"""
    runshell(platform::Symbol)

Open an interactive session inside a Docker environment created for a
particular target `platform`.
"""
function runshell(platform::Symbol)
    runshell(DockerRunner(platform=platform))
end


text: using Test, Bench


text:
n = 200
m = 100
p = 25

info("Testing ", MOD.pairwise)
for kernelobj in additive_kernels
    for T in FloatingPointTypes
        x = rand(T,p)
        y = rand(T,p)
        k = convert(Kernel{T}, (kernelobj)())

        @test MOD.pairwise(k, x[1], y[1]) == MOD.phi(k, x[1], y[1])
        @test MOD.pairwise(k, x, y)       == sum(map((x,y) -> MOD.phi(k,x,y), x, y))
    end
end

info("Testing ", MOD.pairwisematrix!)
steps = length(additive_kernels) + 1
counter = 0
info("    Progress:   0.0%")
for kernelobj in (additive_kernels..., MLKTest.PairwiseTestKernel)
    for T in FloatingPointTypes
        Set_X = [rand(T, p) for i = 1:n]
        Set_Y = [rand(T,p)  for i = 1:m]

        X = transpose(hcat(Set_X...))
        Y = transpose(hcat(Set_Y...))
        k = convert(Kernel{T}, (kernelobj)())
        #println("κ:", k, ", X:", typeof(X), ", Y:", typeof(Y))

        P = [MOD.pairwise(k,x,y) for x in Set_X, y in Set_X]
        #@test_approx_eq MOD.pairwisematrix(Val{:row}, k, X)  P
        #@test_approx_eq MOD.pairwisematrix(Val{:col}, k, X') P
        @test_approx_eq MOD.pairwisematrix!(Val{:row}, Array(T,n,n), k, X,  true)  P
        @test_approx_eq MOD.pairwisematrix!(Val{:col}, Array(T,n,n), k, X', true) P

        P = [MOD.pairwise(k,x,y) for x in Set_X, y in Set_Y]
        #@test_approx_eq MOD.pairwisematrix(Val{:row}, k, X,  Y)  P
        #@test_approx_eq MOD.pairwisematrix(Val{:col}, k, X', Y') P
        @test_approx_eq MOD.pairwisematrix!(Val{:row}, Array(T,n,m), k, X,  Y)  P
        @test_approx_eq MOD.pairwisematrix!(Val{:col}, Array(T,n,m), k, X', Y') P

        @test_throws DimensionMismatch MOD.pairwisematrix!(Val{:row}, Array(T,n+1,n),   k, X, true)
        @test_throws DimensionMismatch MOD.pairwisematrix!(Val{:row}, Array(T,n,n+1),   k, X, true)
        @test_throws DimensionMismatch MOD.pairwisematrix!(Val{:row}, Array(T,n+1,n+1), k, X, true)

        @test_throws DimensionMismatch MOD.pairwisematrix!(Val{:col}, Array(T,n+1,n),   k, X', true)
        @test_throws DimensionMismatch MOD.pairwisematrix!(Val{:col}, Array(T,n,n+1),   k, X', true)
        @test_throws DimensionMismatch MOD.pairwisematrix!(Val{:col}, Array(T,n+1,n+1), k, X', true)

        @test_throws DimensionMismatch MOD.pairwisematrix!(Val{:row}, Array(T,n+1,m), k, X,  Y)
        @test_throws DimensionMismatch MOD.pairwisematrix!(Val{:row}, Array(T,n,m+1), k, X,  Y)
        @test_throws DimensionMismatch MOD.pairwisematrix!(Val{:col}, Array(T,n+1,m), k, X', Y')
        @test_throws DimensionMismatch MOD.pairwisematrix!(Val{:col}, Array(T,n,m+1), k, X', Y')

        @test_throws DimensionMismatch MOD.pairwisematrix!(Val{:row}, Array(T,m+1,n), k, Y,  X)
        @test_throws DimensionMismatch MOD.pairwisematrix!(Val{:row}, Array(T,m,n+1), k, Y,  X)
        @test_throws DimensionMismatch MOD.pairwisematrix!(Val{:col}, Array(T,m+1,n), k, Y', X')
        @test_throws DimensionMismatch MOD.pairwisematrix!(Val{:col}, Array(T,m,n+1), k, Y', X')

        @test_approx_eq MOD.dotvectors(Val{:row}, X) sum((X.*X),2)
        @test_approx_eq MOD.dotvectors(Val{:col}, X) sum((X.*X),1)

        @test_approx_eq MOD.dotvectors!(Val{:row}, Array(T,n), X) sum((X.*X),2)
        @test_approx_eq MOD.dotvectors!(Val{:col}, Array(T,p), X) sum((X.*X),1)

        @test_throws DimensionMismatch MOD.dotvectors!(Val{:row}, Array(T,2), Array(T,3,2))
        @test_throws DimensionMismatch MOD.dotvectors!(Val{:row}, Array(T,4), Array(T,3,4))
        @test_throws DimensionMismatch MOD.dotvectors!(Val{:col}, Array(T,2), Array(T,2,3))
        @test_throws DimensionMismatch MOD.dotvectors!(Val{:col}, Array(T,4), Array(T,4,3))

        @test_throws DimensionMismatch MOD.squared_distance!(Array(T,3,4), Array(T,3), true)
        @test_throws DimensionMismatch MOD.squared_distance!(Array(T,4,3), Array(T,3), true)

        @test_throws DimensionMismatch MOD.squared_distance!(Array(T,3,4), Array(T,2), Array(T,4))
        @test_throws DimensionMismatch MOD.squared_distance!(Array(T,3,4), Array(T,4), Array(T,4))
        @test_throws DimensionMismatch MOD.squared_distance!(Array(T,3,4), Array(T,3), Array(T,3))
        @test_throws DimensionMismatch MOD.squared_distance!(Array(T,3,4), Array(T,3), Array(T,5))
    end
    counter += 1
    info("    Progress: ", @sprintf("%5.1f", counter/steps*100), "%")
end


text:
function ContJitter(n, d, num)

    base_obn = convert(Array{Float32},range(d[1],stop=(n[1]-1)*d[1],length = num)) # spacing = 250m on average
    avg_spc = (n[1]-2)*d[1]/(num-1)
    OBN_loc = zeros(Float32,num)

    for i = 1:num
        OBN_loc[i] = max(min(base_obn[i] + avg_spc*rand(Float32) - avg_spc/2, (n[1]-1)*d[1]),0f0)
    end

    return OBN_loc
end

function Jitt_function(N, gamma, xi)
    # 2020.12.09 Yijun
    # jittered subsampling function for simple test
    # N: the full size
    # gamma: undersampling factor
    # xi: control the maximum gap between the neighbor subsampling indexes
    # gamma=4 --> 75 % subsampling
    #rng(rseed);
    test1 = rand(1,Int64(round(N/gamma))).*(xi-1) .- (xi-1)/2 ;
    test2 = ones(1, Int64(round(N/gamma)));
    test3 = cumsum(test2; dims=2);
    n = Int64.(round.((1-gamma)/2 .+ gamma .* test3 .+ test1));
    for i = 1:length(n)
        if n[i] > N
            n[i] = N
        end
    end
    return n
end


text:
#Requires place_id
function get_geo_id_place_id(auth::TwAuth, options::Dict)
  place_id = options["place_id"]
  Requests.json(get_oauth(auth, "https://api.twitter.com/1.1/geo/id/$place_id.json", options))
end

#Requires lat, long
function get_geo_reverse_geocode(auth::TwAuth, options::Dict)
  Requests.json(get_oauth(auth, "https://api.twitter.com/1.1/geo/reverse_geocode.json", options))
end

function get_geo_search(auth::TwAuth, options::Dict)
  Requests.json(get_oauth(auth, "https://api.twitter.com/1.1/geo/search.json", options))
end

max_stars_repo_name: nick-paul/JuliaTwitter

max_issues_repo_name: nick-paul/JuliaTwitter

max_forks_repo_name: nick-paul/JuliaTwitter


text:
#=##############################################################################
# DESCRIPTION
    Vehicle with vortex-lattice method (VLM) and blade-element momentum (BEM)
    models assuming quasi-steady wakes. This means that at every time step the
    wake of VLM models are represented as rigid, semi-infinite filaments, and
    wakes of BEM models are obtained from a conservation of momentum
    assumptions. In this approach only VLM-on-VLM and VLM-on-BEM interactions
    are accounted for, meaning that the wake of each BEM rotor has no effect
    on other rotors and wings.

# AUTHORSHIP
  * Author    : Eduardo J. Alvarez
  * Email     : Edo.AlvarezR@gmail.com
  * Created   : Apr 2020
  * License   : MIT
=###############################################################################

################################################################################
# QUASI-STEADY VLM VEHICLE TYPE
################################################################################
"""
    `QVLMVehicle(system; optargs...)`

Type handling all geometries and subsystems that define a flight vehicle made
out of VLM (Wing, WingSystem, Rotor) components.

NOTE: For the solver to work correctly, all components in `wake_system` (if any)
need to be also components of `vlm_system`.

NOTE: It is recommended that `wake_system` doesn't include any Rotor object.
Otherwise, blades will generate a wake going straight out of every blade
trailing edge pointing oposite to the direction of rotation instead of
generating a streamtube.

# ARGUMENTS
* `system::vlm.WingSystem`:        System of all FLOWVLM objects. This system
                                    is considered as the entire vehicle. Not all
                                    components in this system will be solved,
                                    but they will all be rotated and translated
                                    during maneuver.
# OPTIONAL ARGUMENTS
* `tilting_systems::Tuple(vlm.WingSystem, ...)`:   Tuple of all FLOWVLM
                                    tilting objects, where `tilting_systems[i]`
                                    contains the i-th FLOWVLM system of lifting
                                    surfaces and rotors that tilt together.
* `rotors_systems::Tuple(Array{vlm.Rotor,1}, ...)`:   Tuple of groups of Rotors
                                    that share a common RPM.
* `vlm_system::vlm.WingSystem`:    System of all FLOWVLM objects to be solved
                                    through the VLM solver.
* `wake_system::vlm.WingSystem`:   System of all FLOWVLM objects that will
                                    have a rigid, semi-infinite VLM wake.
* `grids::Array{gt.GridTypes, 1}`: Array of grids that will be translated and
                                    rotated along with `system`.
"""
struct QVLMVehicle{N, M, R} <: AbstractVLMVehicle{N, M, R}

    # Required inputs
    system::vlm.WingSystem

    # Optional inputs
    tilting_systems::NTuple{N, vlm.WingSystem}
    rotor_systems::NTuple{M, Array{vlm.Rotor, 1}}
    vlm_system::vlm.WingSystem
    wake_system::vlm.WingSystem
    grids::Array{gt.GridTypes, 1}

    # Internal properties
    V::Array{R, 1}                          # Current vehicle velocity
    W::Array{R, 1}                          # Current vehicle angular velocity
    prev_data::Array{Any, 1}                # Information about previous step
    grid_O::Array{Array{R, 1}, 1}           # Origin of every grid

    QVLMVehicle{N, M, R}(
                    system;
                    tilting_systems=NTuple{0, vlm.WingSystem}(),
                    rotor_systems=NTuple{0, Array{vlm.Rotor, 1}}(),
                    vlm_system=vlm.WingSystem(),
                    wake_system=vlm.WingSystem(),
                    grids=Array{gt.GridTypes, 1}(),
                    V=zeros(3), W=zeros(3),
                    prev_data=[deepcopy(vlm_system), deepcopy(wake_system),
                                                    deepcopy(rotor_systems)],
                    grid_O=Array{Array{Float64, 1}, 1}(),
                ) where {N, M, R} = new(
                    system,
                    tilting_systems,
                    rotor_systems,
                    vlm_system,
                    wake_system,
                    grids,
                    V, W,
                    prev_data,
                    grid_O,
                )
end

# Implicit N and M constructor
QVLMVehicle(system::vlm.WingSystem;
        V::Array{R, 1}=zeros(3), W::Array{R, 1}=zeros(3),
        tilting_systems::NTuple{N, vlm.WingSystem}=NTuple{0, vlm.WingSystem}(),
        rotor_systems::NTuple{M, Array{vlm.Rotor, 1}}=NTuple{0, Array{vlm.Rotor, 1}}(),
        grids=Array{gt.GridTypes, 1}(),
        optargs...
        ) where {N, M, R} = QVLMVehicle{N, M, R}( system;
                                V=V, W=W,
                                tilting_systems=tilting_systems,
                                rotor_systems=rotor_systems,
                                grids=Array{gt.GridTypes, 1}(grids),
                                grid_O=[zeros(R, 3) for i in 1:length(grids)],
                                optargs...)

##### FUNCTIONS  ###############################################################
function shed_wake(self::QVLMVehicle, Vinf::Function, pfield::vpm.ParticleField,
                                                dt::Real, nt::Int; optargs...)
    nothing
end

function generate_static_particle_fun(pfield, self::QVLMVehicle, sigma;
                                                                     optargs...)

    function static_particles_function(args...)
        return nothing
    end

    return static_particles_function
end

function save_vtk(self::QVLMVehicle, filename; path=nothing, num=nothing,
                                                save_wopwopin=false, optargs...)

    # Filter optional arguments for wake
    optargs1 = [(sym, val) for (sym, val) in optargs
                         if !(sym in [:infinite_vortex, :save_horseshoes,
                                      :only_horseshoes, :only_infinite_vortex])]

    # Save wake
    strn = vlm.save(self.wake_system, filename; suff="_wake",
                                infinite_vortex=true, save_horseshoes=true,
                                only_horseshoes=true, only_infinite_vortex=true,
                                path=path, num=num, optargs1...)

    # Filter optional arguments of VLM
    optargs2 = [(sym, val) for (sym, val) in optargs if sym!=:infinite_vortex]

    # Save VLM
    strn *= save_vtk_base(self, filename; path=path, num=num,
                                    save_wopwopin=save_wopwopin,
                                    infinite_vortex=false, optargs2...)
    return strn
end
##### INTERNAL FUNCTIONS  ######################################################

##### END OF VEHICLE ###########################################################


text:
using SbmlInterface
using ModelingToolkit
using OrdinaryDiffEq
using Test

@testset "SbmlInterface.jl" begin
    @testset "sbml2modelingtoolkit.jl" begin include("sbml2modelingtoolkit.jl") end
end


text:
# generate examples
using Literate

# Directory where the examples in Literate.jl format are
EXAMPLEDIR = joinpath(@__DIR__, "src", "examples")
# Directory where the examples will be generated
GENERATEDDIR = joinpath(@__DIR__, "src", "examples", "generated")

for example in readdir(EXAMPLEDIR)
    ~endswith(example, ".jl") && continue
    input = abspath(joinpath(EXAMPLEDIR, example))
    script = Literate.script(input, GENERATEDDIR)
    code = strip(read(script, String))
    mdpost(str) = replace(str, "@__CODE__" => code)
    Literate.notebook(input, GENERATEDDIR, execute = true)
    Literate.markdown(input, GENERATEDDIR, postprocess = mdpost)
    # Altenrate comments to avoid running the code when editing the text only!
    #Literate.notebook(input, GENERATEDDIR, execute = false)
    #Literate.markdown(input, GENERATEDDIR, documenter = false)
end


text:
using YaoBlocks.AD
using SparseArrays, LuxurySparse
using YaoArrayRegister
using Test

@testset "outer prod" begin
    D = 8
    T = ComplexF64
    for (l, r) in [(randn(T, D), randn(T, D)), (randn(T, D, 3), randn(T, D, 3))]
        a = outerprod(l, r)
        A = Matrix(a)
        B = randn(T, D, D)
        V = randn(T, D)
        @test A * B ≈ a * B
        @test A * A ≈ a * a
        @test B * A ≈ B * a
        @test V' * A ≈ V' * a
        @test A * V ≈ a * V
        @test size(a) == size(A)
        @test size(a, 2) == size(A, 2)
        @test A[4] == a[4]
        @test A[4, 4] == a[4, 4]
        @test A' == a'
        @test transpose(A) == transpose(a)
        @test conj(A) == conj(a)
        @test conj!(A) == conj!(a)
        outerprod(ArrayReg(l), ArrayReg(r)) == a
    end
end

@testset "projection" begin
    T = ComplexF64
    D = 8
    for y in [pmrand(T, D), sprand(T, D, D, 0.5), Diagonal(randn(T, D))]
        @test projection(y, y) == y
        l, r = randn(T, D), randn(T, D)
        op = outerprod(l, r)
        @test projection(y, op) == projection(y, Matrix(op))
    end
end


text:
using Zygote:@adjoint
using StatsBase

function sort_(x::AbstractArray; by=identity, rev = true)
    p = sortperm(x, by=by, rev = rev)
    return x[p]
end

@adjoint function sort_(x::AbstractArray; by=identity, rev = true)
    p = sortperm(x, by=by, rev = rev)
    return x[p], x̄ -> (x̄[invperm(p)],)
end

"""
    top_k_logits(logits::AbstractArray, k; prob = false)

Masks everything but the k top entries as -infinity (1e10). Incase of `probs=true`, everthing except top-k probabilities are masked to 0.0. `logits` is expected to be a vector.

"""
function top_k_logits(logits::AbstractArray, k; prob = false)
    if k==0
    	return logits
    else
        top_kth = sort_(logits;rev=true)[k]
        mask = logits.>=top_kth
        if prob==false
            return logits_ = map(*, logits, mask) + map(~, mask).*-1e-10
        else
            return logits_ = map(*, logits, mask)
        end
    end
end

function temp_softmax(logits; t=1.2)
  return softmax(logits ./ t)
end

"""
    top_k_sample(probs; k=10)

Sampling function to return index from `top_k` probabilities, based on provided `k`. Function removes all tokens with a probability less than the last token of the top_k before sampling.
"""
function top_k_sample(probs; k=10)
  probs ./=sum(probs)
  sorted = sort(probs, rev = true)
  indexes = partialsortperm(probs, 1:k, rev=true)
  sorted_k = sorted[1:k] ./ sum(sorted[1:k])
  index = sample(indexes, ProbabilityWeights(sorted_k), 1)
  return index[1]
end

"""
    nucleus_sample(probs; p=0.8)

Nuclues sampling function, to return after sampling reverse sorted probabilities `probs` till the index, where cumulative probability remains less than provided `p`. It removes tokens with cumulative probability above the threshold `p` before sampling.
"""
function nucleus_sample(probs; p=0.8)
    probs ./=sum(probs)
    sorted = sort(probs, rev = true)
    indexes = sortperm(probs, rev=true)
    cusum = cumsum(sorted)
    upto_threshold = cusum .<= p
    sorted .*= upto_threshold
    sorted ./= sum(sorted)
    index = sample(indexes, ProbabilityWeights(sorted),1)
    return index[1]
end

# Onecold
function onecold(y)
    onecold_arr = zeros(size(y,2))
    argmax_vec = argmax(y, dims = 1)
    for i in 1:size(y,2)
        onecold_arr[i] = argmax_vec[i].I[1]
    end
    return onecold_arr
end

"""
    binary_accuracy(y_pred, y_true; threshold=0.5)

Calculates Averaged Binary Accuracy based on `y_pred` and `y_true`. Argument `threshold` is used to specify the minimum predicted probability `y_pred` required to be labelled as `1`. Default value set as `0.5`.
"""
function binary_accuracy(y_pred, y_true; threshold=0.5)
    @assert size(y_pred) == size(y_true)
    return sum((y_pred .>= threshold) .== y_true) / size(y_true, 1)
end

"""
    categorical_accuracy(y_pred, y_true)

Calculates Averaged Categorical Accuracy based on `y_pred` and `y_true`.
"""
function categorical_accuracy(y_pred, y_true)
    @assert size(y_pred) == size(y_true)
    return sum(onecold(y_pred) .== onecold(y_true)) / size(y_true, 2)
end

# TODO Implement Kldivergence loss, f1_score


text:
"""
```
corners = imcorner(img; [method])
corners = imcorner(img, threshold, percentile; [method])
```

Performs corner detection using one of the following methods -

    1. harris
    2. shi_tomasi
    3. kitchen_rosenfeld

The parameters of the individual methods are described in their documentation. The
maxima values of the resultant responses are taken as corners. If a threshold is
specified, the values of the responses are thresholded to give the corner pixels.
The threshold is assumed to be a percentile value unless `percentile` is set to false.
"""
function imcorner(img::AbstractArray; method::Function = harris, args...)
    img_gray = convert(Array{Gray}, img)
    responses = method(img_gray; args...)
    corners = falses(size(img))
    maxima = map(CartesianIndex{2}, findlocalmaxima(responses))
    for m in maxima corners[m] = true end
    corners
end

function imcorner(img::AbstractArray, threshold, percentile; method::Function = harris, args...)
    img_gray = convert(Array{Gray}, img)
    responses = method(img_gray; args...)

    if percentile == true
        threshold = StatsBase.percentile(responses[:], threshold * 100)
    end

    corners = map(i -> i > threshold, responses)
    corners
end

"""
```
harris_response = harris(img; [k], [border], [weights])
```

Performs Harris corner detection. The covariances can be taken using either a mean
weighted filter or a gamma kernel.
"""
function harris(img::AbstractArray; k::Float64 = 0.04, args...)
    cov_xx, cov_xy, cov_yy = gradcovs(img, args...)
    corner = map((xx, yy, xy) -> xx * yy - xy ^ 2 - k * (xx + yy) ^ 2, cov_xx, cov_yy, cov_xy)
    corner
end

"""
```
shi_tomasi_response = shi_tomasi(img; [border], [weights])
```

Performs Shi Tomasi corner detection. The covariances can be taken using either a mean
weighted filter or a gamma kernel.
"""
function shi_tomasi(img::AbstractArray; border::AbstractString = "replicate", args...)
    cov_xx, cov_xy, cov_yy = gradcovs(img, border; args...)
    corner = map((xx, yy, xy) -> ((xx + yy) - (sqrt((xx - yy) ^ 2 + 4 * xy ^ 2))) / 2, cov_xx, cov_yy, cov_xy)
    corner
end

"""
```
kitchen_rosenfeld_response = kitchen_rosenfeld(img; [border])
```

Performs Kitchen Rosenfeld corner detection. The covariances can be taken using either a mean
weighted filter or a gamma kernel.
"""
function kitchen_rosenfeld(img::AbstractArray; border::AbstractString = "replicate")
    (grad_x, grad_y) = imgradients(img, "sobel", border)
    (grad_xx, grad_xy) = imgradients(grad_x, "sobel", border)
    (grad_yx, grad_yy) = imgradients(grad_y, "sobel", border)
    numerator = map((x, y, xx, xy, yy) -> xx * (y ^ 2) + yy * (x ^ 2) - 2 * xy * x * y, grad_x, grad_y, grad_xx, grad_xy, grad_yy)
    denominator = map((x, y) -> x ^ 2 + y ^ 2, grad_x, grad_y)
    corner = map((n, d) -> d == 0.0 ? 0.0 : -n / d, numerator, denominator)
    corner
end

"""
```
corners = fastcorners(img, n, threshold)
```

Performs FAST Corner Detection. `n` is the number of contiguous pixels
which need to be greater (lesser) than intensity + threshold (intensity - threshold)
for a pixel to be marked as a corner. The default value for n is 12.
"""
function fastcorners{T}(img::AbstractArray{T}, n::Int = 12, threshold::Float64 = 0.15)
    img_padded = padarray(img, [3,3], [3,3], "value", 0)
    corner = falses(size(img))
    h, w = size(img)
    R = CartesianRange(CartesianIndex((4, 4)), CartesianIndex((h + 3, w + 3)))
    idx = map(CartesianIndex{2}, [(0, 3), (1, 3), (2, 2), (3, 1), (3, 0), (3, -1), (2, -2), (1, -3),
            (0, -3), (-1, -3), (-2, -2), (-3, -1), (-3, 0), (-3, 1), (-2, 2), (-1, 3)])

    idxidx = [1, 5, 9, 13]
    for I in R
        bright_threshold = img_padded[I] + threshold
        dark_threshold = img_padded[I] - threshold
        if n >= 12
            sum_bright = 0
            sum_dark = 0
            for k in idxidx
                pixel = img_padded[I + idx[k]]
                if pixel > bright_threshold
                    sum_bright += 1
                elseif pixel < dark_threshold
                    sum_dark += 1
                end
            end
            if sum_bright < 3 && sum_dark < 3
                continue
            end
        end
        consecutive_bright = 0
        consecutive_dark = 0

        for i in 1:15 + n
            k = mod1(i, 16)
            pixel = img_padded[I + idx[k]]
            if pixel > bright_threshold
                consecutive_dark = 0
                consecutive_bright += 1
            elseif pixel < dark_threshold
                consecutive_bright = 0
                consecutive_dark += 1
            end

            if consecutive_dark == n || consecutive_bright == n
                corner[I - CartesianIndex((3, 3))] = true
                break
            end
        end
    end
    corner
end

function gradcovs(img::AbstractArray, border::AbstractString = "replicate"; weights::Function = meancovs, args...)
    (grad_x, grad_y) = imgradients(img, "sobel", border)

    cov_xx = grad_x .* grad_x
    cov_xy = grad_x .* grad_y
    cov_yy = grad_y .* grad_y

    weights(cov_xx, cov_xy, cov_yy, args...)
end

function meancovs(cov_xx, cov_xy, cov_yy, blockSize::Int = 3)

    box_filter_kernel = (1 / (blockSize * blockSize)) * ones(blockSize, blockSize)

    filt_cov_xx = imfilter(cov_xx, box_filter_kernel)
    filt_cov_xy = imfilter(cov_xy, box_filter_kernel)
    filt_cov_yy = imfilter(cov_yy, box_filter_kernel)

    filt_cov_xx, filt_cov_xy, filt_cov_yy
end

function gammacovs(cov_xx, cov_xy, cov_yy, gamma::Float64 = 1.4)

    filt_cov_xx = imfilter_gaussian(cov_xx, [gamma, gamma])
    filt_cov_xy = imfilter_gaussian(cov_xy, [gamma, gamma])
    filt_cov_yy = imfilter_gaussian(cov_yy, [gamma, gamma])

    filt_cov_xx, filt_cov_xy, filt_cov_yy
end


text:
module SpanningTree

using Kombinator

using DataStructures
using JuMP
using LightGraphs

import Base: copy
import Kombinator:
    solve,
    dimension,
    approximation_ratio,
    approximation_term,
    reward,
    value,
    formulation,
    fastest_exact,
    make_solution

export SpanningTreeInstance,
    SpanningTreeSolution,
    SimpleBudgetedSpanningTreeSolution,
    BudgetedSpanningTreeLagrangianSolution,
    BudgetedSpanningTreeDynamicProgrammingSolution

include("st.jl")

include("helpers.jl")
include("helpers_lp.jl")

include("st_prim.jl")
include("st_dp.jl")
include("st_lp.jl")

include("st_budgeted_lagrangian.jl")
include("st_budgeted_lagrangian_refined.jl")
include("st_budgeted_lagrangian_refined_iterated.jl")
include("st_budgeted_dp.jl")
include("st_budgeted_lp.jl")

end


text:
module TestDistSystem
using Test
using DataFrames, CSV

include("../src/DistSystem.jl")
include("case33_ybus.jl")
import .DistSystem
export runtests

function test_null_der()
    buses = ["a", "2"]
    null_der = DistSystem.null_der(buses)

    @test null_der isa DistSystem.DER
    @test null_der.bus == buses[2]
    @test null_der.S_limit == 0.0
    @test null_der.energy == 0.0
    @test null_der.alpha == 0.0
    @test null_der.beta == 1.0
    @test null_der.P_limit == [0.0, 0.0]
    @test null_der.Q_limit == [0.0, 0.0]
    @test null_der.scenario == [0.0]
    @test null_der.Cost == [0.0, 0.0, 0.0]

end
function test_factory_system_case3_dist()
    VN = 12500
    Sb = 1e6
    Yb = Sb / (VN^2)
    Vsub = 1.0
    P_limt = 3
    Q_limit = 3
    bus_sub = 1
    nbuses = 3

    VL = 0.95
    VH = 1.05

    data = DataFrame(CSV.File("case3_dist.csv"))

    Y_target = [
        20-40im -20+40im 0+0im
        -20+40im 48.5714-97.1429im -28.5714+57.1429im
        0+0im -28.5714+57.1429im 28.5714-56.8429im
    ]

    sub = DistSystem.Substation(VN, bus_sub, Vsub, P_limt, Q_limit, [1.0, 2, 3])
    Bsh = (-(data.Bshunt_MVAr * 1e6) ./ (sub.nominal_voltage^2)) / Yb

    amps = collect(skipmissing(data.Amp_pu))
    fb = collect(skipmissing(data.FB))
    tb = collect(skipmissing(data.TB))

    amp = Dict((fb[k], tb[k]) => amps[k] for k = 1:length(amps))

    sys = DistSystem.factory_system(data, VL, VH, sub)

    @testset "case3_dist" begin

        for i = 1:nbuses, j = 1:nbuses
            @test Y_target[i, j] ≈ sys.Y[i, j] atol = 1e-4
        end
        @test nbuses == sys.nbuses
        @test VL == sys.VL
        @test VH == sys.VH
        @test collect(skipmissing(data.P_MW)) == sys.PL
        @test collect(skipmissing(data.Q_MVAr)) == sys.QL
        @test [1.0] == sys.m_load
        @test [[0.0]] == sys.m_new_dg
        @test sub == sys.substation
        @test amp == sys.amp
        @test collect(skipmissing(data.Bus)) == sys.buses
        @test -collect(skipmissing(data.Bshunt_MVAr)) == sys.Bsh

        # TODO: Find a way to test it better.
        null_der = DistSystem.null_der(sys.buses)
        fields = fieldnames(DistSystem.DER)
        for field in fields
            @test getfield(null_der, field) == getfield(sys.dgs[1], field)
        end
        @test_throws BoundsError sys.dgs[2]
    end
end

function test_factory_system_case33_dist()
    VN = 12660
    Sb = 1e6
    Yb = Sb / (VN^2)
    Vsub = 1.0
    P_limt = 3
    Q_limit = 3
    bus_sub = 1
    nbuses = 33

    VL = 0.95
    VH = 1.05

    data = DataFrame(CSV.File("case33.csv"))

    Y_target = get_ybus()

    sub = DistSystem.Substation(VN, bus_sub, Vsub, P_limt, Q_limit, [1.0, 2, 3])
    Bsh = (-(data.Bshunt_MVAr * 1e6) ./ (sub.nominal_voltage^2)) / Yb

    amps = collect(skipmissing(data.Amp_pu))
    fb = collect(skipmissing(data.FB))
    tb = collect(skipmissing(data.TB))

    amp = Dict((fb[k], tb[k]) => amps[k] for k = 1:length(amps))

    sys = DistSystem.factory_system(data, VL, VH, sub)

    @testset "case33" begin

        for i = 1:nbuses, j = 1:nbuses
            @test Y_target[i, j] ≈ sys.Y[i, j] atol = 1e-3
        end
        @test nbuses == sys.nbuses
        @test VL == sys.VL
        @test VH == sys.VH
        @test collect(skipmissing(data.P_MW)) == sys.PL
        @test collect(skipmissing(data.Q_MVAr)) == sys.QL
        @test [1.0] == sys.m_load
        @test [[0.0]] == sys.m_new_dg
        @test 0.0 == sys.time_curr
        @test sub == sys.substation
        @test amp == sys.amp
        @test collect(skipmissing(data.Bus)) == sys.buses
        @test -collect(skipmissing(data.Bshunt_MVAr)) == sys.Bsh

        # TODO: Find a way to test it better.
        null_der = DistSystem.null_der(sys.buses)
        fields = fieldnames(DistSystem.DER)
        for field in fields
            @test getfield(null_der, field) == getfield(sys.dgs[1], field)
        end
        @test_throws BoundsError sys.dgs[2]
    end
end

function runtests()
    @testset "DistSystem" begin
        test_null_der()
        test_factory_system_case3_dist()
        test_factory_system_case33_dist()
    end
end

end

TestDistSystem.runtests()


text:
macro getoption(name,default)
    quote
        if isdefined(Main,$(QuoteNode(name)))
            const $(esc(name)) = Main.$name
        else
            const $(esc(name)) = $default
        end
    end
end

@osx_only begin
    @getoption MATHEMATICA_HOME "/Applications/Mathematica.app"
    @getoption MATHEMATICA_EXEC joinpath(MATHEMATICA_HOME,"Contents/MacOS/MathKernel")
    @getoption MATHLINK_LIB joinpath(MATHEMATICA_HOME,"SystemFiles/Links/MathLink/DeveloperKit/MacOSX-x86-64/CompilerAdditions/mathlink.framework/mathlink")
end

@linux_only begin
    @getoption MATHEMATICA_HOME "/usr/local/Wolfram/Mathematica"
    @getoption MATHEMATICA_EXEC "math"
    @getoption MATHLINK_LIB if WORD_SIZE == 32
        joinpath(MATHEMATICA_HOME,"SystemFiles/Links/MathLink/DeveloperKit/linux/CompilerAdditions/libML32i3.so")
    else
        joinpath(MATHEMATICA_HOME,"SystemFiles/Links/MathLink/DeveloperKit/linux-x86-64/CompilerAdditions/libML64i3.so")
    end
end

l = dlopen_e(MATHLINK_LIB)
if l != C_NULL
    dlclose(l)
else
    error("Could not find MathLink library: set either MATHEMATICA_HOME or MATHLINK_LIB before loading package.")
end