text:
"""
$(TYPEDEF)

A floating joint type that uses a modified Rodrigues parameter (`Rotations.MRP`,
previously known as `Rotations.SPQuat`) representation for orientation.

Floating joints are 6-degree-of-freedom joints that are in a sense degenerate,
as they impose no constraints on the relative motion between two bodies.

The 6-dimensional configuration vector of a `SPQuatFloating` joint
type consists of a SPQuat representing the orientation that rotates
vectors from the frame 'directly after' the joint to the frame 'directly before'
it, and a 3D position vector representing the origin of the frame after the
joint in the frame before the joint.

The 6-dimensional velocity vector of a `SPQuatFloating` joint is the twist
of the frame after the joint with respect to the frame before it, expressed in
the frame after the joint.
"""
struct SPQuatFloating{T} <: JointType{T} end

Base.show(io::IO, jt::SPQuatFloating) = print(io, "SPQuat floating joint")
Random.rand(::Type{SPQuatFloating{T}}) where {T} = SPQuatFloating{T}()

num_positions(::Type{<:SPQuatFloating}) = 6
num_velocities(::Type{<:SPQuatFloating}) = 6
has_fixed_subspaces(jt::SPQuatFloating) = true
isfloating(::Type{<:SPQuatFloating}) = true

@propagate_inbounds function rotation(jt::SPQuatFloating, q::AbstractVector)
    ModifiedRodriguesParam(q[1], q[2], q[3])
end

@propagate_inbounds function set_rotation!(q::AbstractVector, jt::SPQuatFloating, rot::Rotation{3})
    T = eltype(rot)
    spq = convert(ModifiedRodriguesParam{T}, rot)
    q[1] = spq.x
    q[2] = spq.y
    q[3] = spq.z
    nothing
end

@propagate_inbounds function set_rotation!(q::AbstractVector, jt::SPQuatFloating, rot::AbstractVector)
    q[1] = rot[1]
    q[2] = rot[2]
    q[3] = rot[3]
    nothing
end

@propagate_inbounds translation(jt::SPQuatFloating, q::AbstractVector) = SVector(q[4], q[5], q[6])
@propagate_inbounds set_translation!(q::AbstractVector, jt::SPQuatFloating, trans::AbstractVector) = copyto!(q, 4, trans, 1, 3)

@propagate_inbounds angular_velocity(jt::SPQuatFloating, v::AbstractVector) = SVector(v[1], v[2], v[3])
@propagate_inbounds set_angular_velocity!(v::AbstractVector, jt::SPQuatFloating, ω::AbstractVector) = copyto!(v, 1, ω, 1, 3)

@propagate_inbounds linear_velocity(jt::SPQuatFloating, v::AbstractVector) = SVector(v[4], v[5], v[6])
@propagate_inbounds set_linear_velocity!(v::AbstractVector, jt::SPQuatFloating, ν::AbstractVector) = copyto!(v, 4, ν, 1, 3)

@inline function set_configuration!(q::AbstractVector, joint::Joint{<:Any, <:SPQuatFloating}, config::Transform3D)
    @boundscheck check_num_positions(joint, q)
    @framecheck config.from frame_after(joint)
    @framecheck config.to frame_before(joint)
    @inbounds set_rotation!(q, joint_type(joint), rotation(config))
    @inbounds set_translation!(q, joint_type(joint), translation(config))
    q
end

@propagate_inbounds function set_velocity!(v::AbstractVector, joint::Joint{<:Any, <:SPQuatFloating}, twist::Twist)
    @boundscheck check_num_velocities(joint, v)
    @framecheck twist.base frame_before(joint)
    @framecheck twist.body frame_after(joint)
    @framecheck twist.frame frame_after(joint)
    @inbounds set_angular_velocity!(v, joint_type(joint), angular(twist))
    @inbounds set_linear_velocity!(v, joint_type(joint), linear(twist))
    v
end

@propagate_inbounds function joint_transform(jt::SPQuatFloating, frame_after::CartesianFrame3D, frame_before::CartesianFrame3D,
        q::AbstractVector)
    Transform3D(frame_after, frame_before, rotation(jt, q), translation(jt, q))
end

@inline function motion_subspace(jt::SPQuatFloating, frame_after::CartesianFrame3D, frame_before::CartesianFrame3D,
        q::AbstractVector)
    S = promote_eltype(jt, q)
    angular = hcat(one(SMatrix{3, 3, S}), zero(SMatrix{3, 3, S}))
    linear = hcat(zero(SMatrix{3, 3, S}), one(SMatrix{3, 3, S}))
    GeometricJacobian(frame_after, frame_before, frame_after, angular, linear)
end

@inline function constraint_wrench_subspace(jt::SPQuatFloating, joint_transform::Transform3D)
    S = promote_eltype(jt, joint_transform)
    WrenchMatrix(joint_transform.from, zero(SMatrix{3, 0, S}), zero(SMatrix{3, 0, S}))
end

@propagate_inbounds function bias_acceleration(jt::SPQuatFloating, frame_after::CartesianFrame3D, frame_before::CartesianFrame3D,
        q::AbstractVector, v::AbstractVector)
    S = promote_eltype(jt, q, v)
    zero(SpatialAcceleration{S}, frame_after, frame_before, frame_after)
end

@propagate_inbounds function configuration_derivative_to_velocity!(v::AbstractVector, jt::SPQuatFloating,
        q::AbstractVector, q̇::AbstractVector)
    spq = rotation(jt, q)
    spqdot = SVector(q̇[1], q̇[2], q̇[3])
    ω = angular_velocity_in_body(spq, spqdot)
    posdot = translation(jt, q̇)
    linear = inv(spq) * posdot
    set_angular_velocity!(v, jt, ω)
    set_linear_velocity!(v, jt, linear)
    nothing
end

@propagate_inbounds function configuration_derivative_to_velocity_adjoint!(fq, jt::SPQuatFloating,
        q::AbstractVector, fv)
    spq = ModifiedRodriguesParam(q[1], q[2], q[3])
    rot = velocity_jacobian(angular_velocity_in_body, spq)' * angular_velocity(jt, fv)
    trans = spq * linear_velocity(jt, fv)
    set_rotation!(fq, jt, rot)
    set_translation!(fq, jt, trans)
    nothing
end

@propagate_inbounds function velocity_to_configuration_derivative!(q̇::AbstractVector, jt::SPQuatFloating,
        q::AbstractVector, v::AbstractVector)
    spq = rotation(jt, q)
    ω = angular_velocity(jt, v)
    linear = linear_velocity(jt, v)
    spqdot = spquat_derivative(spq, ω)
    transdot = spq * linear
    set_rotation!(q̇, jt, spqdot)
    set_translation!(q̇, jt, transdot)
    nothing
end

@propagate_inbounds function velocity_to_configuration_derivative_jacobian(jt::SPQuatFloating,
        q::AbstractVector)
    spq = rotation(jt, q)
    vj = velocity_jacobian(spquat_derivative, spq)
    R = RotMatrix(spq)
    # TODO: use hvcat once it's as fast
    @SMatrix(
        [vj[1] vj[4] vj[7]  0    0    0;
         vj[2] vj[5] vj[8]  0    0    0;
         vj[3] vj[6] vj[9]  0    0    0;
         0     0     0      R[1] R[4] R[7];
         0     0     0      R[2] R[5] R[8];
         0     0     0      R[3] R[6] R[9]])
end

@propagate_inbounds function configuration_derivative_to_velocity_jacobian(jt::SPQuatFloating,
        q::AbstractVector)
    spq = rotation(jt, q)
    vj = velocity_jacobian(angular_velocity_in_body, spq)
    R_inv = RotMatrix(inv(spq))
    # TODO: use hvcat once it's as fast
    @SMatrix(
        [vj[1] vj[4] vj[7] 0        0        0;
         vj[2] vj[5] vj[8] 0        0        0;
         vj[3] vj[6] vj[9] 0        0        0;
         0     0     0     R_inv[1] R_inv[4] R_inv[7];
         0     0     0     R_inv[2] R_inv[5] R_inv[8];
         0     0     0     R_inv[3] R_inv[6] R_inv[9]])
end

@propagate_inbounds function zero_configuration!(q::AbstractVector, jt::SPQuatFloating)
    T = eltype(q)
    set_rotation!(q, jt, one(ModifiedRodriguesParam{T}))
    set_translation!(q, jt, zero(SVector{3, T}))
    nothing
end

@propagate_inbounds function rand_configuration!(q::AbstractVector, jt::SPQuatFloating)
    T = eltype(q)
    set_rotation!(q, jt, rand(ModifiedRodriguesParam{T}))
    set_translation!(q, jt, rand(SVector{3, T}) .- 0.5)
    nothing
end

@propagate_inbounds function joint_twist(jt::SPQuatFloating, frame_after::CartesianFrame3D, frame_before::CartesianFrame3D,
        q::AbstractVector, v::AbstractVector)
    S = promote_eltype(jt, q, v)
    angular = convert(SVector{3, S}, angular_velocity(jt, v))
    linear = convert(SVector{3, S}, linear_velocity(jt, v))
    Twist(frame_after, frame_before, frame_after, angular, linear)
end

@propagate_inbounds function joint_spatial_acceleration(jt::SPQuatFloating, frame_after::CartesianFrame3D, frame_before::CartesianFrame3D,
        q::AbstractVector, v::AbstractVector, vd::AbstractVector)
    S = promote_eltype(jt, q, v, vd)
    angular = convert(SVector{3, S}, angular_velocity(jt, vd))
    linear = convert(SVector{3, S}, linear_velocity(jt, vd))
    SpatialAcceleration(frame_after, frame_before, frame_after, angular, linear)
end

@propagate_inbounds function joint_torque!(τ::AbstractVector, jt::SPQuatFloating, q::AbstractVector, joint_wrench::Wrench)
    set_angular_velocity!(τ, jt, angular(joint_wrench))
    set_linear_velocity!(τ, jt, linear(joint_wrench))
    nothing
end

@propagate_inbounds function principal_value!(q::AbstractVector, jt::SPQuatFloating)
    spq = rotation(jt, q)
    set_rotation!(q, jt, principal_value(spq))
    nothing
end


text:
mutable struct _ConflictRefinerData
    confstat::Cint
    rowind::Vector{Cint}
    rowbdstat::Vector{Cint}
    confnumrows::Cint
    colind::Vector{Cint}
    colbdstat::Vector{Cint}
    confnumcols::Cint
end

function MOI.compute_conflict!(model::Optimizer)
    # This function always calls refineconflict first, which starts the conflict
    # refiner. In other words, any call to this function is expensive.
    confnumrows_p = Ref{Cint}(0)
    confnumcols_p = Ref{Cint}(0)
    ret = CPXrefineconflict(model.env, model.lp, confnumrows_p, confnumcols_p)
    _check_ret(model, ret)
    # Then, retrieve it.
    confstat_p = Ref{Cint}()
    rowind = Vector{Cint}(undef, confnumrows_p[])
    rowbdstat = Vector{Cint}(undef, confnumrows_p[])
    confnumrows_p = Ref{Cint}()
    colind = Vector{Cint}(undef, confnumcols_p[])
    colbdstat = Vector{Cint}(undef, confnumcols_p[])
    confnumcols_p  = Ref{Cint}()
    ret = CPXgetconflict(
        model.env,
        model.lp,
        confstat_p,
        rowind,
        rowbdstat,
        confnumrows_p,
        colind,
        colbdstat,
        confnumcols_p,
    )
    if ret == CPXERR_NO_CONFLICT
        model.conflict = _ConflictRefinerData(
            CPX_STAT_CONFLICT_FEASIBLE, Cint[], Cint[], 0, Cint[], Cint[], 0
        )
    else
        _check_ret(model, ret)
        model.conflict = _ConflictRefinerData(
            confstat_p[],
            rowind,
            rowbdstat,
            confnumrows_p[],
            colind,
            colbdstat,
            confnumcols_p[],
        )
    end
    return
end

function _ensure_conflict_computed(model::Optimizer)
    if model.conflict === nothing
        error(
            "Cannot access conflict status. Call " *
            "`MOI.compute_conflict!(model)` first. In case the model is " *
            "modified, the computed conflict will not be purged."
        )
    end
end

"""
    ConflictStatus()

Return the raw status from CPLEX indicating the status of the last computed
conflict.

It returns an integer value defined in
the [CPLEX documentation](https://www.ibm.com/support/knowledgecenter/SSSA5P_12.10.0/ilog.odms.cplex.help/refcallablelibrary/macros/homepagesolutionstatus.html)
(with a name like `CPX_STAT_CONFLICT_*`) or `nothing` if the function
`MOI.compute_conflict!` has not yet been called.
"""
struct ConflictStatus <: MOI.AbstractModelAttribute end

MOI.supports(::Optimizer, ::ConflictStatus) = true

function MOI.get(model::Optimizer, ::ConflictStatus)
    return model.conflict === nothing ? nothing : model.conflict.confstat
end

MOI.supports(::Optimizer, ::MOI.ConflictStatus) = true

function MOI.get(model::Optimizer, ::MOI.ConflictStatus)
    status = MOI.get(model, ConflictStatus())
    if status === nothing
        return MOI.COMPUTE_CONFLICT_NOT_CALLED
    elseif status == CPX_STAT_CONFLICT_MINIMAL
        return MOI.CONFLICT_FOUND
    elseif status == CPX_STAT_CONFLICT_FEASIBLE
        return MOI.NO_CONFLICT_EXISTS
    else
        return MOI.NO_CONFLICT_FOUND
    end
end

function _get_conflict_status(
    model::Optimizer, index::MOI.ConstraintIndex{MOI.SingleVariable, <:Any}
)
    _ensure_conflict_computed(model)
    cindex = findfirst(
        x -> x == Cint(_info(model, index).column - 1), model.conflict.colind
    )
    return cindex === nothing ? nothing : model.conflict.colbdstat[cindex]
end

function MOI.get(
    model::Optimizer,
    ::MOI.ConstraintConflictStatus,
    index::MOI.ConstraintIndex{MOI.SingleVariable, <:MOI.LessThan},
)
    status = _get_conflict_status(model, index)
    if status in (CPX_CONFLICT_MEMBER, CPX_CONFLICT_UB)
        return MOI.IN_CONFLICT
    elseif status in (CPX_CONFLICT_POSSIBLE_MEMBER, CPX_CONFLICT_POSSIBLE_UB)
        return MOI.MAYBE_IN_CONFLICT
    else
        return MOI.NOT_IN_CONFLICT
    end
end

function MOI.get(
    model::Optimizer,
    ::MOI.ConstraintConflictStatus,
    index::MOI.ConstraintIndex{MOI.SingleVariable, <:MOI.GreaterThan},
)
    status = _get_conflict_status(model, index)
    if status in (CPX_CONFLICT_MEMBER, CPX_CONFLICT_LB)
        return MOI.IN_CONFLICT
    elseif status in (CPX_CONFLICT_POSSIBLE_MEMBER, CPX_CONFLICT_POSSIBLE_LB)
        return MOI.MAYBE_IN_CONFLICT
    else
        return MOI.NOT_IN_CONFLICT
    end
end

function MOI.get(
    model::Optimizer,
    ::MOI.ConstraintConflictStatus,
    index::MOI.ConstraintIndex{
        MOI.SingleVariable, <:Union{MOI.EqualTo, MOI.Interval}
    }
)
    status = _get_conflict_status(model, index)
    if status in (CPX_CONFLICT_MEMBER, CPX_CONFLICT_LB, CPX_CONFLICT_UB)
        return MOI.IN_CONFLICT
    elseif status in (CPX_CONFLICT_POSSIBLE_MEMBER, CPX_CONFLICT_POSSIBLE_LB, CPX_CONFLICT_POSSIBLE_UB)
        return MOI.MAYBE_IN_CONFLICT
    else
        return MOI.NOT_IN_CONFLICT
    end
end

function MOI.get(
    model::Optimizer,
    ::MOI.ConstraintConflictStatus,
    index::MOI.ConstraintIndex{
        <:MOI.ScalarAffineFunction,
        <:Union{MOI.LessThan, MOI.GreaterThan, MOI.EqualTo}
    },
)
    _ensure_conflict_computed(model)
    rindex = findfirst(
        x -> x == Cint(_info(model, index).row - 1), model.conflict.rowind
    )
    if rindex === nothing
        return MOI.NOT_IN_CONFLICT
    elseif model.conflict.rowbdstat[rindex] == CPX_CONFLICT_MEMBER
        return MOI.IN_CONFLICT
    else
        return MOI.MAYBE_IN_CONFLICT
    end
end

function MOI.supports(
    ::Optimizer,
    ::MOI.ConstraintConflictStatus,
    ::Type{MOI.ConstraintIndex{<:MOI.SingleVariable, <:_SCALAR_SETS}},
)
    return true
end

function MOI.supports(
    ::Optimizer,
    ::MOI.ConstraintConflictStatus,
    ::Type{MOI.ConstraintIndex{
        <:MOI.ScalarAffineFunction,
        <:Union{MOI.LessThan, MOI.GreaterThan, MOI.EqualTo}
    }},
)
    return true
end

ext: jl


text:
module Opt_PID_work

#
# The main difference between this and the InfDecomp module is:
#  * here, we maximize H(X|YZ) instead of minimizing -H(X|YZ)
#

using Opt_PID_data: My_Eval, condEntropy, ∇f
# using MathProgBase:
#     initialize, features_available,
#     isobjlinear, isobjquadratic, isconstrlinear,
#     eval_f, eval_g, eval_grad_f,
#     jac_structure, eval_jac_g, eval_jac_prod, eval_jac_prod_t,
#     hesslag_structure, eval_hesslag #, eval_hesslag_prod
using MathProgBase

function MathProgBase.initialize(e::My_Eval, requested_features::Vector{Symbol})
    print("initialize")
    for feat in requested_features
        if feat ∉ [:Grad,:Jac,:JacVec,:Hess]
            error("Opt_PID_work.initialize():\n-   JuliaOpt:MathProgBase is asking for a feature ($feat) that I don't have.-   Maybe use another solver?")
        end
    end
    print("END initialize")
    ;
end

# features_available()
MathProgBase.features_available(::My_Eval) = features_list

# Properties:
MathProgBase.isobjlinear(::My_Eval)           = false
MathProgBase.isobjquadratic(::My_Eval)        = false
MathProgBase.isconstrlinear(::My_Eval, ::Integer) = true

#
# ------------------------------------------
# E v a l u a t i o n :   0 t h   o r d e r
# ------------------------------------------
function MathProgBase.eval_f{TFloat_2,TFloat}(e::My_Eval, x::Vector{TFloat_2},dummy::TFloat=Float64(0)) :: TFloat
    # println("eval_f")
    local condent::Float64
    if e.TmpFloat==BigFloat
        prc = precision(BigFloat)
        setprecision(BigFloat,e.bigfloat_nbits)
        condent = condEntropy(e,x,BigFloat(0.))
        setprecision(BigFloat,prc)
    else
        condent = condEntropy(e,x,Float64(0))
    end
    # println("END eval_f")
    return condent
    ;
end #^ eval_f()

# eval_g --- eval of constraint into g
function MathProgBase.eval_g(e::My_Eval, g::Vector{Float64}, x::Vector{Float64})  :: Void
    # println("eval_g")
    g .= reshape( reshape(x,1,e.n)*e.Gt , e.m, ) .- e.rhs
    # println("END eval_g")
    return nothing
    ;
end # eval_g()

# eval_grad_f --- eval gradient of objective function
function MathProgBase.eval_grad_f(e::My_Eval, g::Vector{Float64}, x::Vector{Float64}) :: Void
    # println("eval_grad_f")
    if e.TmpFloat==BigFloat
        prc = precision(BigFloat)
        setprecision(BigFloat,e.bigfloat_nbits)
        ∇f(e,g,x,BigFloat(0.))
        setprecision(BigFloat,prc)
    else
        ∇f(e,g,x,Float64(0))
    end
    # println("END eval_grad_f")
    return nothing
    ;
end # eval_grad_f()

# Constraint Jacobian
# jac_structure() --- zero-nonzero pattern of constraint Jacobian
MathProgBase.jac_structure(e::My_Eval) :: Tuple{Vector{Int64},Vector{Int64}}   =  ( e.Gt_L , e.Gt_K )
# Note: Gt is transposed, so K and L are swapped [K: rows of G^T; L: columns of G^T]

# eval_jac_g() --- constraint Jacobian   -> J
function MathProgBase.eval_jac_g(e::My_Eval, J::Vector{Float64}, x::Vector{Float64}) :: Void
    J .= e.Gt.nzval
    return nothing
    ;
end # eval_jac_g()

# eval_jac_prod() --- constraint_Jacobian * w   -> y
function MathProgBase.eval_jac_prod(e::My_Eval, y::Vector{Float64}, x::Vector{Float64}, w::Vector{Float64}) :: Void
    y .= reshape( reshape(w,1,e.n)*e.Gt , e.m, )
    ;
end # eval_jac_prod()

# eval_jac_prod_t() --- constraint_Jacobian^T * w  -> y
function MathProgBase.eval_jac_prod_t{T<:AbstractFloat}(e::My_Eval, y::Vector{T}, x::Vector{T}, w::Vector{T}) :: Void
    y .= e.Gt*w
    ;
end

# ------------------------------------------
# E v a l u a t i o n :   2 n d   o r d e r
# ------------------------------------------

# Lagrangian:
# L(x, (σ,μ) ) = σ f(x) + μ' G(x)
# Since G(.) is linear, it's Hessian is 0 anyway, so it won't bother us.

# hesslag_structure() --- zero-nonzero pattern of Hessian [wrt x] of the Lagrangian
function MathProgBase.hesslag_structure(e::My_Eval)  :: Tuple{Vector{Int64},Vector{Int64}}
    K = Vector{Int64}()
    L = Vector{Int64}()
    counter = 0
    for y = 1:e.n_y
        for z = 1:e.n_z
            # Start with the diagonal
            for x = 1:e.n_x
                i = e.varidx[x,y,z]
                if i>0
                    counter += 1
                    push!(K,i)
                    push!(L,i)
                end
            end
            # Now off-diagonal.
            # H is  treated as a symmetric matrix, but:
            # if both (k,l) & (l,k) are present, their values will be added!
            for x = 1:e.n_x
                for u = 1:(x-1)
                    i_x = e.varidx[x,y,z]
                    i_u = e.varidx[u,y,z]
                    if i_x>0 && i_u>0
                        counter += 1
                        push!(K,i_x)
                        push!(L,i_u)
                    end
                end
            end
        end# for z
    end #^ for y
    return (K,L)
    ;
end # hesslag_structure()

function Hess{TFloat}(e::My_Eval, H::Vector{Float64}, p::Vector{Float64}, σ::Float64, dummy::TFloat) :: Void
    counter = 0
    for y = 1:e.n_y
        for z = 1:e.n_z
            # make marginal P(*yz)
            P_yz  ::TFloat = TFloat(0.)
            for x = 1:e.n_x
                i = e.varidx[x,y,z]
                if i>0
                    P_yz += p[i]
                end
            end

            # now: for all pairs x,u we have:
            # if x ≠  u:   1/P_yz;
            # if x == u:   -( P_yz - P(xyz) )/(  P_yz * P(xyz) )

            # Start with the diagonal
            for x = 1:e.n_x
                i = e.varidx[x,y,z]
                if i>0
                    counter += 1
                    P_xyz = p[i]
                    H[counter] = Float64(   (P_xyz == 0 ) ?  1.e50  : TFloat(σ)*( -P_yz + P_xyz )/(  P_yz * P_xyz )  ) # flipped sign p_yz - p_xyz
                end
            end
            # Now off-diagonal.
            # H is  treated as a symmetric matrix, but:
            # if both (k,l) & (l,k) are present, their values will be added!
            for x = 1:e.n_x
                for u = 1:(x-1)
                    i_x = e.varidx[x,y,z]
                    i_u = e.varidx[u,y,z]
                    if i_x>0 && i_u>0
                        counter += 1
                        H[counter] = TFloat(σ)/P_yz # flipped sign
                    end
                end
            end
        end#for z
    end# for y
    ;
end # eval_hesslag()

# eval_hesslag() --- Hessian [wrt x] of the Lagrangian
function MathProgBase.eval_hesslag(e::My_Eval, H::Vector{Float64}, x::Vector{Float64}, σ::Float64, μ::Vector{Float64}) :: Void
    if e.TmpFloat==BigFloat
        prc = precision(BigFloat)
        setprecision(BigFloat,e.bigfloat_nbits)
        Hess(e,H,x,σ,BigFloat(0.))
        setprecision(BigFloat,prc)
    else
        Hess(e,H,x,σ,Float64(0))
    end
    return nothing
    ;
end # eval_hesslag()

# eval_hesslag() --- ( Hessian [wrt x] of the Lagrangian ) * v
# function eval_hesslag_prod{T<:AbstractFloat}(e::My_Eval, h::Vector{T}, x::Vector{T}, v::Vector{T}, σ::T, μ::Vector{T}) :: Void
#     h .= σ .* Hf(x)*v
#     return nothing
# end

end #module Opt_PID_work


text:
type LRUNode{K, V}
    k::K
    v::V
    next::LRUNode{K, V}
    prev::LRUNode{K, V}

    # All new created nodes are self referential only
    function LRUNode{K, V}(k::K, v::V)
        x = new(k, v)
        x.next = x
        x.prev = x
        return x
    end
end

type LRUList{K, V}
    first::Nullable{LRUNode{K, V}}
    size::Int

    LRUList() = new(Nullable{LRUNode{K, V}}(), 0)
end

Base.first(l::LRUList) = !isempty(l) ? get(l.first) : error("LRUList is empty")
Base.last(l::LRUList) = !isempty(l) ? get(l.first).prev : error("LRUList is empty")

Base.length(l::LRUList) = l.size
Base.isempty(l::LRUList) = length(l) == 0

function Base.show{K, V}(io::IO, l::LRUNode{K, V})
    print(io, "LRUNode{", K, ", ", V, "}(")
    show(io, l.k)
    print(io, ", ")
    show(io, l.v)
    print(io, ")")
end

function Base.show{K, V}(io::IO, l::LRUList{K, V})
    print(io, "LRUList{", K, ", ", V, "}(")
    if length(l) != 0
        f = first(l)
        show(io, f.k)
        print(io, "=>")
        show(io, f.v)
        n = f.next
        while !is(n, f)
            print(io, ", ")
            show(io, n.k)
            print(io, "=>")
            show(io, n.v)
            n = n.next
        end
    end
    print(io, ")")
end

function Base.push!{K, V}(list::LRUNode{K, V}, new::LRUNode{K, V})
    new.next = list
    new.prev = list.prev
    list.prev.next = new
    list.prev = new
    return list
end

function Base.push!{K, V}(l::LRUList{K, V}, el::LRUNode{K, V})
    if isempty(l)
        l.first = Nullable(el)
    else
        push!(first(l), el)
    end
    l.size += 1
    return l
end

function Base.pop!{K, V}(l::LRUList{K, V}, n::LRUNode{K, V}=last(l))
    if is(n.next, n)
        l.first = Nullable{LRUNode{K, V}}()
    else
        if is(n, first(l))
            l.first = Nullable(n.next)
        end
        n.next.prev = n.prev
        n.prev.next = n.next
    end
    l.size -= 1
    return n
end

function Base.unshift!{K, V}(l::LRUList{K, V}, el::LRUNode{K, V})
    push!(l, el)
    rotate!(l)
end

# Rotate one step forward, so last element is now first
function rotate!(l::LRUList)
    if length(l) > 1
        l.first = Nullable(first(l).prev)
    end
    return l
end

# Move the node n to the front of the list
function move_to_front!{T}(l::LRUList{T}, n::LRUNode{T})
    if !is(first(l), n)
        pop!(l, n)
        unshift!(l, n)
    end
    return l
end

function Base.delete!{K, V}(l::LRUList{K, V}, n::LRUNode{K, V})
    pop!(l, n)
    return l
end

function Base.empty!{K, V}(l::LRUList{K, V})
    l.first = Nullable{LRUNode{K, V}}()
    l.size = 0
    return l
end


text:
# Rows grouping.
# Maps row contents to the indices of all the equal rows.
# Used by groupby(), join(), nonunique()
struct RowGroupDict{T<:AbstractDataFrame}
    "source data table"
    df::T
    "row hashes (optional, can be empty)"
    rhashes::Vector{UInt}
    "hashindex -> index of group-representative row (optional, can be empty)"
    gslots::Vector{Int}
    "group index for each row"
    groups::Vector{Int}
    "permutation of row indices that sorts them by groups"
    rperm::Vector{Int}
    "starts of ranges in rperm for each group"
    starts::Vector{Int}
    "stops of ranges in rperm for each group"
    stops::Vector{Int}
end

# "kernel" functions for hashrows()
# adjust row hashes by the hashes of column elements
function hashrows_col!(h::Vector{UInt},
                       n::Vector{Bool},
                       v::AbstractVector{T},
                       rp::Nothing,
                       firstcol::Bool) where T
    @inbounds for i in eachindex(h)
        el = v[i]
        h[i] = hash(el, h[i])
        if length(n) > 0
            n[i] |= ismissing(el)
        end
    end
    h
end

# should give the same hash as AbstractVector{T}
function hashrows_col!(h::Vector{UInt},
                       n::Vector{Bool},
                       v::AbstractVector,
                       rp::Any,
                       firstcol::Bool)
    # When hashing the first column, no need to take into account previous hash,
    # which is always zero
    # also when the number of values in the pool is more than half the length
    # of the vector avoid using this path. 50% is roughly based on benchmarks
    if firstcol && Int64(2) * length(rp) < length(v)
        hashes = Vector{UInt}(undef, length(rp))
        @inbounds for (i, v) in zip(eachindex(hashes), rp)
            hashes[i] = hash(v)
        end

        fi = firstindex(rp)
        # here we rely on the fact that `DataAPI.refpool` has a continuous
        # block of indices
        @inbounds for (i, ref) in enumerate(DataAPI.refarray(v))
            h[i] = hashes[ref+1-fi]
        end
    else
        @inbounds for (i, x) in enumerate(v)
            h[i] = hash(x, h[i])
        end
    end
    # Doing this step separately is faster, as it would disable SIMD above
    if eltype(v) >: Missing && length(n) > 0
        n .|= ismissing.(v)
    end
    h
end

# Calculate the vector of `df` rows hash values.
function hashrows(cols::Tuple{Vararg{AbstractVector}}, skipmissing::Bool)
    len = length(cols[1])
    rhashes = zeros(UInt, len)
    missings = fill(false, skipmissing ? len : 0)
    for (i, col) in enumerate(cols)
        rp = DataAPI.refpool(col)
        hashrows_col!(rhashes, missings, col, rp, i == 1)
    end
    return (rhashes, missings)
end

# table columns are passed as a tuple of vectors to ensure type specialization
isequal_row(cols::Tuple{AbstractVector}, r1::Int, r2::Int) =
    isequal(cols[1][r1], cols[1][r2])
isequal_row(cols::Tuple{Vararg{AbstractVector}}, r1::Int, r2::Int) =
    isequal(cols[1][r1], cols[1][r2]) && isequal_row(Base.tail(cols), r1, r2)

isequal_row(cols1::Tuple{AbstractVector}, r1::Int, cols2::Tuple{AbstractVector}, r2::Int) =
    isequal(cols1[1][r1], cols2[1][r2])
isequal_row(cols1::Tuple{Vararg{AbstractVector}}, r1::Int,
            cols2::Tuple{Vararg{AbstractVector}}, r2::Int) =
    isequal(cols1[1][r1], cols2[1][r2]) &&
        isequal_row(Base.tail(cols1), r1, Base.tail(cols2), r2)

# IntegerRefarray and IntegerRefPool are two complementary view types that allow
# wrapping arrays with Union{Real, Missing} eltype to satisfy the DataAPI.refpool
# and DataAPI.refarray API when calling row_group_slots.
# IntegerRefarray converts values to Int and replaces missing with an integer
# (set by the caller to the maximum value + 1)
# IntegerRefPool subtracts the minimum value - 1 and replaces back the maximum
# value + 1 to missing. This ensures all values are in 1:length(refpool), while
# row_group_slots knows the number of (potential) groups via length(refpool)
# and is able to skip missing values when skipmissing=true

struct IntegerRefarray{T<:AbstractArray} <: AbstractVector{Int}
    x::T
    offset::Int
    replacement::Int
end

Base.size(x::IntegerRefarray) = size(x.x)
Base.axes(x::IntegerRefarray) = axes(x.x)
Base.IndexStyle(::Type{<:IntegerRefarray{T}}) where {T} = Base.IndexStyle(T)
@inline function Base.getindex(x::IntegerRefarray, i)
    @boundscheck checkbounds(x.x, i)
    @inbounds v = x.x[i]
    if eltype(x.x) >: Missing && v === missing
        return x.replacement
    else
        # Overflow is guaranteed not to happen by checks before calling the constructor
        return Int(v - x.offset)
    end
end

struct IntegerRefpool{T<:Union{Int, Missing}} <: AbstractVector{T}
    max::Int
    function IntegerRefpool{T}(max::Integer) where T<:Union{Int, Missing}
        @assert max < typemax(Int) # to store missing values as max + 1
        new{T}(max)
    end
end

Base.size(x::IntegerRefpool{T}) where {T} = (x.max + (T >: Missing),)
Base.axes(x::IntegerRefpool{T}) where {T} = (Base.OneTo(x.max + (T >: Missing)),)
Base.IndexStyle(::Type{<:IntegerRefpool}) = Base.IndexLinear()
@inline function Base.getindex(x::IntegerRefpool{T}, i::Real) where T
    @boundscheck checkbounds(x, i)
    if T >: Missing && i == x.max + 1
        return missing
    else
        return Int(i)
    end
end
Base.allunique(::IntegerRefpool) = true
Base.issorted(::IntegerRefpool) = true

function refpool_and_array(x::AbstractArray)
    refpool = DataAPI.refpool(x)
    refarray = DataAPI.refarray(x)

    if refpool !== nothing
        return refpool, refarray
    elseif x isa AbstractArray{<:Union{Real, Missing}} &&
        all(v -> ismissing(v) | isinteger(v), x) &&
        !isempty(skipmissing(x))
        minval, maxval = extrema(skipmissing(x))
        ngroups = big(maxval) - big(minval) + 1
        # Threshold chosen with the same rationale as the row_group_slots refpool method:
        # refpool approach is faster but we should not allocate too much memory either
        # We also have to avoid overflow, including with ngroups + 1 for missing values
        # (note that it would be possible to allow minval and maxval to be outside of the
        # range supported by Int by adding a type parameter for minval to IntegerRefarray)
        if typemin(Int) < minval <= maxval < typemax(Int) &&
            ngroups + 1 <= Int64(2) * length(x) <= typemax(Int)
            T = eltype(x) >: Missing ? Union{Int, Missing} : Int
            refpool′ = IntegerRefpool{T}(Int(ngroups))
            refarray′ = IntegerRefarray(x, Int(minval) - 1, Int(ngroups) + 1)
            return refpool′, refarray′
        end
    end
    return nothing, nothing
end

# Helper function for RowGroupDict.
# Returns a tuple:
# 1) the highest group index in the `groups` vector
# 2) vector of row hashes (may be empty if hash=Val(false))
# 3) slot array for a hash map, non-zero values are
#    the indices of the first row in a group
# 4) whether groups are already sorted
# Optional `groups` vector is set to the group indices of each row (starting at 1)
# With skipmissing=true, rows with missing values are attributed index 0.
function row_group_slots(cols::Tuple{Vararg{AbstractVector}},
                         hash::Val = Val(true),
                         groups::Union{Vector{Int}, Nothing} = nothing,
                         skipmissing::Bool = false,
                         sort::Bool = false)::Tuple{Int, Vector{UInt}, Vector{Int}, Bool}
    rpa = refpool_and_array.(cols)
    refpools = first.(rpa)
    refarrays = last.(rpa)
    row_group_slots(cols, refpools, refarrays, hash, groups, skipmissing, sort)
end

# Generic fallback method based on open adressing hash table
function row_group_slots(cols::Tuple{Vararg{AbstractVector}},
                         refpools::Any,  # Ignored
                         refarrays::Any, # Ignored
                         hash::Val = Val(true),
                         groups::Union{Vector{Int}, Nothing} = nothing,
                         skipmissing::Bool = false,
                         sort::Bool = false)::Tuple{Int, Vector{UInt}, Vector{Int}, Bool}
    @assert groups === nothing || length(groups) == length(cols[1])
    rhashes, missings = hashrows(cols, skipmissing)
    # inspired by Dict code from base cf. https://github.com/JuliaData/DataTables.jl/pull/17#discussion_r102481481
    # but using open addressing with a table with at least 5/4 as many slots as rows
    # (rounded up to nearest power of 2) to avoid performance degradation
    # in a corner case of groups having exactly one row
    sz = max(1 + ((5 * length(rhashes)) >> 2), 16)
    sz = 1 << (8 * sizeof(sz) - leading_zeros(sz - 1))
    @assert 4 * sz >= 5 * length(rhashes)
    szm1 = sz-1
    gslots = zeros(Int, sz)
    # If missings are to be skipped, they will all go to group 0,
    # which will be removed by functions down the stream
    ngroups = 0
    @inbounds for i in eachindex(rhashes)
        # find the slot and group index for a row
        slotix = rhashes[i] & szm1 + 1
        # Use -1 for non-missing values to catch bugs if group is not found
        gix = skipmissing && missings[i] ? 0 : -1
        probe = 0
        if !skipmissing || !missings[i]
            while true
                g_row = gslots[slotix]
                if g_row == 0 # unoccupied slot, current row starts a new group
                    gslots[slotix] = i
                    gix = ngroups += 1
                    break
                elseif rhashes[i] == rhashes[g_row] # occupied slot, check if miss or hit
                    if isequal_row(cols, i, g_row) # hit
                        gix = groups !== nothing ? groups[g_row] : 0
                        break
                    end
                end
                slotix = slotix & szm1 + 1 # check the next slot
                probe += 1
                @assert probe < sz
            end
        end
        if groups !== nothing
            groups[i] = gix
        end
    end
    return ngroups, rhashes, gslots, false
end

# Optimized method for arrays for which DataAPI.refpool is defined and returns an AbstractVector
function row_group_slots(cols::NTuple{N, AbstractVector},
                         refpools::NTuple{N, AbstractVector},
                         refarrays::NTuple{N,
                             Union{AbstractVector{<:Real},
                                   Missings.EachReplaceMissing{
                                       <:AbstractVector{<:Union{Real, Missing}}}}},
                         hash::Val{false},
                         groups::Union{Vector{Int}, Nothing} = nothing,
                         skipmissing::Bool = false,
                         sort::Bool = false)::Tuple{Int, Vector{UInt}, Vector{Int}, Bool} where N
    # Computing neither hashes nor groups isn't very useful,
    # and this method needs to allocate a groups vector anyway
    @assert groups !== nothing && all(col -> length(col) == length(groups), cols)

    missinginds = map(refpools) do refpool
        eltype(refpool) >: Missing ?
            something(findfirst(ismissing, refpool), lastindex(refpool)+1) : lastindex(refpool)+1
    end

    # If skipmissing=true, rows with missings all go to group 0,
    # which will be removed by functions down the stream
    ngroupstup = map(refpools, missinginds) do refpool, missingind
        len = length(refpool)
        if skipmissing && missingind <= lastindex(refpool)
            return len - 1
        else
            return len
        end
    end
    ngroups = prod(ngroupstup)

    # Fall back to hashing if there would be too many empty combinations
    # or if the pool does not contain only unique values
    # The first check ensures the computation of ngroups did not overflow.
    # The rationale for the 2 threshold is that while the fallback method is always slower,
    # it allocates a hash table of size length(groups) instead of the remap vector
    # of size ngroups (i.e. the number of possible combinations) in this method:
    # so it makes sense to allocate more memory for better performance,
    # but it needs to remain reasonable compared with the size of the data frame.
    anydups = !all(allunique, refpools)
    if prod(big.(ngroupstup)) > typemax(Int) ||
       ngroups > Int64(2) * length(groups) ||
       anydups
        # In the simplest case, we can work directly with the reference codes
        newcols = (skipmissing && any(refpool -> eltype(refpool) >: Missing, refpools)) ||
                  !(refarrays isa NTuple{<:Any, AbstractVector}) ||
                  sort ||
                  anydups ? cols : refarrays
        return invoke(row_group_slots,
                      Tuple{Tuple{Vararg{AbstractVector}}, Any, Any, Val,
                            Union{Vector{Int}, Nothing}, Bool, Bool},
                      newcols, refpools, refarrays, hash, groups, skipmissing, sort)
    end

    seen = fill(false, ngroups)
    strides = (cumprod(collect(reverse(ngroupstup)))[end-1:-1:1]..., 1)::NTuple{N, Int}
    firstinds = map(firstindex, refpools)
    if sort
        nminds = map(refpools, missinginds) do refpool, missingind
            missingind > lastindex(refpool) ?
                eachindex(refpool) : setdiff(eachindex(refpool), missingind)
        end
        if skipmissing
            sorted = all(issorted(view(refpool, nmind))
                         for (refpool, nmind) in zip(refpools, nminds))
        else
            sorted = all(issorted, refpools)
        end
    else
        sorted = false
    end
    if sort && !sorted
        # Compute vector mapping missing to -1 if skipmissing=true
        refmaps = map(cols, refpools, missinginds, nminds) do col, refpool, missingind, nmind
            refmap = collect(0:length(refpool)-1)
            if skipmissing
                fi = firstindex(refpool)
                if missingind <= lastindex(refpool)
                    refmap[missingind-fi+1] = -1
                    refmap[missingind-fi+2:end] .-= 1
                end
                if sort
                    perm = sortperm(view(refpool, nmind))
                    invpermute!(view(refmap, nmind .- fi .+ 1), perm)
                end
            elseif sort
                # collect is needed for CategoricalRefPool
                invpermute!(refmap, sortperm(collect(refpool)))
            end
            refmap
        end
        tforeach(eachindex(groups), basesize=1_000_000) do i
            @inbounds begin
                local refs_i
                let i=i # Workaround for julia#15276
                    refs_i = map(c -> c[i], refarrays)
                end
                vals = map((m, r, s, fi) -> m[r-fi+1] * s, refmaps, refs_i, strides, firstinds)
                j = sum(vals) + 1
                # x < 0 happens with -1 in refmap, which corresponds to missing
                if skipmissing && any(x -> x < 0, vals)
                    j = 0
                else
                    seen[j] = true
                end
                groups[i] = j
            end
        end
    else
        tforeach(eachindex(groups), basesize=1_000_000) do i
            @inbounds begin
                local refs_i
                let i=i # Workaround for julia#15276
                    refs_i = map(refarrays, missinginds) do ref, missingind
                        r = Int(ref[i])
                        if skipmissing
                            return r == missingind ? -1 : (r > missingind ? r-1 : r)
                        else
                            return r
                        end
                    end
                end
                vals = map((r, s, fi) -> (r-fi) * s, refs_i, strides, firstinds)
                j = sum(vals) + 1
                # x < 0 happens with -1, which corresponds to missing
                if skipmissing && any(x -> x < 0, vals)
                    j = 0
                else
                    seen[j] = true
                end
                groups[i] = j
            end
        end
    end
    # If some groups are unused, compress group indices to drop them
    # sum(seen) is faster than all(seen) when not short-circuiting,
    # and short-circuit would only happen in the slower case anyway
    if sum(seen) < length(seen)
        oldngroups = ngroups
        remap = zeros(Int, ngroups)
        ngroups = 0
        @inbounds for i in eachindex(remap, seen)
            ngroups += seen[i]
            remap[i] = ngroups
        end
        @inbounds for i in eachindex(groups)
            gix = groups[i]
            groups[i] = gix > 0 ? remap[gix] : 0
        end
        # To catch potential bugs inducing unnecessary computations
        @assert oldngroups != ngroups
    end
    return ngroups, UInt[], Int[], sort
end

# Return a 3-tuple of a permutation that sorts rows into groups,
# and the positions of the first and last rows in each group in that permutation
# `groups` must contain group indices in 0:ngroups
# Rows with group index 0 are skipped (used when skipmissing=true)
# Partly uses the code of Wes McKinney's groupsort_indexer in pandas (file: src/groupby.pyx).
function compute_indices(groups::AbstractVector{<:Integer}, ngroups::Integer)
    # count elements in each group
    stops = zeros(Int, ngroups+1)
    @inbounds for gix in groups
        stops[gix+1] += 1
    end

    # group start positions in a sorted table
    starts = Vector{Int}(undef, ngroups+1)
    if length(starts) > 0
        starts[1] = 1
        @inbounds for i in 1:ngroups
            starts[i+1] = starts[i] + stops[i]
        end
    end

    # define row permutation that sorts them into groups
    rperm = Vector{Int}(undef, length(groups))
    copyto!(stops, starts)
    @inbounds for (i, gix) in enumerate(groups)
        rperm[stops[gix+1]] = i
        stops[gix+1] += 1
    end
    stops .-= 1

    # When skipmissing=true was used, group 0 corresponds to missings to drop
    # Otherwise it's empty
    popfirst!(starts)
    popfirst!(stops)

    return rperm, starts, stops
end

# Build RowGroupDict for a given DataFrame, using all of its columns as grouping keys
function group_rows(df::AbstractDataFrame)
    groups = Vector{Int}(undef, nrow(df))
    ngroups, rhashes, gslots, sorted =
        row_group_slots(ntuple(i -> df[!, i], ncol(df)), Val(true), groups, false, false)
    rperm, starts, stops = compute_indices(groups, ngroups)
    return RowGroupDict(df, rhashes, gslots, groups, rperm, starts, stops)
end

# Find index of a row in gd that matches given row by content, 0 if not found
function findrow(gd::RowGroupDict,
                 df::AbstractDataFrame,
                 gd_cols::Tuple{Vararg{AbstractVector}},
                 df_cols::Tuple{Vararg{AbstractVector}},
                 row::Int)
    (gd.df === df) && return row # same table, return itself
    # different tables, content matching required
    rhash = rowhash(df_cols, row)
    szm1 = length(gd.gslots)-1
    slotix = ini_slotix = rhash & szm1 + 1
    while true
        g_row = gd.gslots[slotix]
        if g_row == 0 || # not found
            (rhash == gd.rhashes[g_row] &&
            isequal_row(gd_cols, g_row, df_cols, row)) # found
            return g_row
        end
        slotix = (slotix & szm1) + 1 # miss, try the next slot
        (slotix == ini_slotix) && break
    end
    return 0 # not found
end

# Find indices of rows in 'gd' that match given row by content.
# return empty set if no row matches
function findrows(gd::RowGroupDict,
                  df::AbstractDataFrame,
                  gd_cols::Tuple{Vararg{AbstractVector}},
                  df_cols::Tuple{Vararg{AbstractVector}},
                  row::Int)
    g_row = findrow(gd, df, gd_cols, df_cols, row)
    (g_row == 0) && return view(gd.rperm, 0:-1)
    gix = gd.groups[g_row]
    return view(gd.rperm, gd.starts[gix]:gd.stops[gix])
end

function Base.getindex(gd::RowGroupDict, dfr::DataFrameRow)
    g_row = findrow(gd, parent(dfr), ntuple(i -> gd.df[!, i], ncol(gd.df)),
                    ntuple(i -> parent(dfr)[!, i], ncol(parent(dfr))), row(dfr))
    (g_row == 0) && throw(KeyError(dfr))
    gix = gd.groups[g_row]
    return view(gd.rperm, gd.starts[gix]:gd.stops[gix])
end


text:
@testset "Degree Independent Set" begin

    g0 = SimpleGraph(0)
    for g in testgraphs(g0)
        c = @inferred(independent_set(g, DegreeIndependentSet()))
        @test isempty(c)
    end

    g1 = SimpleGraph(1)
    for g in testgraphs(g1)
        c = @inferred(independent_set(g, DegreeIndependentSet()))
        @test (c == [1,])
    end

    add_edge!(g1, 1, 1)
    for g in testgraphs(g1)
        c = @inferred(independent_set(g, DegreeIndependentSet()))
        @test isempty(c)
    end

    g3 = StarGraph(5)
    for g in testgraphs(g3)
        c = @inferred(independent_set(g, DegreeIndependentSet()))
        @test sort(c) == [2, 3, 4, 5]
    end

    g4 = CompleteGraph(5)
    for g in testgraphs(g4)
        c = @inferred(independent_set(g, DegreeIndependentSet()))
        @test length(c)== 1 #Exactly one vertex
    end

    #PathGraph(5) with additional edge 2-5
    g5 = PathGraph(5)
    add_edge!(g5, 2, 5)
    for g in testgraphs(g5)
        c = @inferred(independent_set(g, DegreeIndependentSet()))
        @test sort(c) == [1, 3, 5]
    end

    add_edge!(g5, 2, 2)
    add_edge!(g5, 3, 3)
    for g in testgraphs(g5)
        c = @inferred(independent_set(g, DegreeIndependentSet()))
        @test sort(c) == [1, 5]
    end
end


text:
### A Pluto.jl notebook ###
# v0.12.15

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 9795c3e0-33f9-11eb-14f0-d12d158539e2
using Plots

# ╔═╡ d6ba3460-33fa-11eb-314e-a102bf58747d
using PlutoUI

# ╔═╡ 9fd82cf0-33f9-11eb-1bcb-afef2c71688b
cardioid(a, t) = a * (1 - cos(t)) .* (cos(t), sin(t))

# ╔═╡ ae072830-33f9-11eb-00dd-87760f6fdbbd
md"""
Vary the $\theta$-slider to draw the cardioid $r = a(1 - \cos \theta),~ 0 \leq \theta \leq 2\pi,~ a = 1$.
"""

# ╔═╡ bceec510-33f9-11eb-3d94-a3d7d660dee5
begin
	θ_slider = @bind θ Slider(0.:1e-3:2π)
	md"""
	θ: $(θ_slider)
	"""
end

# ╔═╡ d508a940-33f9-11eb-2550-1984c73124f8
plot(cardioid.(1, 0:1e-3:θ), label = "Cardioid", xlims=(-2.25,0.75), ylims=(-1.5,1.5))

# ╔═╡ cf4745c0-33f9-11eb-39ce-efbae4b35676
log_spiral(a, b, θ) = a * exp(b * θ) .* (cos(θ), sin(θ))

# ╔═╡ d995f990-33f9-11eb-372c-adc2227982a5
φ = (1 + √5) / 2

# ╔═╡ dfd58dc0-33f9-11eb-2213-014573741740
# Fibonacci constant
b = log(φ) / (π/2)

# ╔═╡ e96ead30-33f9-11eb-29bb-2b64474204b2
md"""
Draw the Fibonacci spiral!
"""

# ╔═╡ fcf32abe-33f9-11eb-0456-cf2e17bc6995
begin
	t_slider = @bind t Slider(0.:1e-2:15π)
	md"""
	θ: $(t_slider)
	"""
end

# ╔═╡ 0da9bb40-33fa-11eb-19df-673cb9106c7a
plot(log_spiral.(1, atan(b), 0:1e-3:t), label = "Fibonacci Spiral")

# ╔═╡ f5f9bc70-33f9-11eb-0b55-5f73e5b1ec96
dickbutt(t) = (((sign(t)*sign(1-t)+1)*(1315*t^3-1905*t^2+840*t+15775)+(sign(t-1)*sign(2-t)+1)*(5455*t^3-36015*t^2+66090*t-19505)+(sign(t-2)*sign(3-t)+1)*(3164*t^3-27165*t^2+78981*t-58723)+(sign(t-3)*sign(4-t)+1)*(-528*t^3+5544*t^2-19008*t+38475)+(sign(t-4)*sign(5-t)+1)*(-273*t^3+222*t^2+17928*t-41156)+(sign(t-5)*sign(6-t)+1)*(3154*t^3-54165*t^2+304623*t-543331)+(sign(t-6)*sign(7-t)+1)*(-1239*t^3+24840*t^2-167043*t+391373)+(sign(t-7)*sign(8-t)+1)*(-1175*t^3+30597*t^2-260517*t+741646)+(sign(t-8)*sign(9-t)+1)*(-354*t^3+8658*t^2-69795*t+199614)+(sign(t-9)*sign(10-t)+1)*(-6562*t^3+195444*t^2-1940280*t+6429945)+(sign(t-10)*sign(11-t)+1)*(-2770*t^3+87867*t^2-926556*t+3258405)+(sign(t-11)*sign(12-t)+1)*(-2673*t^3+96666*t^2-1166910*t+4712751)+(sign(t-12)*sign(13-t)+1)*(578*t^3-20832*t^2+247758*t-960343)+(sign(t-13)*sign(14-t)+1)*(629*t^3-21423*t^2+234534*t-800599)+(sign(t-14)*sign(15-t)+1)*(196*t^3-8433*t^2+122160*t-585251)+(sign(t-15)*sign(16-t)+1)*(588*t^3-28224*t^2+450084*t-2370410)+(sign(t-16)*sign(17-t)+1)*(-1013*t^3+46845*t^2-718551*t+3667782)+(sign(t-17)*sign(18-t)+1)*(1556*t^3-81867*t^2+1434753*t-8362102)+(sign(t-18)*sign(19-t)+1)*(-308*t^3+17067*t^2-316209*t+1972565)+(sign(t-19)*sign(20-t)+1)*(2256*t^3-128730*t^2+2446311*t-15469074)+(sign(t-20)*sign(21-t)+1)*(-1204*t^3+69690*t^2-1339521*t+8567133)+(sign(t-21)*sign(22-t)+1)*(5577*t^3-362229*t^2+7844511*t-56622289)+(sign(t-22)*sign(23-t)+1)*(3443*t^3-234849*t^2+5337006*t-40386267)+(sign(t-23)*sign(24-t)+1)*(138*t^3-8718*t^2+181218*t-1214507)+(sign(t-24)*sign(25-t)+1)*(-1606*t^3+116121*t^2-2797179*t+22468813)+(sign(t-25)*sign(26-t)+1)*(510*t^3-39015*t^2+994500*t-8430184)+(sign(t-26)*sign(27-t)+1)*(874*t^3-68583*t^2+1793244*t-15606096)+(sign(t-27)*sign(28-t)+1)*(4520*t^3-369750*t^2+10079733*t-91558537)+(sign(t-28)*sign(29-t)+1)*(-72*t^3+5118*t^2-115953*t+827743)+(sign(t-29)*sign(30-t)+1)*(-2688*t^3+233913*t^2-6782205*t+65536125)+(sign(t-30)*sign(31-t)+1)*(2466*t^3-221274*t^2+6614814*t-65864145)+(sign(t-31)*sign(32-t)+1)*(-443*t^3+39552*t^2-1173078*t+11568742)+(sign(t-32)*sign(33-t)+1)*(-919*t^3+92427*t^2-3095331*t+34534406)+(sign(t-33)*sign(34-t)+1)*(570*t^3-61758*t^2+2217627*t-26398769)+(sign(t-34)*sign(35-t)+1)*(744*t^3-73845*t^2+2439450*t-26807075)+(sign(t-35)*sign(36-t)+1)*(445*t^3-49371*t^2+1822152*t-22362801)+(sign(t-36)*sign(37-t)+1)*(-1513*t^3+167979*t^2-6214362*t+76618551))/10000, ((sign(t)*sign(1-t)+1)*(552*t^3-1179*t^2+2877*t+12350)+(sign(t-1)*sign(2-t)+1)*(-7*t^3-14604*t^2+50079*t-20868)+(sign(t-2)*sign(3-t)+1)*(-218*t^3+585*t^2+495*t+11511)+(sign(t-3)*sign(4-t)+1)*(836*t^3-8778*t^2+30096*t-21231)+(sign(t-4)*sign(5-t)+1)*(2427*t^3-35316*t^2+167970*t-250507)+(sign(t-5)*sign(6-t)+1)*(2869*t^3-43467*t^2+215052*t-337392)+(sign(t-6)*sign(7-t)+1)*(2650*t^3-49278*t^2+301407*t-599022)+(sign(t-7)*sign(8-t)+1)*(-3964*t^3+84990*t^2-603525*t+1424972)+(sign(t-8)*sign(9-t)+1)*(-390*t^3+9960*t^2-83472*t+236580)+(sign(t-9)*sign(10-t)+1)*(-3567*t^3+103950*t^2-1001430*t+3201045)+(sign(t-10)*sign(11-t)+1)*(-3200*t^3+100104*t^2-1038042*t+3584765)+(sign(t-11)*sign(12-t)+1)*(-183*t^3+6231*t^2-69453*t+270896)+(sign(t-12)*sign(13-t)+1)*(1338*t^3-47853*t^2+566463*t-2210004)+(sign(t-13)*sign(14-t)+1)*(-6596*t^3+262383*t^2-3475191*t+15332612)+(sign(t-14)*sign(15-t)+1)*(-567*t^3+24870*t^2-361902*t+1755538)+(sign(t-15)*sign(16-t)+1)*(-513*t^3+23949*t^2-375684*t+1993260)+(sign(t-16)*sign(17-t)+1)*(1025*t^3-49179*t^2+784419*t-4147268)+(sign(t-17)*sign(18-t)+1)*(1245*t^3-63030*t^2+1061832*t-5940923)+(sign(t-18)*sign(19-t)+1)*(4519*t^3-249483*t^2+4586385*t-28061937)+(sign(t-19)*sign(20-t)+1)*(2191*t^3-128025*t^2+2492208*t-16151160)+(sign(t-20)*sign(21-t)+1)*(3221*t^3-196482*t^2+3992868*t-27022060)+(sign(t-21)*sign(22-t)+1)*(-4713*t^3+300819*t^2-6395349*t+45300018)+(sign(t-22)*sign(23-t)+1)*(1193*t^3-80574*t^2+1810695*t-13525826)+(sign(t-23)*sign(24-t)+1)*(416*t^3-29037*t^2+675039*t-5215052)+(sign(t-24)*sign(25-t)+1)*(-352*t^3+24447*t^2-565107*t+4358500)+(sign(t-25)*sign(26-t)+1)*(582*t^3-44523*t^2+1134900*t-9629639)+(sign(t-26)*sign(27-t)+1)*(811*t^3-64980*t^2+1734936*t-15425356)+(sign(t-27)*sign(28-t)+1)*(-772*t^3+63432*t^2-1738410*t+15910109)+(sign(t-28)*sign(29-t)+1)*(-802*t^3+69090*t^2-1983174*t+18986189)+(sign(t-29)*sign(30-t)+1)*(-3662*t^3+324837*t^2-9602070*t+94602487)+(sign(t-30)*sign(31-t)+1)*(2270*t^3-208773*t^2+6397752*t-65307173)+(sign(t-31)*sign(32-t)+1)*(2160*t^3-203499*t^2+6389070*t-66827661)+(sign(t-32)*sign(33-t)+1)*(-1965*t^3+191052*t^2-6190248*t+66858291)+(sign(t-33)*sign(34-t)+1)*(4280*t^3-428631*t^2+14306211*t-159115377)+(sign(t-34)*sign(35-t)+1)*(-2519*t^3+260031*t^2-8945484*t+102576877)+(sign(t-35)*sign(36-t)+1)*(1757*t^3-185985*t^2+6560850*t-77108863)+(sign(t-36)*sign(37-t)+1)*(-1907*t^3+207681*t^2-7537566*t+91190561))/10000)

# ╔═╡ ccb66c40-33fa-11eb-1089-6921936f1b1f
begin
	butt_slider = @bind butt Slider(0.:1e-3:37)
	md"""
	butt: $(butt_slider)
	"""
end

# ╔═╡ 19671a40-33fa-11eb-2902-6ba2b9d566e2
plot(dickbutt.(0:1e-3:butt), label = "Economist")

# ╔═╡ Cell order:
# ╠═9795c3e0-33f9-11eb-14f0-d12d158539e2
# ╠═d6ba3460-33fa-11eb-314e-a102bf58747d
# ╠═9fd82cf0-33f9-11eb-1bcb-afef2c71688b
# ╟─ae072830-33f9-11eb-00dd-87760f6fdbbd
# ╟─bceec510-33f9-11eb-3d94-a3d7d660dee5
# ╟─d508a940-33f9-11eb-2550-1984c73124f8
# ╠═cf4745c0-33f9-11eb-39ce-efbae4b35676
# ╠═d995f990-33f9-11eb-372c-adc2227982a5
# ╠═dfd58dc0-33f9-11eb-2213-014573741740
# ╟─e96ead30-33f9-11eb-29bb-2b64474204b2
# ╟─fcf32abe-33f9-11eb-0456-cf2e17bc6995
# ╟─0da9bb40-33fa-11eb-19df-673cb9106c7a
# ╟─f5f9bc70-33f9-11eb-0b55-5f73e5b1ec96
# ╟─ccb66c40-33fa-11eb-1089-6921936f1b1f
# ╠═19671a40-33fa-11eb-2902-6ba2b9d566e2

max_stars_repo_name: MathAero/IndEco

max_issues_repo_name: MathAero/IndEco

max_forks_repo_name: MathAero/IndEco


text:
#   Debug.Analysis:
# ===================
# Scoping analysis used to support instrument and graft

module Analysis
using Debug.AST, Debug.Meta
using Compat

export analyze, code_analyzed_only

# ---- Helpers ----------------------------------------------------------------

is_linenumber(ex::LineNumberNode) = true
is_linenumber(ex::Ex)             = is(headof(ex), :line)
is_linenumber(ex)                 = false

is_file_linenumber(ex)            = is_expr(ex, :line, 2)

get_linenumber(ex::Ex)             = argof(ex,1)
get_linenumber(ex::LineNumberNode) = ex.line

# ---- code_analyzed_only: wrap nodes that are only allowed in analyzed code --

analyzed_only(ex, error_message) = error(error_message)
function code_analyzed_only(ex, error_message)
    :( $analyzed_only($(quot(ex)), $(quot(error_message))) )
end

# ---- wrap(): add scoping info to AST ------------------------------------
# Rewrites AST, exchanging Expr:s and leaves for Node:s.
# Adds scope info, classifies nodes, etc.

# TypeEnv: Env used in a type block to throw away non-method definitions
type TypeEnv <: Env
    env::LocalEnv
end
raw(env::TypeEnv) = env.env
raw(env::Env)     = env

raw(s::State) = s
raw{T<:SimpleState}(s::T) = T(raw(s.env))

function wrap(states::Vector, args::Vector)
    [wrap(s, arg) for (s, arg) in zip(states, args)]
end
wrap(s::SimpleState, ex::SymbolNode) = wrap(s, ex.name)
wrap(s::State, ex) = enwrap(s, decorate(s,ex))

# SplitDef: Def with different scopes for left and right side, e.g.
# let x_inner = y_outer   or   type T_outer{S_inner} <: Q_outer
type SplitDef <: State;
    ls::Env;
    rs::Env;
end
wrap(s::SplitDef, ex) = wrap(Def(s.ls), ex)
function wrap(s::SplitDef, ex::Ex)
    head, nargs = headof(ex), nargsof(ex)
    if     head === :(=);   enwrap(s, decorate([Def(s.ls), Rhs(s.rs)], ex))
    elseif head === :(<:);  enwrap(s, decorate([s,         Rhs(s.ls)], ex))
    elseif head === :curly; enwrap(s, decorate([s, fill(Def(s.rs), nargs-1)...], ex))
    else                    wrap(Def(s.ls), ex)
    end
end

enwrap(s::State, value)    = Node(value, raw(s))
enwrap(s::State, loc::Loc) = Node(Location(), raw(s), loc)

## decorate: does most of the work for wrap ##

decorate(states::Vector, ex::Ex) = ExValue(headof(ex), wrap(states,argsof(ex)))

decorate(s::State, ex) = isa(ex, Node) ? valueof(ex) : Plain(ex) # wrap/unwrap
decorate(s::SimpleState, ex::LineNumberNode) = Loc(ex, ex.line)

decorate(s::Def,         ex::Symbol) = (AST.add_defined( s.env, ex); ex)
decorate(s::Lhs,         ex::Symbol) = (AST.add_assigned(s.env, ex); ex)
decorate(s::SimpleState, ex::Symbol) = ex

decorate(::Nonsyntax, ex::Ex) = Plain(ex)

# Sig: (part of) function signature in function f(x) ... / f(x) = ...
type Sig <: State;
    s::SimpleState;  # state with outer Env, to define/assign f
    inner::Env;      # Env inside the method
end
function decorate(s::Sig, ex::Ex)
    @assert headof(ex) in [:call, :curly]
    if is_expr(argof(ex,1), :curly);first = s
    else;                           first = (isa(s.s,Def) ? s.s : Lhs(s.s.env))
    end
    states = [first, fill(Def(s.inner), nargsof(ex)-1)...]
    decorate(states, ex)
end

# SplitState: Expression with different states to be applied for each argument
type SplitState <: State;
    env::Env
    argstates::Vector{State}
end
AST.envof(s::SplitState) = s.env
decorate(s::SplitState, ex::Ex) = decorate(s.argstates, ex)

const updating_ops = @compat Dict(
 :+= => :+,   :-= => :-,  :*= => :*,  :/= => :/,  ://= => ://, :.//= => :.//,
:.*= => :.*, :./= => :./, :\= => :\, :.\= => :.\,  :^= => :^,   :.^= => :.^,
 :%= => :%,   :|= => :|,  :&= => :&,  :$= => :$,  :<<= => :<<,  :>>= => :>>,
 :>>>= => :>>>)

# return a Vector of visit states for each arg of an expr(head, args)
function argstates(state::SimpleState, ex)
    head, args = headof(ex), argsof(ex)
    e, nargs = state.env, length(args)

    if (head in [:function, :(=)]) && is_expr(args[1], :call)
        if isa(e, TypeEnv); return argstates(Def(raw(e)), ex); end
        c = child(ex, e); [Sig(state, c), Rhs(c)]
    elseif head in [:function, :(->)]
        c = child(ex, e); [Def(c),        Rhs(c)]
    elseif head === :kw
        [(isa(state,Def) ? state : Rhs(e)), Rhs(e)]
    elseif head === :parameters
        fill(isa(state,Def) ? state : Rhs(e), nargs)

    elseif head in [:global, :local]; fill(Def(e), nargs)
    elseif head === :while;              [Rhs(e),        Rhs(child(ex, e))]
    elseif head === :try
        cc = child(ex,e); states = [Rhs(child(ex,e)), Def(cc),Rhs(cc)]
        nargs === 4 ? [states; Rhs(child(ex,e))] : states
    elseif head === :for
        if is_expr(args[1], :block)
            last_c = e
            states = State[]
            # Create a nested sequence of environments since each loop
            # variable can depend on the previous for its definition
            for child_ex in argsof(args[1])
                # Is there a better expression to take as ex?
                c = child(ex, last_c)
                push!(states, SplitDef(c,last_c))
                last_c = c
            end
            [SplitState(e,states), Rhs(last_c)]
        else
            c = child(ex, e); [SplitDef(c,e), Rhs(c)]
        end
    elseif head in [:let; untyped_comprehensions]; c = child(ex, e);
        [Rhs(c); fill(SplitDef(c,e), nargs-1)...]
    elseif head in typed_comprehensions; c = child(ex, e)
        [Rhs(e); Rhs(c); fill(SplitDef(c,e), nargs-2)...]

    elseif haskey(updating_ops, head); [Lhs(e), Rhs(e)]
    elseif head === :(=);   [(isa(state,Def) ? state : Lhs(e)), Rhs(e)]
    elseif head === :(<:);  [(isa(state,Def) ? state : Rhs(e)), Rhs(e)]
    elseif head === :tuple; fill(state,  nargs)
    elseif head === :ref;   fill(Rhs(e), nargs)
    elseif head === :(...); [state]
    elseif head === :(::) && nargs == 1; [Rhs(e)]
    elseif head === :(::) && nargs == 2; [state, Rhs(e)]

    # I'm guessing abstract and typealias wrap their args in one scope,
    # except the actual name to be defined
    elseif head === :abstract; c=child(ex, e); [SplitDef(e,c)]
    elseif head === :type;     c=child(ex, e); [Nonsyntax(),
                                                SplitDef(e,c), Rhs(TypeEnv(c))]
    elseif head === :typealias;c=child(ex, e); [SplitDef(e,c), Rhs(c)]

    else fill(Rhs(e), nargs)
    end
end

function decorate(state::SimpleState, ex::Ex)
    head, args  = headof(ex), argsof(ex)
    if head === :line;             return Loc(ex, args...)
    elseif head in [:quote, :top]; return Plain(ex)
    elseif head === :macrocall;    return decorate(state, macroexpand(ex))
    elseif head === :call && length(args) === 3 && args[1] === analyzed_only
        # Detect calls to analyzed_only and substitute with its first argument
        return decorate(state, (args[2]::QuoteNode).value)
    end

    decorate(argstates(state, ex), ex)
end

# ---- post-decoration processing ---------------------------------------------

## set_source!(): propagate source file info ##
function set_source!(ex::LocNode, locex, line, file)
    ex.loc.file = file
end
set_source!(ex::Node, locex, line, file) = (ex.loc = Loc(locex, line, file))
function set_source!(ex::ExNode, locex, line, file)
    ex.loc = Loc(locex, line, file)
    locex  = nothing
    for arg in argsof(ex)
        if isa(arg, LocNode)
            line = arg.loc.line #valueof(arg).line
            if arg.loc.file != empty_symbol;  file = arg.loc.file;  end
        end
        set_source!(arg, locex, line, file)
        locex = isa(arg, LocNode) ? arg.loc : nothing
    end
end

## postprocess_env!: find defined symbols among assigned ##
# todo: Make child links for Env so that this doesn't have to go through the
# nodes?
postprocess_env!(envs::Set{LocalEnv}, ex) = nothing
function postprocess_env!(envs::Set{LocalEnv}, ex::SymNode)
    postprocess_env!(envs,envof(ex))
end
function postprocess_env!(envs::Set{LocalEnv}, ex::Ex)
    if isblocknode(ex); postprocess_env!(envs, envof(ex)); end
    for arg in argsof(ex); postprocess_env!(envs, arg);    end
end

postprocess_env!(envs::Set{LocalEnv}, ::NoEnv) = nothing
function postprocess_env!(envs::Set{LocalEnv}, env::LocalEnv)
    if env in envs; return; end
    push!(envs, env)
    p = env.parent
    p_assigned = isa(p, LocalEnv) ? p.assigned : Set{@compat(Union{})}()
    env.defined  = union(env.defined, setdiff(env.assigned, p_assigned))
    env.assigned = union(env.defined, p_assigned)
end

# ---- analyze(): wrap and then propagate source file info among LocNode's

analyze(ex, process_envs::Bool) = analyze(Rhs(NoEnv()), ex, process_envs)
analyze(env::Env, ex, process_envs::Bool) = analyze(Rhs(env), ex, process_envs)
function analyze(s::State, ex, process_envs::Bool)
    node = wrap(s, ex)
    set_source!(node, nothing, -1, empty_symbol)
    if process_envs; postprocess_env!(Set{LocalEnv}(), node); end
    node
end

end # module


text:
#!/usr/bin/env julia

using neten

adjm = data("butterfly")

netw = enhance(adjm)


text:
import Distances: pairwise
import DelayEmbeddings: AbstractDataset
import StaticArrays: SVector, MVector
import DelayEmbeddings: KDTree
import NearestNeighbors: knn

import Distances: pairwise
import DelayEmbeddings: AbstractDataset
import StaticArrays: SVector, MVector
import DelayEmbeddings: KDTree
import NearestNeighbors: knn

KDTree(D::CustomReconstruction, metric::Metric = Euclidean()) = KDTree(D.reconstructed_pts, metric)

# Define pairwise distances methods for vectors of static vectors, vectors of vectors,
# and Datasets and Customreconstructions.
const VSV = Union{AbstractDataset, Vector{<:SVector}, Vector{<:MVector}, Vector{Vector}}

function pairwise(x::VSV, y::VSV, metric::Metric)
    length(x) == length(y) || error("Lengths of input data does not match")

    dists = zeros(eltype(eltype(x)), length(x), length(y))
    @inbounds for j in 1:length(y)
        for i in 1:length(x)
            dists[i, j] = evaluate(metric, x[i], y[j])
        end
    end
    return dists
end

function pairwise(x::VSV, metric::Metric)
    dists = zeros(eltype(eltype(x)), length(x), length(x))
    @inbounds for j in 1:length(x)
        for i in 1:length(x)
            dists[i, j] = evaluate(metric, x[i], x[j])
        end
    end
    return dists
end

function find_dists(x::VSV, y::VSV, metric::Metric)
    length(x) == length(y) || error("Lengths of input data does not match")
    dists = zeros(eltype(eltype(x)), length(x))

    @inbounds for i = 1:length(x)
        dists[i] = evaluate(metric, x[i], y[i])
    end

    return dists
end

function marginal_NN(points::VSV, dists_to_kth)
    D = pairwise(points, Chebyshev())

    N = zeros(Int, length(points))

    @inbounds for i = 1:length(points)
        N[i] = sum(D[:, i] .<= dists_to_kth[i])
    end

    return N
end

const SVV{N, T} = Union{Vector{<:SVector{N, T}}, Vector{<:MVector{N, T}}} where {N, T}

function get_marginal_pts(pts::SVV{N, T}, vars::TEVars) where {N, T}

     # Create some dummy variable names to avoid cluttering the code too much
    X = vars.target_future
    Y = vars.target_presentpast
    Z = vcat(vars.source_presentpast, vars.conditioned_presentpast)
    XY = vcat(X, Y)
    YZ = vcat(Y, Z)

    nX = length(X)
    nY = length(Y)
    nZ = length(Z)
    nXY = nX + nY
    nYZ = nY + nZ

    pts_X = [SVector{nX, T}(pt[X]) for pt in pts]
    pts_Y = [SVector{nY, T}(pt[Y]) for pt in pts]
    pts_XY = [SVector{nXY, T}(pt[XY]) for pt in pts]
    pts_YZ = [SVector{nYZ, T}(pt[YZ]) for pt in pts]
    pts_XYZ = pts

    return pts_X, pts_Y, pts_XY, pts_YZ, pts_XYZ
end

function get_marginal_pts(pts::Vector{Vector{T}}, vars::TEVars) where T
    get_marginal_pts(Dataset(pts), vars)
end

function get_marginal_pts(pts::AbstractDataset, vars::TEVars)

     # Create some dummy variable names to avoid cluttering the code too much
    X = vars.target_future
    Y = vars.target_presentpast
    Z = vcat(vars.source_presentpast, vars.conditioned_presentpast)
    XY = vcat(X, Y)
    YZ = vcat(Y, Z)

    nX = length(X)
    nY = length(Y)
    nZ = length(Z)
    nXY = nX + nY
    nYZ = nY + nZ

    pts_X = pts[:, X]
    pts_Y = pts[:, Y]
    pts_XY = pts[:, XY]
    pts_YZ = pts[:, YZ]
    pts_XYZ = pts

    return pts_X, pts_Y, pts_XY, pts_YZ, pts_XYZ

end

"""
    transferentropy(pts, vars::TEVars, estimator::NearestNeighbourMI) -> Float64

Compute the transfer entropy for a set of `pts`, computing mutual
information (MI) terms using counting of nearest neighbours.

Uses the MI estimator from [1], which was implemented to compute
transfer entropy in [2].

## Arguments

- **`pts`**: An ordered set of `m`-dimensional points (`pts`) representing
    an appropriate [generalised embedding](@ref) of some data series. Must be
    vector of states, not a vector of variables/time series. Wrap your time
    series in a `DynamicalSystemsBase.Dataset` first if the latter is the case.

- **`vars::TEVars`**: A [`TEVars`](@ref) instance specifying how the `m` different
    variables of `pts` are to be mapped into the marginals required for transfer
    entropy computation.

- **`estimator::NearestNeighbourMI`**: An instance of a [`NearestNeighbourMI`](@ref),
    estimator, which contains information about the number of nearest neighbours to use,
    the distance metric and the base of the logarithm that controls the
    unit of the transfer entropy.

## Returns

A single number that is the transfer entropy estimated over `pts` by assigning marginals
according to `vars`, using nearest neighbour counting method [1] to estimate
the necessary mutual information terms [2].

## References

1. Kraskov, Alexander, Harald Stögbauer, and Peter Grassberger. "Estimating
    mutual information." Physical review E 69.6 (2004): 066138.
2. Diego, David, Kristian Agasøster Haaga, and Bjarte Hannisdal. "Transfer entropy computation
    using the Perron-Frobenius operator." Physical Review E 99.4 (2019): 042212.
"""
function transferentropy(pts, vars::TEVars, estimator::NearestNeighbourMI)

    metric = estimator.metric
    b = estimator.b
    k1 = estimator.k1
    k2 = estimator.k2

    npts = length(pts)

    pts_X, pts_Y, pts_XY, pts_YZ, pts_XYZ = get_marginal_pts(pts, vars)

    # Create trees to search for nearest neighbors
    tree_XYZ = KDTree(pts_XYZ, metric)
    tree_XY = KDTree(pts_XY, metric)
    # Find the k nearest neighbors to all of the points in each of the trees.
    # We sort the points according to their distances, because we want precisely
    # the k-th nearest neighbor.
    knearest_XYZ = [knn(tree_XYZ, pt, k1, true) for pt in pts_XYZ]
    knearest_XY   = [knn(tree_XY, pt, k2, true) for pt in pts_XY]

    # In each of the trees, find the index of the k-th nearest neighbor to all
    # of the points.
    kth_NN_idx_XYZ = [ptnbrs[1][k1] for ptnbrs in knearest_XYZ]
    kth_NN_idx_XY  = [ptnbrs[1][k2] for ptnbrs in knearest_XY]

    # Distances between points in the XYZ and XY spaces and their
    # kth nearest neighbour, along marginals X, YZ (for the XYZ space)
    ϵ_XYZ_X = find_dists(pts_X, pts_X[kth_NN_idx_XYZ], metric)
    ϵ_XYZ_YZ = find_dists(pts_YZ, pts_YZ[kth_NN_idx_XYZ], metric)
    ϵ_XY_X = find_dists(pts_X, pts_X[kth_NN_idx_XY], metric)
    ϵ_XY_Y = find_dists(pts_Y, pts_Y[kth_NN_idx_XY], metric)

    NXYZ_X  = marginal_NN(pts_X,  ϵ_XYZ_X)
    NXYZ_YZ = marginal_NN(pts_YZ, ϵ_XYZ_YZ)
    NXY_X   = marginal_NN(pts_X,  ϵ_XY_X)
    NXY_Y   = marginal_NN(pts_Y,  ϵ_XY_Y)

    # Compute transfer entropy
    te = sum(digamma.(NXY_X) + digamma.(NXY_Y) -
            digamma.(NXYZ_X) - digamma.(NXYZ_YZ)) / npts

    # Convert from nats to to the desired unit (b^x = e^1 => 1/ln(b))
    return te / log(b)
end

function transferentropy(pts::CustomReconstruction, vars::TEVars, estimator::NearestNeighbourMI)
    transferentropy(pts.reconstructed_pts, vars, estimator)
end


text:
module SignedDistanceFunctions

using GeometryBasics
using Rotations
using CoordinateTransformations
using StaticArrays
using LinearAlgebra

using FiniteDifferences
using ForwardDiff
using Zygote

include("types.jl")
include("shapes.jl")

export
    AbstractSignedDistanceFunction,

    normal,
    normalforwarddiff,
    normalzygote,
    normalfinitediff,

    BoundingVolumeSignedDistanceFunction,
    overshootfactor, boundingSDF, interiorSDF,

    UnionSignedDistanceFunction,
    SDFs,

    IntersectSignedDistanceFunction,

    NegatedSignedDistanceFunction,
    sdf,

    TransformedSignedDistanceFunction,
    transform, invtransform, transform_sdf,

    SphereSignedDistanceFunction,
    radius

end


text:
# Pointers with address space information

#
# Address spaces
#

export AS, addrspace

abstract type AddressSpace end

module AS

import ..AddressSpace

struct Generic  <: AddressSpace end
struct Global   <: AddressSpace end
struct Shared   <: AddressSpace end
struct Constant <: AddressSpace end
struct Local    <: AddressSpace end

end

Base.convert(::Type{Int}, ::Type{AS.Generic})  = 0
Base.convert(::Type{Int}, ::Type{AS.Global})   = 1
Base.convert(::Type{Int}, ::Type{AS.Shared})   = 3
Base.convert(::Type{Int}, ::Type{AS.Constant}) = 4
Base.convert(::Type{Int}, ::Type{AS.Local})    = 5

tbaa_addrspace(as::Type{<:AddressSpace}) = tbaa_make_child(lowercase(String(as.name.name)))

#
# Device pointer
#

"""
    DevicePtr{T,A}

A memory address that refers to data of type `T` that is accessible from the GPU. It is the
on-device counterpart of `CuPtr`, additionally keeping track of the address space
`A` where the data resides (shared, global, constant, etc). This information is used to
provide optimized implementations of operations such as `unsafe_load` and `unsafe_store!.`
"""
DevicePtr

if sizeof(Ptr{Cvoid}) == 8
    primitive type DevicePtr{T,A} <: Ref{T} 64 end
else
    primitive type DevicePtr{T,A} <: Ref{T} 32 end
end

# constructors
DevicePtr{T,A}(x::Union{Int,UInt,CuPtr,DevicePtr}) where {T,A<:AddressSpace} = Base.bitcast(DevicePtr{T,A}, x)
DevicePtr{T}(ptr::CuPtr{T}) where {T} = DevicePtr{T,AS.Generic}(ptr)
DevicePtr(ptr::CuPtr{T}) where {T} = DevicePtr{T,AS.Generic}(ptr)

## getters

Base.eltype(::Type{<:DevicePtr{T}}) where {T} = T

addrspace(x::DevicePtr) = addrspace(typeof(x))
addrspace(::Type{DevicePtr{T,A}}) where {T,A} = A

## conversions

# to and from integers
## pointer to integer
Base.convert(::Type{T}, x::DevicePtr) where {T<:Integer} = T(UInt(x))
## integer to pointer
Base.convert(::Type{DevicePtr{T,A}}, x::Union{Int,UInt}) where {T,A<:AddressSpace} = DevicePtr{T,A}(x)
Int(x::DevicePtr)  = Base.bitcast(Int, x)
UInt(x::DevicePtr) = Base.bitcast(UInt, x)

# between host and device pointers
Base.convert(::Type{CuPtr{T}},  p::DevicePtr)  where {T}                 = Base.bitcast(CuPtr{T}, p)
Base.convert(::Type{DevicePtr{T,A}}, p::CuPtr) where {T,A<:AddressSpace} = Base.bitcast(DevicePtr{T,A}, p)
Base.convert(::Type{DevicePtr{T}}, p::CuPtr)   where {T}                 = Base.bitcast(DevicePtr{T,AS.Generic}, p)

# between CPU pointers, for the purpose of working with `ccall`
Base.unsafe_convert(::Type{Ptr{T}}, x::DevicePtr{T}) where {T} = reinterpret(Ptr{T}, x)

# between device pointers
Base.convert(::Type{<:DevicePtr}, p::DevicePtr)                         = throw(ArgumentError("cannot convert between incompatible device pointer types"))
Base.convert(::Type{DevicePtr{T,A}}, p::DevicePtr{T,A})   where {T,A}   = p
Base.unsafe_convert(::Type{DevicePtr{T,A}}, p::DevicePtr) where {T,A}   = Base.bitcast(DevicePtr{T,A}, p)
## identical addrspaces
Base.convert(::Type{DevicePtr{T,A}}, p::DevicePtr{U,A}) where {T,U,A} = Base.unsafe_convert(DevicePtr{T,A}, p)
## convert to & from generic
Base.convert(::Type{DevicePtr{T,AS.Generic}}, p::DevicePtr)               where {T}     = Base.unsafe_convert(DevicePtr{T,AS.Generic}, p)
Base.convert(::Type{DevicePtr{T,A}}, p::DevicePtr{U,AS.Generic})          where {T,U,A} = Base.unsafe_convert(DevicePtr{T,A}, p)
Base.convert(::Type{DevicePtr{T,AS.Generic}}, p::DevicePtr{T,AS.Generic}) where {T}     = p  # avoid ambiguities
## unspecified, preserve source addrspace
Base.convert(::Type{DevicePtr{T}}, p::DevicePtr{U,A}) where {T,U,A} = Base.unsafe_convert(DevicePtr{T,A}, p)

## limited pointer arithmetic & comparison

isequal(x::DevicePtr, y::DevicePtr) = (x === y) && addrspace(x) == addrspace(y)
isless(x::DevicePtr{T,A}, y::DevicePtr{T,A}) where {T,A<:AddressSpace} = x < y

Base.:(==)(x::DevicePtr, y::DevicePtr) = UInt(x) == UInt(y) && addrspace(x) == addrspace(y)
Base.:(<)(x::DevicePtr,  y::DevicePtr) = UInt(x) < UInt(y)
Base.:(-)(x::DevicePtr,  y::DevicePtr) = UInt(x) - UInt(y)

Base.:(+)(x::DevicePtr, y::Integer) = oftype(x, Base.add_ptr(UInt(x), (y % UInt) % UInt))
Base.:(-)(x::DevicePtr, y::Integer) = oftype(x, Base.sub_ptr(UInt(x), (y % UInt) % UInt))
Base.:(+)(x::Integer, y::DevicePtr) = y + x

## memory operations

@generated function pointerref(p::DevicePtr{T,A}, i::Int, ::Val{align}) where {T,A,align}
    sizeof(T) == 0 && return T.instance
    eltyp = convert(LLVMType, T)

    T_int = convert(LLVMType, Int)
    T_ptr = convert(LLVMType, DevicePtr{T,A})

    T_actual_ptr = LLVM.PointerType(eltyp, convert(Int, A))

    # create a function
    param_types = [T_ptr, T_int]
    llvm_f, _ = create_function(eltyp, param_types)

    # generate IR
    Builder(JuliaContext()) do builder
        entry = BasicBlock(llvm_f, "entry", JuliaContext())
        position!(builder, entry)

        ptr = inttoptr!(builder, parameters(llvm_f)[1], T_actual_ptr)
        ptr = gep!(builder, ptr, [parameters(llvm_f)[2]])
        ld = load!(builder, ptr)

        if A != AS.Generic
            metadata(ld)[LLVM.MD_tbaa] = tbaa_addrspace(A)
        end
        alignment!(ld, align)

        ret!(builder, ld)
    end

    call_function(llvm_f, T, Tuple{DevicePtr{T,A}, Int}, :((p, Int(i-one(i)))))
end

@generated function pointerset(p::DevicePtr{T,A}, x::T, i::Int, ::Val{align}) where {T,A,align}
    sizeof(T) == 0 && return
    eltyp = convert(LLVMType, T)

    T_int = convert(LLVMType, Int)
    T_ptr = convert(LLVMType, DevicePtr{T,A})

    T_actual_ptr = LLVM.PointerType(eltyp, convert(Int, A))

    # create a function
    param_types = [T_ptr, eltyp, T_int]
    llvm_f, _ = create_function(LLVM.VoidType(JuliaContext()), param_types)

    # generate IR
    Builder(JuliaContext()) do builder
        entry = BasicBlock(llvm_f, "entry", JuliaContext())
        position!(builder, entry)

        ptr = inttoptr!(builder, parameters(llvm_f)[1], T_actual_ptr)
        ptr = gep!(builder, ptr, [parameters(llvm_f)[3]])
        val = parameters(llvm_f)[2]
        st = store!(builder, val, ptr)

        if A != AS.Generic
            metadata(st)[LLVM.MD_tbaa] = tbaa_addrspace(A)
        end
        alignment!(st, align)

        ret!(builder)
    end

    call_function(llvm_f, Cvoid, Tuple{DevicePtr{T,A}, T, Int},
                  :((p, convert(T,x), Int(i-one(i)))))
end

# operand types supported by llvm.nvvm.ldg.global
# NOTE: CUDA 8.0 supports more caching modifiers, but those aren't supported by LLVM yet
const LDGTypes = Union{UInt8, UInt16, UInt32, UInt64,
                       Int8, Int16, Int32, Int64,
                       Float32, Float64}

# TODO: this functionality should throw <sm_32
@generated function pointerref_ldg(p::DevicePtr{T,AS.Global}, i::Int,
                                   ::Val{align}) where {T<:LDGTypes,align}
    sizeof(T) == 0 && return T.instance
    eltyp = convert(LLVMType, T)

    # TODO: ccall the intrinsic directly with LLVMPtr

    T_int = convert(LLVMType, Int)
    T_int32 = LLVM.Int32Type(JuliaContext())
    T_ptr = convert(LLVMType, DevicePtr{T,AS.Global})

    T_actual_ptr = LLVM.PointerType(eltyp, convert(Int, AS.Global))

    # create a function
    param_types = [T_ptr, T_int]
    llvm_f, _ = create_function(eltyp, param_types)

    # create the intrinsic
    intrinsic_name = let
        class = if isa(eltyp, LLVM.IntegerType)
            :i
        elseif isa(eltyp, LLVM.FloatingPointType)
            :f
        else
            error("Cannot handle $eltyp argument to unsafe_cached_load")
        end
        width = sizeof(T)*8
        typ = Symbol(class, width)
        "llvm.nvvm.ldg.global.$class.$typ.p1$typ"
    end
    mod = LLVM.parent(llvm_f)
    intrinsic_typ = LLVM.FunctionType(eltyp, [T_actual_ptr, T_int32])
    intrinsic = LLVM.Function(mod, intrinsic_name, intrinsic_typ)

    # generate IR
    Builder(JuliaContext()) do builder
        entry = BasicBlock(llvm_f, "entry", JuliaContext())
        position!(builder, entry)

        ptr = inttoptr!(builder, parameters(llvm_f)[1], T_actual_ptr)
        ptr = gep!(builder, ptr, [parameters(llvm_f)[2]])
        ld = call!(builder, intrinsic,
                   [ptr, ConstantInt(Int32(align), JuliaContext())])

        metadata(ld)[LLVM.MD_tbaa] = tbaa_addrspace(AS.Global)

        ret!(builder, ld)
    end

    call_function(llvm_f, T, Tuple{DevicePtr{T,AS.Global}, Int}, :((p, Int(i-one(i)))))
end

# interface

export unsafe_cached_load

Base.unsafe_load(p::DevicePtr{T}, i::Integer=1, align::Val=Val(1)) where {T} =
    pointerref(p, Int(i), align)

Base.unsafe_store!(p::DevicePtr{T}, x, i::Integer=1, align::Val=Val(1)) where {T} =
    pointerset(p, convert(T, x), Int(i), align)

unsafe_cached_load(p::DevicePtr{<:LDGTypes,AS.Global}, i::Integer=1, align::Val=Val(1)) =
    pointerref_ldg(p, Int(i), align)
# NOTE: fall back to normal pointerref for unsupported types. we could be smarter here,
#       e.g. destruct/load/reconstruct, but that's too complicated for what it's worth.
unsafe_cached_load(p::DevicePtr, i::Integer=1, align::Val=Val(1)) =
    pointerref(p, Int(i), align)


text:
# This file is part of Kpax3. License is MIT.

function logmarglik(y::Vector{Float64},
                    n::Real,
                    α::Vector{Float64},
                    β::Vector{Float64})
  SpecialFunctions.logbeta.(α .+ y, β .- y .+ n) .-
  SpecialFunctions.logbeta.(α, β)
end

function logmarglik(y::Real,
                    n::Real,
                    α::Real,
                    β::Real)
  SpecialFunctions.logbeta(α + y, β + n - y) - SpecialFunctions.logbeta(α, β)
end

function logcondmarglik(x::Vector{UInt8},
                        y::Vector{Float64},
                        n::Real,
                        α::Vector{Float64},
                        β::Vector{Float64})
  float(x) .* log.(α .+ y) .+
  (one(Float64) .- float(x)) .* log.(β .- y .+ n) .-
  log.(α .+ β .+ n)
end

function logcondmarglik(x::Real,
                        y::Real,
                        n::Real,
                        α::Real,
                        β::Real)
  (x * log(α + y) + (one(x) - x) * log(β + n - y) - log(α + β + n))::Float64
end

function marglik(y::Vector{Float64},
                 n::Real,
                 α::Vector{Float64},
                 β::Vector{Float64})
  exp.(logmarglik(y, n, α, β))
end

function marglik(y::Real,
                 n::Real,
                 α::Real,
                 β::Real)
  exp(logmarglik(y, n, α, β))
end

function condmarglik(x::Vector{UInt8},
                     y::Vector{Float64},
                     n::Real,
                     α::Vector{Float64},
                     β::Vector{Float64})
  exp.(logcondmarglik(x, y, n, α, β))
end

function condmarglik(x::Real,
                     y::Real,
                     n::Real,
                     α::Real,
                     β::Real)
  exp(logcondmarglik(x, y, n, α, β))
end


text:
const TypeMASK = Union{Nothing,AbstractArray{Bool,2}}

function st_mask!(ga::AbstractGeoArray, mask::AbstractGeoArray)
  # arr = zeros(T, size(mask)..., ntime)
  ind = findall(.!(mask.A)) # 只保留为true的部分
  ntime = size(ga, 3)

  for i = 1:ntime
    x = @views ga.A[:, :, i]
    x[ind] .= NaN
  end
end

export st_mask!


text:
function get_positional_encoding(song::AbstractArray)

    # Positional data
    word_size = size(song)[1]
    sequence_length = size(song)[2]
    num_ins = size(song)[3]

    pos_data = [k%2==0 ? sin(1/(10000^(2*k/word_size))*t) : cos(1/(10000^(2*k/word_size))*t)
                    for k in 1:word_size,
                        t in (-sequence_length+1):(sequence_length-1),
                        i in (-num_ins+1):(num_ins-1)]

    return pos_data
end


text:
using Test
using Pitchjx
using DataFrames

# test for no game day.
empty = DataFrame(
  date=String[],
  pitcherid=String[],
  pitcher_teamid=String[],
  pitcher_firstname=String[],
  pitcher_lastname=String[],
  pitcher_teamname=String[],
  pitcherthrow=String[],
  batterid=String[],
  batter_teamid=String[],
  batter_firstname=String[],
  batter_lastname=String[],
  batter_teamname=String[],
  batterstand=String[],
  eventdesc=String[],
  pitchresult=String[],
  x=String[],
  y=String[],
  px=String[],
  pz=String[],
  sztop=String[],
  szbottom=String[],
  pitchtype=String[],
  startspeed=String[],
  endspeed=String[],
  spindir=String[],
  spinrate=String[]
)

@test pitchjx("2018-01-01") == empty


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

module PeriodsTest

import Dates
using ExtendedDates
using Test

@testset "basic arithmetic" begin
    @test -Year(1) == Year(-1)
    @test Year(1) > Year(0)
    @test (Year(1) < Year(0)) == false
    @test Year(1) == Year(1)
    @test Year(1) != 1
    @test Year(1) + Year(1) == Year(2)
    @test Year(1) - Year(1) == zero(Year)
    @test 1 == one(Year)
    @test_throws MethodError Year(1) * Year(1) == Year(1)
    t = Year(1)
    t2 = Year(2)
    @test ([t, t, t, t, t] .+ Year(1)) == ([t2, t2, t2, t2, t2])
    @test (Year(1) .+ [t, t, t, t, t]) == ([t2, t2, t2, t2, t2])
    @test ([t2, t2, t2, t2, t2] .- Year(1)) == ([t, t, t, t, t])
    @test_throws MethodError ([t, t, t, t, t] .* Year(1)) == ([t, t, t, t, t])
    @test ([t, t, t, t, t] * 1) == ([t, t, t, t, t])
    @test ([t, t, t, t, t] .% t2) == ([t, t, t, t, t])
    @test div.([t, t, t, t, t], Year(1)) == ([1, 1, 1, 1, 1])
    @test mod.([t, t, t, t, t], Year(2)) == ([t, t, t, t, t])
    @test [t, t, t] / t2 == [0.5, 0.5, 0.5]
    @test abs(-t) == t
    @test sign(t) == sign(t2) == 1
    @test sign(-t) == sign(-t2) == -1
    @test sign(Year(0)) == 0
end
@testset "div/mod/gcd/lcm/rem" begin
    @test Year(10) % Year(4) == Year(2)
    @test gcd(Year(10), Year(4)) == Year(2)
    @test lcm(Year(10), Year(4)) == Year(20)
    @test div(Year(10), Year(3)) == 3
    @test div(Year(10), Year(4)) == 2
    @test div(Year(10), 4) == Year(2)
    @test Year(10) / Year(4) == 2.5

    @test mod(Year(10), Year(4)) == Year(2)
    @test mod(Year(-10), Year(4)) == Year(2)
    @test mod(Year(10), 4) == Year(2)
    @test mod(Year(-10), 4) == Year(2)

    @test rem(Year(10), Year(4)) == Year(2)
    @test rem(Year(-10), Year(4)) == Year(-2)
    @test rem(Year(10), 4) == Year(2)
    @test rem(Year(-10), 4) == Year(-2)
end

y = Year(1)
q = Quarter(1)
m = Month(1)
w = Week(1)
d = Day(1)
emptyperiod = ((y + d) - d) - y
@testset "Period arithmetic" begin
    @test Year(y) == y
    @test Quarter(q) == q
    @test Month(m) == m
    @test Week(w) == w
    @test Day(d) == d
    @test Year(convert(Int8, 1)) == y
    @test Year(convert(UInt8, 1)) == y
    @test Year(convert(Int16, 1)) == y
    @test Year(convert(UInt16, 1)) == y
    @test Year(convert(Int32, 1)) == y
    @test Year(convert(UInt32, 1)) == y
    @test Year(convert(Int64, 1)) == y
    @test Year(convert(UInt64, 1)) == y
    @test Year(convert(Int128, 1)) == y
    @test Year(convert(UInt128, 1)) == y
    @test Year(convert(BigInt, 1)) == y
    @test Year(convert(BigFloat, 1)) == y
    @test Year(convert(Complex, 1)) == y
    @test Year(convert(Rational, 1)) == y
    @test Year(convert(Float16, 1)) == y
    @test Year(convert(Float32, 1)) == y
    @test Year(convert(Float64, 1)) == y
    @test y == y
    @test m == m
    @test w == w
    @test d == d
    y2 = Year(2)
    @test y < y2
    @test y2 > y
    @test y != y2

    @test Year(Int8(1)) == y
    @test Year(UInt8(1)) == y
    @test Year(Int16(1)) == y
    @test Year(UInt16(1)) == y
    @test Year(Int(1)) == y
    @test Year(UInt(1)) == y
    @test Year(Int64(1)) == y
    @test Year(UInt64(1)) == y
    @test Year(UInt128(1)) == y
    @test Year(UInt128(1)) == y
    @test Year(big(1)) == y
    @test Year(BigFloat(1)) == y
    @test Year(float(1)) == y
    @test Year(Float32(1)) == y
    @test Year(Rational(1)) == y
    @test Year(complex(1)) == y
    @test_throws InexactError Year(BigFloat(1.2)) == y
    @test_throws InexactError Year(1.2) == y
    @test_throws InexactError Year(Float32(1.2)) == y
    @test_throws InexactError Year(3 // 4) == y
    @test_throws InexactError Year(complex(1.2)) == y
    @test_throws InexactError Year(Float16(1.2)) == y
    @test Year(true) == y
    @test Year(false) != y
    @test_throws MethodError Year(:hey) == y
    @test Year(real(1)) == y
    @test_throws InexactError Year(m) == y
    @test_throws MethodError Year(w) == y
    @test_throws MethodError Year(d) == y
    @test Year(Date(2013, 1, 1)) == Year(2013)
    @test Year(DateTime(2013, 1, 1)) == Year(2013)
    @test typeof(y + m) <: Dates.CompoundPeriod
    @test typeof(m + y) <: Dates.CompoundPeriod
    @test typeof(y + w) <: Dates.CompoundPeriod
    @test typeof(y + d) <: Dates.CompoundPeriod
    @test y > m
    @test d < w
    @test typemax(Year) == Year(typemax(Int64))
    @test typemax(Year) + y == Year(-9223372036854775808)
    @test typemin(Year) == Year(-9223372036854775808)
end
@testset "Period-Real arithmetic" begin
    @test_throws MethodError y + 1 == Year(2)
    @test_throws MethodError y + true == Year(2)
    @test_throws InexactError y + Year(1.2)
    @test y + Year(1.0f0) == Year(2)
    @test y * 4 == Year(4)
    @test y * 4.0f0 == Year(4)
    @test Year(2) * 0.5 == y
    @test Year(2) * 3 // 2 == Year(3)
    @test_throws InexactError y * 0.5
    @test_throws InexactError y * 3 // 4
    @test Year(4) / 2 == Year(2)
    @test Year(4) / 2.0f0 == Year(2)
    @test Year(4) / 0.5 == Year(8)
    @test Year(4) / 2 // 3 == Year(6)
    @test_throws InexactError Year(4) / 3.0
    @test_throws InexactError Year(4) / 3 // 2
    @test div(y, 2) == Year(0)
    @test_throws MethodError div(2, y) == Year(2)
    @test div(y, y) == 1
    @test y * 10 % Year(5) == Year(0)
    @test_throws MethodError (y > 3) == false
    @test_throws MethodError (4 < y) == false
    @test 1 != y
    t = [y, y, y, y, y]
    @test t .+ Year(2) == [Year(3), Year(3), Year(3), Year(3), Year(3)]

    let x = Year(5), y = Year(2)
        @test div(x, y) * y + rem(x, y) == x
        @test fld(x, y) * y + mod(x, y) == x
    end
end
@testset "Associativity" begin
    dt = DayDate(2012, 12, 21)
    test = ((((dt + y) - m) + w) - d)
    @test test == dt + y - m + w - d
    @test test == y - m + w - d + dt
    @test test == dt - m + y - d + w
    @test test == dt + (y - m + w - d)
    @test test == dt + y - m + w - d
    @test (dt + Year(4)) + Day(1) == dt + (Year(4) + Day(1))
    @test DayDate(2014, 1, 29) + Month(1) + Day(1) + Month(1) + Day(1) ==
          DayDate(2014, 1, 29) + Day(1) + Month(1) + Month(1) + Day(1)
    @test DayDate(2014, 1, 29) + Month(1) + Day(1) ==
          DayDate(2014, 1, 29) + Day(1) + Month(1)
end
@testset "traits" begin
    @test ExtendedDates._units(Year(0)) == " years"
    @test ExtendedDates._units(Year(1)) == " year"
    @test ExtendedDates._units(Year(-1)) == " year"
    @test ExtendedDates._units(Year(2)) == " years"
    @test string(Year(0)) == "0 years"
    @test string(Year(1)) == "1 year"
    @test string(Year(-1)) == "-1 year"
    @test string(Year(2)) == "2 years"
    @test isfinite(Year)
    @test isfinite(Year(0))
    @test zero(Year) == Year(0)
    @test zero(Year(10)) == Year(0)
    @test zero(Month) == Month(0)
    @test zero(Month(10)) == Month(0)
    @test zero(Day) == Day(0)
    @test zero(Day(10)) == Day(0)
    @test Year(-1) < Year(1)
    @test !(Year(-1) > Year(1))
    @test Year(1) == Year(1)
    @test Year(1) != 1
    @test 1 != Year(1)
    @test Month(-1) < Month(1)
    @test !(Month(-1) > Month(1))
    @test Month(1) == Month(1)
    @test Day(-1) < Day(1)
    @test !(Day(-1) > Day(1))
    @test Day(1) == Day(1)

    # issue #27076
    @test Year(1) != Dates.Millisecond(1)
    @test Dates.Millisecond(1) != Year(1)
end

@testset "basic properties" begin
    @test Year("1") == y
    @test Quarter("1") == q
    @test Month("1") == m
    @test Week("1") == w
    @test Day("1") == d
    @test_throws ArgumentError Year("1.0")
    @test Year(parse(Float64, "1.0")) == y

    dt = DayDate(2014)
    @test typeof(Year(dt)) <: Year
    @test typeof(Quarter(dt)) <: Quarter
    @test typeof(Month(dt)) <: Month
    @test typeof(Week(dt)) <: Week
    @test typeof(Day(dt)) <: Day
end
@testset "Default values" begin
    @test ExtendedDates.default(Year) == y
    @test ExtendedDates.default(Quarter) == q
    @test ExtendedDates.default(Month) == m
    @test ExtendedDates.default(Week) == w
    @test ExtendedDates.default(Day) == d
end
@testset "Conversions" begin
    @test ExtendedDates.days(d) == 1
    @test ExtendedDates.days(w) == 7
end
@testset "issue #9214" begin
    @test emptyperiod == ((d + y) - y) - d == ((d + y) - d) - y
    @test string(emptyperiod) == "empty period"
    @test y - m == 11m
end
@testset "compound periods and types" begin
    # compound periods should avoid automatically converting period types
    @test sprint(show, y + m) == string(y + m)
    @test convert(Dates.CompoundPeriod, y) + m == y + m
    if VERSION >= v"1.7"
        @test Dates.periods(convert(Dates.CompoundPeriod, y)) ==
            convert(Dates.CompoundPeriod, y).periods
    end
end
@testset "compound period simplification" begin
    # reduce compound periods into the most basic form
    @test Dates.canonicalize(-y + d).periods == Period[-y, d]
    @test Dates.canonicalize(-y + m - w + d).periods == Period[-11m, -6d]
    @test DayDate(2009, 2, 1) - (Month(1) + Day(1)) == DayDate(2008, 12, 31)
    @test_throws MethodError (Month(1) + Day(1)) - DayDate(2009, 2, 1)
end

@testset "Dates.canonicalize Period" begin
    # reduce individual Period into most basic Dates.CompoundPeriod
    @test Dates.canonicalize(Day(7)) == Dates.canonicalize(Week(1))
    @test Dates.canonicalize(Month(12)) == Dates.canonicalize(Year(1))
    Dates.canonicalize(Dates.CompoundPeriod([Day(1), Hour(12)]))
end
@testset "unary ops and vectorized period arithmetic" begin
    pa = [1y 1m 1w 1d]
    cpa = [1y 1m 1w 1d]

    @test +pa == pa == -(-pa)
    @test -pa == map(-, pa)
    @test 1y .+ pa == [2y 1y + 1m 1y + 1w 1y + 1d]
    @test (1y + 1m) .+ pa == [2y + 1m 1y + 2m 1y + 1m + 1w 1y + 1m + 1d]
    @test pa .+ 1y == [2y 1y + 1m 1y + 1w 1y + 1d]
    @test pa .+ (1y + 1m) == [2y + 1m 1y + 2m 1y + 1m + 1w 1y + 1m + 1d]

    @test 1y .+ cpa == [2y 1y + 1m 1y + 1w 1y + 1d]
    @test (1y + 1m) .+ cpa == [2y + 1m 1y + 2m 1y + 1m + 1w 1y + 1m + 1d]
    @test cpa .+ 1y == [2y 1y + 1m 1y + 1w 1y + 1d]
    @test cpa .+ (1y + 1m) == [2y + 1m 1y + 2m 1y + 1m + 1w 1y + 1m + 1d]

    @test 1y .+ pa == [2y 1y + 1m 1y + 1w 1y + 1d]
    @test (1y + 1m) .+ pa == [2y + 1m 1y + 2m 1y + 1m + 1w 1y + 1m + 1d]
    @test pa .+ 1y == [2y 1y + 1m 1y + 1w 1y + 1d]
    @test pa .+ (1y + 1m) == [2y + 1m 1y + 2m 1y + 1m + 1w 1y + 1m + 1d]

    @test 1y .+ cpa == [2y 1y + 1m 1y + 1w 1y + 1d]
    @test (1y + 1m) .+ cpa == [2y + 1m 1y + 2m 1y + 1m + 1w 1y + 1m + 1d]
    @test cpa .+ 1y == [2y 1y + 1m 1y + 1w 1y + 1d]
    @test cpa .+ (1y + 1m) == [2y + 1m 1y + 2m 1y + 1m + 1w 1y + 1m + 1d]

    @test 1y .- pa == [0y 1y - 1m 1y - 1w 1y - 1d]
    @test (1y + 1m) .- pa == [1m 1y 1y + 1m - 1w 1y + 1m - 1d]
    @test pa .- (1y + 1m) == [-1m -1y -1y - 1m + 1w -1y - 1m + 1d]
    @test pa .- 1y == [0y 1m - 1y -1y + 1w -1y + 1d]

    @test 1y .- cpa == [0y 1y - 1m 1y - 1w 1y - 1d]
    @test (1y + 1m) .- cpa == [1m 1y 1y + 1m - 1w 1y + 1m - 1d]
    @test cpa .- 1y == [0y -1y + 1m -1y + 1w -1y + 1d]
    @test cpa .- (1y + 1m) == [-1m -1y -1y - 1m + 1w -1y - 1m + 1d]

    @test [1y 1m; 1w 1d] + [1y 1y; 1y 1y] == [2y 1y+1m; 1y+1w 1y+1d]
    @test [1y 1m; 1w 1d] - [1y 1y; 1y 1y] == [0y 1m-1y; 1w-1y 1d-1y]
end
@testset "Equality and hashing between OtherPeriod types" begin
    for x in (0, 1, 235, -4677, 15250)
        local x, y, z
        y = Year(x)
        z = convert(Month, y)
        @test y == z
        @test hash(y) == hash(z)

        y = Quarter(x)
        z = convert(Month, y)
        @test y == z
        @test hash(y) == hash(z)

        y = Year(x)
        z = convert(Quarter, y)
        @test y == z
        @test hash(y) == hash(z)
    end
end
@testset "Equality and hashing between FixedPeriod/OtherPeriod/Dates.CompoundPeriod (#37459)" begin
    function test_hash_equality(x, y)
        @test x == y
        @test y == x
        @test isequal(x, y)
        @test isequal(y, x)
        @test hash(x) == hash(y)
    end
    for FP in (Week, Day)
        for OP in (Year, Quarter, Month)
            test_hash_equality(FP(0), OP(0))
        end
    end
end

@testset "CompoundPeriod and Period isless()" begin
    #tests for allowed comparisons
    #OtherPeriod
    @test (2y - m < 25m + 1y) == true
    @test (2y < 25m + 1y) == true
    @test (25m + 1y < 2y) == false
end

if VERSION >= v"1.7"
    @testset "Convert CompoundPeriod to Period" begin
        @test convert(Month, Year(1) + Month(1)) === Month(13)
    end
end

end


text:
code_dic = Dict()
code = []
open("day09/day09.in") do io
    line_num = 1
    global instructions
    for line in readlines(io)
        num = parse(Int64, line)
        code_dic[num] = line_num
        line_num += 1
        push!(code, num)
    end
end

check = 25
target_num = 0
for num_pos in (check+1):length(code)
    global target_num
    num = code[num_pos]
    num_check = false
    for pos_back in 1:check
        pair_num = num - code[num_pos-pos_back]
        if get(code_dic, pair_num, 0) != 0 && (num_pos-pos_back <= code_dic[pair_num] <= num_pos-1)
            num_check = true
        end
    end
    if !num_check
        target_num = num
    end
end

code_cumsum = vcat([0], cumsum(code))
set_found = false
start_pos, end_pos = 2, 3
while !set_found
    global start_pos, end_pos, set_found
    set_sum = code_cumsum[end_pos] - code_cumsum[start_pos-1]
    while set_sum < target_num
        end_pos += 1
        set_sum = code_cumsum[end_pos] - code_cumsum[start_pos-1]
    end
    if set_sum == target_num
        println(minimum(code[start_pos-1:end_pos-1])+maximum(code[start_pos-1:end_pos-1]))
        set_found = true
    end
    start_pos += 1
    end_pos = start_pos+1
end

max_stars_repo_name: detrin/advent-of-code

max_issues_repo_name: detrin/advent-of-code

max_forks_repo_name: detrin/advent-of-code


text:
#=
Created on Wednesday 18 March 2020
Last update: Tuesday 23 March 2020

@author: Michiel Stock
michielfmstock@gmail.com

This is an illustration of solving the Krontrovich formulation of the optimal
transport problem. It also demonstrates the use of Convex.jl, a package for
convex optimization problems.
=#

using Convex, SCS

# first give the problem parameters, we use the dessert problem

# portions per person
a = [3.0, 3, 3, 4, 2, 2, 2, 1]
# quantities of each dessert
b = [4.0, 2, 6, 4, 4]

# should both have the same sum
@assert sum(a) ≈ sum(b)

preferences = [2 2 1 0 0;
              0 -2 -2 -2 2;
              1 2 2 2 -1;
              2 1 0 1 -1;
              0.5 2 2 1 0;
              0 1 1 1 -1;
             -2 2 2 1 1;
              2 1 2 1 -1]

C = -preferences

n, m = size(C)

P = Variable(n, m)

problem = minimize(sum(P .* C),  # objective
                    P >= 0,  # non-neg constraints
                    P * ones(m) == a,  # row marginals
                    P' * ones(n) == b)  # column marginals

solve!(problem, () -> SCS.Optimizer(verbose=false))

pstar = problem.optval  # best objective
Pstar = P.value  # distribution

using StatsPlots

groupedbar(Pstar, bar_position = :stack, title="Optimal distribution desserts\nobjective = $pstar")

savefig("figures/kantovorisc.png")


text:
@testset "Misc/MPolyAbsFact" begin

  f = Hecke.MPolyFact.example(wildanger_field(3, 13)[1], 3, 5)
  @test length(factor_absolute(f)) >= 2

  f = Hecke.MPolyFact.example(cyclotomic_field(4)[1], 3, 5)
  @test length(factor_absolute(f)) >= 2

  @test length(factor(f, ComplexField(40)).fac) >= 2
  @test length(factor(f, ArbField(40)).fac) >= 1
end