text:
export NystromApprox, RandomizedNystromPreconditioner

struct NystromApprox{T}
    U::Matrix{T}
    Λ::Diagonal{T, Vector{T}}
end

# Constructs Â_nys in factored form
# Â_nys = (AΩ)(ΩᵀAΩ)^†(AΩ)^ᵀ = UΛUᵀ
# [Martinsson & Tropp, Algorithm 16]
function NystromApprox(A::Matrix{T}, k::Int, r::Int; check=false) where {T <: Real}
    check && check_psd(A)
    n = size(A, 1)

    Ω = randn(n, r)
    Ω .= Array(qr(Ω).Q)
    Y = A * Ω
    ν = sqrt(n)*eps(norm(Y))
    @. Y += ν * Ω
    B = Y / cholesky(Symmetric(Ω' * Y)).U
    U, Σ, _ = svd(B)
    Λ = Diagonal(max.(0, Σ.^2 .- ν)[1:k])

    return NystromApprox(U[:, 1:k], Λ)
end

LinearAlgebra.eigvals(Anys::NystromApprox) = Anys.Λ.diag

function LinearAlgebra.mul!(y, Anys::NystromApprox, x; cache=zeros(size(Anys.U, 2)))
    length(y) != length(x) || length(y) != size(Anys.U, 1) && error(DimensionMismatch())
    r = size(Anys.U, 2)
    @views mul!(cache[1:r], Anys.U', x)
    @views cache[1:r] .*= Anys.Λ.diag
    @views mul!(y, Anys.U, cache[1:r])
    return nothing
end

function LinearAlgebra.:*(Anys::NystromApprox, x::AbstractVector)
    y = Anys.U*x
    y .*= Anys.Λ
    return Anys.U'*y
end

# Doubles rank until the approximation is sufficiently good
function adaptive_nystrom_approx(A::Matrix{T}, r0::Int; tol=1e-6, check=false, q=10) where {T <: Real}
    check && check_psd(A)
    n = size(A, 1)
    cache = (
        v0=zeros(n),
        v=zeros(n),
        Anys_mul=zeros(n)
    )
    r = r0
    Enorm = Inf
    Anys = nothing
    while Enorm > tol && r < n
        k = Int(round(.9*r))
        Anys = NystromApprox(A, k, r; check=false)
        Enorm = estimate_norm_E(A, Anys; q=q, cache=cache)
        r = 2r
    end
    return Anys
end

# By prop 5.3, have that κ(P^{-1/2} * A * P^{-1/2}) ≤ (λᵣ + μ + ||E||)/μ
function estimate_norm_E(A, Anys; q=10, cache=nothing)
    n = size(A, 1)
    if !isnothing(cache)
        v0, v = cache.v0, cache.v
    else
        v0, v = zeros(n), zeros(n)
        cache = (Anys_mul=zeros(n),)
    end

    v0 .= randn(n)
    normalize!(v0)

    Ehat = Inf
    for _ in 1:q
        mul!(v, Anys, v0; cache=cache.Anys_mul)
        mul!(v, A, v0, 1.0, -1.0)
        Ehat = dot(v0, v)
        normalize!(v)
        v0 .= v
    end
    return Ehat
end

# theoretically, choose r = 2⌈1.5deff(µ)⌉ + 1.
function deff(A::AbstractMatrix, μ; check=false)
    check && check_psd(A)
    λ = eigvals(A)
    return sum(x->x/(x+μ), λ)
end

function check_psd(A)
    n = size(A, 1)
    psd_tol = sqrt(n)*eps(norm(A))
    !isposdef(A + psd_tol*I) && error(ArgumentError("A must be PSD"))
end

struct RandomizedNystromPreconditioner{T <: Real}
    A_nys::NystromApprox{T}
    λ::MVector{1, T}
    μ::MVector{1, T}
    cache::Vector{T}
    function RandomizedNystromPreconditioner(A_nys::NystromApprox{T}, μ::T) where {T <: Real}
        return new{T}(A_nys, SA[A_nys.Λ.diag[end]], SA[μ], zeros(size(A_nys.U, 2)))
    end
end
Base.eltype(P::RandomizedNystromPreconditioner{T}) where {T} = T

function LinearAlgebra.:\(P::RandomizedNystromPreconditioner{T}, x::Vector{T}) where {T <: Real}
    λ = P.λ[1]
    μ = P.μ[1]
    # P⁻x = x - U*((λ + μ)*(Λ + μI)^-1 + I)*Uᵀ*x
    mul!(P.A_nys.cache, P.A_nys.U', x)
    @. P.A_nys.cache *= (P.λ + P.μ) * 1 / (Λ.diag + μ) + 1
    return x - U*P.A_nys.cache
end

function LinearAlgebra.ldiv!(y::Vector{T}, P::RandomizedNystromPreconditioner{T}, x::Vector{T}) where {T <: Real}
    λ = P.λ[1]
    μ = P.μ[1]
    length(y) != length(x) && error(DimensionMismatch())
    mul!(P.cache, P.A_nys.U', x)
    @. P.cache *= (P.λ + P.μ) * 1 / (P.A_nys.Λ.diag + P.μ) - 1
    mul!(y, P.A_nys.U, P.cache)
    @. y = x + y
    return nothing
end

function LinearAlgebra.ldiv!(P::RandomizedNystromPreconditioner{T}, x::Vector{T}) where {T <: Real}
    λ = P.λ[1]
    μ = P.μ[1]
    mul!(P.cache, P.A_nys.U', x)
    @. P.cache *= (P.λ + P.μ) * 1 / (P.A_nys.Λ.diag + P.μ) + 1
    x .-= P.A_nys.U*P.cache
    return nothing
end

# Used for Krylov method
struct RandomizedNystromPreconditionerInverse{T <: Real}
    A_nys::NystromApprox{T}
    λ::MVector{1, T}
    μ::MVector{1, T}
    cache::Vector{T}
    function RandomizedNystromPreconditionerInverse(A_nys::NystromApprox{T}, μ::T) where {T <: Real}
        return new{T}(A_nys, SA[A_nys.Λ.diag[end]], SA[μ], zeros(size(A_nys.U, 2)))
    end
end
Base.eltype(P::RandomizedNystromPreconditionerInverse{T}) where {T} = T

function LinearAlgebra.mul!(y, P::RandomizedNystromPreconditionerInverse, x)
    length(y) != length(x) && error(DimensionMismatch())

    λ = P.λ[1]
    μ = P.μ[1]

    mul!(P.cache, P.A_nys.U', x)
    @. P.cache *= (P.λ + P.μ) / (P.A_nys.Λ.diag + P.μ) - 1
    mul!(y, P.A_nys.U, P.cache)
    @. y = x + y
end


text:
"""
    get_surface_velocities!(Vcp, Vh, Vv, Vte, current_surface, previous_surface, dt)

Calculate the velocities experienced by the surface at the control points, horizontal
bound vortex centers, vertical bound vortex centers and trailing edge vertices due
to surface motion.
"""
function get_surface_velocities!(Vcp, Vh, Vv, Vte, current_surface, previous_surface, dt)

    nc, ns = size(current_surface)

    # velocity at the control points
    for j = 1:ns, i = 1:nc
        Vcp[i,j] = (controlpoint(previous_surface[i,j]) - controlpoint(current_surface[i,j]))/dt
    end

    # velocity at the horizontal bound vortices
    for j = 1:ns
        for i = 1:nc
            Vh[i,j] = (top_center(previous_surface[i,j]) - top_center(current_surface[i,j]))/dt
        end
        Vh[end,j] = (bottom_center(previous_surface[end,j]) - bottom_center(current_surface[end,j]))/dt
    end

    # velocity at the vertical bound vortices
    for i = 1:nc
        for j = 1:ns
            Vv[i,j] = (left_center(previous_surface[i,j]) - left_center(current_surface[i,j]))/dt
        end
        Vv[i,end] = (right_center(previous_surface[i,end]) - right_center(current_surface[i,end]))/dt
    end

    # velocity at the trailing edge vertices
    for j = 1:ns
        Vte[j] = (bottom_left(previous_surface[end,j]) - bottom_left(current_surface[end,j]))/dt
    end
    Vte[end] = (bottom_right(previous_surface[end,end]) - bottom_right(current_surface[end,end]))/dt

end

"""
    normal_velocity!(w, surfaces, wakes, ref, fs; additional_velocity,
        Vcp, symmetric, nwake, surface_id, wake_finite_core, trailing_vortices, xhat)

Compute the downwash at the control points on `surfaces` due to the freestream
velocity, rotational velocity, additional velocity field, surface motion, and
induced velocity from the wake panels.

This forms the right hand side of the circulation linear system solve.
"""
@inline function normal_velocity!(w, surfaces, wakes, ref, fs; additional_velocity,
    Vcp, symmetric, nwake, surface_id, wake_finite_core, trailing_vortices, xhat)

    nsurf = length(surfaces)

    # index for keeping track of where we are in the w vector
    iw = 0

    # loop through receiving surfaces
    for isurf = 1:nsurf

        # current surface
        surface = surfaces[isurf]

        # iterate through all control points on this surface
        for i = 1:length(surface)

            # control point
            rcp = controlpoint(surface[i])

            # normal vector
            nhat = normal(surface[i])

            # freestream velocity
            V = freestream_velocity(fs)

            # rotational velocity
            V += rotational_velocity(rcp, fs, ref)

            # additional velocity field
            if !isnothing(additional_velocity)
                V += additional_velocity(rcp)
            end

            # velocity due to surface motion
            if !isnothing(Vcp)
                V += Vcp[isurf][i]
            end

            # velocity due to the wake panels
            for jsurf = 1:length(wakes)
                if nwake[jsurf] > 0
                    V += induced_velocity(rcp, wakes[jsurf];
                        finite_core = wake_finite_core[isurf] || (surface_id[isurf] != surface_id[jsurf]),
                        symmetric = symmetric[jsurf],
                        nc = nwake[jsurf],
                        trailing_vortices = trailing_vortices[jsurf],
                        xhat = xhat)
                end
            end

            # get normal component
            w[iw+i] = -dot(V, nhat)
        end

        # increment position in AIC matrix
        iw += length(surface)
    end

    return w
end

"""
    normal_velocity_derivatives!(w, dw, surfaces, wakes, ref, fs;
        additional_velocity, Vcp, symmetric, nwake, surface_id,
        wake_finite_core, trailing_vortices, xhat)

Compute the downwash at the control points on `surfaces` due to the freestream
velocity, rotational velocity, additional velocity field, surface motion, and
induced velocity from the wake panels.  Also calculate its derivatives with
respect to the freestream parameters.

This forms the right hand side of the circulation linear system solve (and its derivatives).
"""
@inline function normal_velocity_derivatives!(w, dw, surfaces, wakes, ref, fs;
    additional_velocity, Vcp, symmetric, nwake, surface_id, wake_finite_core,
    trailing_vortices, xhat)

    nsurf = length(surfaces)

    # unpack derivatives
    (w_a, w_b, w_p, w_q, w_r) = dw

    # index for keeping track of where we are in the w vector
    iw = 0

    # loop through receiving surfaces
    for isurf = 1:nsurf

        # current surface
        surface = surfaces[isurf]

        # iterate through all control points on this surface
        for i = 1:length(surface)

            # control point
            rcp = controlpoint(surface[i])

            # normal vector
            nhat = normal(surface[i])

            # freestream velocity and its derivatives
            V, dVf = freestream_velocity_derivatives(fs)

            V_a, V_b = dVf

            # rotational velocity and its derivatives
            Vrot, dVrot = rotational_velocity_derivatives(rcp, fs, ref)

            V += Vrot

            V_p, V_q, V_r = dVrot

            # additional velocity field
            if !isnothing(additional_velocity)
                V += additional_velocity(rcp)
            end

            # velocity due to surface motion
            if !isnothing(Vcp)
                V += Vcp[isurf][i]
            end

            # velocity due to the wake panels
            for jsurf = 1:length(wakes)
                if nwake[jsurf] > 0
                    V += induced_velocity(rcp, wakes[jsurf];
                        finite_core = wake_finite_core[jsurf] || (surface_id[isurf] != surface_id[jsurf]),
                        symmetric = symmetric[jsurf],
                        nc = nwake[jsurf],
                        trailing_vortices = trailing_vortices[jsurf],
                        xhat = xhat)
                end
            end

            # right hand side vector
            w[iw+i] = -dot(V, nhat)

            # associated derivatives
            w_a[iw+i] = -dot(V_a, nhat)
            w_b[iw+i] = -dot(V_b, nhat)
            w_p[iw+i] = -dot(V_p, nhat)
            w_q[iw+i] = -dot(V_q, nhat)
            w_r[iw+i] = -dot(V_r, nhat)

        end

        # increment position in AIC matrix
        iw += length(surface)
    end

    # pack up derivatives
    dw = (w_a, w_b, w_p, w_q, w_r)

    return w, dw
end

# --- circulation solve --- #

"""
    circulation(AIC, w)

Solve for the circulation distribution.
"""
circulation(AIC, w) = AIC\w

"""
    circulation!(Γ, AIC, w)

Pre-allocated version of `circulation`
"""
circulation!(Γ, AIC, w) = ldiv!(Γ, lu(AIC), w)

"""
    circulation_derivatives(AIC, w, dw)

Solve for the circulation distribution and its derivatives with respect to
the freestream parameters.
"""
function circulation_derivatives(AIC, w, dw)

    # unpack derivatives
    (w_a, w_b, w_p, w_q, w_r) = dw

    # factorize AIC matrix (since we'll be reusing it)
    fAIC = factorize(AIC)

    # solve for circulation and its derivatives
    Γ = fAIC\w

    Γ_a = fAIC\w_a
    Γ_b = fAIC\w_b
    Γ_pb = fAIC\w_p
    Γ_qb = fAIC\w_q
    Γ_rb = fAIC\w_r

    # pack up derivatives
    dΓ = (Γ_a, Γ_b, Γ_pb, Γ_qb, Γ_rb)

    return Γ, dΓ
end

"""
    circulation_derivatives!(Γ, dΓ, AIC, w, dw)

Pre-allocated version of `circulation_derivatives`
"""
function circulation_derivatives!(Γ, dΓ, AIC, w, dw)

    # unpack derivatives
    (Γ_a, Γ_b, Γ_p, Γ_q, Γ_r) = dΓ
    (w_a, w_b, w_p, w_q, w_r) = dw

    # factorize AIC matrix (since we'll be reusing it)
    fAIC = lu(AIC)

    # solve for circulation and its derivatives
    ldiv!(Γ, fAIC, w)

    ldiv!(Γ_a, fAIC, w_a)
    ldiv!(Γ_b, fAIC, w_b)
    ldiv!(Γ_p, fAIC, w_p)
    ldiv!(Γ_q, fAIC, w_q)
    ldiv!(Γ_r, fAIC, w_r)

    return Γ, dΓ
end

ext: jl


text:
"""
	BafflingLaneChanger
A ramdom lane change behavior that changes lanes:
 - whenever the lead car is going slower than our desired speed.
 - whenever randomly decided (based on a uniform distribution).

[For collision avoidance] Lane changes are made when:
 - there is an available lane and
 - fore/rear gaps exceed our thresholds.
[For practicality] However, lane changes can be also made when:
 - we are slower than a rear vehicle in the target lane,
 - any lead vehicle in the target lane is slower than we can currently go.

# Constructors
	BafflingLaneChanger(timestep::Float64,
                    v_des::Float64=29.0,
                    rec::SceneRecord=SceneRecord(2,timestep),
                    threshold_fore::Float64 = 50.0,
                    threshold_lane_change_gap_fore::Float64 = 10.0,
                    threshold_lane_change_gap_rear::Float64 = 10.0,
                    dir::Int=DIR_MIDDLE)

# Fields
- `dir::Int = DIR_MIDDLE` the desired lane to go to eg: left,middle (i.e. stay in same lane) or right
- `rec::SceneRecord` TODO
- `v_des::Float64 = 29.0` desired velocity
- `threshold_fore::Float64 = 50.0` Distance from lead vehicle
- `threshold_lane_change_gap_fore::Float64 = 10.0` Space in frontLaneChangeChoice
- `threshold_lane_change_gap_rear::Float64 = 10.0` Space rear
"""
mutable struct BafflingLaneChanger <: LaneChangeModel
    dir::Int
    rec::SceneRecord

    v_des::Float64
    threshold_fore::Float64
    threshold_lane_change_gap_fore::Float64
    threshold_lane_change_gap_rear::Float64
    threshold_lane_change_rand::Float64

    function BafflingLaneChanger(
        timestep::Float64;
        v_des::Float64=29.0,
        rec::SceneRecord=SceneRecord(2,timestep),
        threshold_fore::Float64 = 50.0,
        threshold_lane_change_gap_fore::Float64 = 5.0,
        threshold_lane_change_gap_rear::Float64 = 5.0,
        threshold_lane_change_rand::Float64 = 0.00,
        dir::Int=DIR_MIDDLE,
        )

        retval = new()
        retval.dir = dir
        retval.rec = rec
        retval.v_des = v_des
        retval.threshold_fore = threshold_fore
        retval.threshold_lane_change_gap_fore = threshold_lane_change_gap_fore
        retval.threshold_lane_change_rand = threshold_lane_change_rand
        retval.threshold_lane_change_gap_rear = threshold_lane_change_gap_rear
        retval
    end
end
get_name(::BafflingLaneChanger) = "BafflingLaneChanger"
function AutomotiveDrivingModels.set_desired_speed!(model::BafflingLaneChanger, v_des::Float64)
    model.v_des = v_des
    model
end
function AutomotiveDrivingModels.observe!(model::BafflingLaneChanger, scene::Scene, roadway::Roadway, egoid::Int)

    rec = model.rec
    update!(rec, scene)
    vehicle_index = findfirst(egoid, scene)

    veh_ego = scene[vehicle_index]
    v = veh_ego.state.v

    left_lane_exists = convert(Float64, get(N_LANE_LEFT, rec, roadway, vehicle_index)) > 0
    right_lane_exists = convert(Float64, get(N_LANE_RIGHT, rec, roadway, vehicle_index)) > 0
    fore_M = get_neighbor_fore_along_lane(scene, vehicle_index, roadway, VehicleTargetPointFront(), VehicleTargetPointRear(), VehicleTargetPointFront(), max_distance_fore=model.threshold_fore)
    fore_L = get_neighbor_fore_along_left_lane(scene, vehicle_index, roadway, VehicleTargetPointRear(), VehicleTargetPointRear(), VehicleTargetPointFront())
    fore_R = get_neighbor_fore_along_right_lane(scene, vehicle_index, roadway, VehicleTargetPointRear(), VehicleTargetPointRear(), VehicleTargetPointFront())
    rear_L = get_neighbor_rear_along_left_lane(scene, vehicle_index, roadway, VehicleTargetPointFront(), VehicleTargetPointFront(), VehicleTargetPointRear())
    rear_R = get_neighbor_rear_along_right_lane(scene, vehicle_index, roadway, VehicleTargetPointFront(), VehicleTargetPointFront(), VehicleTargetPointRear())

    model.dir = DIR_MIDDLE
    r = Base.rand(Float64)
    if fore_M.Δs < model.threshold_fore # there is a lead vehicle
        veh_M = scene[fore_M.ind]
        speed_M = veh_M.state.v
        if speed_M ≤ min(model.v_des, v) # they are driving slower than we want

            speed_ahead = speed_M

            # consider changing to a different lane
            if right_lane_exists &&
               fore_R.Δs > model.threshold_lane_change_gap_rear && # there is space rear
               rear_R.Δs > model.threshold_lane_change_gap_fore && # there is space fore
               (rear_R.ind == nothing || scene[rear_R.ind].state.v ≤ v) && # we are faster than any follower
               (fore_R.ind == nothing || scene[fore_R.ind].state.v > speed_ahead) # lead is faster than current speed

                speed_ahead = fore_R.ind != nothing ? scene[fore_R.ind].state.v : Inf
                model.dir = DIR_RIGHT
            end
            if left_lane_exists &&
               fore_L.Δs > model.threshold_lane_change_gap_rear && # there is space rear
               rear_L.Δs > model.threshold_lane_change_gap_fore && # there is space fore
               (rear_L.ind == nothing || scene[rear_L.ind].state.v ≤ v) && # we are faster than any follower
               (fore_L.ind == nothing || scene[fore_L.ind].state.v > speed_ahead) # lead is faster than current speed

                speed_ahead = fore_L.ind != nothing ? scene[fore_L.ind].state.v : Inf
                model.dir = DIR_LEFT
            end
        end
    elseif r <= model.threshold_lane_change_rand
        # consider changing to a different lane
        if right_lane_exists &&
           fore_R.Δs > model.threshold_lane_change_gap_rear && # there is space rear
           rear_R.Δs > model.threshold_lane_change_gap_fore # there is space fore
            model.dir = DIR_RIGHT
        end
        if left_lane_exists &&
           fore_L.Δs > model.threshold_lane_change_gap_rear && # there is space rear
           rear_L.Δs > model.threshold_lane_change_gap_fore # there is space fore
            model.dir = DIR_LEFT
        end
    end

    model
end

Base.rand(model::BafflingLaneChanger) = LaneChangeChoice(model.dir)


text:
# NegativeBinomial is the distribution of the number of failures
# before the r-th success in a sequence of Bernoulli trials.
# We do not enforce integer size, as the distribution is well defined
# for non-integers, and this can be useful for e.g. overdispersed
# discrete survival times.

immutable NegativeBinomial <: DiscreteUnivariateDistribution
    r::Float64
    prob::Float64

    function NegativeBinomial(r::Real, p::Real)
        zero(p) < p <= one(p) || error("prob must be in (0, 1].")
        zero(r) < r || error("r must be positive.")
        new(float64(r), float64(p))
    end

    NegativeBinomial() = new(1.0, 0.5)
end

@_jl_dist_2p NegativeBinomial nbinom

insupport(::NegativeBinomial, x::Real) = isinteger(x) && zero(x) <= x
insupport(::Type{NegativeBinomial}, x::Real) = isinteger(x) && zero(x) <= x

function mgf(d::NegativeBinomial, t::Real)
    r, p = d.size, d.prob
    return ((1.0 - p) * exp(t))^r / (1.0 - p * exp(t))^r
end

function cf(d::NegativeBinomial, t::Real)
    r, p = d.size, d.prob
    return ((1.0 - p) * exp(im * t))^r / (1.0 - p * exp(im * t))^r
end

function mean(d::NegativeBinomial)
    p = d.prob
    (1.0 - p) * d.r / p
end

function var(d::NegativeBinomial)
    p = d.prob
    (1.0 - p) * d.r / (p * p)
end

function std(d::NegativeBinomial)
    p = d.prob
    sqrt((1.0 - p) * d.r) / p
end

function skewness(d::NegativeBinomial)
    p = d.prob
    (2.0 - p) / sqrt((1.0 - p) * r)
end

function kurtosis(d::NegativeBinomial)
    p = d.prob
    6.0 / r + (p * p) / ((1.0 - p) * r)
end

function mode(d::NegativeBinomial)
    p = d.prob
    ifloor((1.0 - p) * (d.r - 1.) / p)
end

modes(d::NegativeBinomial) = [mode(d)]


text:
function solve_polymerflood(Si::T, Sj::T, kr1::RelPerms, kr2::RelPerms, D::T, μw::T, μo::T) where T <: Float64

    f1(s) = fractional_flow(s, kr1, μw, μo)
    df1(s) = fw_derivative(s, k1, μw, μo)
    f2(s) = fractional_flow(s, kr2, μw, μo)
    df2(s) = fw_derivative(s, kr2, μw, μo)

    # Velocity of polymer shock
    vΔc(s) = f2(s) ./ (s + D)

    # * Step 1 - Calculate Ŝ1 and vΔc
    Ŝ = fzero(s -> vΔc.(s) - df2(s), 0.6)
    VΔc = vΔc(Ŝ)

    # * Step2 - Calculate S1B and v2b
    Sb = fzero(s -> (s + D) * vΔc(Ŝ) - f1(s), 0.4)
    Vb = (f1(Sb) - f1(Si)) / (Sb - Si)

    PolymerFlooding(Ŝ, Si, Sj, kr1, kr2, μw, μo, Sb, D, Vb, VΔc)

end

function saturation_profile(pf::PolymerFlooding, time::Real)
    si = pf.Si
    sj = pf.Sj
    sb = pf.Sb
    ŝ = pf.S̃

    df(s) = fw_derivative(s, pf.kr2, pf.μw, pf.μo)

    s = vcat([si, si, sb ,sb, ŝ], collect(ŝ:0.01:sj))
    x = zeros(length(s))
    x[1] = 1.0
    x[2:3] .= pf.Vb * time
    x[4:5] .= pf.VΔc * time
    x[6:end] = df(s[6:end]) * time

    return x, s
end

function polymer_profile(pf::PolymerFlooding, time::Real)
    sj = pf.Sj
    ŝ = pf.S̃

    df(s) = fw_derivative(s, pf.kr2, pf.μw, pf.μo)

    s = collect(ŝ:0.01:sj)
    x = zeros(length(s))
    x = df(s) * time

    return x, s
end


text:
module TestSamples

using GeneratorsX
using Test
using Transducers: Map

@generator noone() = nothing
@generator oneone() = @yield 1
@generator function onetwothree()
    @yield 1
    @yield 2
    @yield 3
end

@generator function iid(xs)
    for x in xs
        @yield x
    end
end

@generator function iflatten(xs)
    for y in xs
        for x in y
            @yield x
        end
    end
end

struct Count
    start::Int
    stop::Int
end

Base.length(itr::Count) = max(0, itr.stop - itr.start + 1)
Base.eltype(::Type{<:Count}) = Int

@generator function (itr::Count)
    i = itr.start
    i > itr.stop && return
    while true
        @yield i
        i == itr.stop && break
        i += 1
    end
end

raw_testdata = """
noone() == []
oneone() == [1]
onetwothree() == [1, 2, 3]
iid(1:1) == [1]
iid(noone()) == []
iid(oneone()) == [1]
iid(onetwothree()) == [1, 2, 3]
iflatten([[1], [2, 3], [4]]) == [1, 2, 3, 4]
iflatten([[1], (2, 3), 4]) == [1, 2, 3, 4]
Count(0, 2) == [0, 1, 2]
Count(0, -1) == Int[]
"""

args_and_kwargs(args...; kwargs...) = args, (; kwargs...)

# An array of `(label, (f, args, kwargs, comparison, desired))`
testdata = map(split(raw_testdata, "\n", keepempty = false)) do x
    comp_ex = Meta.parse(x)
    @assert comp_ex.head == :call
    @assert length(comp_ex.args) == 3
    comparison, ex, desired = comp_ex.args
    f = ex.args[1]
    ex.args[1] = args_and_kwargs

    label = strip(x[1:prevind(x, first(findlast(String(comparison), x)))])
    Meta.parse(label)  # validation

    @eval ($label, ($(Symbol(f)), $ex..., $comparison, $desired))
end

@testset "$label" for (label, (f, args, kwargs, comparison, desired)) in testdata
    ==′ = comparison
    @test collect(f(args...; kwargs...)) ==′ desired
    @test collect(Map(identity), f(args...; kwargs...)) ==′ desired
end

end  # module


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Random
using Base: Experimental

@testset "single-threaded Condition usage" begin
    a = Condition()
    t = @async begin
        Base.notify(a, "success")
        "finished"
    end
    @test wait(a) == "success"
    @test fetch(t) == "finished"
end

@testset "various constructors" begin
    c = Channel()
    @test eltype(c) == Any
    @test c.sz_max == 0

    c = Channel(1)
    @test eltype(c) == Any
    @test put!(c, 1) == 1
    @test isready(c) == true
    @test take!(c) == 1
    @test isready(c) == false
    @test eltype(Channel(1.0)) == Any

    c = Channel{Int}(1)
    @test eltype(c) == Int
    @test_throws MethodError put!(c, "Hello")

    c = Channel{Int}()
    @test eltype(c) == Int
    @test c.sz_max == 0

    c = Channel{Int}(Inf)
    @test eltype(c) == Int
    pvals = map(i->put!(c,i), 1:10^6)
    tvals = Int[take!(c) for i in 1:10^6]
    @test pvals == tvals

    @test_throws ArgumentError Channel(-1)
    @test_throws InexactError Channel(1.5)
end

@testset "Task constructors" begin
    c = Channel{Int}() do c; map(i->put!(c,i), 1:100); end
    @test eltype(c) == Int
    @test c.sz_max == 0
    @test collect(c) == 1:100

    c = Channel() do c; put!(c, 1); put!(c, "hi") end
    @test c.sz_max == 0
    @test collect(c) == [1, "hi"]

    c = Channel(Inf) do c; put!(c,1); end
    @test eltype(c) == Any
    @test c.sz_max == typemax(Int)
    c = Channel{Int}(Inf) do c; put!(c,1); end
    @test eltype(c) == Int
    @test c.sz_max == typemax(Int)

    taskref = Ref{Task}()
    c = Channel{Int}(0, taskref=taskref) do c; put!(c, 0); end
    @test eltype(c) == Int
    @test c.sz_max == 0
    @test istaskstarted(taskref[])
    @test !istaskdone(taskref[])
    take!(c); wait(taskref[])
    @test istaskdone(taskref[])

    # Legacy constructor
    c = Channel(ctype=Float32, csize=2) do c; map(i->put!(c,i), 1:100); end
    @test eltype(c) == Float32
    @test c.sz_max == 2
    @test isopen(c)
    @test collect(c) == 1:100
end
@testset "Multithreaded task constructors" begin
    taskref = Ref{Task}()
    c = Channel(spawn=true, taskref=taskref) do c; put!(c, 0); end
    # Test that the task is using the multithreaded scheduler
    @test taskref[].sticky == false
    @test collect(c) == [0]
end

@testset "multiple concurrent put!/take! on a channel for different sizes" begin
    function testcpt(sz)
        c = Channel{Int}(sz)
        size = 0
        inc() = size += 1
        dec() = size -= 1
        @sync for i = 1:10^4
            @async (sleep(rand()); put!(c, i); inc())
            @async (sleep(rand()); take!(c); dec())
        end
        @test size == 0
    end
    testcpt(0)
    testcpt(1)
    testcpt(32)
    testcpt(Inf)
end

@testset "type conversion in put!" begin
    c = Channel{Int64}(0)
    @async put!(c, Int32(1))
    wait(c)
    @test isa(take!(c), Int64)
    @test_throws MethodError put!(c, "")
    @assert !islocked(c.cond_take)
end

@testset "multiple for loops waiting on the same channel" begin
    # Test multiple "for" loops waiting on the same channel which
    # is closed after adding a few elements.
    c = Channel(32)
    results = []
    @sync begin
        for i in 1:20
            @async for ii in c
                push!(results, ii)
            end
        end
        sleep(1.0)
        for i in 1:5
            put!(c,i)
        end
        close(c)
    end
    @test sum(results) == 15
end

# Tests for channels bound to tasks.
using Distributed
@testset "channels bound to tasks" for N in [0, 10]
    # Normal exit of task
    c = Channel(N)
    bind(c, @async (GC.gc(); yield(); nothing))
    @test_throws InvalidStateException take!(c)
    @test !isopen(c)

    # Error exception in task
    c = Channel(N)
    task = @async (GC.gc(); yield(); error("foo"))
    bind(c, task)
    @test_throws TaskFailedException(task) take!(c)
    @test !isopen(c)

    # Multiple channels closed by the same bound task
    cs = [Channel(N) for i in 1:5]
    tf2() = begin
        GC.gc()
        if N > 0
            foreach(c -> (@assert take!(c) === 2), cs)
        end
        yield()
        error("foo")
    end
    task = Task(tf2)
    foreach(c -> bind(c, task), cs)
    schedule(task)

    if N > 0
        for i in 1:5
            @test put!(cs[i], 2) === 2
        end
    end
    for i in 1:5
        while isopen(cs[i])
            yield()
        end
        @test_throws TaskFailedException(task) wait(cs[i])
        @test_throws TaskFailedException(task) take!(cs[i])
        @test_throws TaskFailedException(task) put!(cs[i], 1)
        N == 0 || @test_throws TaskFailedException(task) fetch(cs[i])
        N == 0 && @test_throws ErrorException fetch(cs[i])
    end

    # Multiple tasks, first one to terminate closes the channel
    nth = rand(1:5)
    ref = Ref(0)
    tf3(i) = begin
        GC.gc()
        if i == nth
            ref[] = i
        else
            sleep(2.0)
        end
    end

    tasks = [Task(() -> tf3(i)) for i in 1:5]
    c = Channel(N)
    foreach(t -> bind(c, t), tasks)
    foreach(schedule, tasks)
    @test_throws InvalidStateException wait(c)
    @test !isopen(c)
    @test ref[] == nth
    @assert !islocked(c.cond_take)

    # channeled_tasks
    for T in [Any, Int]
        tf_chnls1 = (c1, c2) -> (@assert take!(c1) == 1; put!(c2, 2))
        chnls, tasks = Base.channeled_tasks(2, tf_chnls1; ctypes=[T,T], csizes=[N,N])
        put!(chnls[1], 1)
        @test take!(chnls[2]) === 2
        @test_throws InvalidStateException wait(chnls[1])
        @test_throws InvalidStateException wait(chnls[2])
        @test istaskdone(tasks[1])
        @test !isopen(chnls[1])
        @test !isopen(chnls[2])

        f = Future()
        tf4 = (c1, c2) -> begin
            @assert take!(c1) === 1
            wait(f)
        end

        tf5 = (c1, c2) -> begin
            put!(c2, 2)
            wait(f)
        end

        chnls, tasks = Base.channeled_tasks(2, tf4, tf5; ctypes=[T,T], csizes=[N,N])
        put!(chnls[1], 1)
        @test take!(chnls[2]) === 2
        yield()
        put!(f, 1) # allow tf4 and tf5 to exit after now, eventually closing the channel

        @test_throws InvalidStateException wait(chnls[1])
        @test_throws InvalidStateException wait(chnls[2])
        @test istaskdone(tasks[1])
        @test istaskdone(tasks[2])
        @test !isopen(chnls[1])
        @test !isopen(chnls[2])
    end

    # channel
    tf6 = c -> begin
        @assert take!(c) === 2
        error("foo")
    end

    for T in [Any, Int]
        taskref = Ref{Task}()
        chnl = Channel{T}(tf6, N, taskref=taskref)
        put!(chnl, 2)
        yield()
        @test_throws TaskFailedException(taskref[]) wait(chnl)
        @test istaskdone(taskref[])
        @test !isopen(chnl)
        @test_throws TaskFailedException(taskref[]) take!(chnl)
    end
end

@testset "timedwait" begin
    alwaystrue() = true
    alwaysfalse() = false
    @test timedwait(alwaystrue, 0) === :ok
    @test timedwait(alwaysfalse, 0) === :timed_out
    @test_throws ArgumentError timedwait(alwaystrue, 0; pollint=0)

    # Allowing a smaller positive `pollint` results in `timewait` hanging
    @test_throws ArgumentError timedwait(alwaystrue, 0, pollint=1e-4)

    # Callback passed in raises an exception
    failure_cb = function (fail_on_call=1)
        i = 0
        function ()
            i += 1
            i >= fail_on_call && error("callback failed")
            return false
        end
    end

    try
        timedwait(failure_cb(1), 0)
        @test false
    catch e
        @test e isa CapturedException
        @test e.ex isa ErrorException
    end

    try
        timedwait(failure_cb(2), 0)
        @test false
    catch e
        @test e isa CapturedException
        @test e.ex isa ErrorException
    end

    duration = @elapsed timedwait(alwaysfalse, 1)  # Using default pollint of 0.1
    @test duration ≈ 1 atol=0.4

    duration = @elapsed timedwait(alwaysfalse, 0; pollint=1)
    @test duration ≈ 1 atol=0.4
end

@testset "timedwait on multiple channels" begin
    @Experimental.sync begin
        rr1 = Channel(1)
        rr2 = Channel(1)
        rr3 = Channel(1)

        callback() = all(map(isready, [rr1, rr2, rr3]))
        # precompile functions which will be tested for execution time
        @test !callback()
        @test timedwait(callback, 0) === :timed_out

        @async begin sleep(0.5); put!(rr1, :ok) end
        @async begin sleep(1.0); put!(rr2, :ok) end
        @async begin sleep(2.0); put!(rr3, :ok) end

        et = @elapsed timedwait(callback, 1)

        # assuming that 0.5 seconds is a good enough buffer on a typical modern CPU
        try
            @assert (et >= 1.0) && (et <= 1.5)
            @assert !isready(rr3)
        catch
            @warn "`timedwait` tests delayed. et=$et, isready(rr3)=$(isready(rr3))"
        end
        @test isready(rr1)
    end
end

@testset "yield/wait/event failures" begin
    # garbage_finalizer returns `nothing` rather than the garbage object so
    # that the interpreter doesn't accidentally root the garbage when
    # interpreting the calling function.
    @noinline garbage_finalizer(f) = (finalizer(f, "gar" * "bage"); nothing)
    run = Ref(0)
    garbage_finalizer(x -> nothing) # warmup
    @test GC.enable(false)
    # test for finalizers trying to yield leading to failed attempts to context switch
    garbage_finalizer((x) -> (run[] += 1; sleep(1)))
    garbage_finalizer((x) -> (run[] += 1; yield()))
    garbage_finalizer((x) -> (run[] += 1; yieldto(@task () -> ())))
    t = @task begin
        @test !GC.enable(true)
        GC.gc()
        true
    end
    oldstderr = stderr
    newstderr = redirect_stderr()
    local errstream
    try
        errstream = @async read(newstderr[1], String)
        yield(t)
    finally
        redirect_stderr(oldstderr)
        close(newstderr[2])
    end
    @test istaskdone(t)
    @test fetch(t)
    @test run[] == 3
    output = fetch(errstream)
    @test 3 == length(findall(
        """error in running finalizer: ErrorException("task switch not allowed from inside gc finalizer")""", output))
    # test for invalid state in Workqueue during yield
    t = @async nothing
    t._state = 66
    newstderr = redirect_stderr()
    try
        errstream = @async read(newstderr[1], String)
        yield()
    finally
        redirect_stderr(oldstderr)
        close(newstderr[2])
    end
    @test fetch(errstream) == "\nWARNING: Workqueue inconsistency detected: popfirst!(Workqueue).state != :runnable\n"
end

@testset "throwto" begin
    t = @task(nothing)
    ct = current_task()
    testerr = ErrorException("expected")
    @async Base.throwto(t, testerr)
    @test (try
        Base.wait(t)
        false
    catch ex
        ex
    end).task.exception === testerr
end

@testset "Timer / AsyncCondition triggering and race #12719" begin
    let tc = Ref(0)
        t = Timer(0) do t
            tc[] += 1
        end
        cb = first(t.cond.waitq)
        Libc.systemsleep(0.005)
        @test isopen(t)
        Base.process_events()
        @test !isopen(t)
        @test tc[] == 0
        yield()
        @test tc[] == 1
        @test istaskdone(cb)
    end

    let tc = Ref(0)
        t = Timer(0) do t
            tc[] += 1
        end
        cb = first(t.cond.waitq)
        Libc.systemsleep(0.005)
        @test isopen(t)
        close(t)
        @test !isopen(t)
        wait(cb)
        @test tc[] == 0
        @test t.handle === C_NULL
    end

    let tc = Ref(0)
        async = Base.AsyncCondition() do async
            tc[] += 1
        end
        cb = first(async.cond.waitq)
        @test isopen(async)
        ccall(:uv_async_send, Cvoid, (Ptr{Cvoid},), async)
        ccall(:uv_async_send, Cvoid, (Ptr{Cvoid},), async)
        @test isempty(Base.Workqueue)
        Base.process_events() # schedule event
        Sys.iswindows() && Base.process_events() # schedule event (windows?)
        @test length(Base.Workqueue) == 1
        ccall(:uv_async_send, Cvoid, (Ptr{Cvoid},), async)
        @test tc[] == 0
        yield() # consume event
        @test tc[] == 1
        Sys.iswindows() && Base.process_events() # schedule event (windows?)
        yield() # consume event
        @test tc[] == 2
        sleep(0.1) # no further events
        @test tc[] == 2
        ccall(:uv_async_send, Cvoid, (Ptr{Cvoid},), async)
        ccall(:uv_async_send, Cvoid, (Ptr{Cvoid},), async)
        Base.process_events() # schedule event
        Sys.iswindows() && Base.process_events() # schedule event (windows?)
        close(async) # and close
        @test !isopen(async)
        @test tc[] == 2
        @test tc[] == 2
        yield() # consume event & then close
        @test tc[] == 3
        sleep(0.1) # no further events
        wait(cb)
        @test tc[] == 3
        @test async.handle === C_NULL
    end

    let tc = Ref(0)
        async = Base.AsyncCondition() do async
            tc[] += 1
        end
        cb = first(async.cond.waitq)
        @test isopen(async)
        ccall(:uv_async_send, Cvoid, (Ptr{Cvoid},), async)
        Base.process_events() # schedule event
        Sys.iswindows() && Base.process_events() # schedule event (windows)
        close(async)
        @test !isopen(async)
        Base.process_events() # and close
        @test tc[] == 0
        yield() # consume event & then close
        @test tc[] == 1
        sleep(0.1) # no further events
        wait(cb)
        @test tc[] == 1
        @test async.handle === C_NULL
    end
end

@testset "check_channel_state" begin
    c = Channel(1)
    close(c)
    @test !isopen(c)
    c.excp == nothing # to trigger the branch
    @test_throws InvalidStateException Base.check_channel_state(c)
end

# PR #36641
# Ensure that `isempty()` does not mutate a Channel's state:
@testset "isempty(::Channel) mutation" begin
    function isempty_timeout(c::Channel)
        inner_c = Channel{Union{Bool,Nothing}}()
        @async put!(inner_c, isempty(c))
        @async begin
            sleep(0.01)
            if isopen(inner_c)
                put!(inner_c, nothing)
            end
        end
        result = take!(inner_c)
        if result === nothing
            error("isempty() timed out!")
        end
        return result
    end
    # First, with a non-buffered channel
    c = Channel()
    @test isempty_timeout(c)
    t_put = @async put!(c, 1)
    @test !isempty_timeout(c)
    # check a second time to ensure `isempty(c)` didn't just consume the element.
    @test !isempty_timeout(c)
    @test take!(c) == 1
    @test isempty_timeout(c)
    wait(t_put)

    # Next, with a buffered channel:
    c = Channel(2)
    @test isempty_timeout(c)
    t_put = put!(c, 1)
    @test !isempty_timeout(c)
    @test !isempty_timeout(c)
    @test take!(c) == 1
    @test isempty_timeout(c)
end

# issue #12473
# make sure 1-shot timers work
let a = []
    Timer(t -> push!(a, 1), 0.01, interval = 0)
    sleep(0.2)
    @test a == [1]
end

# make sure that we don't accidentally create a one-shot timer
let
    t = Timer(t->nothing, 10, interval=0.00001)
    @test ccall(:uv_timer_get_repeat, UInt64, (Ptr{Cvoid},), t) == 1
    close(t)
end

# make sure repeating timers work
@noinline function make_unrooted_timer(a)
    t = Timer(0.0, interval = 0.1)
    finalizer(t -> a[] += 1, t)
    wait(t)
    e = @elapsed for i = 1:5
        wait(t)
    end
    @test 1.5 > e >= 0.4
    @test a[] == 0
    nothing
end
let a = Ref(0)
    make_unrooted_timer(a)
    GC.gc()
    @test a[] == 1
end

# trying to `schedule` a finished task
let t = @async nothing
    wait(t)
    @test_throws ErrorException("schedule: Task not runnable") schedule(t, nothing)
end

@testset "push!(c, v) -> c" begin
    c = Channel(Inf)
    @test push!(c, nothing) === c
end

# Channel `show`
let c = Channel(3)
    @test repr(c) == "Channel{Any}(3)"
    @test repr(MIME("text/plain"), c) == "Channel{Any}(3) (empty)"
    put!(c, 0)
    @test repr(MIME("text/plain"), c) == "Channel{Any}(3) (1 item available)"
    put!(c, 1)
    @test repr(MIME("text/plain"), c) == "Channel{Any}(3) (2 items available)"
    close(c)
    @test repr(MIME("text/plain"), c) == "Channel{Any}(3) (closed)"
end

max_stars_repo_name: udesou/julia

max_issues_repo_name: udesou/julia

max_forks_repo_name: udesou/julia


text:
module OpenGLBackend

using ModernGL
using CSyntax
using CSyntax.CEnum
using ..LibCImGui
using ..CImGui: GetIO, GetTexDataAsRGBA32, SetTexID

include("impl.jl")
export ImGui_ImplOpenGL3_Init, ImGui_ImplOpenGL3_Shutdown
export ImGui_ImplOpenGL3_NewFrame
export ImGui_ImplOpenGL3_RenderDrawData
export ImGui_ImplOpenGL3_CreateFontsTexture, ImGui_ImplOpenGL3_DestroyFontsTexture
export ImGui_ImplOpenGL3_CreateDeviceObjects, ImGui_ImplOpenGL3_DestroyDeviceObjects
export ImGui_ImplOpenGL3_CreateImageTexture, ImGui_ImplOpenGL3_UpdateImageTexture, ImGui_ImplOpenGL3_DestroyImageTexture

function __init__()
    global g_GlslVersion = 130
    global g_FontTextures = GLuint[]
    global g_ShaderHandle = GLuint(0)
    global g_VertHandle = GLuint(0)
    global g_FragHandle = GLuint(0)
    global g_AttribLocationTex = GLint(0)
    global g_AttribLocationProjMtx = GLint(0)
    global g_AttribLocationPosition = GLint(0)
    global g_AttribLocationUV = GLint(0)
    global g_AttribLocationColor = GLint(0)
    global g_VboHandle = GLuint(0)
    global g_ElementsHandle = GLuint(0)
    global g_ImageTexture = Dict{Int,GLuint}()
end

end # module


text:
module Stabilizers

## I am abondoning this experiment. This is not a useful way to implement stabilizers

import Random
import ..AbstractStabilizers: AbstractStabilizer
import ..AbstractOps: op_index
import ..AbstractPaulis
import ..Paulis
import ..Paulis: Pauli

export Stabilizer
export pI, mI, pX, mX, pY, mY, pZ, mZ

struct Stabilizer <: AbstractStabilizer
    x::Bool
    z::Bool
    p::Bool
end

const pI = Stabilizer(0, 0, 0)
const mI = Stabilizer(0, 0, 1)
const pX = Stabilizer(1, 0, 0)
const mX = Stabilizer(1, 0, 1)
const pZ = Stabilizer(0, 1, 0)
const mZ = Stabilizer(0, 1, 1)
const pY = Stabilizer(1, 1, 0)
const mY = Stabilizer(1, 1, 1)

function Base.show(io::IO, s::Stabilizer)
    if s.p
        print(io, "-")
    else
        print(io, "+")
    end
    if s.x && s.z
        print(io, "𝒴")
    elseif s.x
        print(io, "𝒳")
    elseif s.z
        print(io, "𝒵")
    else
        print(io, "ℐ")
    end
end

function Random.rand(rng::Random.AbstractRNG, ::Random.SamplerType{Stabilizer})
    return Stabilizer(rand(rng, Bool), rand(rng, Bool), rand(rng, Bool))
end

function Base.:*(o1::Stabilizer, o2::Stabilizer)
    (x1, z1) = (o1.x, o1.z)
    (x2, z2) = (o2.x, o2.z)
    _sign = (x1 & z2 & (x2 | z1)) | (~x1 & x2 & z1 & ~z2)
    phase = (o1.p + o2.p + _sign) % 2
    return Stabilizer(x1 ⊻ x2, z1 ⊻ z2, phase)
end

## The symplectic order is X, Z, Y if we interpret bits as a binary index.
## So, we have to use another method
"""
    op_index(s::Stabilizer)

Return the index in `0--7` of `s`. The order is
`pI, pX, pY, pZ, mI, mX, mY, mZ`.
"""
function op_index(s::Stabilizer)
    if s.x
        if s.z
            res = 2
        else
            res = 1
        end
    else
        if s.z
            res = 3
        else
            res = 0
        end
    end
    return res + 4 * s.p
end

"""
    Pauli(s::Stabilizer)

Return the `Pauli` corresponding to `s`, disregarding
the sign and the factor of `i` in the `Y` stabilizer.
"""
function Paulis.Pauli(s::Stabilizer)
    if s.x
        if s.z
            return(Paulis.Y)
        else
            return(Paulis.X)
        end
    else
        if s.z
            return(Paulis.Z)
        else
            return(Paulis.I)
        end
    end
end

const _pI_mat = Matrix(Paulis.I)
const _pX_mat = Matrix(Paulis.X)
const _pZ_mat = Matrix(Paulis.Z)
const _pY_mat = [0.0 1.0; -1.0 0.0]
const _mI_mat = - _pI_mat
const _mX_mat = - _pX_mat
const _mY_mat = - _pY_mat
const _mZ_mat = - _pZ_mat

# function Base.Matrix(s::Stabilizer)
# end

end # module Stabilizers


text:
module ViZDoom

using Artifacts
using CxxWrap
const vizdoom_artifact_dir = artifact"ViZDoom"
@wrapmodule(joinpath(vizdoom_artifact_dir, "ViZDoom-1.1.8", "bin", "libvizdoomjl.so"))

function __init__()
    @initcxx
end

include("util.jl")
export get_scenario_path, set_game

include("games/games.jl")
export basic_game

end # module


text:
using Manifolds, Documenter

makedocs(
    # for development, we disable prettyurls
    format = Documenter.HTML(prettyurls = false),
    modules = [Manifolds],
    sitename = "Manifolds",
    pages = [
        "Home" => "index.md",
        "Manifolds" => [
            "Sphere" => "manifolds/sphere.md",
            "Rotations" => "manifolds/rotations.md",
        ],
        "Distributions" => "distributions.md"
    ]
)


text:
using Test
using CLIMAParameters: AbstractEarthParameterSet
using CLIMAParameters.Planet
using CLIMAParameters.Atmos.SubgridScale

@testset "Atmos SubgridScale" begin
  struct EarthParameterSet <: AbstractEarthParameterSet end

  ps = EarthParameterSet()
  # Atmos.SubgridScale
  @test c_1_KASM(ps)       ≈ c_a_KASM(ps)*0.76^2
  @test c_2_KASM(ps)       ≈ c_e2_KASM(ps)+2*c_1_KASM(ps)
  @test c_3_KASM(ps)       ≈ c_a_KASM(ps)^(3/2)

end


text:
"""
    SelectTransform(dims::AbstractVector{Int})

Select the dimensions `dims` that the kernel is applied to.
```
    dims = [1,3,5,6,7]
    tr = SelectTransform(dims)
    X = rand(100,10)
    transform(tr,X,obsdim=2) == X[dims,:]
```
"""
struct SelectTransform{T<:AbstractVector{Int}} <: Transform
    select::T
    function SelectTransform{V}(dims::V) where {V<:AbstractVector{Int}}
        @assert all(dims .> 0) "Selective dimensions should all be positive integers"
        return new{V}(dims)
    end
end

SelectTransform(x::T) where {T<:AbstractVector{Int}} = SelectTransform{T}(x)

set!(t::SelectTransform{<:AbstractVector{T}}, dims::AbstractVector{T}) where {T<:Int} = t.select .= dims

duplicate(t::SelectTransform,θ) = t

(t::SelectTransform)(x::AbstractVector) = view(x, t.select)

_map(t::SelectTransform, x::ColVecs) = ColVecs(view(x.X, t.select, :))
_map(t::SelectTransform, x::RowVecs) = RowVecs(view(x.X, :, t.select))

Base.show(io::IO, t::SelectTransform) = print(io, "Select Transform (dims: ", t.select, ")")


text:
# RAPTOR: Round-Based Public Transit Optimized Routing
# Described in Delling, D., Pajor, T., & Werneck, R. (2012). Round-Based Public Transit Routing. http://research.microsoft.com/pubs/156567/raptor_alenex.pdf

using Dates

const MAX_TIME = 2147483646
const XFER_ROUTE = -1
const INT_MISSING = -2
const ORIGIN = -3
const BOARD_SLACK_SECONDS = 60
const EMPTY_SET = BitSet()

struct StopAndTime
    stop::Int64
    time::Int32
end

struct RaptorRequest
    origins::Vector{StopAndTime}
    max_rides::Int64
    date::Date
    walk_speed_meters_per_second::Float64
end

struct RaptorResult
    times_at_stops_each_round::Array{Int32, 2}
    prev_stop::Array{Int64, 2}
    prev_route::Array{Int64, 2}
    prev_boardtime::Array{Int32, 2}
end

function empty_no_resize!(s::BitSet)
    for set_bit in s
        delete!(s, set_bit)
    end
end

function raptor(net::TransitNetwork, req::RaptorRequest)
    nstops = length(net.stops)
    # * 2 for transfer rounds, + 1 for origin times, -1 for skipped final transfer round
    nrounds = req.max_rides * 2 + 1 - 1

    # these get allocated here, and the core RAPTOR algorithm should have 0 allocations
    times_at_stops::Array{Int32,2} = fill(MAX_TIME, (nrounds, nstops))
    prev_stop::Array{Int64,2} = fill(INT_MISSING, (nrounds, nstops))
    prev_route::Array{Int64,2} = fill(INT_MISSING, (nrounds, nstops))
    prev_boardtime::Array{Int32,2} = fill(INT_MISSING, (nrounds, nstops))
    # set bit 0 so that offset is forced to zero and there aren't allocations later
    prev_touched_stops::BitSet = BitSet([0])
    touched_stops::BitSet = BitSet([0])
    sizehint!(prev_touched_stops, nstops)
    sizehint!(touched_stops, nstops)

    # unset the 0 set bits that were to force offset
    delete!(prev_touched_stops, 0)
    delete!(touched_stops, 0)

    # initialize times at stops
    for sat in req.origins
        times_at_stops[1, sat.stop] = sat.time
        push!(prev_touched_stops, sat.stop)
    end

    @assert prev_touched_stops.offset == 0
    @assert touched_stops.offset == 0

    # get which service idxes are running
    services_running = BitSet(map(t -> t[1], filter(t -> is_service_running(t[2], req.date), collect(enumerate(net.services)))))

    @info "$(length(services_running)) services running on requested date"

    # ideally this would have no allocations, although it does have a few due to empty!ing and push!ing to the bitsets - would be nice to have
    # a bounded bitset implementation that did not dynamically resize.
    @time run_raptor!(net, times_at_stops, prev_stop, prev_route, prev_boardtime, req, services_running, prev_touched_stops, touched_stops)

    return RaptorResult(
        times_at_stops,
        prev_stop,
        prev_route,
        prev_boardtime
    )
end

function run_raptor!(net::TransitNetwork, times_at_stops::Array{Int32, 2}, prev_stop::Array{Int64, 2},
    prev_route::Array{Int64, 2}, prev_boardtime::Array{Int32, 2}, req::RaptorRequest, services_running::BitSet, prev_touched_stops::BitSet, touched_stops::BitSet)
    for round in 1:req.max_rides
        # where the results of this round will be recorded
        target = round * 2

        for stop in prev_touched_stops
            # find all patterns that touch this stop
            for patidx in net.patterns_for_stop[stop]
                tp = net.patterns[patidx]

                if !in(tp.service, services_running)
                    continue
                end

                # TODO handle loop trips
                stoppos = INT_MISSING

                # not using a vectorized findfirst to avoid allocations
                for (i, tpstop) in enumerate(tp.stops)
                    if tpstop == stop
                        stoppos = i
                        break
                    end
                end

                @assert stoppos != INT_MISSING

                # find the trip that departs at or after the earliest possible departure
                earliest_departure::Int32 = times_at_stops[target - 1, stop] + BOARD_SLACK_SECONDS
                best_departure::Int32 = MAX_TIME
                best_trip_idx::Int64 = INT_MISSING

                for tripidx in net.trips_for_pattern[patidx]
                    trip = net.trips[tripidx]
                    time_at_stop = trip.stop_times[stoppos].departure_time
                    if (time_at_stop >= earliest_departure && time_at_stop < best_departure)
                        best_trip_idx = tripidx
                        best_departure = time_at_stop
                        # pre-sorting trips by departure time would allow us to break here, but I think we found in R5 that it didn't really help
                    end
                end

                if best_trip_idx == INT_MISSING
                    continue  # no possible trip to board
                end

                best_trip = net.trips[best_trip_idx]

                for stidx in stoppos + 1:length(best_trip.stop_times)
                    stop_time = best_trip.stop_times[stidx]
                    if stop_time.arrival_time < times_at_stops[target, stop_time.stop]
                        # we have found a new fastest way to get to this stop!
                        times_at_stops[target, stop_time.stop] = stop_time.arrival_time
                        prev_stop[target, stop_time.stop] = stop
                        prev_route[target, stop_time.stop] = best_trip.route
                        prev_boardtime[target, stop_time.stop] = best_trip.stop_times[stoppos].departure_time
                        push!(touched_stops, stop_time.stop)
                    end
                end
            end
        end # loop over prev_touched_stop

        #@info "round $round found $(length(touched_stops)) stops accessible by transit"

        # clear prev_touched_stops and reuse as next_touched_stops, avoid allocation
        empty!(prev_touched_stops)
        next_touched_stops = prev_touched_stops

        # do transfers, but skip after last iteration
        if round < req.max_rides
            times_at_stops[target + 1, :] = times_at_stops[target, :]
            # leave other things missing if there were no transfers
            for stop in touched_stops
                push!(next_touched_stops, stop)  # this stop was touched by this round

                for xfer in net.transfers[stop]
                    xfer_walk_time = Base.round(xfer.distance_meters / req.walk_speed_meters_per_second)
                    pre_xfer_time = times_at_stops[target, stop]
                    time_after_xfer = pre_xfer_time + xfer_walk_time
                    if time_after_xfer < times_at_stops[target + 1, xfer.target_stop]
                        # transferring to this stop is optimal!
                        times_at_stops[target + 1, xfer.target_stop] = time_after_xfer
                        prev_stop[target + 1, xfer.target_stop] = stop
                        prev_route[target + 1, xfer.target_stop] = XFER_ROUTE
                        prev_boardtime[target + 1, stop] = pre_xfer_time
                        push!(next_touched_stops, xfer.target_stop)
                    end
                end
            end
        end

        # prepare for next iteration
        prev_touched_stops = next_touched_stops
        empty!(touched_stops)
    end
end


text:
module JuliaStudies

# Write your package code here.

end


text:
using Test
import CryptoGroups: PGroup, validate, order, Enc, Dec, modulus, value, specialize, Specs, generator, <|

q = 11
p = 2*q + 1

G = PGroup(p, q; name=:G)

@test validate(G(3)) == true
@test validate(G(11)) == false

n = let
    n = 0
    for i in 2:p-1
        validate(G(i)) && (n+=1)
    end
    n
end
@test n == q - 1

g = G(3)

@test g*g^2 == g^3
@test (g^2)^2 == g^4
@test g^(order(G) + 1) == g

h = g^7

@test h*h^2 == h^3
@test (h^2)^2 == h^4
@test h^(order(G) + 1) == h

@test inv(g)*g^2 == g
@test (g^7)^6 == g^(7*6) # This is only true for a cyclic group
@test g*g*g == g^3 # Checking multiplication
@test g^2/g == g

################################ Legacy ##################################

import CryptoGroups
import CryptoGroups: <|, PGroup, specialize, value, order, ECGroup, Group, modulus, generator, Specs

using Test

### Testing prime groups
using Primes

function testgroup(g::G) where G <: Group
    @test g*g^2 == g^3
    @test (g^2)^2 == g^4
    @test g^(order(G) + 1) == g
end

G = PGroup(23, 11)
g = G(3)

@test value(g) == 3
@test order(G) == 11
@test modulus(G) == 23

testgroup(g)

### Testing ECGroup

let
    spec = Specs.Curve_P_256

    G = specialize(ECGroup, spec; name = :P_256)
    g = G <| generator(spec)

    testgroup(g)
end

###  RFC standart group test with PGroup

let
    spec = Specs.MODP_1024
    G = specialize(PGroup, spec)
    g = G <| generator(spec)

    testgroup(g)
end
### PGroup generation algorithms for

using Random: MersenneTwister
rng = MersenneTwister(0)

# Broken
let
modp_spec = Specs.sophie_germain_group(rng, 2, 100)
G = specialize(PGroup, modp_spec)
g = G <| generator(modp_spec)

#testgroup(g)
end

### Seems to have an issue with large numbers
let
modp_spec = Specs.dsa_standart_group(rng, 10, 10)
G = specialize(PGroup, modp_spec)
g = G <| generator(modp_spec)

testgroup(g)
end


text:
using DemeNet: ID, DemeID
using Base: UUID, Dict
using Sockets

import Base.Dict
import PeaceVote: BraiderConfig, RecorderConfig, BraidChainConfig
import Recruiters: CertifierConfig

struct Port
    port::Int
    ip::Union{Nothing,IPv4,IPv6}
end

Port(port::Int) = Port(port,nothing)

function Dict(port::Port)
    dict = Dict{String,Union{String,Int}}("port"=>port.port)
    if !isnothing(port.ip)
        dict["ip"] = string(port.ip)
        dict["type"] = "$(typeof(port.ip))"
    else
        dict["type"] = "Int"
    end

    return dict
end

function Port(dict::Dict)
    port = dict["port"]
    type = dict["type"]
    if type == "Int"
        ip = nothing
    elseif type == "IPv4"
        ip = IPv4(dict["ip"])
    elseif type == "IPv6"
        ip = IPv6(dict["ip"])
    end
    return Port(port,ip)
end

import Sockets.listen
function listen(port::Port)
    if port.ip isa Nothing
        listen(port.port)
    elseif port.ip isa IPv4
        listen(IPv4(0),port.port)
    else
        listen(IPv6(0),port.port)
    end
end

import Sockets.connect
function connect(port::Port)
    if port.ip isa Nothing
        connect(port.port)
    else
        connect(port.ip,port.port)
    end
end

# I could add also a hash of demefile here!
struct AddressRecord
    id::Union{ID,DemeID}
    hash::Union{Nothing,BigInt}
    ip::Union{IPv4,IPv6}
end

function ip(machines::Vector{AddressRecord},id::Union{ID,DemeID})
    for m in machines
        if m.id==id
            return m.ip
        end
    end
end

function Dict(record::AddressRecord)
    dict = Dict()

    if typeof(record.id)==DemeID
        dict["uuid"] = string(record.id.uuid)
    end

    if record.hash!=nothing
        dict["hash"] = string(record.hash,base=16)
    end

    dict["id"] = string(record.id.id,base=16)
    dict["ip"] = string(record.ip)
    return dict
end

Dict(arecords::Vector{AddressRecord}) = Dict[Dict(i) for i in arecords]

function AddressRecord(config::Dict)
    localid = parse(BigInt,config["id"],base=16)

    if "uuid" in config.keys
        uuid = UUID(config["uuid"])
        id = DemeID(uuid,localid)
    else
        id = ID(localid)
    end

    if "hash" in config.keys
        hash = parse(BigInt,config["hash"],base=16)
    else
        hash = nothing
    end

    ip = IPv4(config["ip"])

    return AddressRecord(id,hash,ip)
end

function Dict(config::CertifierConfig{Port})
    dict = Dict()
    dict["ca"] = string(config.tookenca.id,base=16)
    dict["server"] = string(config.serverid.id,base=16)
    dict["tport"] = config.tookenport.port
    dict["cport"] = config.certifierport.port
    return dict
end

function CertifierConfig{Port}(dict::Dict,arecords::Vector{AddressRecord})
    ca = ID(parse(BigInt,dict["ca"],base=16))
    server = ID(parse(BigInt,dict["server"],base=16))

    addr = ip(arecords,server)
    tport = Port(dict["tport"],addr)
    cport = Port(dict["cport"],addr)

    CertifierConfig(ca,server,tport,cport)
end

function Dict(config::BraiderConfig{Port})
    dict = Dict()
    dict["N"] = Int(config.N)
    dict["M"] = Int(config.M)
    dict["server"] = string(config.gateid.id,base=16)
    dict["port"] = config.port.port
    mixer = Dict()
    mixer["uuid"] = string(config.mixerid.uuid)
    mixer["server"] = string(config.mixerid.id,base=16)
    mixer["port"] = config.ballotport.port
    dict["mixer"] = mixer
    return dict
end

function BraiderConfig{Port}(dict::Dict,arecords::Vector{AddressRecord})
    N = UInt8(dict["N"])
    M = UInt8(dict["M"])

    server = ID(parse(BigInt,dict["server"],base=16))
    port = Port(dict["port"],ip(arecords,server))

    muuid = UUID(dict["mixer"]["uuid"])
    mserver = parse(BigInt,dict["mixer"]["server"],base=16)

    mid = DemeID(muuid,mserver)
    bport = Port(dict["mixer"]["port"],ip(arecords,mid)) ### Braider needs to know the ip address of the mixer

    BraiderConfig(port,bport,N,M,server,mid)
end

function Dict(config::RecorderConfig{Port})
    ca = Dict[]
    for mca in config.membersca
        dict = Dict()
        dict["id"] = string(mca.id,base=16)
        if typeof(mca)==DemeID
            dict["uuid"] = string(mca.uuid)
        end
        push!(ca,dict)
    end

    dict = Dict()
    dict["ca"] = ca
    dict["server"] = string(config.serverid.id,base=16)
    dict["rport"] = config.registratorport.port
    dict["vport"] = config.votingport.port
    dict["pport"] = config.proposalport.port

    return dict
end

function RecorderConfig{Port}(dict::Dict,arecords::Vector{AddressRecord})

    ca = Union{ID,DemeID}[]
    for i in dict["ca"]
        localid = parse(BigInt,i["id"],base=16)
        if haskey(i,"uuid")
            uuid = UUID(i["uuid"])
            push!(ca,DemeID(uuid,localid))
        else
            push!(ca,ID(localid))
        end
    end

    server = ID(parse(BigInt,dict["server"],base=16))

    addr = ip(arecords,server)
    rport = Port(dict["rport"],addr)
    vport = Port(dict["vport"],addr)
    pport = Port(dict["pport"],addr)

    RecorderConfig(ca,server,rport,vport,pport)
end

function Dict(config::BraidChainConfig{Port})

    dict = Dict()
    dict["server"] = string(config.server,base=16)
    dict["mport"] = config.mixerport.port
    dict["sport"] = config.syncport.port
    dict["braider"] = Dict(config.braider)
    dict["recorder"] = Dict(config.recorder)

    return dict
end

function BraidChainConfig{Port}(dict::Dict,arecords::Vector{AddressRecord})

    server = ID(parse(BigInt,dict["server"],base=16))
    addr = ip(arecords,server)
    mport = Port(dict["mport"],addr)
    sport = Port(dict["sport"],addr)

    braider = BraiderConfig{Port}(dict["braider"],arecords)
    recorder = RecorderConfig{Port}(dict["recorder"],arecords)

    BraidChainConfig(server,mport,sport,braider,recorder)
end

struct PeaceFounderConfig
    braidchain::BraidChainConfig{Port}
    certifier::CertifierConfig{Port}
    arecords::Vector{AddressRecord}
end

function Dict(config::PeaceFounderConfig)
    dict = Dict("braidchain"=>Dict(config.braidchain),"certifier"=>Dict(config.certifier))
    length(config.arecords)>0 && (dict["arecords"]=Dict(config.arecords))
    return dict
end

function PeaceFounderConfig(dict::Dict)

    if haskey(dict,"arecords")
        arecords = AddressRecord[AddressRecord(i) for i in dict["arecords"]]
    else
        arecords = AddressRecord[]
    end

    certifier = CertifierConfig{Port}(dict["certifier"],arecords)
    braidchain = BraidChainConfig{Port}(dict["braidchain"],arecords)

    return PeaceFounderConfig(braidchain,certifier,arecords)
end

max_stars_repo_name: akels/DemocracyCredit

max_issues_repo_name: akels/DemocracyCredit

max_forks_repo_name: akels/DemocracyCredit


text:
"""
    abstract type Encoding

Transformation of `Block`s. Can encode some `Block`s ([`encode`]), and optionally
decode them [`decode`]

## Interface

- `encode(::E, ::Context, block::Block, data)` encodes `block` of `data`.
    The default is to do nothing. This should be overloaded for an encoding `E`,
    concrete `Block` types and possibly a context.
- `decode(::E, ::Context, block::Block, data)` decodes `block` of `data`. This
    should correspond as closely as possible to the inverse of `encode(::E, ...)`.
    The default is to do nothing, as not all encodings can be reversed. This should
    be overloaded for an encoding `E`, concrete `Block` types and possibly a context.
- `encodedblock(::E, block::Block) -> block'` returns the block that is obtained by
    encoding `block` with encoding `E`. This needs to be constant for an instance of `E`,
    so it cannot depend on the sample or on randomness. The default is to return `nothing`,
    meaning the same block is returned and not changed. Encodings that return the same
    block but change the data (e.g. `ProjectiveTransforms`) should return `block`.
- `decodedblock(::E, block::Block) -> block'` returns the block that is obtained by
    decoding `block` with encoding `E`. This needs to be constant for an instance of `E`,
    so it cannot depend on the sample or on randomness. The default is to return `nothing`,
    meaning the same block is returned and not changed.
- `encode!(buf, ::E, ::Context, block::Block, data)` encodes `data` inplace.
- `decode!(buf, ::E, ::Context, block::Block, data)` decodes `data` inplace.

"""
abstract type Encoding end

"""

    fillblock(inblocks, outblocks)

Replaces all `nothing`s in outblocks with the corresponding block in `inblocks`.
`outblocks` may be obtained by
"""
fillblock(inblocks::Tuple, outblocks::Tuple) = map(fillblock, inblocks, outblocks)
fillblock(inblock::AbstractBlock, ::Nothing) = inblock
fillblock(::AbstractBlock, outblock::AbstractBlock) = outblock

function encodedblock(enc, block, fill::Bool)
    outblock = encodedblock(enc, block)
    return fill ? fillblock(block, outblock) : outblock
end

function decodedblock(enc, block, fill::Bool)
    outblock = decodedblock(enc, block)
    return fill ? fillblock(block, outblock) : outblock
end

# ## `encode` methods

# By default an encoding doesn't change the data
encode(encoding::Encoding, ctx, block::Block, data; kwargs...) =
    isempty(kwargs) ? data : encode(encoding, ctx, block, data)

# By default, a tuple of encodings encodes by encoding the data one encoding
# after the other
"""
    encode(encoding, context, block, data)
    encode(encoding, context, blocks, data)
    encode(encodings, context, blocks, data)
"""
function encode(encodings::NTuple{N, Encoding}, context, blocks, data) where N
    for encoding in encodings
        data = encode(encoding, context, blocks, data)
        blocks = encodedblock(encoding, blocks, true)
    end
    return data
end

# By default, an encoding encodes every element in a tuple separately
function encode(encoding::Encoding, context, blocks::Tuple, datas::Tuple)
    @assert length(blocks) == length(datas)
    return Tuple(encode(encoding, context, block, data)
                    for (block, data) in zip(blocks, datas))
end

# Named tuples of data are handled like tuples, but the keys are preserved
function encode(encoding::Encoding, context, blocks::NamedTuple, datas::NamedTuple)
    @assert length(blocks) == length(datas)
    return NamedTuple(zip(
        keys(datas),
        encode(encoding, context, values(blocks), values(datas))
    ))
end

# ## `decode` methods

# By default an encoding doesn't change the data
decode(encoding::Encoding, ctx, block::Block, data; kwargs...) =
    isempty(kwargs) ? data : decode(encoding, ctx, block, data)

# By default, a tuple of encodings decodes by decoding the data one encoding
# after the other, with encodings iterated in reverse order
function decode(encodings::NTuple{N, Encoding}, context, blocks, data) where N
    for encoding in Iterators.reverse(encodings)
        data = decode(encoding, context, blocks, data)
        blocks = decodedblock(encoding, blocks, true)
    end
    return data
end

# By default, an encoding decodes every element in a tuple separately
function decode(encoding::Encoding, context, blocks::Tuple, datas::Tuple)
    @assert length(blocks) == length(datas)
    return Tuple(decode(encoding, context, block, data)
                    for (block, data) in zip(blocks, datas))
end

# Named tuples of data are handled like tuples, but the keys are preserved
function decode(encoding::Encoding, context, blocks::NamedTuple, datas::NamedTuple)
    @assert length(blocks) == length(datas)
    return NamedTuple(zip(
        keys(datas),
        decode(encoding, context, values(blocks), values(datas))
    ))
end

"""
    encodedblock(encoding, block)
    encodedblock(encoding, blocks)
    encodedblock(encodings, blocks)

Return the block that is obtained by encoding `block` with encoding `E`.
This needs to be constant for an instance of `E`, so it cannot depend on the
sample or on randomness. The default is to return `nothing`,
meaning the same block is returned and not changed. Encodings that return the same
block but change the data (e.g. `ProjectiveTransforms`) should return `block`.
"""
encodedblock(enc, block) = nothing
function encodedblock(encoding::Encoding, blocks::Tuple)
    Tuple(encodedblock(encoding, block) for block in blocks)
end
function encodedblock(encodings::NTuple{N, Encoding}, blocks) where N
    encoded = false
    for encoding in encodings
        encoded = encoded || !isnothing(encodedblock(encoding, blocks))
        blocks = encodedblock(encoding, blocks, true)
    end
    return encoded ? blocks : nothing
end

"""
    decodedblock(encoding, block)
    decodedblock(encoding, blocks)

Return the block that is obtained by decoding `block` with encoding `E`.
This needs to be constant for an instance of `E`, so it cannot depend on the
sample or on randomness. The default is to return `nothing`,
meaning the same block is returned and not changed. Encodings that return the same
block but change the data when decoding should return `block`.
"""
decodedblock(::Encoding, ::Block) = nothing
function decodedblock(encoding::Encoding, blocks::Tuple)
    Tuple(decodedblock(encoding, block) for block in blocks)
end
function decodedblock(encodings::NTuple{N, Encoding}, blocks) where N
    decoded = false
    for encoding in Iterators.reverse(encodings)
        decoded = decoded || !isnothing(decodedblock(encoding, blocks))
        blocks = decodedblock(encoding, blocks, true)
    end
    return decoded ? blocks : nothing
end

"""
    abstract type StatefulEncoding <: Encoding

Encoding that needs to compute some state from the whole sample, even
if it only transforms some of the blocks. This could be random state
for stochastic augmentations that needs to be the same for every block
that is encoded.

The state is created by calling `encodestate(encoding, context, blocks, sample)`
and passed to recursive calls with the keyword argument `state`.
As a result, you need to implement `encode`, `decode`, `encode!`, `decode!` with a
keyword argument `state` that defaults to the above call.

Same goes for `decode`, which should accept a `state` keyword argument defaulting
to `decodestate(encoding, context, blocks, sample)`
"""
abstract type StatefulEncoding <: Encoding end

encodestate(encoding, context, blocks, data) = nothing
decodestate(encoding, context, blocks, data) = nothing

function encode(
        encoding::StatefulEncoding,
        context,
        blocks::Tuple,
        datas::Tuple;
        state = encodestate(encoding, context, blocks, datas))

    @assert length(blocks) == length(datas)
    return Tuple(encode(encoding, context, block, data; state = state)
                    for (block, data) in zip(blocks, datas))
end

function decode(
        encoding::StatefulEncoding,
        context,
        blocks::Tuple,
        datas::Tuple;
        state = decodestate(encoding, context, blocks, datas))

    @assert length(blocks) == length(datas)
    return Tuple(decode(encoding, context, block, data; state = state)
                    for (block, data) in zip(blocks, datas))
end

"""
    testencoding(encoding, block, data)

Performs some tests that the encoding interface is set up properly for
`encoding` and `block`. Tests that

- `data` is a valid `block`
- `encode` returns a valid `encodedblock(encoding, block)`
- `decode` returns a valid `decodedblock(encoding, encodedblock(encoding, block))`
    and that the block is identical to `block`
"""
function testencoding(encoding, block, data = mockblock(block))
    Test.@testset "Encoding `$(typeof(encoding))` for block `$block`" begin
        # Test that `data` is a valid instance of `block`
        Test.@test checkblock(block, data)
        Test.@test !isnothing(encodedblock(encoding, block))
        outblock = encodedblock(encoding, block, true)
        outdata = encode(encoding, Training(), block, data)
        # The encoded data should be a valid instance for the `encodedblock`
        Test.@test checkblock(outblock, outdata)

        # Test decoding (if supported) works correctly
        if (outblock isa Tuple)
            for idx in length(outblock)
                inblock = decodedblock(encoding, outblock[idx])
                if !isnothing(inblock)
                    Test.@test wrapped(block[idx]) == wrapped(inblock)
                    indata = decode(encoding, Training(), outblock[idx], outdata[idx])
                    Test.@test checkblock(inblock, indata)
                end
            end
        else
            inblock = decodedblock(encoding, outblock)
            if !isnothing(inblock)
                Test.@test wrapped(block) == wrapped(inblock)
                indata = decode(encoding, Training(), outblock, outdata)
                Test.@test checkblock(inblock, indata)
            end
        end
    end
end


text:
# nice idea but with the current need to build js ourselves
# it makes more sense to include it.

# function getfile(fname)
# 	url = "http://cdn.pydata.org/" * fname
# 	path = joinpath(dirname(@__FILE__), "bokehjs", fname))
# 	println("downloading from $url to $path")
# 	download(url, path)
# end

# bokehjs = "bokeh-0.7.1"
# getfile("$bokehjs.js")
# getfile("$bokehjs.min.js")
# getfile("$bokehjs.css")
# getfile("$bokehjs.min.css")


text:
using ..CaimCore

abstract type Feature <: CaimCore.Feature end

function Base.iterate(f::Feature)
    b = box(f)
    next = iterate(b)
    while !isnothing(next)
        p, state = next
        if p in f
            return p, (b, state)
        end
        next = iterate(b, state)
    end
    nothing
end

function Base.iterate(f::Feature, (b, state))
    next = iterate(b, state)
    while !isnothing(next)
        p, state = next
        if p in f
            return p, (b, state)
        end
        next = iterate(b, state)
    end
    nothing
end

function CaimCore.series(d::VideoGraphicDataset{C}, f::Feature;
						 Acc::Type{A}=SqrtColorAccumulator{C},
						 kwargs...) where {C, D, A <: ColorAccumulator{C,D}}
    width, height, duration = size(d)
    series = Array{D}(undef, duration)
    acc = Acc()
    for t in 1:duration
        for p in f
            if all([1,1] .<= p .<= [width,height])
            	add!(acc, d.frames[p[1], p[2], t])
            end
        end
        series[t] = color(acc)
        reset!(acc)
    end
    series
end

const Point = AbstractVector{<:Integer}

macro ensure2d(ps...)
    Expr(:block, map(ps) do p
        :(length($(esc(p))) == 2 || throw(ArgumentError("point must have length 2")))
    end...)
end

struct Box{P <: Point} <: Feature
    tl::P
    br::P
    function Box(a::P, b::P) where {P <: Point}
        @ensure2d a b
        new{P}(min.(a, b), max.(a, b))
    end
end
width(b::Box) = 1 + b.br[1] - b.tl[1]
height(b::Box) = 1 + b.br[2] - b.tl[2]
area(b::Box) = width(b) * height(b)

Base.iterate(b::Box) = (b.tl, (b.tl, 1))
function Base.iterate(b::Box, state)
    p, i = state
    if i >= area(b)
        nothing
    else
        q = (p[1] < b.br[1]) ? [p[1] + 1, p[2]] : [b.tl[1], p[2] + 1]
        q, (q, i + 1)
    end
end
Base.length(b::Box) = area(b)
Base.eltype(b::Box{P}) where P = P

Base.:(==)(b::Box{P}, c::Box{P}) where P = b.tl == c.tl && b.br == c.br
Base.:(==)(b::Box{P}, c::Box{Q}) where {P, Q} = false

Base.in(p::P, b::Box{P}) where {P <: Point} = all(b.tl .<= p .<= b.br)
Base.in(p, b::Box) = false

struct Circle{P <: Point} <: Feature
    center::P
    radius::Float64
    radius²::Float64
    function Circle(c::Point, r::Real)
        @ensure2d c
        new{typeof(c)}(c, abs(r), r^2)
    end
end

function Circle(c::Point, b::Point)
    @ensure2d c b
    r = norm(b - c)
    Circle(c, r)
end

function box(c::Circle)
    tl = round.(Int, c.center .- c.radius)
    br = round.(Int, c.center .+ c.radius)

    Box(tl, br)
end

function Base.in(p::P, c::Circle{P}) where {P <: Point}
    if length(p) != 2
        false
    else
        rx, ry = p .- c.center
        rx^2 + ry^2 ≤ c.radius²
    end
end
Base.in(p, c::Circle) = false
Base.IteratorSize(::Circle) = Base.SizeUnknown()

struct FreeForm{P <: Point} <: Feature
    boundary::Vector{P}
    function FreeForm(b::Vector{P}) where {P <: Point}
        if isempty(b)
            throw(ArgumentError("boundary must be non-empty"))
        end
        for p in b
            @ensure2d p
        end
        new{P}(b)
    end
end

function box(f::FreeForm)
    tl = f.boundary[1]
    br = f.boundary[1]
    for p in f.boundary
        tl = min.(tl, p)
        br = max.(br, p)
    end
    Box(tl, br)
end

function Base.in(p::P, f::FreeForm{P}) where {P <: Point}
    if length(p) != 2
        false
    else
        windingnumber = 0
        len = length(f.boundary)
        for i in 1:len
            r, s = f.boundary[i], f.boundary[(i + 1 <= len) ? i + 1 : 1]
            isleft = (s[1] - r[1]) * (p[2] - r[2])  - (p[1] -  r[1]) * (s[2] - r[2])

            if all(((r .<= p) .& (p .<= s)) .| ((s .<= p) .& (p .<= r)))
                if isleft == 0
                    windingnumber = 1
                    break
                end
            end

            if r[2] <= p[2]
                if s[2] > p[2] && isleft > 0
                    windingnumber += 1
                else
                end
            elseif s[2] <= p[2] && isleft < 0
                windingnumber -= 1
            end
        end

        windingnumber != 0
    end
end
Base.in(p, f::FreeForm) = false
Base.IteratorSize(::FreeForm) = Base.SizeUnknown()


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    PlainGraph, ImmutablePlainGraph,
    plaingraph, immutableplaingraph

struct PlainGraph <: OrderedGraph
    neighbormap::Vector{Dict{Int,Int}}
    edges::Vector{Tuple{Int,Int}}
    cache::Dict{Symbol,Any}
end

struct ImmutablePlainGraph <: OrderedGraph
    neighbormap::Vector{Dict{Int,Int}}
    edges::Vector{Tuple{Int,Int}}
    cache::Dict{Symbol,Any}
end

"""
    plaingraph() -> PlainGraph

Generate empty `PlainGraph`.
"""
plaingraph() = PlainGraph([], [], Dict())

"""
    plaingraph(nodes, edges) -> PlainGraph{PlainNode,Edge}

Generate vector graph that have given nodes and edges represented by the list of
node indices in integer and the list of pairs of node indices, respectively.
"""
function plaingraph(size::Int, edges)
    nbrmap = [Dict{Int,Int}() for i in 1:size]
    edges = collect(edges)
    for (i, (u, v)) in enumerate(edges)
        nbrmap[u][i] = v
        nbrmap[v][i] = u
    end
    return PlainGraph(nbrmap, edges, Dict{Symbol,Any}())
end

"""
    plaingraph(graph::OrderedGraph) -> PlainGraph

Convert an arbitrary `OrderedGraph` into a `PlainGraph`.
"""
function plaingraph(graph::OrderedGraph)
    newg = plaingraph()
    for nbr in graph.neighbormap
        push!(newg.neighbormap, copy(nbr))
    end
    append!(newg.edges, graph.edges)
    return newg
end

"""
    plaingraph(graph::UndirectedGraph) -> PlainGraph

Convert an arbitrary `UndirectedGraph` into a `PlainGraph`.

If the given graph has non-sequential indices (ex. subgraph view), node and
edge indices are sorted in ascending order and are re-indexed.
"""
function plaingraph(graph::UndirectedGraph)
    newg = plaingraph()
    nkeys = sort(collect(nodeset(graph)))
    ekeys = sort(collect(edgeset(graph)))
    nmap = Dict{Int,Int}()
    for (i, n) in enumerate(nkeys)
        nmap[n] = i
        push!(newg.neighbormap, Dict())
    end
    for (i, e) in enumerate(ekeys)
        (oldu, oldv) = getedge(graph, e)
        u = nmap[oldu]
        v = nmap[oldv]
        push!(newg.edges, (u, v))
        newg.neighbormap[u][i] = v
        newg.neighbormap[v][i] = u
    end
    return newg
end

"""
    immutableplaingraph() -> ImmutablePlainGraph
    immutableplaingraph(size::Int, edges) -> ImmutablePlainGraph
    immutableplaingraph(graph::UndirectedGraph) -> ImmutablePlainGraph

Convert an arbitrary `UndirectedGraph` into a `ImmutablePlainGraph`.

`ImmutablePlainGraph` is same as `PlainGraph` except for availability of
`addnode!` and `addedge!`, and no performance difference exists between these
graphs. `ImmutablePlainGraph` is mainly used for unit testing to check adverse
effect of some graph algorithm implementations.
"""
immutableplaingraph() = ImmutablePlainGraph([], [], Dict())
function immutableplaingraph(size::Int, edges)
    g = plaingraph(size, edges)
    return ImmutablePlainGraph(g.neighbormap, g.edges, g.cache)
end
function immutableplaingraph(graph::UndirectedGraph)
    g = plaingraph(graph)
    return ImmutablePlainGraph(g.neighbormap, g.edges, g.cache)
end

function addnode!(graph::ImmutablePlainGraph)
    throw(ErrorException("Immutable"))
end

function addedge!(graph::ImmutablePlainGraph, u::Int, v::Int)
    throw(ErrorException("Immutable"))
end

function unlinknodes(graph::PlainGraph, nodes)
    # TODO: costful
    subg = nodesubgraph(graph, setdiff(nodeset(graph), nodes))
    return plaingraph(subg)
end

function unlinkedges(graph::PlainGraph, edges)
    # TODO: costful
    subg = SubgraphView(graph, nodeset(graph), setdiff(edgeset(graph), edges))
    return plaingraph(subg)
end


text:
### ============== ### ============== ###
##    Inertial Spin Model             ##
##    Martin Zumaya Hernandez          ##
##    EXAMPLE SIMULATION SCRIPT        ##
### ============== ### ============== ###

### ============ INCLUDE PACKAGES ============ ###

using CollectiveDynamics.InertialSpin

### =============== ### =============== ###
###   DEFINITION OF INITIAL PARAMETERS  ###
### =============== ### =============== ###

# N   = 128
# η   = 1.0
# τ   = 6
# rep = 1
# T   = 8*exp10(-5)
# δ   = 0.35

N    = parse(Int, ARGS[1]) # Number of particles
η    = parse(Float64, ARGS[2]) # Dissipation term
T    = parse(Float64, ARGS[3]) # temperature, noise
τ    = parse(Int, ARGS[4]) # number of iterations
rep  = parse(Int, ARGS[5]) # repetition number

χ   = 1.25
J   = 0.8
ρ   = 0.3
v0  = 0.1
n_t = 6

pars = InertialParameters(χ, J, η, v0*sqrt(J/χ), ρ, v0, N, T, n_t, 0)

σ = sqrt((2pars.d) * η * T) # noise std deviation ( square root of variance )
L = cbrt(N / pars.ρ) # 3D

### =============== ### =============== ###
### SET UP SYSTEM AND OUTPUT STRUCTURE  ###
### =============== ### =============== ###

flock = InertialFlock(N, L, pars.v0)

output_path = set_output_data_structure("INERTIAL_SPIN", N, η, T)

pos_file  = open(output_path * "/pos_$(rep).dat", "w+")
vel_file  = open(output_path * "/vel_$(rep).dat", "w+")
spin_file = open(output_path * "/spin_$(rep).dat", "w+")

### ============== ### ============== ###
###          SYSTEM EVOLUTION         ###
### ============== ### ============== ###

times = [convert(Int, exp10(i)) for i in 0:τ]

  for i in 1:(length(times) - 1)

      if i > 1

          for t in (times[i]+1):times[i+1]

              evolve_system(flock.pos, flock.vel, flock.v_t, flock.spin, flock.Rij, pars, σ)

              if t % times[i] == 0 || t % times[i-1] == 0
                  println("//////// ", t)
                  write(pos_file, vcat(flock.pos...))
                  write(vel_file, vcat(flock.vel...))
                  write(spin_file, vcat(flock.spin...))
              end
          end

      else

          for t in (times[i]+1):times[i+1]

              evolve_system(flock.pos, flock.vel, flock.v_t, flock.spin, flock.Rij, pars, σ)

              if t % times[i] == 0
                  println("//////// ", t)
                  write(pos_file, vcat(flock.pos...))
                  write(vel_file, vcat(flock.vel...))
                  write(spin_file, vcat(flock.spin...))
              end
          end

      end

  end

close(pos_file)
close(vel_file)

println("Done all")


text:
using Compat

immutable IRef{T}
    value::T
    @compat (::Type{IRef{T}}){T}() = new{T}()
    @compat (::Type{IRef{T}}){T}(value) = new{T}(value)
end

function Base.show(io::IO, r::IRef)
    if isdefined(r, :value)
        print(io, r.value)
    else
        print(io, "∘")
    end
end

immutable Result{T,S}
    issuccess::Bool
    value::IRef{T}
    error::IRef{S}
    function (::Type{Result{T,S}}){T,S}(issuccess, val)
        issuccess ?
            new{T,S}(issuccess, IRef{T}(val), IRef{S}()) :
            new{T,S}(issuccess, IRef{T}(), IRef{S}(val))
    end
end

immutable FailedResult end
immutable SuccessResult end

issuccess(e::Result) = e.issuccess
value(e::Result) = e.issuccess ? e.value.value : throw(FailedResult())
geterror(e::Result) = e.issuccess ? throw(SuccessResult()) : e.error.value

function Base.show(io::IO, r::Result)
    if issuccess(r)
        print(io, "✓ $(value(r))")
    else
        print(io, "✗ $(geterror(r))")
    end
end


text:
using DDF

using LinearAlgebra
using SparseArrays
using StaticArrays
using Test

@testset "Sparse matrix iteration" begin
    nrows = 11
    ncols = 12
    A = sprand(Complex{Int8}, nrows, ncols, 0.25)
    nvals = 0
    for col in 1:ncols
        for (row, val) in sparse_column(A, col)
            @test val == A[row, col]
            nvals += 1
        end
    end
    @test nvals == nnz(A)

    for col in 1:ncols
        Ac = sparse_column(A, col)
        Acr = sparse_column_rows(A, col)
        Acv = sparse_column_values(A, col)
        colvals = [row_val for row_val in Ac]
        colvals::Vector{Tuple{Int,Complex{Int8}}}
        @test length(Ac) == length(colvals)
        @test isempty(Ac) == (length(Ac) == 0)
        if !isempty(Ac)
            @test first(Ac) == colvals[1]
            @test last(Ac) == colvals[end]
        end
        for val in colvals
            @test val ∈ Ac
            @test .-val ∉ Ac
        end
        @test collect(Ac) == colvals
        @test collect(Acr) == first.(colvals)
        @test collect(Acv) == last.(colvals)
    end
end

@testset "SparseOp vector space operations" begin
    T = Rational{Int64}
    R1 = 1
    R2 = 2

    nelts1 = rand(5:10)
    nelts2 = rand(5:10)

    z = zero(SparseOp{R1,R2,T}, nelts1, nelts2)
    A = rand(SparseOp{R1,R2,T}, nelts1, nelts2)
    B = rand(SparseOp{R1,R2,T}, nelts1, nelts2)
    C = rand(SparseOp{R1,R2,T}, nelts1, nelts2)
    a = rand(T)
    b = rand(T)

    # Vector space

    @test +A == A
    @test (A + B) + C == A + (B + C)

    @test z + A == A
    @test A + z == A

    @test A + (-A) == z
    @test (-A) + A == z
    @test A - B == A + (-B)

    @test (a * b) * A == a * (b * A)
    @test A * (a * b) == (A * a) * b

    @test one(T) * a == a
    @test one(T) * A == A
    @test A * one(T) == A

    @test zero(T) * a == zero(T)
    @test zero(T) * A == z
    @test A * zero(T) == z

    if a != zero(T)
        @test a * inv(a) == one(T)

        @test inv(a) * (a * A) == A
        @test (A * a) * inv(a) == A

        @test inv(a) * A == a \ A
        @test A * inv(a) == A / a
    end

    @test (a + b) * A == a * A + b * A
    @test A * (a + b) == A * a + A * b

    @test a * (A + B) == a * A + a * B
    @test (A + B) * a == A * a + B * a

    # Adjoints

    @test A'' == A
    @test (a * A)' == A' * a'
    @test (A + B)' == A' + B'
end

@testset "SparseOp category operations" begin
    T = Rational{Int64}
    R1 = 1
    R2 = 2
    R3 = 3
    R4 = 4

    nelts1 = rand(5:10)
    nelts2 = rand(5:10)
    nelts3 = rand(5:10)
    nelts4 = rand(5:10)

    e = one(SparseOp{R1,R1,T}, nelts1, nelts1)
    e2 = one(SparseOp{R2,R2,T}, nelts2, nelts2)
    z = zero(SparseOp{R1,R2,T}, nelts1, nelts2)
    z2 = zero(SparseOp{R2,R3,T}, nelts2, nelts3)
    z3 = zero(SparseOp{R1,R3,T}, nelts1, nelts3)
    A = rand(SparseOp{R1,R2,T}, nelts1, nelts2)
    B = rand(SparseOp{R2,R3,T}, nelts2, nelts3)
    C = rand(SparseOp{R3,R4,T}, nelts3, nelts4)
    a = rand(T)

    @test e * A == A
    @test A * e2 == A

    @test z * B == z3
    @test A * z2 == z3

    @test nnz((z * B).op) == 0
    @test nnz((A * z2).op) == 0

    @test (A * B) * C == A * (B * C)

    @test a * (A * B) == (a * A) * B
    @test (A * B) * a == A * (B * a)

    # Adjoints

    @test (A * B)' == B' * A'
end


text:
using JuliaDB, UUIDs
using CSV, Tables
import IntervalSets: width, (..), AbstractInterval, leftendpoint
function DungBase.AbstractTimeLine(x)
    files = Vector{VideoFile}(undef, length(x))
    for i in rows(x)
        files[i.index] = VideoFile(i.file_name, i.date_time, Millisecond(Nanosecond(i.duration)))
    end
    y = x[1]
    comment = ismissing(y.comment) ? "" : y.comment
    AbstractTimeLine(files, comment)
end
DungBase.AbstractPeriod(x, y::Missing) = AbstractPeriod(x)
function DungBase.Temporal(xs)
    x = xs[]
    video = x.AbstractTimeLine
    p1 = Millisecond(Nanosecond(x.start))
    p = ismissing(x.stop) ? p1 : p1..Millisecond(Nanosecond(x.stop))
    time = AbstractPeriod(p)
    try
        Temporal(video, time, ismissing(x.comment) ? "" : x.comment)
    catch
        Temporal(video, Prolonged(), ismissing(x.comment) ? "the time is made up because it was missing" : x.comment*" the time is made up because it was missing")
    end
end
function DungBase.Calibration(_x)
    x = _x[]
    board = Board(x.board, x.checker_width_cm, (x.checker_per_width, x.checker_per_height), x.board_description)
    comment = ismissing(x.comment) ? "" : x.comment
    Calibration(x.intrinsic, x.extrinsic, board, comment)
end
DungBase.POI(x) = POI(x.Calibration, x.Temporal)
function DungBase.Run(x)
    @assert allunique(x.type) "POI types are not unique: $(x.type)"
    pois = Dict(Symbol(i.type) => i.POI for i in x)
    tmp = select(table(x), Not(:date, :comment, :POI, :type))
    setup = Dict(k => v for (k, v) in pairs(tmp[1]) if !isempty(v))
    _x = x[1]
    date = _x.date
    comment = ismissing(_x.comment) ? "" : _x.comment
    metadata = Metadata(setup, comment, date)
    Run(pois, metadata)
end

function loaddeomcsv()

    video = loadtable(datadep"coffeebeetle/video.csv", indexcols = :video)
    videofile = loadtable(datadep"coffeebeetle/videofile.csv", indexcols = :file_name)
    x = join(video, videofile, rkey = :video)
    timeline = groupby(AbstractTimeLine,  x)
    interval = dropmissing(loadtable(datadep"coffeebeetle/interval.csv", indexcols = :interval), :start)
    x = join(interval, timeline, lkey = :video)
    temporal = dropmissing(groupby(Temporal, x, :interval), :Temporal)

    board = loadtable(datadep"coffeebeetle/board.csv", indexcols = :designation)
    calibration = loadtable(datadep"coffeebeetle/calibration.csv", indexcols = :calibration)
    x = join(calibration, board, lkey = :board, rkey = :designation)
    x = join(x, temporal, rkey = :interval, lkey = :extrinsic)
    x = select(x, Not(:extrinsic))
    x = JuliaDB.rename(x, :Temporal => :extrinsic)
    x = join(x, temporal, rkey = :interval, lkey = :intrinsic, how = :left)
    x = select(x, Not(:intrinsic))
    x = JuliaDB.rename(x, :Temporal => :intrinsic)
    calibration = groupby(Calibration, x, :calibration)

    poi = filter(!isempty, loadtable(datadep"coffeebeetle/poi.csv", indexcols = :poi), select = :calibration)
    coord = first.(splitext.(readdir(datadep"coffeebeetle/pixel")))
    poi = filter(in(coord), poi, select = :interval)
    x = join(poi, temporal, lkey = :interval)
    x = join(x, calibration, lkey = :calibration)
    x = select(x, Not(:calibration))
    x = JuliaDB.transform(x, :POI => POI.(x))

    #=a = []
    b = []
    j = 0
    for i in x
        y = Pair(Symbol(i.type), i.POI)
        if y∈ a
            global j
            j += 1
            push!(b, y)
        else
            push!(a,  y)
        end
    end=#

    temporal2pixel = Dict(Pair(Symbol(i.type), i.POI) => joinpath(datadep"coffeebeetle/pixel", string(i.interval, ".csv")) for i in x)

    JLSO.save(joinpath(datadep"coffeebeetle", "temporal2pixel.jlso"), temporal2pixel)
    poi = select(x, Not(:Temporal, :Calibration, :poi, :interval))

    runs = loadtable(datadep"coffeebeetle/run.csv", indexcols = :run)
    x = JuliaDB.transform(runs, :date => :date => Date)
    x = join(x, poi, rkey = :run)
    experiment = loadtable(datadep"coffeebeetle/experiment.csv", indexcols = :experiment)
    x = join(x, experiment, lkey = :experiment)
    data = groupby(x, :experiment) do xx
        runs = groupby(Run, select(table(xx), Not(:experiment_folder, :experiment_description)), :run)
        (Experiment = Experiment(select(runs, :Run), xx[1].experiment_description), name = xx[1].experiment_folder)
    end
    data = Dict(i.name => i.Experiment for i in data)

    JLSO.save(joinpath(datadep"coffeebeetle", "data.jlso"), data)
end


text:
struct RBDFAST <: GSAMethod
    num_harmonics::Int
end
RBDFAST(;num_harmonics = 6) = RBDFAST(num_harmonics)

"""
Code based on the theory presented in:
    Saltelli, A. (2008). Global sensitivity analysis: The primer. Chichester: Wiley, pp. 167-169.
and
    S. Tarantola, D. Gatelli and T. Mara (2006)
    "Random Balance Designs for the Estimation of First Order Global Sensitivity Indices",
    Reliability Engineering and System Safety, 91:6, 717-727
"""

using FFTW, Random, Statistics, StatsBase, Distributions
allsame(x) = all(y -> y == first(x), x)

function gsa(f, method::RBDFAST; num_params, N, rng::AbstractRNG = Random.default_rng(), batch = false, kwargs...)
    # Initalize matrix containing range of values of the parametric variable
    # along each column (factor).
    s0 = range(-π, stop = π, length = N)

    # Compute inputs
    s = [s0[randperm(rng, N)] for i in 1:num_params]
    x = hcat([0.5.+asin.(sin.(s[i]))./pi for i in 1:num_params]...)

    # Compute outputs
    if batch
        Y = f(x')
    else
        Y = [f(@view x[i,:]) for i in axes(x, 1)]
    end
    # Iterate over factors

    sensitivites = zeros(num_params)
    for i in 1:num_params

        s_order = sortperm(s[i])
        # Order Ys by how they would occur if they were
        # monotonically increasing as the
        # parametric variable s (not its permutation) increased.
        y_reordered = @view Y[s_order]

        ft = fft(y_reordered)
        ys = abs2.(ft) ./ N
        V = sum(ys[2:N])
        Vi = 2*sum(ys[2:method.num_harmonics+1])
        Si = Vi/V
        # println(ys)
        # unskew the sensitivies
        # lambda = 2*method.num_harmonics/N
        # sensitivites[i] = Si - (lambda / (1 - lambda)) * (1-Si)
        sensitivites[i] = Si
    end

    return sensitivites
end


text:
function getPath(G, idStart, idEnd)
	path = Int[]
	temp = idEnd
	while (temp != idStart)
		path = [temp; path]
		temp = G.V[temp].pred
	end
	return [idStart; path]
end

function DFS!(G, pocetak, kraj)
	v = 1:size(G.AdjMatrix, 1);
	for u = v
		G.V[u].color = 'W'
		G.V[u].pred = -1
	end
	DFS_Visit!(G, pocetak, kraj, pocetak)
end

function DFS_Visit!(G, u, kraj, pocetak)
        if u==kraj
            path = getPath(G, pocetak, kraj)
            println("Put: $path")
        else
   	   G.V[u].color = 'G';
           for v in findall(G.AdjMatrix[u,:] .==1)
		if G.V[v].color == 'W'
                        G.V[v].pred = u
			DFS_Visit!(G, v, kraj, pocetak)
		end
	   end
	   G.V[u].color = 'W';
        end
end


text:
using Flux
using Transformers
using Transformers.Basic

struct Translator
    transformers::TransformersModel
    beam_search::BeamSearch
end

Flux.@functor Translator
Flux.trainable(translator::Translator) = (translator.transformers,)

function (translator::Translator)(
    inputs::AbstractVector{T},
    vocabulary::Vocabulary{T};
    start_token::T="[CLS]",
    end_token::T="[SEP]"
)::AbstractVector{T} where T
    predict(outputs) = translator.transformers(vocabulary, inputs, outputs)
    return translator.beam_search(
        vocabulary,
        predict,
        start_token=start_token,
        end_token=end_token,
    )
end

function Translator(
    transformers::TransformersModel,
    beam_width::Integer,
    max_length::Integer;
    length_normalization::AbstractFloat=0.0
)
    return Translator(
        transformers,
        BeamSearch(beam_width, max_length, length_normalization)
    )
end


text:
using HydrogenDimerMRCC
using Test

@testset "HydrogenDimerMRCC" begin
	@testset "Invalid Inputs" begin
		# negative bond length
		@test_throws ErrorException p = PolarCoords(-1.0, pi/2, pi/2)
		# zero bond length
		@test_throws ErrorException p = PolarCoords(0.0, pi/2, pi/4)
		# polar angle below bounds
		@test_throws ErrorException p = PolarCoords(1.0, -0.5*pi, pi/2)
		# polar angle above bounds
		@test_throws ErrorException p = PolarCoords(1.0, 1.5*pi, pi/2)
		# azimuthal angle below bounds
		@test_throws ErrorException p = PolarCoords(1.0, pi/2, -0.5*pi)
		# azimuthal angle above bounds
		@test_throws ErrorException p = PolarCoords(1.0, pi/2, 2.5*pi)
		# invalid distance directly into DimerParameters
		@test_throws ErrorException d = DimerParameters(-5.0, 1.0, pi/2, pi/2, 1.0, pi/2, pi/2)
		@test_throws ErrorException d = DimerParameters(0.0, 1.0, pi/2, pi/2, 1.0, pi/2, pi/2)
	end

	@testset "Ideal Case 1" begin
		import HydrogenDimerMRCC: atom_positions, Cartesian3D, distance

		# rotate dimer (at origin) of length 2 onto x-y plane
		# one atom should be at  ( 1, 0, 0)
		# the other should be at (-1, 0, 0)
		polar = PolarCoords(2.0, pi/2, 0.0)
		com   = Cartesian3D(0.0, 0.0, 0.0)
		atom1, atom2 = atom_positions(polar, com)

		ideal1 = Cartesian3D( 1.0, 0.0, 0.0)
		ideal2 = Cartesian3D(-1.0, 0.0, 0.0)

		tolerance = 1.0e-15
		@test distance(atom1, ideal1) < tolerance
		@test distance(atom2, ideal2) < tolerance
	end

	@testset "Ideal Case 2" begin
		import HydrogenDimerMRCC: atom_positions, Cartesian3D, distance

		# rotate dimer (at origin) of length 2 onto x-y plane
		# then rotate 45 degrees in x-y place
		# one atom should be at  ( 1/sqrt(2),  1/sqrt(2), 0)
		# the other should be at (-1/sqrt(2), -1/sqrt(2), 0)
		polar = PolarCoords(2.0, pi/2, pi/4)
		com   = Cartesian3D(0.0, 0.0, 0.0)
		atom1, atom2 = atom_positions(polar, com)

		ideal1 = Cartesian3D( 1.0/sqrt(2.0),  1.0/sqrt(2.0), 0.0)
		ideal2 = Cartesian3D(-1.0/sqrt(2.0), -1.0/sqrt(2.0), 0.0)

		tolerance = 1.0e-15
		@test distance(atom1, ideal1) < tolerance
		@test distance(atom2, ideal2) < tolerance
	end
end