text:
function parentind(ind::Symbol, op::Operation)
    for (id,opp) ∈ enumerate(parents(op))
        name(opp) === ind && return id
    end
    -1
end
function symbolind(ind::Symbol, op::Operation, td::UnrollArgs)
    id = parentind(ind, op)
    id == -1 && return ind
    @unpack u₁, u₁loopsym, u₂loopsym, suffix = td
    parent = parents(op)[id]
    pvar = if u₂loopsym ∈ loopdependencies(parent)
        variable_name(parent, suffix)
    else
        mangledvar(parent)
    end
    ex = u₁loopsym ∈ loopdependencies(parent) ? Symbol(pvar, u₁) : pvar
    Expr(:call, lv(:staticm1), ex)
end

staticexpr(x::Integer) = Expr(:call, Expr(:curly, lv(:Static), convert(Int, x)))
staticexpr(x) = Expr(:call, lv(:Static), x)
maybestatic(x::Integer) = staticexpr(x)
maybestatic(x) = x
_MMind(ind) = Expr(:call, lv(:MM), VECTORWIDTHSYMBOL, ind)
_MMind(ind::Integer) = Expr(:call, lv(:MM), VECTORWIDTHSYMBOL, staticexpr(ind))
function addoffset!(ret::Expr, ex, offset::Integer, _mm::Bool = false)
    if iszero(offset)
        if _mm
            ind = _MMind(ex)
        else
            push!(ret.args, ex)
            return
        end
    elseif _mm
        ind = _MMind(Expr(:call, lv(:vadd_fast), ex, staticexpr(offset)))
    else
        ind = Expr(:call, lv(:vadd_fast), ex, staticexpr(offset))
    end
    push!(ret.args, ind)
    nothing
end
function addoffset!(ret::Expr, offset::Int, _mm::Bool = false)
    push!(ret.args, _mm ? _MMind(offset) : staticexpr(offset))
    nothing
end

"""
unrolled loads are calculated as offsets with respect to an initial gesp. This has proven important to helping LLVM generate efficient code in some cases.
Therefore, unrolled === true results in inds being ignored.
_mm means to insert `mm`s.
"""
function mem_offset(op::Operation, td::UnrollArgs, inds_calc_by_ptr_offset::Vector{Bool}, _mm::Bool = true)
    # @assert accesses_memory(op) "Computing memory offset only makes sense for operations that access memory."
    ret = Expr(:tuple)
    indices = getindicesonly(op)
    offsets = getoffsets(op)
    loopedindex = op.ref.loopedindex
    # inds_calc_by_ptr_offset = indices_calculated_by_pointer_offsets(ls, op.ref)
    @unpack vectorized = td
    for (n,ind) ∈ enumerate(indices)
        offset = offsets[n] % Int
        # if ind isa Int # impossible
            # push!(ret.args, ind + offset)
        # else
        if loopedindex[n]
            if inds_calc_by_ptr_offset[n]
                addoffset!(ret, offset, _mm & (ind === vectorized))
            else
                addoffset!(ret, ind, offset, _mm & (ind === vectorized))
            end
        else
            addoffset!(ret, symbolind(ind, op, td), offset)
        end
    end
    ret
end

function add_vectorized_offset!(ret::Expr, ind, offset, incr)
    if isone(incr)
        if iszero(offset)
            push!(ret.args, _MMind(Expr(:call, lv(:vadd_fast), VECTORWIDTHSYMBOL, maybestatic(ind))))
        else
            push!(ret.args, _MMind(Expr(:call, lv(:vadd_fast), ind, Expr(:call, lv(:vadd_fast), VECTORWIDTHSYMBOL, staticexpr(offset)))))
        end
    elseif iszero(incr)
        if iszero(offset)
            push!(ret.args, _MMind(ind))
        else
            addoffset!(ret, ind, offset, true)
        end
    elseif iszero(offset)
        push!(ret.args, _MMind(Expr(:call, lv(:vadd_fast), Expr(:call, lv(:vmul_fast), VECTORWIDTHSYMBOL, maybestatic(incr)), maybestatic(ind))))
    else
        push!(ret.args, _MMind(Expr(:call, lv(:vadd_fast), ind, Expr(:call, lv(:vadd_fast), Expr(:call, lv(:vmul_fast), VECTORWIDTHSYMBOL, maybestatic(incr)), staticexpr(offset)))))
    end
end
function add_vectorized_offset_unrolled!(ret::Expr, offset, incr)
    if isone(incr)
        if iszero(offset)
            push!(ret.args, _MMind(Expr(:call, lv(:Static), VECTORWIDTHSYMBOL)))
        else
            push!(ret.args, _MMind(Expr(:call, lv(:vadd_fast), VECTORWIDTHSYMBOL, staticexpr(offset))))
        end
    elseif iszero(incr)
        if iszero(offset)
            push!(ret.args, _MMind(Expr(:call, lv(:Zero))))
        else
            push!(ret.args, _MMind(staticexpr(offset)))
        end
    elseif iszero(offset)
        push!(ret.args, _MMind(Expr(:call, lv(:*), VECTORWIDTHSYMBOL, maybestatic(incr))))
    else
        push!(ret.args, _MMind(Expr(:call, lv(:vadd_fast), Expr(:call, lv(:vmul_fast), VECTORWIDTHSYMBOL, maybestatic(incr)), staticexpr(offset))))
    end
end
function add_vectorized_offset!(ret::Expr, ind, offset, incr, unrolled)
    if unrolled
        add_vectorized_offset_unrolled!(ret, offset, incr)
    else
        add_vectorized_offset!(ret, ind, offset, incr)
    end
end
function mem_offset_u(op::Operation, td::UnrollArgs, inds_calc_by_ptr_offset::Vector{Bool})
    @assert accesses_memory(op) "Computing memory offset only makes sense for operations that access memory."
    @unpack u₁, u₁loopsym, u₂loopsym, vectorized, suffix = td
    incr₁ = u₁
    incr₂ = isnothing(suffix) ? 0 : suffix
    ret = Expr(:tuple)
    indices = getindicesonly(op)
    offsets = getoffsets(op)
    # allbasezero = all(inds_calc_by_ptr_offset) && all(iszero, offsets)
    loopedindex = op.ref.loopedindex
    if iszero(incr₁) & iszero(incr₂)
        return mem_offset(op, td, inds_calc_by_ptr_offset)
        # append_inds!(ret, indices, loopedindex)
    else
        for (n,ind) ∈ enumerate(indices)
            ind_by_offset = inds_calc_by_ptr_offset[n]
            offset = convert(Int, offsets[n])
            indvectorized = ind === vectorized
            # if n == zero_offset
                # addoffset!(ret, 0)
            # elseif ind === u₁loopsym
            if ind === u₁loopsym
                if indvectorized
                    add_vectorized_offset!(ret, ind, offset, incr₁, ind_by_offset)
                elseif ind_by_offset
                    addoffset!(ret, incr₁ + offset)
                else
                    addoffset!(ret, ind, incr₁ + offset)
                end
            elseif ind === u₂loopsym
                if indvectorized
                    add_vectorized_offset!(ret, ind, offset, incr₂, ind_by_offset)
                elseif ind_by_offset
                    addoffset!(ret, incr₂ + offset)
                else
                    addoffset!(ret, ind, incr₂ + offset)
                end
            elseif loopedindex[n]
                if ind_by_offset
                    addoffset!(ret, offset, indvectorized)
                else
                    addoffset!(ret, ind, offset, indvectorized)
                end
            else
                addoffset!(ret, symbolind(ind, op, td), offset)
            end
        end
    end
    ret
end

function varassignname(var::Symbol, u::Int, isunrolled::Bool)
    isunrolled ? Symbol(var, u) : var
end
# name_memoffset only gets called when vectorized
function name_memoffset(var::Symbol, op::Operation, td::UnrollArgs, u₁unrolled::Bool, inds_calc_by_ptr_offset::Vector{Bool})
    @unpack u₁, u₁loopsym, u₂loopsym, suffix = td
    if isnothing(suffix) && u₁ < 0 # sentinel value meaning not unrolled
        name = var
        mo = mem_offset(op, td, inds_calc_by_ptr_offset)
    else
        name = u₁unrolled ? Symbol(var, u₁) : var
        mo = mem_offset_u(op, td, inds_calc_by_ptr_offset)
    end
    name, mo
end

function condvarname_and_unroll(cond, u₁loop, u₂loop, suffix, opu₂)
    if isnothing(suffix) || !opu₂
        condvar, condu₁, condu₂ = variable_name_and_unrolled(cond, u₁loop, u₂loop, nothing)
    else
        condvar, condu₁, condu₂ = variable_name_and_unrolled(cond, u₁loop, u₂loop, suffix)
    end
    condvar, condu₁
end


text:
"""
    create_cell_grid(img, tracks; win) -> grid

Given a xyt image and set of tracks, creates a grid of all tracked cells with
the cells centered in their grid location. Helps to quickly diagnose weird cell
segmentation behavior.
"""
function create_cell_grid(img::AbstractArray, tracks::DataFrame; particle_col=:particle, win=30)
    create_cell_grid(img, [particle for particle in groupby(tracks, particle_col)]; win=win)
end

function create_cell_grid(img::AbstractArray{T, 3}, tracks::AbstractArray{S, 1}; win=30) where {T, S <: SubDataFrame}
    n_tracks = length(tracks)
    nx, ny, nt = size(img)
    cols = ceil(Int, sqrt(n_tracks))
    fontface = newface(normpath(@__DIR__, "..", "assets", "ubuntu-font-family-0.83", "Ubuntu-R.ttf"))
    w = cols*(win*2+2)
    h = ceil(Int, n_tracks / cols)*(win*2+2)
    # the default element should be different for binary masks vs greyscales
    # this picks something intelligent for both
    default_el = T <: Colorant ? maximum(img[:, :, 1]) : oneunit(UInt8)
    grid = fill(default_el, h, w, nt)

    for (pidx, subdf) in enumerate(tracks)
        # calculate the x,y center of this grid location
        xcent, ycent = (pidx-1)÷cols*(win*2+2)+win+1, mod((pidx-1)*(win*2+2)+win+1, w)
        particle = subdf[1, :particle]

        # for each frame in the movie
        for (tidx, xpos, ypos, id) in zip(1:nt, subdf[!, :x], subdf[!, :y], subdf[!, :id])

            # find approximate index of the center of the cell
            xidx, yidx = round(Int, ypos), round(Int, xpos)
            # handle edges intelligently
            xrange, yrange = max(xidx-win, 1):min(xidx+win, nx), max(yidx-win, 1):min(yidx+win, ny)
            # range - idx centers the view on the cell to be +/- the window and then we offset it by
            # our current location in the grid
            xloc, yloc = xrange .- xidx .+ xcent,yrange .- yidx .+ ycent

            if T <: Colorant
                cell = img[xrange, yrange, tidx]
            else # if a color wasn't provided then we need to select only cells that are labeled
                cell = UInt8.(img[xrange, yrange, tidx] .== id)
            end

            txt_offset = win > 20 ? 2 : 0
            renderstring!(cell, "$particle", fontface, (10, 10), txt_offset,
                          txt_offset, halign=:hleft, valign=:vtop, fcolor=default_el)

            grid[xloc, yloc, tidx] .= cell
        end
    end
    grid
end

ext: jl


text:
@inline function getVelocityDelta(joint::Translational1, body1::AbstractBody, body2::Body{T}, v::SVector{2,T}) where T
    Δv = joint.V12' * v # in body1 frame
    return Δv
end

@inline function getPositionDelta(joint::Translational1, body1::AbstractBody, body2::Body{T}, x::SVector{2,T}) where T
    Δx = joint.V12' * x # in body1 frame
    return Δx
end

@inline function setForce!(joint::Translational1, body1::Body, body2::Body{T}, F::SVector{2,T}, No) where T
    q1 = body1.q[No]
    q2 = body2.q[No]

    F1 = vrotate(joint.V12' * -F, q1)
    F2 = -F1

    τ1 = torqueFromForce(F1, vrotate(joint.vertices[1], q1))
    τ2 = torqueFromForce(F2, vrotate(joint.vertices[2], q2))

    body1.F[No] = F1
    body1.τ[No] = τ1

    body2.F[No] = F2
    body2.τ[No] = τ2
    return
end

@inline function setForce!(joint::Translational1, body1::Origin, body2::Body{T}, F::SVector{2,T}, No) where T
    F2 = joint.V12' * F
    τ2 = torqueFromForce(F2, vrotate(joint.vertices[2], body2.q[No]))

    body2.F[No] = F2
    body2.τ[No] = τ2
    return
end

@inline function minimalCoordinates(joint::Translational1, body1::Body, body2::Body, No)
    vertices = joint.vertices
    q1 = body1.q[No]
    joint.V12 * vrotate(body2.x[No] + vrotate(vertices[2], body2.q[No]) - (body1.x[No] + vrotate(vertices[1], q1)), inv(q1))
end

@inline function minimalCoordinates(joint::Translational1, body1::Origin, body2::Body, No)
    vertices = joint.vertices
    joint.V12 * (body2.x[No] + vrotate(vertices[2], body2.q[No]) - vertices[1])
end

@inline function g(joint::Translational1, body1::Body, body2::Body, Δt, No)
    vertices = joint.vertices
    q1 = getq3(body1, Δt)
    joint.V3 * vrotate(getx3(body2, Δt) + vrotate(vertices[2], getq3(body2, Δt)) - (getx3(body1, Δt) + vrotate(vertices[1], q1)), inv(q1))
end

@inline function g(joint::Translational1, body1::Origin, body2::Body, Δt, No)
    vertices = joint.vertices
    joint.V3 * (getx3(body2, Δt) + vrotate(vertices[2], getq3(body2, Δt)) - vertices[1])
end

@inline function ∂g∂posa(joint::Translational1{T}, body1::Body, body2::Body, No) where T
    if body2.id == joint.cid
        q1 = body1.q[No]
        point2 = body2.x[No] + vrotate(joint.vertices[2], body2.q[No])

        X = -joint.V3 * VLᵀmat(q1) * RVᵀmat(q1)
        R = joint.V3 * 2 * VLᵀmat(q1) * (Lmat(Quaternion(point2)) - Lmat(Quaternion(body1.x[No]))) * LVᵀmat(q1)

        return [X R]
    else
        return ∂g∂posa(joint)
    end
end

@inline function ∂g∂posb(joint::Translational1{T}, body1::Body, body2::Body, No) where T
    if body2.id == joint.cid
        q1 = body1.q[No]
        q2 = body2.q[No]

        X = joint.V3 * VLᵀmat(q1)RVᵀmat(q1)
        R = joint.V3 * 2 * VLᵀmat(q1) * Rmat(q1) * Rᵀmat(q2) * Rmat(Quaternion(joint.vertices[2])) * LVᵀmat(q2)

        return [X R]
    else
        return ∂g∂posb(joint)
    end
end

@inline function ∂g∂vela(joint::Translational1{T}, body1::Body, body2::Body, Δt, No) where T
    if body2.id == joint.cid
        q1 = body1.q[No]
        ωbar1 = ωbar(body1, Δt)
        point2 = body2.x[No] + Δt * getvnew(body2) + vrotate(vrotate(joint.vertices[2], ωbar(body2, Δt)), body2.q[No])

        V = -Δt * joint.V3 * VLᵀmat(ωbar1)Lᵀmat(q1)Rmat(ωbar1)RVᵀmat(q1)
        Ω = 2 * joint.V3 * VLᵀmat(ωbar1) * Lᵀmat(q1) * (Lmat(Quaternion(point2)) - Lmat(Quaternion(body1.x[No] + Δt * getvnew(body1)))) * Lmat(q1) * derivωbar(body1, Δt)

        return [V Ω]
    else
        return ∂g∂vela(joint)
    end
end

@inline function ∂g∂velb(joint::Translational1{T}, body1::Body, body2::Body, Δt, No) where T
    if body2.id == joint.cid
        q1 = body1.q[No]
        q2 = body2.q[No]
        ωbar1 = ωbar(body1, Δt)

        V = Δt * joint.V3 * VLᵀmat(ωbar1)Lᵀmat(q1)Rmat(ωbar1)RVᵀmat(q1)
        Ω = 2 * joint.V3 * VLᵀmat(ωbar1) * Lᵀmat(q1) * Lmat(q2) * Rmat(ωbar1) * Rmat(q1) * Rᵀmat(q2) * Rᵀmat(ωbar(body2, Δt)) * Rmat(Quaternion(joint.vertices[2])) * derivωbar(body2, Δt)

        return [V Ω]
    else
        return ∂g∂velb(joint)
    end
end

@inline function ∂g∂posb(joint::Translational1{T}, body1::Origin, body2::Body, No) where T
    if body2.id == joint.cid
        q2 = body2.q[No]

        X = joint.V3
        R = joint.V3 * 2 * VRᵀmat(q2) * Rmat(Quaternion(joint.vertices[2])) * LVᵀmat(q2)

        return [X R]
    else
        return ∂g∂posb(joint)
    end
end

@inline function ∂g∂velb(joint::Translational1{T}, body1::Origin, body2::Body, Δt, No) where T
    if body2.id == joint.cid
        q2 = body2.q[No]

        V = Δt * joint.V3
        Ω = 2 * joint.V3 * VLmat(q2) * Rᵀmat(q2) * Rᵀmat(ωbar(body2, Δt)) * Rmat(Quaternion(joint.vertices[2])) * derivωbar(body2, Δt)

        return [V Ω]
    else
        return ∂g∂velb(joint)
    end
end


text:
mutable struct CCSAState{T<:AbstractFloat}
    n::Integer # number of variables > 0
    m::Integer # number of inequality constraints ≥ 0
    lb::AbstractVector{T} # n lower bounds
    ub::AbstractVector{T} # n upper bounds
    f_and_∇f::Function # f(x) = (m+1, (m+1) x n linear operator)
    ρ::AbstractVector{T} # m + 1 penality weight
    σ::AbstractVector{T} # n radius of trust region
    x::AbstractVector{T} # current best feasible point
    xtol_rel::T # relative tolerence
    max_iters::Integer # max number of iterations
    iters::Integer # interation count

    fx::AbstractVector{T} # (m+1) x 1 function values at x
    ∇fx::AbstractVecOrMat{T} # (m+1) x n linear operator of Jacobian at x
    a::AbstractVector{T} # n
    b::AbstractVector{T} # n
    Δx::AbstractVector{T} # n xᵏ⁺¹ - xᵏ
    Δx_last::AbstractVector{T} # n xᵏ - xᵏ⁻¹
    gλ::T # Lagrange dual function value
    ∇gλ::AbstractVector{T} # m Lagrange dual function gradient
    dual::CCSAState{T} # Lagrange dual problem if the primal problem has constraints
    function CCSAState(
        n::Integer, # number of variables
        m::Integer, # number of inequality constraints
        f_and_∇f::Function,
        x₀::AbstractVector{T}; # initial feasible point
        ρ::AbstractVector{T}=ones(T, m + 1), # (m + 1) penality weight ρ
        σ::AbstractVector{T}=ones(T, n), # n radius of trust region σ
        lb::AbstractVector{T}=fill(typemin(T), n), # lower bounds, default -Inf
        ub::AbstractVector{T}=fill(typemax(T), n), # upper bounds, default Inf
        xtol_rel::T=T(1e-5), # relative tolerence
        max_iters::Integer=typemax(Int64) # max number of iterations, defalut 2⁶³ - 1
    ) where {T<:AbstractFloat}
        fx, ∇fx = f_and_∇f(x₀)
        opt = new{T}(
            n,
            m,
            lb,
            ub,
            f_and_∇f,
            ρ,
            σ,
            x₀,
            xtol_rel,
            max_iters,
            zero(max_iters),
            fx,
            ∇fx,
            Vector{T}(undef, n),
            Vector{T}(undef, n),
            Vector{T}(undef, n),
            Vector{T}(undef, n),
            T(0),
            Vector{T}(undef, m + 1)
        )
        if opt.m > 0
            opt.dual = CCSAState( # Lagrange dual problem
                opt.m, # number of Lagrange multipliers
                0, # no inequality constraints
                λ -> (T[], T[]),
                ones(T, opt.m), # initial feasible point for Lagrange multipliers
                lb=zeros(T, opt.m) # Lagrange multipliers ≥ 0 for inequality constraints
            )
        end
        return opt
    end
end
# Returns the dual function g(λ) and ∇g(λ)
function dual_func!(λ::AbstractVector{T}, st::CCSAState{T}) where {T}
    λ_all = CatView([one(T)], λ)
    st.a .= dot(st.ρ, λ_all) ./ (2 .* st.σ .^ 2)
    mul!(st.b, st.∇fx', λ_all)
    for j in 1:st.n
        st.Δx[j] = clamp(-st.b[j] / (2 * st.a[j]), -st.σ[j], st.σ[j])
        st.Δx[j] = clamp(st.Δx[j], st.lb[j] - st.x[j], st.ub[j] - st.x[j])
    end
    st.gλ = dot(λ_all, st.fx) + sum(@. st.a * st.Δx^2 + st.b * st.Δx)
    st.∇gλ .= st.fx .+ sum(abs2, st.Δx ./ st.σ) / 2 .* st.ρ
    mul!(st.∇gλ, st.∇fx, st.Δx, true, true)
    return [st.gλ], (@view st.∇gλ[2:end, :])'
end

# optimize problem with no constraint
function optimize_simple(opt::CCSAState{T}) where {T}
    monotonic = BitVector(undef, opt.n)
    while opt.iters < opt.max_iters
        opt.fx, opt.∇fx = opt.f_and_∇f(opt.x)
        opt.a .= opt.ρ[1] / 2 ./ (opt.σ) .^ 2
        opt.b .= opt.∇fx'
        while true
            @. opt.Δx = clamp(-opt.b / (2 * opt.a), -opt.σ, opt.σ)
            @. opt.Δx = clamp(opt.Δx, opt.lb - opt.x, opt.ub - opt.x)
            opt.gλ = opt.fx[1] + sum(@. opt.a * opt.Δx^2 + opt.b * opt.Δx)
            if opt.gλ ≥ opt.f_and_∇f(opt.x + opt.Δx)[1][1] # check conservative
                break
            end
            opt.ρ *= 2
            opt.a *= 2
        end
        opt.ρ /= 2
        monotonic .= signbit.(opt.Δx_last) .== signbit.(opt.Δx) # signbit avoid multiplication
        opt.σ[monotonic] *= 2 # double trust region if xⱼ moves monotomically
        opt.σ[.!monotonic] /= 2 # shrink trust region if xⱼ oscillates
        opt.x .+= opt.Δx
        opt.Δx_last .= opt.Δx
        if norm(opt.Δx, Inf) < opt.xtol_rel
            return
        end
    end
end

function inner_iterations(opt::CCSAState{T}) where {T}
    gᵢxᵏ⁺¹ = Vector{T}(undef, opt.m + 1)
    conservative = BitVector(undef, opt.m + 1)
    while true
        opt.fx, opt.∇fx = opt.f_and_∇f(opt.x)
        opt.dual.f_and_∇f = function (λ) # negative Lagrange dual function and gradient
            gλ, ∇gλ = dual_func!(λ, opt)
            -gλ, -∇gλ # minus signs used to change max problem to min problem
        end
        opt.dual.fx, opt.dual.∇fx = opt.dual.f_and_∇f(opt.dual.x)
        optimize_simple(opt.dual)
        #dual_func!(opt.dual.x , opt)# optimal solution of Lagrange dual problem
        gᵢxᵏ⁺¹ .= opt.fx .+ sum(abs2, opt.Δx ./ opt.σ) / 2 .* opt.ρ
        mul!(gᵢxᵏ⁺¹, opt.∇fx, opt.Δx, true, true)
        fᵢxᵏ⁺¹ = opt.f_and_∇f(opt.x + opt.Δx)[1]
        conservative .= gᵢxᵏ⁺¹ .≥ fᵢxᵏ⁺¹
        if all(conservative)
            return
        end
        opt.ρ[.!conservative] *= 2 # increase ρ until achieving conservative approxmation
        opt.dual.ρ .= one(T) # reinitialize penality weights
        opt.dual.σ .= one(T) # reinitialize radii of trust region
        opt.dual.x .= one(T) # reinitialize starting point of Lagrange multipliers
    end
end

function optimize(opt::CCSAState{T}; callback=nothing) where {T}
    if opt.m == 0
        return optimize_simple(opt)
    end
    monotonic = BitVector(undef, opt.n)
    while opt.iters < opt.max_iters
        inner_iterations(opt)
        monotonic .= signbit.(opt.Δx_last) .== signbit.(opt.Δx)
        opt.σ[monotonic] *= 2 # double trust region if xⱼ moves monotomically
        opt.σ[.!monotonic] /= 2 # shrink trust region if xⱼ oscillates
        opt.ρ /= 2
        opt.x .+= opt.Δx
        if norm(opt.Δx, Inf) < opt.xtol_rel
            return
        end
        opt.Δx_last .= opt.Δx
        opt.iters += 1
        if callback !== nothing
            callback()
        end
    end
    return iters
end


text:
using PtFEM

data = Dict(
  # Plane(ndim, nst, nxe, nye, nip, direction, finite_element(nod, nodof), axisymmetric)
  :struc_el => Plane(2, 3, 2, 1, 12, :y, Triangle(15, 2), false),
  :properties => [1.0e5 0.2;],
  :x_coords => [0.0,  1.0,  6.0],
  :y_coords => [0.0,  -2.0],
  :support => [
    (1, [0 1]),
    (2, [0 1]),
    (3, [0 1]),
    (4, [0 1]),
    (5, [0 0]),
    (10, [0 0]),
    (15, [0 0]),
    (20, [0 0]),
    (25, [0 0]),
    (30, [0 0]),
    (35, [0 0]),
    (40, [0 0]),
    (41, [0 1]),
    (42, [0 1]),
    (43, [0 1]),
    (44, [0 1]),
    (45, [0 0])
    ],
  :loaded_nodes => [
    (1, [0.0 -7/90]),
    (6, [0.0 -16/45]),
    (11, [0.0 -2/15]),
    (16, [0.0 -16/45]),
    (21, [0.0 -7/90])
    ]
)

data |> display
println()

@time m, dis_df, fm_df = p51(data)
println()

dis_df |> display
println()

println("\nThe integration point (nip = $(data[:struc_el].nip)) stresses are:")
fm_df |> display
println()


text:
#############################################################
# Running Simulations
#############################################################

#Include KG_Module to run simulation
include("../src/KG_Module.jl")

p0 = Param(

    #Temporal variables
    tmax           = 1000.0     ,
    out_every_t    = 0.1        ,
    deltat         = 0.05       ,

    #Spacial variables
    xmin           =  -30.0     ,
    xmax           =  30.0      ,
    xnodes         =  128       ,

    #Spacial variables
    ymin           =  -30.0     ,
    ymax           =  30.0      ,
    ynodes         =   128      ,

    #Initial Conditions
    A0              =   3.5     ,
    σ               =   1.5     ,
    r0              =   17.0    ,
    ω               =   0.5     ,
    m               =   2       ,

    #Klein Gordon Mass
    μ               =   0.0     ,

    #Dissipation terms
    dissipation     =  :true    ,

    R_orbit         = 0.0       ,
    R1              = 1.5       ,
    R2              = 0.0       ,

    Ω               = 0.45      , #(3)^(-3/2) ,      ,
    alpha1          = 10.0      ,
    alpha2          = 0.0       ,

    #Periodic BC?
    Boundaries      = :radial   , # :square :periodic

    #Data variables
    folder          = "/test_runs/"  ,       #must have a bar before the folder
    fname           = "run_004"

)

solve_wave_equation_2D(p0)

#############################################################
# Analysis
#############################################################

#Include plot utils
include("../src/KG_Plot_Module.jl")

#name of file to analyse
#filename = "./test_runs/run_003.h5"

#Slices Observer (also works with filename)
slices_observer(p0)

#Generate animation
#generate_animation_heatmap_makie2( filename , "test_animation.mp4")


text:
@testset "Multivariate Basic Structural Model" begin
    # Multivariate
    front_rear_df = CSV.read(StateSpaceModels.FRONT_REAR_SEAT_KSI, DataFrame)
    front_rear = [front_rear_df.front_seat_KSI front_rear_df.rear_seats_KSI]

    model = MultivariateBasicStructural(front_rear, 12)
    fit!(model)
    forec = forecast(model, 10)
    @test monotone_forecast_variance(forec)
end


text:
# Kmers.jl
# ========
#
# Module for kmers and short sequence fragments.
#
# This file is a part of the Kmers.jl, a package in the BioJulia ecosystem.
# License is MIT: https://github.com/BioJulia/Kmers.jl/blob/master/LICENSE

__precompile__()

module Kmers

export
    # BioSymbols re-exports.
    NucleicAcid,
    DNA,
    RNA,
    DNA_A,
    DNA_C,
    DNA_G,
    DNA_T,
    DNA_M,
    DNA_R,
    DNA_W,
    DNA_S,
    DNA_Y,
    DNA_K,
    DNA_V,
    DNA_H,
    DNA_D,
    DNA_B,
    DNA_N,
    DNA_Gap,
    ACGT,
    ACGTN,
    RNA_A,
    RNA_C,
    RNA_G,
    RNA_U,
    RNA_M,
    RNA_R,
    RNA_W,
    RNA_S,
    RNA_Y,
    RNA_K,
    RNA_V,
    RNA_H,
    RNA_D,
    RNA_B,
    RNA_N,
    RNA_Gap,
    ACGU,
    ACGUN,
    AminoAcid,
    AA_A,
    AA_R,
    AA_N,
    AA_D,
    AA_C,
    AA_Q,
    AA_E,
    AA_G,
    AA_H,
    AA_I,
    AA_L,
    AA_K,
    AA_M,
    AA_F,
    AA_P,
    AA_S,
    AA_T,
    AA_W,
    AA_Y,
    AA_V,
    AA_O,
    AA_U,
    AA_B,
    AA_J,
    AA_Z,
    AA_X,
    AA_Term,
    AA_Gap,

    # BioSequences re-exports
    Alphabet,
    BioSequence,
    NucleicAcidAlphabet,
    AminoAcidAlphabet,
    DNAAlphabet,
    RNAAlphabet,

    ###
    ### Mers
    ###

    # Type & aliases
    Kmer,
    DNAKmer,
    DNAKmer27,
    DNAKmer31,
    DNAKmer63,
    RNAKmer,
    RNAKmer27,
    RNAKmer31,
    RNAKmer63,
    AAKmer,
    DNACodon,
    RNACodon,

    # Iteration
    fw_neighbors,
    bw_neighbors,

    ###
    ### Sequence literals
    ###

    @mer_str,
    @bigmer_str

using BioSequences

include("kmer.jl")

end # module


text:
#Run sample code
#-------------------------------------------------------------------------------

module CMDimCircuits_Sample

#==Helper functions
===============================================================================#
function printheader(title)
	println("\n", title, "\n", repeat("-", 80))
end

filelist = ["demo_cktanalysis.jl",
	"demo_psf_r.jl", "demo_snp_rw.jl", "demo_tr0_r.jl",
]

for file in filelist
	printheader("Executing $file...")
	evalfile(file)
end

end

:SampleCode_Executed


text:
#
# Query composition.
#

function ComposeQuery(qs::Vector{Query})
    if length(qs) >= 2
        lq = qs[1]
        for rq in qs[2:end]
            @assert fits(output(lq), input(rq)) "($lq) >> ($rq)"
            lq = rq
        end
    end
    ity =
        if isempty(qs)
            Input(Any)
        else
            Input(
                domain(input(qs[1])),
                ibound([mode(input(q)) for q in qs]))
        end
    oty =
        if isempty(qs)
            Output(Any)
        else
            Output(
                domain(output(qs[end])),
                obound([mode(output(q)) for q in qs]))
        end
    return Query(ComposeSig(), qs, ity, oty)
end

>>(q1::Query, q2::Query) = ComposeQuery([q1, q2])

immutable ComposeSig <: AbstractSignature
end

function ev(::ComposeSig, args::Vector{Query}, ::Input, oty::Output, iflow::InputFlow)
    oflow = OutputFlow(
            domain(iflow),
            PlainColumn(values(iflow)))
    first = true
    for arg in args
        iflow′ =
            if first
                narrow(iflow, input(arg))
            else
                distribute(narrow(iflow, input(arg)), oflow)
            end
        oflow′ = ev(arg, iflow′)
        oflow =
            if first
                oflow′
            else
                optional = isoptional(oflow) || isoptional(oflow′)
                plural = isplural(oflow) || isplural(oflow′)
                OutputFlow(
                    output(oflow′),
                    Column{optional,plural}(
                        compose_impl(offsets(oflow), offsets(oflow′)),
                        values(oflow′)))
            end
        first = false
    end
    return OutputFlow(oty, column(oflow))
end

compose_impl(offs1::AbstractVector{Int}, offs2::AbstractVector{Int}) =
    Int[offs2[off] for off in offs1]

compose_impl(offs1::OneTo, offs2::OneTo) = offs1

compose_impl(offs1::OneTo, offs2::AbstractVector{Int}) = offs2

compose_impl(offs1::AbstractVector{Int}, offs2::OneTo) = offs1

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
export move_particle!, equilibrate_diagonal!

"""
    move_particle!(m::Model, e::Ensemble, c::Configuration) :: Tuple{Float64, Step}
Update that removes occupation in one kinkfree Orbital and adds it to another.
"""
function move_particle!(m::Model, e::Ensemble, c::Configuration) :: Tuple{Float64, Step}
    free_orbitals = filter(x -> isunaffected(c.kinks, x), c.occupations)
    if isempty(free_orbitals)
        return 1.0, Step()
    else
        x = rand( free_orbitals )
    end
    oe = filter(x -> isunaffected(c.kinks,x), setdiff!(sphere_with_same_spin(x, dk = ex_radius), c.occupations))

    #if there are no empty non interacting orbitals in neighbourhood make no change
    if isempty(oe)
        return 1.0, Step()
    end
    y = rand(oe)
    @assert x != y "same Configuration proposed."

    delta_di = ΔWdiag_element(m, e, c, y, x, ImgTime(0), ImgTime(1))
    @assert delta_di != Inf
    # weight change
    dw = exp(-(e.β*(energy(m,y)-energy(m,x)) + e.β*delta_di))

    # MC Step generated by this update

    Δ = Step((x,), (y,))
    apply_step!(c,Δ)

    # get orbitals for reverse update
    oe2 = filter(x -> isunaffected(c.kinks,x), setdiff!(sphere_with_same_spin(y, dk = ex_radius), c.occupations) )

    # quotient of proposal probabilities
    dv = length(oe)/length(oe2)
    @assert dv >= 0
    @assert dw >= 0 "dw=$(dw), x=$(x), y=$(y), delta_di=$(delta_di)"
    @assert !isinf(abs(dv*dw))
    return dv*dw, Δ
end

"""
    equilibrate_diagonal(m::Model, e::Ensemble, c::Configuration)
This function will propose 100 type-a-updates for each particle in c::Configuration.
As move_particle! can only be called on kinkfreee orbitals, calling this before starting the true equilibration
might reduce the error of the kinetic energy in runs with lots of kinks.
"""
function equilibrate_diagonal!(m::Model, e::Ensemble, c::Configuration)
    for _ in 1:e.N*1e2
        dv, Δ = move_particle!(m, e, c)
        if rand() > dv
            #deny step
            apply_back_step!(c, Δ)
        end
    end
end


text:
include("model.jl")
include("statistics.jl")
include("transforms.jl")
include("tools/plotting.jl")
include("validation/validate.jl")

#used for testing derivatives
include("computation/derivatives.jl")
include("computation/buffers.jl")

using DataFrames
using CSV
#using StatsBase
using Plots
#using Profile
#using ProfileView

df = DataFrame(CSV.File("datasets/abalone.csv"))
data = convert(Matrix, df[:,2:8]) #length, diameter, height, whole weight, shucked weight, viscera weight, shell weight
target = convert(Array, df[:, 9]) #age
normalizing_constant = maximum(target)
target = target/normalizing_constant #normalization

#pick training points
#ind = 1:30
ind = 1:5
s = data[ind, :]
#choose a subset of variables to be regressors for the mean
X = data[ind, 1:1]
z = float(target[ind])
#normalizing_constant = maximum(z)
#z = z ./ normalizing_constant

#select range_lambda and range_theta
if false
    histogram(boxCoxObj.f.(target, .5), bins=:scott)
end
#define ranges for theta and lambda
rangeθ = [200.0, 700.0]
rangeλ = [.4, .6]

if false #look at eigenspectrum of kernel matrix
    θ = 200
    gm = K(s, s, θ, rbf)
    display(plot(sort(eigvals(gm))))
end

#initialize train and test data structures
i = 120
s0 = data[i:i,:] #covariates and coordinates
X0  = data[i:i, 1:1]
#example = setting(s, s0, X, X0, z)#abalone data
train = trainingData(s, X, z)
test = testingData(s0, X0)
#example2 = getExample(1, 10, 1, 1, 2)
println("Data loaded.")

#experiment to see how large n has to be for Cholesky to have the most expensive computational cost
if false
    data  = repeat(data, 10, 1)
    times2 = zeros(2, 4)
    for i = 10000:10000:40000
        println("iteration: ", i/10000)
        ind = 1:i
        s = data[ind, :]
        tK = @time K = fastK(s, s, 1.2)
        tC = @time cholesky(K)
        times[1, Int64(i/10000)] = tK;
        times[2, Int64(i/10000)] = tC;
    end
end

  #check derivative of p(z0|theta, lambda, z)
if false
    theta_params = funcθ(1.5, train, test, "Gaussian")
    (f, df) = partial_z0(1.5, 2.0, train, test, boxCoxObj, theta_params, "Gaussian")
    (h, A) = checkDerivative(f, df, .6, 9, 15, 10)
    plt1 = Plots.plot(h, A, title = "Finite Difference Derivative Checker", xlabel = "log of h", ylabel = "log of error",fontfamily=font(48, "Courier") , reuse = false)
    #plot(polyfit(h, A, 1), reuse = true)
    println("partial theta of p(z0|theta, lambda, z)")
    println(polyfit(h, A, 1))
    nums = collect(.05:.05:1)
    g = x -> f(x)[1]
    plt2 = Plots.plot(nums, g.(nums),xlabel = "theta", ylabel = "p(z0|theta, lambda, z)", fontfamily=font(48, "Courier") ,title = "theta vs p(z0| theta, lambda, z)")
    display(Plots.plot(plt1, plt2, fontfamily=font(48, "CoSurier")))
    gui()
end

if false # use this blockc to get pdf and cdf, sanity check them, and optionally plot them
    reset_timer!()
    choleskytime = 0

    #(f, g, df) = getBtgDensity(train, test, rangeθ, rangeλ, boxCoxObj, "Gaussian", "Uniform")

    i=43
    ind = [collect(1:i-1); collect(i+1:length(z))]
    train_cur = trainingData(s[ind, :], X[ind, :], z[ind])
    test_cur = testingData(s[i:i, :], X[i:i, :])
    (f, g, df) = getBtgDensity(train_cur, test_cur, rangeθ, rangeλ, boxCoxObj, "Gaussian", "Uniform")

    if true
    println("Plotting...")
    #plt(df, 0.1, 2, 100, "dpdf")
    plt(f, 0.1, 2, 300, "pdf")
    plt!(g, 0.1, 2, 300, "cdf")
    println(median(f, g))
    #plt(f, 0.1, 0.9, 500, "pdf")
    #plt!(g, 0.1, 0.9, 500, "cdf")
    end
end

if true #check derivative of reference function
    f = x->x[1]^2+x[2]
    df = x -> [2*x[1]; 1]
    #f = x -> x[1]^2*x[2]^2 + x[2]^4
    #df = x -> [2*x[1]*x[2]^2 ; 2*x[1]^2*x[2] + 4*x[2]^3]
    (h, A) = checkDerivative(f, df, [1.2, 2.1])
    println(polyfit(h, A, 1))
end

if true #check derivatives of CDF(z, s) w.r.t s
    function phi(s0_new)
        #println("here1")
        test.s0 = s0_new
        #println("here2")
        theta_params = funcθ(1.2, train, test, "Gaussian")
        #println("here3")
        (r1, r2, r3, r4, r5) = partial_s(1.2, 10.5, train, test, boxCoxObj, theta_params, "Gaussian")
    end
    s0_initial = similar(test.s0)
    for i = 1: max(size(s0, 1), size(s0, 2))
        s0_initial[i] = test.s0[i] + 2*(rand()-.5)
    end
    (h, A) = checkDerivative(s0 -> phi(s0)[3], s0 -> phi(s0)[4], s0_initial, nothing,7, 20)
    println(polyfit(h, A, 1))
    plt1 = Plots.plot(h, A, title = "Finite Difference Derivative Checker", xlabel = "log of h", ylabel = "log of error",fontfamily=font(48, "Courier") , reuse = false)
    Plots.display(plt1)
    gui()
end

if false #check derivatives of df, f, or g
    (h, A) = checkDerivative(f, df, .3, nothing, 7, 15, 10)
    plt1 = Plots.plot(h, A, title = "Finite Difference Derivative Checker", xlabel = "log of h", ylabel = "log of error",fontfamily=font(48, "Courier") , reuse = false)
    #plot(polyfit(h, A, 1), reuse = true)
    println("derivative of p(z0|z)")
    println(polyfit(h, A, 1))
    nums = collect(.05:.05:1)
    gg = x -> f(x)[1]
    plt2 = Plots.plot(nums, gg.(nums),xlabel = "theta", ylabel = "p(z0|theta, lambda, z)", fontfamily=font(48, "Courier") ,title = "theta vs p(z0| theta, lambda, z)")
    display(Plots.plot(plt1, plt2, fontfamily=font(48, "CoSurier")))
    gui()
end

if false #cross validation on training set
    Xs, Ys = cross_validate(train, rangeθ, rangeλ, boxCoxObj, "Gaussian", "Uniform")
    z = train.z
   # _, Xs, Ys = cross_validate(X, s, boxCox, boxCoxPrime, pθ, pλ, z, range_theta, range_lambda, 500, 2, 24)
    display(Plots.plot(Xs, Ys,
    layout = length(z),
    legend=nothing,
    xtickfont = Plots.font(4, "Courier"),
    ytickfont = Plots.font(4, "Courier"),
    lw=0.5))
    display(plot!(z', layout = length(z),
    legend=nothing,
    seriestype = :vline,
    xtickfont = Plots.font(4, "Courier"),
    ytickfont = Plots.font(4, "Courier"),
    lw=0.5))
    Plots.savefig("abalone_1.pdf")
end

if false #delete one group cross validation
    Xs, Ys = cross_validate_groups(X, s, boxCox, boxCoxPrime, pθ, pλ, z, rangeθ, rangeλ, 5, 500, 2, 24)
    display(plot(Xs, Ys,
    layout = length(z),
    legend=nothing,
    xtickfont = Plots.font(4, "Courier"),
    ytickfont = Plots.font(4, "Courier"),
    lw=0.5))
    display(plot!(z', layout = length(z),
    legend=nothing,
    seriestype = :vline,
    xtickfont = Plots.font(4, "Courier"),
    ytickfont = Plots.font(4, "Courier"),
    lw=0.5))
    savefig("results//abalone//abalone_group_cross_validation4.pdf")
end

if false #compute loss for Abalone dataset
    graph = false #generate plots for debugging

    println("Computing loss...")
    (med,  _, _) = compute_loss(train, rangeθ, rangeλ, boxCoxObj, "Gaussian", "Uniform", graph, 100, 0, 2.5)
    println("maximum computed median: ", maximum(med))
    println("minimum computed median: ", minimum(med))
    SE = norm(med .- z)^2 * normalizing_constant^2
    println("SE: ", SE)

    if graph == true
        display(Plots.plot(Xs, Ys,
        layout = length(z),
        legend=nothing,
        xtickfont = Plots.font(4, "Courier"),
        ytickfont = Plots.font(4, "Courier"),
        lw=0.5))
        display(Plots.plot!(z', layout = length(z),
        legend=nothing,
        seriestype = :vline,
        xtickfont = Plots.font(4, "Courier"),
        ytickfont = Plots.font(4, "Courier"),
        lw=0.5))
        display(Plots.plot!(med', layout = length(z),
        legend=nothing,
        seriestype = :vline,
        xtickfont = Plots.font(4, "Courier"),
        ytickfont = Plots.font(4, "Courier"),
        lw=0.5))
    end
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
export size, drop_monomorphic, drop_monomorphic!

## experimental and not exported or documented!
function adjacency_matrix(data::PopData)
    data_loci = groupby(data.loci, :locus)
    out_vec = Vector{Array{Int8,2}}(undef, length(data_loci))
    for (j,i) in enumerate(data_loci)
        uniq = unique(skipmissing(i.genotype))
        adj_mat = fill(Int8(0), length(samples(data)), length(uniq))
        for (j,k) in zip(i.genotype, eachrow(adj_mat))
            k .= Ref(j) .=== uniq
        end
        out_vec[j] = adj_mat
    end
    return out_vec
end

#TODO change location in API docs and rename allele_pool?
#TODO replace alleles with universal Symbol type?
"""
    alleles(locus::T) where T<:GenoArray
Return an array of all the non-missing alleles of a locus.
"""
@inline function alleles(locus::T) where T<:GenoArray
    if all(ismissing.(locus))
        int_type = eltype(typeof(locus)) |> nonmissingtype |> eltype
        return Vector{Union{Missing, int_type}}(undef, length(locus))
    end
    alle_out = Base.Iterators.flatten(skipmissing(locus)) |> collect
end

"""
    allele_count(locus::T) where T<:GenoArray
Return the number of unique alleles present at a locus.
"""
@inline function allele_count(locus::T) where T<:GenoArray
    Base.Iterators.flatten(skipmissing(locus)) |> collect |> unique |> length
end

"""
    alleles(locus::T, miss::Bool = false) where T<:GenoArray
Return an array of all the non-missing alleles of a locus. Use the second positional
argument as `true` to include missing values.
"""
@inline function alleles(locus::T, miss::Bool) where T<:GenoArray
    int_type = eltype(typeof(locus)) |> nonmissingtype |> eltype
    if all(ismissing.(locus))
        return Vector{Union{Missing, int_type}}(undef, length(locus))
    end
    alle_out = Vector{Union{Missing, int_type}}(Base.Iterators.flatten(skipmissing(locus)) |> collect)
    if miss == true
        append!(alle_out, locus[locus .=== missing])
    end
    return alle_out
end

function Base.size(data::PopData)
    return (samples = size(data.meta)[1], loci = length(loci(data)))
end

"""
    unique_alleles(locus::T) where T<:GenoArray
Return an array of all the unique non-missing alleles of a locus.
"""
@inline function unique_alleles(locus::GenoArray)
    unique(alleles(locus))
end

"""
    convert_coord(coordinate::String)
Takes non-decimal-degree format as a `String` and returns it as a decimal degree
`Float32`. Can be broadcasted over an array of coordinate strings to convert them.
## Formatting requirements
- Coordinates as a `String` separated by spaces (`"11 43 41"`) or colons (`"11:43:41"`)
- Must use negative sign (`"-11 43.52"`) or single-letter cardinal direction (`"11 43.52W"`)
- Missing data should be coded as the string `"missing"` (can be accomplished with `replace!()`)
- Can mix colons and spaces (although it's bad practice)
### Example
```
julia> convert_coord("-41 31.52")
-41.5253f0

julia> convert_coord.(["-41 31.52", "25 11:54S"])
2-element Array{Float32,1}:
-41.5253
-25.1983
```
"""
function convert_coord(coordinate::String)
    lowercase(coordinate) == "missing" && return missing
    coord_strip = replace(uppercase(coordinate), r"[NSEW]" => "")
    split_coord = parse.(Float32, split(coord_strip, r"\s|:"))
    split_coord[2] /=60.0
    if length(split_coord) == 3
        split_coord[3] /=3600.0
    end
    conv = mapreduce(abs, +, split_coord)
    # N + E are positive | S + W are negative
    if split_coord[1] < 0 || occursin(r"[SW]", uppercase(coordinate))
        # negative
        return round(conv * -1, digits = 4)
    else
        # positive
        return round(conv, digits = 4)
    end
end

function Base.copy(data::PopData)
    PopData(copy.([data.meta,data.loci])...)
end

"""
    count_nonzeros(x::AbstractVector{T}) where T<:Real
Return the number of non-zero values in a vector
"""
function count_nonzeros(x::AbstractVector{T}) where T<:Real
    mapreduce(!iszero, +, x)
end

"""
    drop_monomorphic(data::PopData)
Return a `PopData` object omitting any monomorphic loci. Will inform you which
loci were removed.
"""
function drop_monomorphic(data::PopData)
    rm_loci = Vector{String}()
    for (loc, loc_sdf) in pairs(groupby(data.loci, :locus))
        length(unique(skipmissing(loc_sdf[:, :genotype]))) == 1 && push!(rm_loci, loc.locus)
    end

    if length(rm_loci) == 0
        return data
    elseif length(rm_loci) == 1
        @info "Dropping monomorphic locus " * rm_loci[1]
    else
        @info "Dropping $(length(rm_loci)) monomorphic loci" * "\n $rm_loci"
    end
    exclude(data, locus = rm_loci)
end

"""
    drop_monomorphic!(data::PopData)
Edit a `PopData` object in place by omitting any monomorphic loci. Will inform you which
loci were removed.
"""
function drop_monomorphic!(data::PopData)
    rm_loci = Vector{String}()
    for (loc, loc_sdf) in pairs(groupby(data.loci, :locus))
        length(unique(skipmissing(loc_sdf[:, :genotype]))) == 1 && push!(rm_loci, loc.locus)
    end
    if length(rm_loci) == 0
        return data
    elseif length(rm_loci) == 1
        @info "Dropping monomorphic locus " * rm_loci[1]
    else
        @info "Dropping $(length(rm_loci)) monomorphic loci" * "\n $rm_loci"
    end
    exclude!(data, locus = rm_loci)
end

"""
    loci_dataframe(data::PopData)
Return a wide `DataFrame` of samples as columns, ommitting population information.

**Example**
```
julia> loci_dataframe(@nancycats)
9×237 DataFrame. Omitted printing of 232 columns
│ Row │ N215       │ N216       │ N217       │ N218       │ N219       │
│     │ Tuple…?    │ Tuple…?    │ Tuple…?    │ Tuple…?    │ Tuple…?    │
├─────┼────────────┼────────────┼────────────┼────────────┼────────────┤
│ 1   │ missing    │ missing    │ (135, 143) │ (133, 135) │ (133, 135) │
│ 2   │ (136, 146) │ (146, 146) │ (136, 146) │ (138, 138) │ (140, 146) │
│ 3   │ (139, 139) │ (139, 145) │ (141, 141) │ (139, 141) │ (141, 145) │
│ 4   │ (116, 120) │ (120, 126) │ (116, 116) │ (116, 126) │ (126, 126) │
│ 5   │ (156, 156) │ (156, 156) │ (152, 156) │ (150, 150) │ (152, 152) │
│ 6   │ (142, 148) │ (142, 148) │ (142, 142) │ (142, 148) │ (142, 148) │
│ 7   │ (199, 199) │ (185, 199) │ (197, 197) │ (199, 199) │ (193, 199) │
│ 8   │ (113, 113) │ (113, 113) │ (113, 113) │ (91, 105)  │ (113, 113) │
│ 9   │ (208, 208) │ (208, 208) │ (210, 210) │ (208, 208) │ (208, 208) │
```
"""
function loci_dataframe(data::PopData)
    unstack(select(data.loci, Not(:population)), :name, :genotype)[:, Not(:locus)]
end

#TODO make a SMatrix instead?
"""
    loci_matrix(data::PopData)
Return a matrix of genotypes with dimensions `samples × loci`.
Rows are samples and columns are loci. Will return an error if ploidy varies between samples.

**Example**
```
julia> loci_matrix(@nancycats)
237×9 Array{Union{Missing, Tuple{Int16,Int16}},2}:
 missing     (136, 146)  (139, 139)  …  (199, 199)  (113, 113)  (208, 208)
 missing     (146, 146)  (139, 145)     (185, 199)  (113, 113)  (208, 208)
 (135, 143)  (136, 146)  (141, 141)     (197, 197)  (113, 113)  (210, 210)
 (133, 135)  (138, 138)  (139, 141)     (199, 199)  (91, 105)   (208, 208)
 (133, 135)  (140, 146)  (141, 145)     (193, 199)  (113, 113)  (208, 208)
 (135, 143)  (136, 146)  (145, 149)  …  (193, 195)  (91, 113)   (208, 208)
 (135, 135)  (136, 146)  (139, 145)     (199, 199)  (105, 113)  (208, 208)
 (135, 143)  (136, 146)  (135, 149)     (193, 197)  (91, 91)    (208, 212)
 (137, 143)  (136, 146)  (139, 139)     (197, 197)  (105, 113)  (208, 212)
 (135, 135)  (132, 132)  (141, 145)     (197, 197)  (91, 105)   (208, 208)
 (137, 141)  (130, 136)  (137, 145)  …  (193, 199)  (91, 91)    (182, 182)
 (129, 133)  (130, 136)  (135, 145)     (193, 199)  (91, 113)   (182, 208)
 ⋮                                   ⋱
 (133, 135)  (136, 136)  (135, 139)  …  (199, 199)  (113, 113)  (182, 182)
 (133, 141)  (136, 136)  (135, 139)     (197, 197)  (113, 113)  (182, 208)
 (133, 141)  (130, 146)  (141, 141)     (191, 199)  missing     (208, 208)
 (123, 133)  (138, 138)  (141, 145)     (191, 197)  missing     (208, 208)
 (123, 133)  (138, 138)  (139, 139)     (197, 199)  missing     (208, 208)
 (133, 141)  (136, 146)  (139, 139)  …  (197, 197)  missing     (208, 208)
 (133, 141)  (130, 136)  (139, 145)     (191, 199)  missing     (208, 208)
 (133, 141)  (136, 146)  (139, 145)     (199, 199)  missing     (208, 220)
 (133, 143)  (130, 130)  (135, 145)     (197, 197)  missing     (208, 208)
 (135, 141)  (136, 144)  (143, 143)     (191, 197)  (113, 117)  (208, 208)
 (137, 143)  (130, 136)  (135, 145)  …  (193, 199)  (113, 117)  (208, 208)
 (135, 141)  (130, 146)  (135, 139)     (197, 197)  missing     (208, 208)
 ```
"""
function loci_matrix(data::PopData)
    dims = size(data)
    sort_df = issorted(data.loci, [:name, :locus]) ? sort(data.loci, [:name, :locus]) : data.loci
    reshape(sort_df.genotype, (dims.samples, dims.loci)) |> collect
end

"""
    multitest_missing(pvals::Vector{T}, method::String) where T <: Union{Missing, <:AbstractFloat}
Modification to `MultipleTesting.adjust` to include `missing` values in the
returned array. See MultipleTesting.jl docs for full more detailed information.

**Example**
```
julia> multitest_missing([0.1, 0.01, 0.005, 0.3], "bh")

```

### `correction` methods (case insensitive)
- `"bonferroni"` : Bonferroni adjustment
- `"holm"` : Holm adjustment
- `"hochberg"` : Hochberg adjustment
- `"bh"` : Benjamini-Hochberg adjustment
- `"by"` : Benjamini-Yekutieli adjustment
- `"bl"` : Benjamini-Liu adjustment
- `"hommel"` : Hommel adjustment
- `"sidak"` : Šidák adjustment
- `"forwardstop"` or `"fs"` : Forward-Stop adjustment
- `"bc"` : Barber-Candès adjustment
"""
@inline function multitest_missing(pvals::Vector{T}, method::String) where T <: Union{Missing, <:AbstractFloat}
    # make a dict of all possible tests and their respective functions
    d = Dict(
        "bonferroni" => Bonferroni(),
        "holm" => Holm(),
        "hochberg" => Hochberg(),
        "bh" => BenjaminiHochberg(),
        "by" => BenjaminiYekutieli(),
        "bl" => BenjaminiLiu(),
        "hommel" => Hommel(),
        "sidak" => Sidak(),
        "forwardstop" => ForwardStop(),
        "fs" => ForwardStop(),
        "bc" => BarberCandes(),
    )
    p_copy = copy(pvals)
    p_copy[.!ismissing.(p_copy)] .= adjust(p_copy[.!ismissing.(p_copy)] |> Vector{Float64}, d[lowercase(method)])
    return p_copy
end

"""
    nonmissing(vec::T) where T<:AbstractArray
Convenience function to count the number of non-`missing` values
in a vector.
"""
@inline function nonmissing(vec::T) where T<:AbstractArray
    mapreduce(!ismissing, +, vec)
end

"""
    nonmissing(data::PopData, locus::String)
Convenience function to count the number of non-`missing` samples
at a locus.
"""
@inline function nonmissing(data::PopData, locus::String)
    data.loci[data.loci.locus .== locus, :genotype] |> nonmissing
end

"""
    nonmissings(vec1::AbstractVector, vec2::AbstractVector)
Return a vector of indices where neither input vectors have a `missing` value, i.e. an
intersection of the indices of their non-missing elements.
"""
@inline function nonmissings(vec1::T, vec2::T) where T <: AbstractVector
    intersect(map(i -> findall(!ismissing, i), (vec1, vec2))...)
end

"""
    pairwise_pairs(smp_names::Vector{T}) where T
Given a vector of some iterable, returns a vector of tuples of unique all x
all combinations of pairs, excluding self-comparisons.

**Example**
```
julia> samps = ["red_1", "red_2", "blue_1", "blue_2"] ;

julia> pairwise_pairs(samps)
6-element Array{Tuple{String,String},1}:
 ("red_1", "red_2")
 ("red_1", "blue_1")
 ("red_1", "blue_2")
 ("red_2", "blue_1")
 ("red_2", "blue_2")
 ("blue_1", "blue_2")
```
"""
@inline function pairwise_pairs(smp_names::AbstractVector{T}) where T
    [tuple(smp_names[i], smp_names[j]) for i in 1:length(smp_names)-1 for j in i+1:length(smp_names)]
end

"""
    partitionarray(array::AbstractArray, steps::AbstractVector{<:Integer})
Like Base.Iterators.Partition, except you can apply arbitrary sizes to
partition the array by. The `steps` must add up to the total row length
of the array.

**Example**
```
julia> partitionmatrix(rand(20,5), [10,3,4,3]) .|> size
((10, 5), (3, 5), (4, 5), (3, 5))
```
"""
# solution brilliantly provided by @stevengj and @mcabbott on Slack and Discourse (https://discourse.julialang.org/t/is-there-a-simple-intuitive-way-to-partition-a-matrix-by-arbitrary-strides-like-i/55863)
function partitionarray(array::AbstractArray, steps::AbstractVector{<:Integer})
    v = axes(array,1)
    v == 1:sum(steps) || error("Steps provided do not sum to length of the first dimension")
    i = firstindex(v)
    tmp = (view(v, i:(i+=s)-1) for s in steps)
    [view(array,r,:) for r in tmp]
end

"""
    phase(data::PopData)
Return a `Vector` of length `ploidy` composed of allele matrices with dimensions `samples × loci`.
Rows are samples and columns are loci. Will return an error if ploidy varies between samples.

**Example**
```
julia> mtx = phase(@nancycats)
2-element Array{Array{Union{Missing, Int16},2},1}:
 [missing 136 … 113 208; missing 146 … 113 208; … ; 137 130 … 113 208; 135 130 … missing 208]
 [missing 146 … 113 208; missing 146 … 113 208; … ; 143 136 … 117 208; 141 146 … missing 208]

julia> mtx[1]
237×9 Array{Union{Missing, Int16},2}:
    missing  136  139  116         156  142  199  113         208
    missing  146  139  120         156  142  185  113         208
 135         136  141  116         152  142  197  113         210
 133         138  139  116         150  142  199   91         208
 133         140  141  126         152  142  193  113         208
 135         136  145  120         150  148  193   91         208
 135         136  139  116         152  142  199  105         208
 135         136  135  120         154  142  193   91         208
 137         136  139  116         150  142  197  105         208
 135         132  141  120         150  148  197   91         208
 137         130  137  128         152  142  193   91         182
 129         130  135  126         144  140  193   91         182
   ⋮                                      ⋮
 133         136  135     missing  146  142  199  113         182
 133         136  135     missing  150  142  197  113         182
 133         130  141     missing  148  142  191     missing  208
 123         138  141     missing  148  142  191     missing  208
 123         138  139     missing  150  142  197     missing  208
 133         136  139     missing  150  142  197     missing  208
 133         130  139     missing  152  142  191     missing  208
 133         136  139     missing  150  142  199     missing  208
 133         130  135     missing  148  142  197     missing  208
 135         136  143     missing  144  142  191  113         208
 137         130  135     missing  150  142  193  113         208
 135         130  135     missing  150  142  197     missing  208
```
"""
function phase(data::PopData)
    dims = size(data)
    ploidy = unique(data.meta.ploidy)
    ploidy = length(ploidy) != 1 ? error("Phasing will not work on mixed-ploidy samples") : ploidy[1]
    sort_df = issorted(data.loci, [:name, :locus]) ? sort(data.loci, [:name, :locus]) : data.loci
    matrices = map(j -> map(i -> ismissing(i) ? missing : i[j] , sort_df.genotype), 1:ploidy)
    map(i -> collect(reshape(i, (dims.samples, dims.loci))), matrices)
end

"""
    reciprocal(num::T) where T <: Signed
Returns the reciprocal (1/number) of a number. Will return `0` when
the number is `0` instead of returning `Inf`.
"""
function reciprocal(num::T) where T <: Real
    !iszero(num) ? 1.0/float(num) : 0.0
end

"""
    reciprocal_sum(x::AbstractVector{T}) where T<:Real
Return the sum of the reciprocal values of `x`, skipping the `Inf` values
resulting from divide-by-zero errors.
"""
function reciprocal_sum(x::AbstractVector{T}) where T<:Real
    mapreduce(reciprocal, +, x)
end

"""
    safemean(::AbstractVector{T}) where T<:Real
A wrapper for StatsBase.mean to calculate a mean after skipping `Inf`, `-Inf`, and `NaN` values.
"""
function safemean(x::AbstractVector{T}) where T<:Real
    y = x
    if any(isinf.(x))
        y =  x[.!isinf.(x)]
    end
    if any(isnan.(x))
        y = y[.!isnan.(y)]
    end
    mean(y)
end

"""
    sim_pairs(data::Vector{String})
Takes a Vector of sample names and returns a Tuple of sample pairs, grouped by simulation
number. This is an internal function used for isolating sibship pairs from simulated shipship
pairs (via `PopGenSims.jl`) to perform `relatedness` estimates only on those pairs.

**Example**
```julia
julia> a = ["sim1_off1", "sim1_off2", "sim2_off1", "sim2_off2"] ;

julia> sim_pairs(a)
("sim1_off1", "sim1_off2")
("sim2_off1", "sim2_off2")
```
"""
function sim_pairs(data::Vector{String})
    n = length(data)
    isodd(n) && error("Expected an even number of samples, but got $n")
    Tuple.(Base.Iterators.partition(sort(data), 2))
end

"""
    Base.sort(x::NTuple{N,T}) where N where T <: Signed
Sort the integers within a Tuple and return the sorted Tuple.
"""
function Base.sort(x::NTuple{N,T}) where N where T <: Signed
    Tuple(sort(SVector(x)))
end

"""
    strict_shuffle(x::T) where T <: AbstractArray
Shuffle only the non-missing values of a Vector and return a copy of the vector,
keeping the `missing` values at their original locations.
Use `strict_shuffle!` to edit in-place instead of returning a copy.
"""
@inline function strict_shuffle(x::T) where T <: AbstractArray
    # get indices of where original missing are
    miss_idx = findall(i -> i === missing, x)
    out_vec = shuffle(Xoroshiro128Star(), x[.!ismissing.(x)])

    insert!.(Ref(out_vec), miss_idx, missing)
    return out_vec
end

"""
    strict_shuffle!(x::T)! where T <: AbstractArray
Shuffle only the non-missing values of a Vector, keeping the
`missing` values at their original locations. Use `strict_shuffle`
to return a copy instead of editing in-place.
"""
function strict_shuffle!(x::T) where T <: AbstractArray
    @inbounds shuffle!(Xoroshiro128Star(), @view x[.!ismissing.(x)])
    return x
end

"""
    generate_meta(data::DataFrame)
Given a genotype DataFrame formatted like `PopData.loci`, generates a corresponding
`meta` DataFrame. In other words, it creates the `.meta` part of `PopData` from the `.loci` part.

**Example**
```
julia> cats = @nancycats ;

julia> cats_nometa = cats.loci ;

julia> cats_meta = generate_meta(cats_nometa)
237×5 DataFrame
 Row │ name    population  ploidy  longitude  latitude
     │ String  String      Int8    Float32?   Float32?
─────┼─────────────────────────────────────────────────
   1 │ N215    1                2   missing   missing
   2 │ N216    1                2   missing   missing
   3 │ N217    1                2   missing   missing
   4 │ N218    1                2   missing   missing
   5 │ N219    1                2   missing   missing
   6 │ N220    1                2   missing   missing
   7 │ N221    1                2   missing   missing
  ⋮  │   ⋮         ⋮         ⋮         ⋮         ⋮
 232 │ N295    17               2   missing   missing
 233 │ N296    17               2   missing   missing
 234 │ N297    17               2   missing   missing
 235 │ N281    17               2   missing   missing
 236 │ N289    17               2   missing   missing
 237 │ N290    17               2   missing   missing
                                       224 rows omitted
```
"""
function generate_meta(data::DataFrame)
    grp = groupby(data, :name)
    nms = map(z -> z.name, keys(grp))
    pops = map(z -> first(z.population), grp)
    ploids = map(z -> find_ploidy(z.genotype), grp)
    DataFrame(
        :name => nms,
        :population => pops,
        :ploidy => ploids,
        :longitude => Vector{Union{Missing, Float32}}(undef, (length(nms))),
        :latitude => Vector{Union{Missing, Float32}}(undef, (length(nms)))
    )
end


text:
# Abstract type: whether the expert is zero-inflated
abstract type ZeroInflation end
struct ZI <: ZeroInflation end
struct NonZI <: ZeroInflation end

abstract type ExpertSupport end
struct RealValued <: ExpertSupport end
struct NonNegative <: ExpertSupport end

# Abstract type: AnyExpert
abstract type AnyExpert{s<:ExpertSupport, z<:ZeroInflation, d<:UnivariateDistribution} end
# Discrete or Continuous
const DiscreteExpert{s<:ExpertSupport, z<:ZeroInflation} = AnyExpert{s, z, DiscreteUnivariateDistribution}
const ContinuousExpert{s<:ExpertSupport, z<:ZeroInflation} = AnyExpert{s, z, ContinuousUnivariateDistribution}
# Real-valued expert distributions
const RealDiscreteExpert = DiscreteExpert{RealValued, NonZI}
const RealContinuousExpert = ContinuousExpert{RealValued, NonZI}
# Nonnegative-valued expert distributions: actuarial-specific
const NonNegDiscreteExpert{z<:ZeroInflation} = DiscreteExpert{NonNegative, z}
const NonNegContinuousExpert{z<:ZeroInflation} = ContinuousExpert{NonNegative, z}
# Zero-inflated
const ZIDiscreteExpert = NonNegDiscreteExpert{ZI}
const ZIContinuousExpert = NonNegContinuousExpert{ZI}
# Non zero-inflated
const NonZIDiscreteExpert = NonNegDiscreteExpert{NonZI}
const NonZIContinuousExpert = NonNegContinuousExpert{NonZI}

## Broadcast issues
Broadcast.broadcastable(d::RealDiscreteExpert) = Ref(d)
Broadcast.broadcastable(d::RealContinuousExpert) = Ref(d)
Broadcast.broadcastable(d::ZIDiscreteExpert) = Ref(d)
Broadcast.broadcastable(d::ZIContinuousExpert) = Ref(d)
Broadcast.broadcastable(d::NonZIDiscreteExpert) = Ref(d)
Broadcast.broadcastable(d::NonZIContinuousExpert) = Ref(d)

##### specific distributions #####

const add_experts = [
    "burr",
    "gammacount"
]

const discrete_experts = [
    "binomial", "zibinomial",
    "gammacount", "zigammacount",
    "negativebinomial", "zinegativebinomial",
    "poisson", "zipoisson"
]

const continuous_experts = [
    "burr", "ziburr",
    "gamma", "zigamma",
    "inversegaussian", "ziinversegaussian",
    "lognormal", "zilognormal",
    "weibull", "ziweibull"

]

const calcs = [
    # "expert_ll_pos",
    # "expert_tn_pos",
    # "expert_tn_bar_pos",
    # "expert_ll",
    # "expert_tn",
    # "expert_tn_bar",
    "exposurize"
]

for dname in add_experts
    include(joinpath("experts", "add_dist", "$(dname).jl"))
end

for dname in discrete_experts
    include(joinpath("experts", "discrete", "$(dname).jl"))
end

for dname in continuous_experts
    include(joinpath("experts", "continuous", "$(dname).jl"))
end

for dname in calcs
    include(joinpath("experts", "ll", "$(dname).jl"))
end

max_stars_repo_name: sparktseung/JuliaLRMoE

max_issues_repo_name: sparktseung/JuliaLRMoE

max_forks_repo_name: sparktseung/JuliaLRMoE


text:
# Code derived from metabolicEP (https://github.com/anna-pa-m/Metabolic-EP)

function scaleepfield!(scalefact, epfields, vals...)
    @extract epfields : μ s av va
    rmul!(μ, scalefact)
    rmul!(s, scalefact^2)
    rmul!(av, scalefact)
    rmul!(va, scalefact^2)

    for val in vals
        rmul!(val, scalefact)
    end

end


text:
# vanderPol winding number arnold tongues
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("style.jl"))
using DynamicalSystems, PyPlot

function winding_number(ds, ω, a, u0, Ttr, T)
    period = 2π/ω
    ds.p[2:3] .= (a, period)
    # T = period*ceil(T/period)   # no. of periods required for transient
    # Ttr = period*ceil(Ttr/period)     # no. of periods used for averaging
    tr = trajectory(ds, T, u0; Δt = period/20, Ttr)
    tvec = 0:period/20:T
    u1, u2 = columns(tr)
    θ = 0.0
    θ_old = atan(u2[1], u1[1])
    for (x, y) ∈ zip(u1, u2)
        θ_new = atan(y, x)
        θ += mod(θ_new - θ_old + π, 2π) - π
        θ_old = θ_new
    end
    W = abs(θ / (tvec[end] * 2π/period))
    return tr[end], W
end

using ProgressMeter

d = 5.0

a_min = 0.0
a_max = 1.0
na = 100
a_vec = LinRange(a_min, a_max, na)

ω_min = 0.2
ω_max = 2.0
nω =  1000
ω_vec = LinRange(ω_min, ω_max, nω)

ds = Systems.vanderpol(; μ = d)

wmat = zeros(nω,na)

Ttr =  5000 # transient time
T =   4000 #  time on the attractor used for estimating winding no.
u0 = ds.u0

W_vec = zeros(nω)
@showprogress 1 "Computing..." for ia = 1:na
    a = a_vec[ia]
    for iω = 1:length(ω_vec)
        ω = ω_vec[iω]
        u0, W = winding_number(ds, ω, a, u0, Ttr, T)
        wmat[iω,ia] = 1/W
    end # of frequency loop
end  # of amplitude loop

# %% Plot with all numbers
fig = figure(figsize=(0.6*figx,figy))
im1 = plt.pcolormesh(ω_vec,a_vec,wmat',cmap = "gnuplot" )

xlabel(L"\omega")
ylabel(L"a")

ax = gca()
ax.set_xlim([ω_min,ω_max])

cbar = plt.colorbar(im1)  # , ticks=[-0.2, 0.2, 0.6, 1.], orientation="horizontal")
cbar.set_label(L"W")

fig.tight_layout(pad=0.3)

# %% Plot version with only wanted numbers
figure()
qmat = 4*ones(nω,na)   # array with winding numbers to be displayed

rat_num = [1/2, 2/3,   1, 4/3, 3/2, 5/3, 2,  7/3, 5/2, 8/3,  3,  10/3,  7/2]
nrn = length(rat_num)
for ia = 1:na
    for iomega = 1:nω
        wind_no = wmat[iomega,ia]
        # check whether equal p/q
        for rat in rat_num
            if abs(rat - wind_no) < 0.008  # this threshold should be as small as possible
               qmat[iomega,ia] = rat
            end
        end
    end
end

im1 = plt.pcolormesh(ω_vec,a_vec,qmat',cmap = "CMRmap") # a color map which ends with white
plt.clim(0.,4.)
ylabel(L"a")
xlabel(L"\omega"; labelpad = -10)
cbar = plt.colorbar(im1)
cbar.set_label(L"W")
fig.tight_layout(pad=0.3)
wsave(plotsdir("9", "vanderPol_tongues"), fig)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
module MetaDataGraphs

using Graphs
using SparseArrays

export AbstractDataGraph, DataDiGraph
export get_vertex, get_label
export get_data, set_data!

include("abstract/abstractdatagraph.jl")
include("abstract/data.jl")
include("abstract/graphs_interface.jl")
include("abstract/weights.jl")
include("abstract/dict.jl")

include("datadigraph/datadigraph.jl")
include("datadigraph/data.jl")
include("datadigraph/graphs_interface.jl")

end


text:
module CauseMap
using Base.LinAlg.BLAS
using PyCall

include("plotting.jl")
include("coord_descent_tuning.jl")
include("calc_manifolds.jl")

export makesingleplot, makeoptimizationplots, optandcalcCCM, precalc_manif_dists, calcCCM, load_example_data

function load_example_data(dsname)
    if dsname == "ParaDidi"
        ds   = readdlm(Pkg.dir("CauseMap") * "/examples/vr.raw_fixed.txt", '\t', Float64, header=false)  # read in data
        para = ds[:, 2]
        didi = ds[:, 3]
        return para, didi
    else
        throw(ArgumentError("Did not recognize $dsname dataset"))
    end
end

function processpredvals_simple(predvals::Array{Float64,2}, targvals::Vector{Float64}, lib_size::Int)
    nanval = NaN

    if all(isnan(predvals))
        println("All predvals were nanval in process predvals of lib size of $lib_size")
        return nanval
    end
    rhos = fill(nanval, size(predvals,2))

    for xx in 1:size(predvals,2)
        touse = !isnan(predvals[:,xx]) & !isnan(targvals)

        if sum(touse) < 10
            rhos[xx] = nanval
        else
            rhos[xx] = cor(predvals[touse,xx], targvals[touse])
            if isnan(rhos[xx])
                warn("in process predvals, why is rho nan?")
            end
        end
    end
    toret = !isnan(rhos)
    if sum(toret) > 0
        return median(rhos[toret])
    else
        println("no non-nan values")
        return nanval
    end

end

function weightfunc(distances::Array{Float64}; kernelargs...)
    w = [exp(-d) for d in distances]
    return w / sum(w)
end

function getdist!(dist_top::Vector{Float64}, inds_touse::AbstractVector{Int}, num_neighbors::Int)
    if dist_top[1] == 0
        for xx in 1:num_neighbors
            if xx == 1
                dist_top[xx] = 1
            else
                dist_top[xx] = 1e200
            end
        end
        return 1e-100
    else
        min_dist::Float64 = dist_top[1]
        for xx in num_neighbors:-1:1
            dist_top[xx] = dist_top[xx] / min_dist
        end
        return min_dist
    end
end

function getpredstartstop(nobs::Int, ll::Int, lib_size::Int, npred::Int, pred_start_min::Int)
    left::Int = iceil(npred / 2)
    right::Int = npred - left - 1
    midpoint::Int = iceil(ll + lib_size / 2)
    lp::Int = midpoint - left
    rp::Int = midpoint + right

    if lp < pred_start_min
        rp += (pred_start_min - lp)
        lp  = pred_start_min
    end

    if rp > nobs
        lp -= rp % nobs
        rp  = nobs
    end

    if (rp - lp + 1) != npred
        error("start and stop not equal to npred")
    end

    return lp, rp
end

### this function accounts for ~70% of algorithm run time
function calcdistslice!(source_dists::Array{Float64, 2}, dist_top::Vector{Float64}, slice_inds::AbstractVector{Int}, topred::Int, nn::Range1{Int})
    if in(topred, slice_inds)
        slice_inds = convert(Vector{Int}, slice_inds)
        splice!(slice_inds, findfirst([yy == topred for yy in slice_inds])) # remove topred, but don't move any of the preceding elements in the array
    end

    # this line accounts for ~60% of algorithm run time
    inds_touse::Vector{Int} = slice_inds[sortperm(source_dists[slice_inds, topred])[nn]]
    for xx in nn
        dist_top[xx] = source_dists[inds_touse[xx], topred]  # this is not returned because it is written in-place
    end
    return inds_touse
end

function prepgetpred(source_manifold, lib_start, lib_size, tau_p, npred, num_neighbors; nboots=0)
    nobs::Int = size(source_manifold, 2)
    lib_stop::Int = nobs - lib_size - tau_p
    nlib::Int = lib_stop - lib_start + 1
    if nlib < 1
        println("No libs of this size in dataset. Lib_start is $lib_start, lib_stop is $lib_stop")
        return fill(NaN, 1, 1), fill(NaN, 1,1)
    end

        if nboots > 0
        min_distances = fill(NaN, npred, nboots)
        predvals = fill(NaN, npred, nboots)
        targvals = fill(NaN, npred)
    else
        min_distances = fill(NaN, npred, nlib)
        predvals = fill(NaN, npred, nlib)
        targvals = fill(NaN, npred)
    end
    lib_end::Int = lib_start + lib_size - 1
    dist_top = fill(NaN, num_neighbors)
    return nobs, lib_stop, nlib, predvals, targvals, min_distances, lib_end, dist_top
end

function getpredvals_boot(source_manifold::Array{Float64,2}, source_dists::Array{Float64,2},
                                             target_series::Vector{Float64}, lib_size::Int,
                                             lib_start::Int, num_neighbors::Int,
                                             tau_p::Int, npred::Int, pred_start_min::Int, nboots::Int)
    # Use bootstrap samples to generate libraries

    nobs, lib_stop, nlib, predvals, targvals, min_distances, lib_end, dist_top = prepgetpred(source_manifold,
                                                                                lib_start, lib_size,
                                                                                tau_p, npred, num_neighbors;
                                                                                nboots=nboots)
    # lib_stop is the last possible sliding window START index (stop of the start...confusing...)
    # lib_end is the last in the current window
    nn = 1:num_neighbors
    for (pp, topred) in enumerate(
                                                    pred_start_min:min(pred_start_min + npred, length(target_series))
                                                    )
        targvals[pp] = target_series[topred]
        for xx in 1: nboots
            slice_inds::Array{Int, 1} = [rand(lib_start:lib_stop) for xx in 1:lib_size]

            inds_touse = calcdistslice!(source_dists, dist_top, slice_inds, topred, nn)
            min_distances[pp, xx] = getdist!(dist_top, inds_touse, num_neighbors)
            weights = weightfunc(dist_top)
            predvals[pp, xx] = dot(weights, target_series[inds_touse + tau_p])
        end
    end

    rhos = processpredvals_simple(predvals, targvals, lib_size)
    return predvals, min_distances, rhos
end

function getpredvals_sw(source_manifold::Array{Float64,2}, source_dists::Array{Float64,2},
    target_series::Vector{Float64}, lib_size::Int,
    lib_start::Int, num_neighbors::Int,
    tau_p::Int, npred::Int, pred_start_min::Int)
    # Use a sliding window for your library. This is the traditional form of CCM

    nobs, lib_stop, nlib, predvals, targvals, min_distances, lib_end, dist_top = prepgetpred(source_manifold,
                                                                                lib_start, lib_size,
                                                                                tau_p, npred, num_neighbors)

    # lib_stop is the last possible sliding window START index (stop of the start..confusing...)
    # lib_end is the last in the current window
    nn = 1:num_neighbors

    start_count::Int = 1
    for ll in lib_start:lib_stop
        ll_end::Int     = ll + lib_size - 1
        slice_inds = ll:ll_end
        predstart, predstop = getpredstartstop(nobs, ll, lib_size, npred, pred_start_min)

        pred_count::Int=1
        for topred in predstart:predstop
            slice_inds = ll:ll_end
            inds_touse = calcdistslice!(source_dists, dist_top, slice_inds, topred, nn)

            min_distances[pred_count, start_count] = getdist!(dist_top, inds_touse, num_neighbors)

            weights = weightfunc(dist_top)
            predvals[pred_count, start_count] = dot(weights, target_series[inds_touse+tau_p])
            targvals[pred_count] = target_series[topred]

            pred_count += 1
        end
        start_count += 1
        lib_end += 1
    end
    rhos = processpredvals_simple(predvals, targvals, lib_size)
    return predvals, min_distances, rhos
end

function cross_mapping(source_manif_dict::Dict, source_dist_dict::Dict,
                                        target_series::AbstractVector,
                                        nobs::Int,
                                        libsizemin::Int, libsizemax::Int, E::Int,
                                        tau_s::Int, tau_p::Int, npred::Int,
                                        pred_start_min::Int, num_neighbors, args...;
                                        lib_start::Int=0, nboots=0)

    res12 = fill(NaN, libsizemax - libsizemin + 1)
    count = 1

    nlibpluslibsize = nobs - tau_p - lib_start + 1

    for lib_size in libsizemin:libsizemax
        nlib = nlibpluslibsize - lib_size
        if nboots > 0
            predvals, min_distances, rhos = getpredvals_boot(source_manif_dict[tau_s][E], source_dist_dict[tau_s][E],
                                target_series, lib_size, lib_start, num_neighbors, tau_p, npred, pred_start_min, nboots)
        else
            predvals, min_distances, rhos = getpredvals_sw(source_manif_dict[tau_s][E], source_dist_dict[tau_s][E],
                                target_series, lib_size, lib_start, num_neighbors, tau_p, npred, pred_start_min)
        end
        res12[count] = rhos

        count += 1
    end
    if count == 1
        warn("Why was there no loop?")
    end
    return res12
end

function calclibstart(shadowmat_dict::Dict, E::Int, tau_s::Int)

    nanrows::Array{Bool, 1} = [in(NaN, shadowmat_dict[tau_s][E][:,xx]) for xx in 1:size(shadowmat_dict[tau_s][E], 2)]
    if in(true, nanrows)
        return (maximum(find(nanrows)) + 1)
    else
        return 1
    end
end

function calclibsizemax(nobs::Int, E::Int, tau_s::Int, tau_p::Int)  ## REMOVE
    return ifloor(nobs/tau_s) - E - tau_p
end

function check_lib_pred_defaults(libsizemin, libsizemax, lib_start,
                                                        npred, pred_start,
                                                        nobs, E, b_offset, nlag,
                                                        tau_p, tau_s, shadowmat_dict;
                                                        quick=false, lmin_def=8, lmax_offset=2)
    if libsizemin == 0
        libsizemin = lmin_def
    end

    if libsizemax == 0
        libsizemax = calclibsizemax(nobs, E, tau_s, tau_p)
    end
    libsizemax = min(libsizemax, nobs - tau_p - lib_start)

    if quick
        libsizemin = max(E + b_offset + 1, libsizemax - nlag)
    else
        libsizemin = max(libsizemin, E + b_offset + 1)
    end

    if lib_start == 0
        lib_start = calclibstart(shadowmat_dict, E, tau_s)
    end

    if npred == 0
        npred = nobs
    end

    if pred_start == 0
        pred_start = 1
    end
    pred_start = max(lib_start, pred_start)

    if pred_start + npred - 1 > nobs
        npred -= (pred_start + npred - 1) - nobs
    end

    return libsizemin, libsizemax, lib_start, npred, pred_start
end

function calcCCM(var1::AbstractVector, var2::AbstractVector,
                                shadowmat_dict::Dict, distmat_dict::Dict,
                                E::Int, tau_s::Int, tau_p::Int;
                                libsizemin::Int=0, libsizemax::Int=0,
                                npred::Int=0, pred_start::Int=0,
                                lib_start::Int=0, b_offset::Int=1,
                                quick=false, nlag=10, nboots=0)

    nobs = length(var1)
    num_neighbors = E + b_offset

    libsizemin, libsizemax, lib_start, npred, pred_start = check_lib_pred_defaults(libsizemin, libsizemax, lib_start,
                                                                                npred, pred_start,
                                                                                nobs, E, b_offset, nlag,
                                                                                tau_p, tau_s, shadowmat_dict; quick=quick)

    ################## clean input parameters

    #########
    ##### Raise error or warning messages as needed.
    if libsizemin > libsizemax
        println("Libsizemin, libsizemax: ($libsizemin, $libsizemax)")
        println("E, tau_s, tau_p, nobs, lib_start : ($E, $tau_s, $tau_p, $nobs, $lib_start)")
        warn("why is libsizemin less than libsizemax? returning NaN")
        return 0:0, fill(NaN, 1)
    end

    ###################
    if nboots > 0
        warn("Bootstrapped libraries are still experimental")
    end
    res12 = cross_mapping(shadowmat_dict, distmat_dict, var2, nobs, libsizemin, libsizemax, E,
                            tau_s, tau_p, npred, pred_start, num_neighbors; lib_start = lib_start, nboots=nboots)

    return libsizemin:libsizemax, res12
end

function optandcalcCCM(vec1::AbstractVector, vec2::AbstractVector,
                                        E_vals::AbstractVector,  tau_s_vals::AbstractVector, tau_p_vals::AbstractVector;
                                        nreps=5, b_offset=1, nboots=0,
                                        npred::Int=0, pred_start::Int=0,
                                        libsizemin::Int=0, libsizemax::Int=0)
    """
    vec1: Time series 1
    vec2: Time series 2
    libsizemin: Minimum library size
    libsizemax: Maximum library size
    E_vals: A vector of dimensions to try
    tau_s_vals: A vector of lag lengths to use for manifold reconstruction
    tau_p_vals: A vector of lag lengths to try for the causal effect
    npred: Number of points to predict
    pred_start: First point to predict
    ## kwargs
    nreps: Number of coordinate descent runs
    """

    shadowmat_dict_vec1, distmat_dict_vec1 = precalc_manif_dists(E_vals, tau_s_vals, vec1)
    shadowmat_dict_vec2, distmat_dict_vec2 = precalc_manif_dists(E_vals, tau_s_vals, vec2)

    println("CoordDescentOpt opt1")
    res12 = CoordDescentOpt(vec1, vec2,
                                                shadowmat_dict_vec1, distmat_dict_vec1,
                                                E_vals, tau_s_vals, tau_p_vals;
                                                libsizemin=libsizemin, libsizemax=libsizemax,
                                                npred=npred, pred_start=pred_start,
                                                nreps=nreps, nboots=nboots)

    println("CoordDescentOpt opt2")
    res21 = CoordDescentOpt(vec2, vec1,
                                                shadowmat_dict_vec2, distmat_dict_vec2,
                                                E_vals, tau_s_vals, tau_p_vals;
                                                libsizemin=libsizemin, libsizemax=libsizemax,
                                                npred=npred, pred_start=pred_start,
                                                nreps=nreps, nboots=nboots)

    libsizemin_12 = max(res12["E"] + b_offset + 1, libsizemin)
    libsizemin_21 = max(res21["E"] + b_offset + 1, libsizemin)

    println("starting calcCCM1")
    librange12, yval_12 = calcCCM(vec1, vec2,
                                                        shadowmat_dict_vec1, distmat_dict_vec1,
                                                        res12["E"], res12["tau_s"], res12["tau_p"];
                                                        libsizemin=libsizemin_12, libsizemax=libsizemax,
                                                        npred=npred, pred_start=pred_start,
                                                        nboots=nboots)
    println("starting calcCCM2")
    librange21, yval_21 = calcCCM(vec2, vec1,
                                                        shadowmat_dict_vec2, distmat_dict_vec2,
                                                        res21["E"], res21["tau_s"], res21["tau_p"], ;
                                                        libsizemin=libsizemin_21, libsizemax=libsizemax,
                                                        npred=npred, pred_start=npred,
                                                        nboots=nboots)

    return (librange12, yval_12), (librange21, yval_21)
end
########## end CCM functions
end # module