text:
abstract type Simulation{D} end
abstract type UserParams end
abstract type TransformLibrary end
abstract type Field end

struct XField{D} <: Field
    psiX::Array{Complex{Float64},D}
    X::NTuple{D}
    K::NTuple{D}
    K2::Array{Float64,D}
end

struct KField{D} <: Field
    psiK::Array{Complex{Float64},D}
    X::NTuple{D}
    K::NTuple{D}
    K2::Array{Float64,D}
end

@with_kw mutable struct Transforms{D,N} <: TransformLibrary
    Txk::AbstractFFTs.ScaledPlan{Complex{Float64},FFTW.cFFTWPlan{Complex{Float64},-1,false,D,UnitRange{Int64}},Float64} = 0.1*plan_fft(crandn_array(D))
    Txk!::AbstractFFTs.ScaledPlan{Complex{Float64},FFTW.cFFTWPlan{Complex{Float64},-1,true,D,UnitRange{Int64}},Float64} = 0.1*plan_fft!(crandn_array(D))
    Tkx::AbstractFFTs.ScaledPlan{Complex{Float64},FFTW.cFFTWPlan{Complex{Float64},1,false,D,UnitRange{Int64}},Float64} = 0.1*plan_ifft(crandn_array(D))
    Tkx!::AbstractFFTs.ScaledPlan{Complex{Float64},FFTW.cFFTWPlan{Complex{Float64},1,true,D,UnitRange{Int64}},Float64} = 0.1*plan_ifft!(crandn_array(D))
    psi::ArrayPartition = crandnpartition(D,N)
end

# ==== define user parameters =======
# rename and make this whatevery you require and pass it to sim
@with_kw mutable struct Params <: UserParams @deftype Float64
    κ = 0.0 # a placeholder
end

@with_kw mutable struct Sim{D} <: Simulation{D} @deftype Float64
    # Add more parameters as necessary, or add to params (see examples)
    L::NTuple{D,Float64} # box length scales
    N::NTuple{D,Int64}  # grid points in each dimensions
    μ = 15.0    # chemical potential
    g = 0.1     # interaction parameter
    γ = 0.5; @assert γ >= 0.0 # damping parameter
    ti = 0.0    # initial time
    tf = 2/γ    # final time
    Nt::Int64 = 200     # number of saves over (ti,tf)
    params::UserParams = Params() # optional user parameters
    V0::Array{Float64,D} = zeros(N)
    t::LinRange{Float64} = LinRange(ti,tf,Nt) # time of saves
    ϕi::Array{Complex{Float64},D} = zeros(N) |> complex # initial condition
    alg::OrdinaryDiffEq.OrdinaryDiffEqAdaptiveAlgorithm = Tsit5() # default solver
    reltol::Float64 = 1e-6 # default tolerance; may need to use 1e-7 for corner cases
    flags::UInt32 = FFTW.MEASURE # choose a plan. PATIENT, NO_TIMELIMIT, EXHAUSTIVE
    # === saving
    nfiles::Bool = false
    path::String = nfiles ? joinpath(@__DIR__,"data") : @__DIR__
    filename::String = "save"
    # === arrays, transforms, spectral operators
    X::NTuple{D,Array{Float64,1}} = xvecs(L,N)
    K::NTuple{D,Array{Float64,1}} = kvecs(L,N)
    espec::Array{Complex{Float64},D} = 0.5*k2(K)
    T::TransformLibrary = makeT(X,K,flags=flags)
end

Sim(L,N,par) = Sim{length(L)}(L=L,N=N,params=par)
Sim(L,N) = Sim{length(L)}(L=L,N=N,params=Params())

FGPETypes = [Transforms Sim]
for type in FGPETypes
    Base.show(io::IO, ::MIME"application/prs.juno.inline",x::type) = x
end


text:
ENV["PYTHON"] = "/opt/conda/bin/python"

using Pkg

Pkg.activate()
Pkg.resolve()
Pkg.instantiate()
Pkg.precompile()

using IJulia
installkernel(
    "Julia",
    julia=`julia -t auto --project=@base -i --color=yes`,
)

# Force Gen to precompile
using Gen

ext: jl

max_stars_repo_name: jmuchovej/containers

max_issues_repo_name: jmuchovej/docker-bases

max_forks_repo_name: jmuchovej/docker-bases


text:
# traits for space dimension
abstract type SpaceDimension end
struct D1 <: SpaceDimension end
struct D2 <: SpaceDimension end
struct D3 <: SpaceDimension end

# traits for the type of confinement
abstract type ConfinementType end
struct BulkConfinement <: ConfinementType end
struct BoxConfinement <: ConfinementType end
struct SlabConfinement <: ConfinementType end
struct DiskConfinement <: ConfinementType end
struct SphereConfinement <: ConfinementType end
struct CylinderConfinement <: ConfinementType end
isconfined(::BulkConfinement) = false
isconfined(::ConfinementType) = true

# traits for the type of polymer system
abstract type PolymerSystemType end
struct NeatPolymer <: PolymerSystemType end
struct PolymerBlend <: PolymerSystemType end
struct PolymerSolution <: PolymerSystemType end

# traits for the type of charge distribution
abstract type ChargedType end
struct Neutral <: ChargedType end
struct SmearedCharge <: ChargedType end
struct DiscreteCharge <: ChargedType end
ischarged(::Neutral) = false
ischarged(::ChargedType) = true

abstract type AbstractSpecie end
struct KuhnSegment <: AbstractSpecie
    label::Symbol
    b::Real # length
    M::Real # molecular weight in g/mol
end
KuhnSegment(label; b=1.0, M=1.0) = KuhnSegment(label, b, M)

abstract type BlockEnd end

struct FreeEnd <: BlockEnd
    label::Symbol
end
FreeEnd(; label=:EF) = FreeEnd(label)

struct BranchPoint <: BlockEnd
    label::Symbol
end
isfreeblockend(::BlockEnd) = false
isfreeblockend(::FreeEnd) = true

abstract type AbstractBlock end
struct PolymerBlock <: AbstractBlock
    label::Symbol
    segment::KuhnSegment
    f::Real # = N_b / N, N is the total number of segments in a whole polymer chain
    E1::BlockEnd
    E2::BlockEnd
end

function PolymerBlock(; label=:A, specie=:A, f=1.0, E1=FreeEnd(), E2=FreeEnd())
    return PolymerBlock(label, KuhnSegment(specie), f, E1, E2)
end

"""
Check if the length of all blocks in a chain sum to 1.0.
"""
function _isachain(blocks)
    mapreduce(x->x.f, +, blocks) ≈ 1.0 ? true : false
end

abstract type AbstractMolecule end

struct SmallMolecule <: AbstractMolecule
    label::Symbol
    b::Real # length
    M::Real # molecular weight in g/mol
end
SmallMolecule(label; b=1.0, M=1.0) = SmallMolecule(label, b, M)

abstract type AbstractPolymer <: AbstractMolecule end

struct BlockCopolymer{T<:AbstractBlock} <: AbstractPolymer
    label::Symbol
    blocks::Vector{T}

    function BlockCopolymer(label, blocks::Vector{T}) where {T<:PolymerBlock}
        @argcheck _isachain(blocks)
        new{T}(label, blocks)
    end
end

nblocks(bc::BlockCopolymer) = length(bc.blocks)

struct RandomCopolymer <: AbstractPolymer end
struct AlternatingCopolymer <: AbstractPolymer end

struct Particle <: AbstractMolecule end
struct GiantMolecule <: AbstractMolecule end

abstract type AbstractComponent end

struct Component{T<:AbstractMolecule} <: AbstractComponent
    molecule::T
    α::Real # N / N_ref, N_ref is the total number of segments in a reference polymer chain
    ϕ::Real # = n*N/V, number density of this component in the system.

    function Component(molecule::T, α, ϕ) where {T<:AbstractMolecule}
        new{T}(molecule, α, ϕ)
    end
end
Component(molecule::T; α=1.0, ϕ=1.0) where {T<:AbstractMolecule} = Component(molecule, α, ϕ)

abstract type AbstractSystem end

"""
The key of `χN_map` should be a two-element `Set`. Each element is the unique symbol for a specie. For example, the `χN_map` of an AB diblock copolymer is a `Dict` with one entry `Set([:A, :B]) => χN`.

Note: For Edwards Model A (homopolymer + implicit solvent), one has to add a dummy solvent component in the `components` array to make the function `multicomponent` return correct result.
"""
struct PolymerSystem <: AbstractSystem
    components::Vector{AbstractComponent}
    confinement::ConfinementType
    χNmatrix::χNMatrix
    C::Real # = \rho_0 R_g^3 / N, dimensionless chain density

    function PolymerSystem(components::Vector{T}, conf, χNmatrix::χNMatrix, C) where {T<:AbstractComponent}
        @argcheck _isasystem(components)
        @argcheck _isasystem(components, χNmatrix)
        new(components, conf, χNmatrix, C)
    end
end

PolymerSystem(components::Vector{T}, χNmatrix::χNMatrix; conf=BulkConfinement(), C=1.0) where {T<:AbstractComponent} = PolymerSystem(components, conf, χNmatrix, C)

PolymerSystem(components::Vector{T}, χNmap; conf=BulkConfinement(), C=1.0) where {T<:AbstractComponent} = PolymerSystem(components, χNMatrix(χNmap); conf=conf, C=C)

"""
Check if the volume fraction of all compnents sums to 1.0.
"""
function _isasystem(components)
    return mapreduce(x->x.ϕ, +, components) ≈ 1.0 ? true : false
end

isconfined(s::PolymerSystem) = isconfined(s.confinement)

multicomponent(s::PolymerSystem) = length(s.components) == 1 ? false : true
ncomponents(s::PolymerSystem) = length(s.components)

specie(s::AbstractSpecie) = s.label
specie(m::SmallMolecule) = m.label
specie(b::PolymerBlock) = specie(b.segment)

species(c::BlockCopolymer) = [specie(b) for b in c.blocks] |> unique |> sort
nspecies(c::BlockCopolymer) = species(c) |> length
species(m::SmallMolecule) = [specie(m)]
nspecies(m::SmallMolecule) = 1
function _species(components)
    sps = Symbol[]
    for c in components
        append!(sps, species(c))
    end
    return unique(sps) |> sort
end

species(c::Component) = species(c.molecule)
nspecies(c::Component) = species(c) |> length
species(s::PolymerSystem) = _species(s.components)
nspecies(s::PolymerSystem) = species(s) |> length

function _isasystem(components, χNmatrix::χNMatrix)
    return _species(components) == species(χNmatrix)
end

function systemtype(s::PolymerSystem)
    if !multicomponent(s)
        return NeatPolymer()
    end
    for c in s.components
        if c.molecule isa SmallMolecule
            return PolymerSolution()
        end
    end
    return PolymerBlend()
end


text:
@testset "filter.jl" begin
    @testset "IdentityFilter" begin
        f = IdentityFilter()
        g = IdentityFilter()
        @test f === g
    end

    @testset "TopHatFilter" begin
        h(x) = 1 - 1 / 2 * cos(x)
        f = TopHatFilter(h)
        @test f.width == h
    end

    @testset "ConvolutionalFilter" begin
        h = x -> 1
        G = x -> x^2
        f = ConvolutionalFilter(h, G)
        @test f.width == h
        @test f.kernel == G
    end

    @testset "GaussianFilter" begin
        σ = 0.1
        G = GaussianFilter(σ)
        @test G.kernel(0.0) ≈ 1 / √(2π * σ^2)
    end
end


text:
# algorithm 3.6 in N&W
function zoom(ϕ, ϕ0, α0)
    function interpolate(a_i1::Real, a_i::Real,
                     ϕ_a_i1::Real, ϕ_a_i::Real,
                     dϕ_a_i1::Real, dϕ_a_i::Real)
        d1 = dϕ_a_i1 + dϕ_a_i - 3 * (ϕ_a_i1 - ϕ_a_i) / (a_i1 - a_i)
        d2 = sqrt(d1 * d1 - dϕ_a_i1 * dϕ_a_i)
        return a_i - (a_i - a_i1) *((dϕ_a_i + d2 - d1) /(dϕ_a_i - dϕ_a_i1 + 2 * d2))
    end
    ϕj = ϕ0
    αj = α0
    for i = 1:10
        interpolate()
        if

        else

        end

    end
    return αj
end


text:
module Transactions

export transaction, with_transaction, begin_transaction, commit_transaction, rollback_transaction

function begin_transaction end

function commit_transaction end

function rollback_transaction end

function transaction(f::Function)
  begin_transaction()

  try
    f()
  catch ex
    @error ex
    @info "Exception encountered during transaction, rolling back"

    rollback_transaction()
    rethrow(ex)
  end

  try
    commit_transaction()
  catch ex
    @error ex
    @info "Exception encountered when commiting transaction, rolling back"

    rollback_transaction()
    rethrow(ex)
  end
end
const with_transaction = transaction

end


text:
function rankfractional{T<:Real}(a::Array{T,1})
    r = float64(collect(1:length(a)))
    1.0 < r[end] || return r
    for i in ties(a)
        r[a.==i] = mean(r[a.==i])
    end
    return r
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using .Base.Order
import Base.Sort.QuickSortAlg
import Base.Sort.InsertionSortAlg
import Base.Sort.MergeSort
import Base.Order.lt

using .Base: sub_with_overflow, add_with_overflow

trunc2int(x) = unsafe_trunc(Int, x)
trunc2int(::Missing) = missing
_is_intable(x) = (typemin(Int) <= x <= typemax(Int)) && (round(x, RoundToZero) == x)
_is_intable(::Missing) = true
# x is sorted based on o
function _fill_starts!(ranges, x, rangescpy, last_valid_range, o::Ordering, ::Val{T}) where T

    cnt = 1
    st = 1
    lo::T = 0
    hi::T = 0
    @inbounds for j in 1:last_valid_range
        lo = rangescpy[j]
        j == last_valid_range ? hi = length(x) : hi = rangescpy[j+1] - 1
        cnt = _find_blocks_sorted!(ranges, x, lo, hi, cnt, o, Val(T))
    end
    @inbounds for j in 1:(cnt - 1)
        rangescpy[j] = ranges[j]
    end
    return cnt - 1
end

function _find_blocks_sorted!(ranges, x, lo, hi, cnt, o::Ordering, ::Val{T}) where T
    n = hi - lo + 1
    counter = 0
    st::T = lo
    ranges[cnt] = st
    cnt += 1
    @inbounds while true
        stopval::T = searchsortedlast(x, x[st], st, hi, o)
        # # the last obs in the current group
        st = stopval + 1
        st > hi && return cnt
        ranges[cnt] = st
        cnt += 1
        counter += 1
        #
        # # if too many levels switch strategy
        # #TODO supplied by should be take into account
        # # the decision is to always assume by = identity
        if counter > div(n, 2)
            # ranges[cnt] = st
            # cnt += 1
            for i in st:hi - 1
                if !isequal(x[i], x[i+1])
                # if lt(o, x[i], x[i+1])
                    ranges[cnt] = i + 1
                    cnt += 1
                end
            end
            return cnt
        end
    end
end

# inbits is zeros(Bool, length(x))
function _fill_starts_v2!(ranges, inbits, x, last_valid_range, o::Ordering, ::Val{T}) where T
    # first split x to chunks
    # if last_valid_range == 1
    #     @error "not yet implemented"
    # else
    #
    # inbit = Vector{}
    fill!(inbits, false)
    Threads.@threads for j in 1:last_valid_range
        inbits[ranges[j]] = true
    end
    Threads.@threads for j in 1:last_valid_range
        lo::T = 0
        hi::T = 0
        lo = ranges[j]
        j == last_valid_range ? hi = length(x) : hi = ranges[j+1] - 1
        _mark_start_of_groups_sorted!(inbits, x, lo, hi, o, Val(T))
    end
    cnt = 1
    @inbounds for i in 1:length(inbits)
        if inbits[i]
            ranges[cnt] = i
            cnt += 1
        end
    end
    cnt - 1
end

function _mark_start_of_groups_sorted!(inbits, x, lo, hi, o, ::Val{T}) where T
    n = hi - lo + 1
    n == 1 && return
    cp = ceil(Int, n/log2(n))
    # cp = div(n,2)
    counter = 0
    st::T = lo
    @inbounds while true
        stopval::T = searchsortedlast(x, x[st], st, hi, o)
        # # the last obs in the current group
        st = stopval + 1
        st > hi && break
        inbits[st] = true
        counter += 1
        #
        # # if too many levels switch strategy
        # #TODO supplied by should be take into account
        # # the decision is to always assume by = identity
        if counter > cp
            # ranges[cnt] = st
            # cnt += 1
            for i in st:hi - 1
                !isequal(x[i], x[i+1]) ? inbits[i + 1] = true : nothing
            end
            break
        end
    end
end


text:
using POMDPs
using ContinuousPOMDPTreeSearchExperiments
using POMDPToolbox
using POMCP
using POMCPOW
using ProgressMeter
using PmapProgressMeter
using ParticleFilters
using JLD
using CPUTime
using VDPTag
using Plots

N = 1000

@everywhere begin
    using VDPTag
    using POMDPs
    using ContinuousPOMDPTreeSearchExperiments
    using POMDPToolbox
    using POMCP
    using POMCPOW
    using ParticleFilters
    using CPUTime

    pomdp = VDPTagPOMDP()
    p_rng = MersenneTwister(4)
    solvers = Dict{Tuple, Union{Solver, Policy}}()
    for m in [0.01, 0.1, 1, 2, 5, 10]

        solvers = merge(solvers, Dict{Tuple, Union{Solver,Policy}}(

            ("pomcpow", m) => begin
                rollout_policy = ToNextML(mdp(pomdp))
                n = 10_000*m
                solver = POMCPOWSolver(tree_queries=n,
                                       criterion=MaxUCB(40.0),
                                       final_criterion=MaxTries(),
                                       max_depth=20,
                                       k_action=8.0,
                                       alpha_action=1/20,
                                       k_observation=4.0,
                                       alpha_observation=1/20,
                                       estimate_value=FORollout(rollout_policy),
                                       check_repeat_act=false,
                                       check_repeat_obs=false,
                                       rng=MersenneTwister(13)
                                      )
            end,

            ("bt_100", m) => begin
                n = 1000*m
                rng3 = copy(p_rng)
                rollout_policy = ToNextML(mdp(pomdp))
                node_updater = ObsAdaptiveParticleFilter(pomdp, LowVarianceResampler(100), 0.05, rng3)
                solver = POMCPDPWSolver(tree_queries=n,
                                        c=40.0,
                                        max_depth=20,
                                        k_action=8.0,
                                        alpha_action=1/20,
                                        k_observation=4.0,
                                        alpha_observation=1/20,
                                        estimate_value=FORollout(rollout_policy),
                                        node_belief_updater=node_updater,
                                        default_action=TagAction(false, 0.0),
                                        rng=rng3
                                       )
            end,

            # ("dpw", m) => begin
            #     rng3 = copy(p_rng)
            #     rollout_policy = ToNextML(mdp(pomdp))
            #     n = 10_000*m
            #     solver = POMCPDPWSolver(tree_queries=n,
            #                             c=40.0,
            #                             max_depth=20,
            #                             k_action=8.0,
            #                             alpha_action=1/20,
            #                             k_observation=4.0,
            #                             alpha_observation=1/20,
            #                             estimate_value=FORollout(rollout_policy),
            #                             default_action=TagAction(false, 0.0),
            #                             rng=rng3
            #                            )
            # end,
        ))
    end
end

solver_keys = collect(keys(solvers))
rewards = Dict{Tuple, AbstractVector{Float64}}()
# state_hists = Dict{String, AbstractVector{AbstractVector{state_type(pomdp)}}}()
times = Dict{Tuple, AbstractVector{Float64}}()
steps = Dict{Tuple, AbstractVector{Int}}()

for (j, sk) in enumerate(solver_keys)
    s_rewards = SharedArray(Float64, N)
    # s_hists = SharedArray(AbstractVector{state_type(pomdp)}, N)
    s_times = SharedArray(Float64, N)
    s_steps = SharedArray(Float64, N)
    prog = Progress(N, desc="$sk ($j of $(length(solver_keys)))...")
    pmap(prog, 1:N) do i
        sim_rng = MersenneTwister(i)
        up_rng = MersenneTwister(i+100_000)
        policy = solve(solvers[sk], pomdp)
        up = ObsAdaptiveParticleFilter(deepcopy(pomdp), LowVarianceResampler(10_000), 0.05, up_rng)
        sim = HistoryRecorder(max_steps=40, rng=sim_rng)
        s_times[i] = @CPUelapsed hist = simulate(sim, deepcopy(pomdp), policy, up)
        s_rewards[i] = discounted_reward(hist)
        s_steps[i] = n_steps(hist)
    end
    rewards[sk] = sdata(s_rewards)
    @show mean(rewards[sk])
    # state_hists[sk] = sdata(s_hists)
    times[sk] = sdata(s_times)
    steps[sk] = sdata(s_steps)
    @show sum(times[sk])/sum(steps[sk])
end

for k in solver_keys
    println("$k mean: $(mean(rewards[k])) sem: $(std(rewards[k])/sqrt(N))")
    println("$k time per step: $(sum(times[k])/sum(steps[k]))")
end

filename = Pkg.dir("ContinuousPOMDPTreeSearchExperiments", "data", "pomdp_trends_$(Dates.format(now(), "E_d_u_HH_MM")).jld")
println("saving to $filename...")
@save(filename, solver_keys, rewards, times, steps)
println("done.")

solver_types = unique([k[1] for k in solver_keys])
series = Dict{String, Any}()
for t in solver_types
    x = []
    y = []
    for (sk, rew) in rewards
        if sk[1] == t
            push!(x, sk[2])
            push!(y, mean(rew))
        end
    end
    p = sortperm(x)
    series[t] = (x[p], y[p])
end
unicodeplots()
for (t, p) in series
    plot(p[1], p[2], ylim=(-30, 100), title=t)
    gui()
end


text:
if !(@isdefined LinkSpherePairDistance)
    @computed struct LinkSpherePairDistance{Rm,R1,S,JC} <: TaskMap{Rm,R1,S}
        state::MechanismState{S,S,S,JC}
        state_cache::StateCache{S,JC}
        point1_cache::fulltype(FramePointCache{dim(Rm),S})
        point2_cache::fulltype(FramePointCache{dim(Rm),S})
        radius1::S
        radius2::S
    end
end
struct LinkSpherePairDistanceT{Rm,R1,S,JC} <: TaskMapT{Rm,R1,S}
    base_map::LinkSpherePairDistance{Rm,R1,S,JC}
end
function LinkSpherePairDistance{ℝ{m},R1}(state::MechanismState{S,S,S,JC}, state_cache,
        point1_cache::FramePointCache{m,S}, point2_cache::FramePointCache{m,S},
        radius1, radius2) where {m,R1,S,JC}
    LinkSpherePairDistance{ℝ{m},R1,S,JC}(state, state_cache, point1_cache, point2_cache,
        radius1, radius2)
end

function task_map_emb(::EmbRep, ::EmbRep, xme, task_map::LinkSpherePairDistance{ℝ{m},R1}) where m
    rframe = root_frame(task_map.state.mechanism)
    x1 = transform(task_map.state, task_map.point1_cache.point, rframe).v
    x2 = transform(task_map.state, task_map.point2_cache.point, rframe).v
    SA[norm(x2 - x1) - task_map.radius1 - task_map.radius2]
end

function task_jacobian_emb(xme, task_map::LinkSpherePairDistance{ℝ{m},R1}, arg...) where m
    state_cache = task_map.state_cache
    point1_cache, point2_cache = task_map.point1_cache, task_map.point2_cache
    link1_cache, link2_cache = point1_cache.link_cache, point2_cache.link_cache
    link1, link2 = link1_cache.link, link2_cache.link
    point1_frame_position, point2_frame_position = point1_cache.point.v, point2_cache.point.v
    radius1, radius2 = task_map.radius1, task_map.radius2

    Jl1c, Jl2c = link1_cache.transform_jacobian, link2_cache.transform_jacobian
    update_link_transform_jacobian!(Jl1c, xme, state_cache, link1)
    update_link_transform_jacobian!(Jl2c, xme, state_cache, link2)
    Jl1, Jl2 = Jl1c.data[1], Jl2c.data[1]

    Jp1c, Jp2c = point1_cache.world_position_jacobian, point2_cache.world_position_jacobian
    Rt1 = link_frame_transform(xme, state_cache, link1)
    Rt2 = link_frame_transform(xme, state_cache, link2)
    update_point_world_position_jacobian!(Jp1c, Rt1, point1_frame_position)
    update_point_world_position_jacobian!(Jp2c, Rt2, point2_frame_position)
    Jp1, Jp2 = Jp1c.data[1], Jp2c.data[1]

    Jlp1, Jlp2 = Jp1*Jl1, Jp2*Jl2
    Jlp = [Jlp1; Jlp2]

    point1_world_position = point_world_position(Rt1, point1_frame_position)
    point2_world_position = point_world_position(Rt2, point2_frame_position)
    world_positions = [point1_world_position; point2_world_position]
    Jn = stacked_position_sphere_distance_jacobian(world_positions, radius1, radius2)

    Jf = Jn*Jlp
end
task_jacobian_chart(xm, task_map::LinkSpherePairDistance{ℝ{m},R1},
    CM::Chart{1,ℝ{m}}, CN::Chart{1,R1}) where m = task_jacobian_emb(xm, task_map)

function task_jacobian_emb_dot(xme, vme, task_map::LinkSpherePairDistance{ℝ{m},R1}, arg...) where m
    state_cache = task_map.state_cache
    point1_cache, point2_cache = task_map.point1_cache, task_map.point2_cache
    link1_cache, link2_cache = point1_cache.link_cache, point2_cache.link_cache
    link1, link2 = link1_cache.link, link2_cache.link
    point1_frame_position, point2_frame_position = point1_cache.point.v, point2_cache.point.v
    radius1, radius2 = task_map.radius1, task_map.radius2

    Jl1c, Jl2c = link1_cache.transform_jacobian, link2_cache.transform_jacobian
    Jl1_dotc, Jl2_dotc = link1_cache.transform_jacobian_dot, link2_cache.transform_jacobian_dot
    update_link_transform_jacobian!(Jl1c, xme, state_cache, link1)
    update_link_transform_jacobian!(Jl2c, xme, state_cache, link2)
    update_link_transform_jacobian_dot!(Jl1_dotc, xme, vme, state_cache, link1)
    update_link_transform_jacobian_dot!(Jl2_dotc, xme, vme, state_cache, link2)
    Jl1, Jl2 = Jl1c.data[1], Jl2c.data[1]
    Jl1_dot, Jl2_dot = Jl1_dotc.data[1], Jl2_dotc.data[1]

    Jp1c, Jp2c = point1_cache.world_position_jacobian, point2_cache.world_position_jacobian
    Jp1_dotc = point1_cache.world_position_jacobian_dot
    Jp2_dotc = point2_cache.world_position_jacobian_dot
    Rt1 = link_frame_transform(xme, state_cache, link1)
    Rt2 = link_frame_transform(xme, state_cache, link2)
    Rt1_dot, Rt2_dot = Jl1*vme, Jl2*vme
    update_point_world_position_jacobian!(Jp1c, Rt1, point1_frame_position)
    update_point_world_position_jacobian!(Jp2c, Rt2, point2_frame_position)
    update_point_world_position_jacobian_dot!(Jp1_dotc, Rt1, Rt1_dot, point1_frame_position)
    update_point_world_position_jacobian_dot!(Jp2_dotc, Rt2, Rt2_dot, point2_frame_position)
    Jp1, Jp2 = Jp1c.data[1], Jp2c.data[1]
    Jp1_dot, Jp2_dot = Jp1_dotc.data[1], Jp2_dotc.data[1]

    Jlp1, Jlp2 = Jp1*Jl1, Jp2*Jl2
    Jlp1_dot = Jp1_dot*Jl1 + Jp1*Jl1_dot
    Jlp2_dot = Jp2_dot*Jl2 + Jp2*Jl2_dot

    Jlp = [Jlp1; Jlp2]
    Jlp_dot = [Jlp1_dot; Jlp2_dot]

    point1_world_position = point_world_position(Rt1, point1_frame_position)
    point2_world_position = point_world_position(Rt2, point2_frame_position)
    point1_world_velocity, point2_world_velocity = Jlp1*vme, Jlp2*vme
    world_positions = [point1_world_position; point2_world_position]
    world_velocities = [point1_world_velocity; point2_world_velocity]
    Jn = stacked_position_sphere_distance_jacobian(world_positions, radius1, radius2)
    Jn_dot = stacked_position_sphere_distance_jacobian_dot(world_positions, world_velocities,
        radius1, radius2)

    Jf_dot = Jn_dot*Jlp + Jn*Jlp_dot
end
task_jacobian_chart_dot(xm, task_map::LinkSpherePairDistance{ℝ{m},R1},
    CM::Chart{1,ℝ{m}}, CN::Chart{1,R1}) where m = task_jacobian_emb_dot(xm, task_map)

# Distances between points, with positions input as a single stacked vector
function stacked_position_distance(positions)
    SA[norm(positions[static(1):static(3)] - positions[static(4):static(6)])]
end

function stacked_position_distance_jacobian(positions)
    f = (positions) -> stacked_position_distance(positions)
    ForwardDiff.jacobian(f, positions)
end

function stacked_position_distance_jacobian_dot(xs, xs_dot)
    m, n = 6, 1
    ∇xf = ForwardDiff.jacobian(xs -> stacked_position_distance_jacobian(xs),
        xs)::SArray{Tuple{m,m},eltype(xs),2,m*m}
    Jfdot = reshape(∇xf*xs_dot, Size(n,m))
end

# Distances between spheres, with positions input as a single stacked vector
function stacked_position_sphere_distance(positions, radius1, radius2)
    SA[norm(positions[static(1):static(3)] - positions[static(4):static(6)]) - radius1 - radius2]
end

function stacked_position_sphere_distance_jacobian(positions, radius1, radius2)
    f = (positions) -> stacked_position_sphere_distance(positions, radius1, radius2)
    ForwardDiff.jacobian(f, positions)
end

function stacked_position_sphere_distance_jacobian_dot(xs, xs_dot, radius1, radius2)
    m, n = 6, 1
    ∇xf = ForwardDiff.jacobian(xs -> stacked_position_sphere_distance_jacobian(xs, radius1,
        radius2), xs)::SArray{Tuple{m,m},eltype(xs),2,m*m}
    Jfdot = reshape(∇xf*xs_dot, Size(n,m))
end

# For verification
# A bit slower than custom function after precomputation and caching
function task_jacobian_emb2(xme, task_map::LinkSpherePairDistance{ℝ{m},R1}, arg...) where m
    S = eltype(xme)
    state = task_map.state
    rframe = root_frame(state.mechanism)
    point1_cache, point2_cache = task_map.point1_cache, task_map.point2_cache
    link1_cache, link2_cache = point1_cache.link_cache, point2_cache.link_cache
    point1, point2 = point1_cache.point, point2_cache.point
    point1_world = transform(state, point1, rframe)
    point2_world = transform(state, point2, rframe)
    x1, x2 = point1_world.v, point2_world.v

    Jf1 = SMatrix{3,m,S}(point_jacobian(state, link1_cache.path, point1_world).J)
    Jf2 = SMatrix{3,m,S}(point_jacobian(state, link2_cache.path, point2_world).J)

    SMatrix{1,m,S}((Jf2' - Jf1')*(x2 - x1)/norm(x2-x1))
end

function task_map_emb_for_jdot(xme, task_map::LinkSpherePairDistance{ℝ{m},R1}) where m
    state = task_map.state_cache[eltype(xme)]
    rframe = root_frame(state.mechanism)
    set_configuration!(state, xme)
    x1 = transform(state, task_map.point1_cache.point, rframe).v
    x2 = transform(state, task_map.point2_cache.point, rframe).v
    SA[norm(x2 - x1) - task_map.radius1 - task_map.radius2]
end
task_jacobian_emb_for_jdot(xme, task_map::LinkSpherePairDistance{ℝ{m},R1}) where m =
    ForwardDiff.jacobian(xme -> task_map_emb_for_jdot(xme, task_map), xme)

function task_jacobian_emb_dot2(xme, vme, task_map::LinkSpherePairDistance{ℝ{m},R1}, arg...) where m
    n, S = embdim(codomain_manifold(task_map)), eltype(xme)
    ∇xf = ForwardDiff.jacobian(xme -> SVector{m*n,S}(
        reshape(task_jacobian_emb_for_jdot(xme, task_map), m*n)),
        xme)::SArray{Tuple{m*n,m},S,2,m*m*n}
    Jfdot_emb = reshape(∇xf*vme, Size(n,m))
end

task_jacobian_chart_dot(xm, vm, task_map::LinkSpherePairDistance{ℝ{m},R1},
    CM::Chart{1,ℝ{m}}, CN::Chart{1,R1}) where m = task_jacobian_emb_dot(xm, vm, task_map)


text:
#UTILITY PROGRESSMETER, REPORTS WORKER NUMBER AND CURRENT ITERATE
mutable struct ProgressNS{T<:Real} <: AbstractProgress
    interval::T
    dt::AbstractFloat
    start_it::Integer
    counter::Integer
    triggered::Bool
    tfirst::AbstractFloat
    tlast::AbstractFloat
    tstp::AbstractFloat
    printed::Bool        # true if we have issued at least one status update
    desc::AbstractString # prefix to the percentage, e.g.  "Computing..."
    color::Symbol        # default to green
    output::IO           # output stream into which the progress is written
    numprintedvalues::Integer   # num values printed below progress in last iteration
    offset::Integer             # position offset of progress bar (default is 0)
    total_step::AbstractFloat
    information::AbstractFloat
    etc::AbstractFloat
    contour::AbstractFloat
    max_lh::AbstractFloat
    naive::AbstractFloat
    li_dist::Vector{AbstractFloat}
    stepworker::Integer
    workers::Vector{Integer}
    wk_totals::Vector{Integer}
    total_li_delta::Vector{AbstractFloat}
    wk_li_delta::Matrix{AbstractFloat}
    wk_instruction::Vector{Integer}
    model_exhaust::Vector{Integer}
    wk_eff::Vector{Vector{AbstractFloat}}
    wk_jobs::Vector{String}
    inst_counters::Vector{Integer}
    mean_stp_time::AbstractFloat
    eff_iterates::Integer
    no_displayed_srcs::Integer
    sources::Vector{Tuple{Matrix{AbstractFloat},Integer}}
    mix::BitMatrix
    job_times::Vector{AbstractFloat}
    job_i1_ct::Vector{Integer}

    function ProgressNS{T}(    naive::AbstractFloat,
                               interval::T,
                               workers::Vector{Integer};
                               dt::Real=0.1,
                               eff_iterates::Integer,
                               desc::AbstractString="Nested Sampling::",
                               color::Symbol=:green,
                               output::IO=stderr,
                               offset::Int=0,
                               start_it::Int=1,
                               no_displayed_srcs::Integer=0) where T
        tfirst = tlast = time()
        printed = false
        new{T}(interval,
         dt,
         start_it,
         start_it,
         false,
         tfirst,
         tlast,
         0.,
         printed,
         desc,
         color,
         output,
         0,
         offset,
         0.,0.,0.,0.,0.,
         naive,
         [0.],
         0,
         workers,
         zeros(Integer,length(workers)),
         zeros(length(workers)),
         zeros(length(workers), 65),
         zeros(Integer,length(workers)),
         zeros(Integer,length(workers)),
         [[0.] for i in 1:length(workers)],
         ["none" for worker in 1:length(workers)],
         zeros(Integer,9),
         0.,
         eff_iterates,
         no_displayed_srcs,
         [(zeros(0,0),0)],
         falses(0,0),
         zeros(9),
         zeros(Integer,9))
    end
end

ProgressNS(naive::AbstractFloat, interval::Real, workers::Vector{Integer}, dt::Real=0.1, desc::AbstractString="Nested Sampling::",
         color::Symbol=:green, output::IO=stderr, offset::Integer=0, start_it::Integer=1, eff_iterates=250) =
            ProgressNS{typeof(interval)}(naive, interval, workers, dt=dt, eff_iterates=eff_iterates, desc=desc, color=color, output=output, offset=offset, start_it=start_it)

ProgressNS(naive::AbstractFloat, interval::Real, workers::Vector{Integer}, desc::AbstractString, offset::Integer=0, start_it::Integer=1; eff_iterates=250, no_displayed_srcs=1) = ProgressNS{typeof(interval)}(naive, interval, workers, desc=desc, offset=offset, start_it=start_it, eff_iterates=eff_iterates, no_displayed_srcs=no_displayed_srcs)

function update!(p::ProgressNS, contour, max, val, thresh, info, li_dist, worker, wk_time, job, model, old_li, new_li, instruction; sources=[(zeros(0,0),0)], bitmatrix=falses(0,0), options...)

    instruction == "PS" && (p.inst_counters[1]+=1)
    instruction == "PS" && job==1 && (p.job_i1_ct[1]+=1; p.job_times[1]+=wk_time)
    instruction == "PM" && (p.inst_counters[2]+=1)
    instruction == "PM"  && job==1 && (p.job_i1_ct[2]+=1; p.job_times[2]+=wk_time)
    instruction == "PSFM" && (p.inst_counters[3]+=1)
    instruction == "PSFM" && job==1 && (p.job_i1_ct[3]+=1; p.job_times[3]+=wk_time)
    instruction == "FM" && (p.inst_counters[4]+=1)
    instruction == "FM" && job==1 && (p.job_i1_ct[4]+=1; p.job_times[4]+=wk_time)
    instruction == "DM" && (p.inst_counters[5]+=1)
    instruction == "DM" && job==1 && (p.job_i1_ct[5]+=1; p.job_times[5]+=wk_time)
    instruction == "SM" && (p.inst_counters[6]+=1)
    instruction == "SM" && job==1 && (p.job_i1_ct[6]+=1; p.job_times[6]+=wk_time)
    instruction == "RD" && (p.inst_counters[7]+=1)
    instruction == "RD" && job==1 && (p.job_i1_ct[7]+=1; p.job_times[7]+=wk_time)
    instruction == "RI" && (p.inst_counters[8]+=1)
    instruction == "RI" && job==1 && (p.job_i1_ct[8]+=1; p.job_times[8]+=wk_time)
    instruction == "EM" && (p.inst_counters[9]+=1)
    instruction == "EM" && job==1 && (p.job_i1_ct[9]+=1; p.job_times[9]+=wk_time)

    p.counter += 1
    stps_elapsed=p.counter-p.start_it
    p.tstp=time()-p.tlast
    p.mean_stp_time=(p.tlast-p.tfirst)/stps_elapsed

    widx=findfirst(isequal(worker), p.workers)
    p.wk_totals[widx]+=1
    li_delta=new_li-old_li
    p.total_li_delta[widx]+=new_li-old_li

    p.wk_li_delta[widx,1:end-1]=p.wk_li_delta[widx,2:end]
    p.wk_li_delta[widx,end]=li_delta

    push!(p.wk_eff[widx],li_delta/wk_time)
    length(p.wk_eff[widx])>p.eff_iterates && (p.wk_eff[widx]=p.wk_eff[widx][end-p.eff_iterates+1:end]) #keep worker efficiency array under iterate size limit
    p.wk_instruction[widx]=job; p.model_exhaust[widx]=model
    p.wk_jobs[widx]=instruction

    p.contour = contour
    p.max_lh = max

    interval = val - thresh
    step = p.interval - interval
    !isinf(step) && step>0 && (p.total_step+=step)

    p.etc= (p.interval/(p.total_step/stps_elapsed))*p.mean_stp_time

    p.interval=interval
    p.information = info
    p.li_dist=li_dist
    p.stepworker = worker
    p.sources=sources
    p.mix=bitmatrix
    updateProgress!(p; options...)
end

function updateProgress!(p::ProgressNS; showvalues = Any[], valuecolor = :blue, offset::Integer = p.offset, keep = (offset == 0))
    p.offset = offset
    t = time()
    if p.interval <= 0 && !p.triggered
        p.triggered = true
        if p.printed
            p.triggered = true
            wk_msgs = [@sprintf "Wk:%g:: I:%i, %s M:%i S:%.2f" p.workers[widx] p.wk_instruction[widx] p.wk_jobs[widx] p.model_exhaust[widx] (p.wk_totals[widx]/(p.counter-p.start_it)) for widx in 1:length(p.workers)]
            wk_inst=UnicodePlots.boxplot(wk_msgs, p.wk_eff, title="Worker Diagnostics", xlabel="Likelihood surface velocity (sec^-1)", color=:magenta)

            lh_heatmap=UnicodePlots.heatmap(p.wk_li_delta[end:-1:1,:], xoffset=-size(p.wk_li_delta,2)-1, colormap=:viridis, title="Worker lhΔ history", xlabel="Lh stride/step")

            dur = ProgressMeter.durationstring(t-p.tfirst)
            msg1 = @sprintf "%s Converged. Time: %s (%d iterations). H: %s" p.desc dur p.counter p.information
            print(p.output, "\n" ^ (p.offset + p.numprintedvalues))
            msg2 = @sprintf "Ensemble Stats:: Contour: %g MaxLH: %g Max/Naive: %g H: %g" p.contour p.max_lh (p.max_lh-p.naive) p.information
            hist=UnicodePlots.histogram(p.li_dist, title="Ensemble Likelihood Distribution", color=:green)
            srclines=p.no_displayed_srcs+1
            p.numprintedvalues=nrows(wk_inst.graphics)+nrows(lh_heatmap.graphics)+nrows(hist.graphics)+16+srclines
            print(p.output, "\n" ^ (p.offset + p.numprintedvalues))
            ProgressMeter.move_cursor_up_while_clearing_lines(p.output, p.numprintedvalues)
            show(p.output, wk_inst)
            print(p.output, "\n")
            show(p.output, lh_heatmap)
            print(p.output, "\n")
            p.no_displayed_srcs>0 && printsources(p)
            ProgressMeter.printover(p.output, msg1, :magenta)
            print(p.output, "\n")
            ProgressMeter.printover(p.output, msg2, p.color)
            print(p.output, "\n")
            show(p.output, hist)
            print(p.output, "\n")

            if keep
                println(p.output)
            else
                print(p.output, "\r\u1b[A" ^ (p.offset + p.numprintedvalues))
            end
        end
        return
    end

    if t > p.tlast+p.dt && !p.triggered
        wk_msgs = [@sprintf "Wk:%g:: I:%i, %s M:%i S:%.2f" p.workers[widx] p.wk_instruction[widx] p.wk_jobs[widx] p.model_exhaust[widx] (p.wk_totals[widx]/(p.counter-p.start_it)) for widx in 1:length(p.workers)]
        wk_inst=UnicodePlots.boxplot(wk_msgs, p.wk_eff, title="Worker Diagnostics", xlabel="Likelihood surface velocity (sec^-1)", color=:magenta)

        jobtime_msg= @sprintf "PS:%s|PM:%s|PSFM:%s|FM:%s|DM:%s|SM:%s|RD:%s|RI:%s|EM:%s" hmss(p.job_times[1]/p.job_i1_ct[1]) hmss(p.job_times[2]/p.job_i1_ct[2]) hmss(p.job_times[3]/p.job_i1_ct[3]) hmss(p.job_times[4]/p.job_i1_ct[4]) hmss(p.job_times[5]/p.job_i1_ct[5]) hmss(p.job_times[6]/p.job_i1_ct[6]) hmss(p.job_times[7]/p.job_i1_ct[7]) hmss(p.job_times[8]/p.job_i1_ct[8]) hmss(p.job_times[9]/p.job_i1_ct[9])

        lh_heatmap=UnicodePlots.heatmap(p.wk_li_delta[end:-1:1,:], xoffset=-size(p.wk_li_delta,2)-1, colormap=:viridis, title="Worker lhΔ history", xlabel="Lh stride/step")

        msg1 = @sprintf "%s Step %i::Wk:%g: PS:%s|PM:%s|PSFM:%s|FM:%s|DM:%s|SM:%s|RD:%s|RI:%s|EM:%s" p.desc p.counter p.stepworker p.inst_counters[1] p.inst_counters[2] p.inst_counters[3] p.inst_counters[4] p.inst_counters[5] p.inst_counters[6] p.inst_counters[7] p.inst_counters[8] p.inst_counters[9]
        msg2 = @sprintf "Step time μ, last Δ: %s,%s Convergence Interval: %g ETC: %s" hmss(p.mean_stp_time) hmss(p.tstp-p.mean_stp_time) p.interval hmss(p.etc)
        msg3 = @sprintf "Ensemble Stats:: Contour: %g MaxLH: %g Max/Naive: %g H: %g" p.contour p.max_lh (p.max_lh-p.naive) p.information

        hist=UnicodePlots.histogram(p.li_dist, title="Ensemble Likelihood Distribution", color=:green)

        #p.numprintedvalues=nrows(wk_inst.graphics)+nrows(hist.graphics)+nrows(lh_heatmap.graphics)+1
        srclines=p.no_displayed_srcs+1
        p.numprintedvalues=nrows(wk_inst.graphics)+nrows(lh_heatmap.graphics)+nrows(hist.graphics)+18+srclines
        print(p.output, "\n" ^ (p.offset + p.numprintedvalues))
        ProgressMeter.move_cursor_up_while_clearing_lines(p.output, p.numprintedvalues)
        show(p.output, wk_inst)
        print(p.output, "\n")
        ProgressMeter.printover(p.output, jobtime_msg, :magenta)
        print(p.output, "\n")
        show(p.output, lh_heatmap)
        print(p.output, "\n")
        p.no_displayed_srcs>0 && printsources(p)
        ProgressMeter.printover(p.output, msg1, :magenta)
        print(p.output, "\n")
        ProgressMeter.printover(p.output, msg2, :cyan)
        print(p.output, "\n")
        ProgressMeter.printover(p.output, msg3, p.color)
        print(p.output, "\n")
        show(p.output, hist)
        print(p.output, "\n")

        print(p.output, "\r\u1b[A" ^ (p.offset + p.numprintedvalues))

        # Compensate for any overhead of printing. This can be
        # especially important if you're running over a slow network
        # connection.
        p.tlast = t + 2*(time()-t)
        p.printed = true
    end
end

                function hmss(dt)
                    isnan(dt) && return "NaN"
                    (h,r) = divrem(dt,60*60)
                    (m,r) = divrem(r, 60)
                    (isnan(h)||isnan(m)||isnan(r)) && return "NaN"
                    string(Int(h),":",Int(m),":",Int(ceil(r)))
                end

                function printsources(p; freqsort=true)
                    printidxs=Vector{Integer}()
                    printsrcs=Vector{Matrix{AbstractFloat}}()
                    printfreqs=Vector{AbstractFloat}()

                    if freqsort
                        freqs=vec(sum(p.mix,dims=1)); total=size(p.mix,1)
                        sortfreqs=sort(freqs,rev=true)
                        sortidxs=sortperm(freqs)
                        for srcidx in 1:p.no_displayed_srcs
                            push!(printidxs, sortidxs[srcidx])
                            push!(printsrcs, p.sources[sortidxs[srcidx]][1])
                            push!(printfreqs, sortfreqs[srcidx]/total)
                        end
                    else
                        freqs=vec(sum(p.mix,dims=1)); total=size(p.mix,1)
                        for srcidx in 1:p.no_displayed_srcs
                            push!(printidxs, srcidx)
                            push!(printsrcs, p.sources[srcidx][1])
                            push!(printfreqs, freqs[srcidx]/total)
                        end
                    end

                    printstyled(p.output, "MLE Top Sources\n", bold=true)

                    for src in 1:p.no_displayed_srcs
                        print(p.output, "S$(printidxs[src]), $(printfreqs[src]*100)%: ")
                        pwmstr_to_io(p.output, printsrcs[src])
                        print(p.output, "\n")
                    end
                end


text:
using NorahPhd
using Test

@testset "NorahPhd.jl" begin
    # Write your tests here.
end


text:
using ..DBUtils
using SQLite
using MySQL
using DataFrames
using Dates

"""
    create_tables!(conn)

Create and initialize tables to save results from an Entrez/PubMed search or a medline file load.
Caution, all related tables are dropped if they exist
"""
function create_tables!(conn)

    # Determine engine
    sql_engine = (typeof(conn)== MySQL.Connection) ? MySQL : SQLite
    AUTOINCREMENT = (sql_engine == MySQL) ? "AUTO_INCREMENT" : "AUTOINCREMENT"
    engine_info = (sql_engine == MySQL) ? "ENGINE=InnoDB DEFAULT CHARSET=utf8mb4" : ""

   #purge related tables
   DBUtils.disable_foreign_checks(conn)
   sql_engine.execute!(conn, "DROP TABLE IF EXISTS basic")
   sql_engine.execute!(conn, "DROP TABLE IF EXISTS author2article")
   sql_engine.execute!(conn, "DROP TABLE IF EXISTS author_ref")
   sql_engine.execute!(conn, "DROP TABLE IF EXISTS pub_type")
   sql_engine.execute!(conn, "DROP TABLE IF EXISTS abstract_full")
   sql_engine.execute!(conn, "DROP TABLE IF EXISTS abstract_structured")
   sql_engine.execute!(conn, "DROP TABLE IF EXISTS mesh_heading")
   sql_engine.execute!(conn, "DROP TABLE IF EXISTS mesh_desc")
   sql_engine.execute!(conn, "DROP TABLE IF EXISTS mesh_qual")
   sql_engine.execute!(conn, "DROP TABLE IF EXISTS file_meta")
   sql_engine.execute!(conn, "DROP TABLE IF EXISTS mesh2umls")
   DBUtils.enable_foreign_checks(conn)

   #Create tables to store
   sql_engine.execute!(conn, "CREATE TABLE IF NOT EXISTS `basic` (
         `pmid` int(9) NOT NULL,
         `pub_year` smallint(4) DEFAULT NULL,
         `pub_month` tinyint(2) DEFAULT NULL,
         `pub_dt_desc` varchar(50) DEFAULT NULL,
         `title` text DEFAULT NULL,
         `authors` text DEFAULT NULL,
         `journal_title` varchar(500) DEFAULT NULL,
         `journal_ISSN` varchar(9) DEFAULT NULL,
         `journal_volume` varchar(30) DEFAULT NULL,
         `journal_issue` varchar(30) DEFAULT NULL,
         `journal_pages` varchar(50) DEFAULT NULL,
         `journal_iso_abbreviation` varchar(255) DEFAULT NULL,
         `url` varchar(100) DEFAULT NULL,
         `ins_dt_time` timestamp DEFAULT CURRENT_TIMESTAMP,
         PRIMARY KEY (`pmid`)
       ) $engine_info;"
       )

   sql_engine.execute!(conn, "CREATE TABLE `author_ref` (
         `pmid` int(9) NOT NULL,
         `last_name` varchar(60) DEFAULT NULL,
         `first_name` varchar(60) DEFAULT NULL,
         `initials` varchar(10) DEFAULT NULL,
         `suffix` varchar(10) DEFAULT NULL,
         `orcid` varchar(19) DEFAULT NULL,
         `collective` varchar(200) DEFAULT NULL,
         `affiliation` text DEFAULT NULL,
         `ins_dt_time` timestamp DEFAULT CURRENT_TIMESTAMP,
         FOREIGN KEY(`pmid`) REFERENCES basic(`pmid`)
       ) $engine_info;"
       )

   # sql_engine.execute!(conn, "CREATE TABLE IF NOT EXISTS `author2article` (
   #       `pmid` int(9) NOT NULL,
   #       `author_id` int(10) NOT NULL,
   #       `ins_dt_time` timestamp DEFAULT CURRENT_TIMESTAMP,
   #       PRIMARY KEY (`pmid`, `author_id`),
   #       FOREIGN KEY(`pmid`) REFERENCES basic(`pmid`),
   #       FOREIGN KEY(`author_id`) REFERENCES author_ref(`author_id`)
   #     ) ENGINE=InnoDB DEFAULT CHARSET=utf8mb4;"
   #     )

   sql_engine.execute!(conn, "CREATE TABLE IF NOT EXISTS `pub_type` (
         `pmid` int(9) NOT NULL,
         `uid` int(6) NOT NULL,
         `name` varchar(100) NOT NULL,
         `ins_dt_time` timestamp DEFAULT CURRENT_TIMESTAMP,
         PRIMARY KEY (`pmid`,`uid`),
         FOREIGN KEY(`pmid`) REFERENCES basic(`pmid`)
       ) $engine_info;"
       )

   sql_engine.execute!(conn, "CREATE TABLE IF NOT EXISTS `abstract_full` (
         `pmid` int(9) NOT NULL,
         `abstract_text` text DEFAULT NULL,
         `ins_dt_time` timestamp DEFAULT CURRENT_TIMESTAMP,
         PRIMARY KEY (`pmid`),
         FOREIGN KEY(`pmid`) REFERENCES basic(`pmid`)
       ) $engine_info;"
       )

   sql_engine.execute!(conn, "CREATE TABLE IF NOT EXISTS `abstract_structured` (
         `abstract_structured_id` INTEGER PRIMARY KEY $AUTOINCREMENT,
         `pmid` int(9) NOT NULL,
         `nlm_category` varchar(20) DEFAULT NULL,
         `label` varchar(40) DEFAULT NULL,
         `abstract_text` text DEFAULT NULL,
         `ins_dt_time` timestamp DEFAULT CURRENT_TIMESTAMP,
         FOREIGN KEY(`pmid`) REFERENCES basic(`pmid`)
       ) $engine_info;"
       )

   sql_engine.execute!(conn, "CREATE TABLE IF NOT EXISTS `file_meta` (
         `file_name` varchar(30) NOT NULL,
         `ins_start_time` timestamp NOT NULL DEFAULT CURRENT_TIMESTAMP,
         `ins_end_time` timestamp NULL,
         PRIMARY KEY (`file_name`)
       ) $engine_info;"
       )

   # # --------------------------
   # #  MeshHeading Tables
   # # --------------------------
   # # --
   # # Descriptor
   # # The id corresponds to the DUI of mesh library
   # # Adding a "D" at the beginning of the id, allows for
   # # lookup in the mesh browerser
   # #  https://www.nlm.nih.gov/mesh/MBrowser.html

   # Qualifier
   sql_engine.execute!(conn, "CREATE TABLE IF NOT EXISTS `mesh_desc` (
         `uid` int(6) NOT NULL,
         `name` varchar(255) NOT NULL,
         `ins_dt_time` timestamp DEFAULT CURRENT_TIMESTAMP,
         PRIMARY KEY (`uid`)
       ) $engine_info;"
       )

   # Heading
   sql_engine.execute!(conn, "CREATE TABLE `mesh_qual` (
         `uid` int(6) NOT NULL,
         `name` varchar(255) NOT NULL,
         `ins_dt_time` timestamp DEFAULT CURRENT_TIMESTAMP,
         PRIMARY KEY (`uid`)
       ) $engine_info;"
       )

   sql_engine.execute!(conn, "CREATE TABLE IF NOT EXISTS `mesh_heading` (
         `pmid` int(9) NOT NULL,
         `desc_uid` int(6) NOT NULL,
         `desc_maj_status` tinyint(1) NOT NULL,
         `qual_uid` int DEFAULT NULL,
         `qual_maj_status` tinyint DEFAULT NULL,
         `ins_dt_time` timestamp DEFAULT CURRENT_TIMESTAMP,
         FOREIGN KEY(`pmid`) REFERENCES basic(`pmid`),
         FOREIGN KEY(`desc_uid`) REFERENCES mesh_desc(`uid`),
         FOREIGN KEY(`qual_uid`) REFERENCES mesh_qual(`uid`)
       ) $engine_info"
       )

end

function create_post_tables!(conn::MySQL.Connection)

    # engine
    engine_info = "ENGINE=InnoDB DEFAULT CHARSET=utf8mb4"

   #purge related tables
   DBUtils.disable_foreign_checks(conn)
   MySQL.execute!(conn, "DROP TABLE IF EXISTS author2article")
   MySQL.execute!(conn, "DROP TABLE IF EXISTS author")
   DBUtils.enable_foreign_checks(conn)

   MySQL.execute!(conn, "CREATE TABLE `author` (
         `auth_id` int(12) NOT NULL AUTO_INCREMENT,
         `last_name` varchar(60) DEFAULT NULL,
         `first_name` varchar(60) DEFAULT NULL,
         `initials` varchar(10) DEFAULT NULL,
         `suffix` varchar(10) DEFAULT NULL,
         `orcid` varchar(19) DEFAULT NULL,
         `collective` varchar(200) DEFAULT NULL,
         `affiliation` text DEFAULT NULL,
         `ins_dt_time` timestamp DEFAULT CURRENT_TIMESTAMP,
         PRIMARY KEY (`auth_id`)
       ) $engine_info;"
       )

   #Create tables to store
   MySQL.execute!(conn, "CREATE TABLE IF NOT EXISTS `author2article` (
         `pmid` int(9) NOT NULL,
         `auth_id` int(12) NOT NULL,
         PRIMARY KEY (`pmid`,`auth_id`),
         FOREIGN KEY (`pmid`) REFERENCES basic(`pmid`),
         FOREIGN KEY (`auth_id`) REFERENCES author(`auth_id`)
       ) $engine_info;"
       )

end

"""
    add_mysql_keys!(conn)

Adds indices/keys to MySQL PubMed tables.
"""
function add_mysql_keys!(conn::MySQL.Connection)

    res = MySQL.Query(conn, "SHOW INDEX FROM basic WHERE key_name = 'pub_year'") |> DataFrame
    size(res)[1] == 1 && return nothing

    MySQL.execute!(conn, "ALTER TABLE `basic`
          ADD KEY `pub_year` (`pub_year`),
          ADD KEY `pub_month` (`pub_month`),
          ADD KEY `pub_year_month` (`pub_year`,`pub_month`),
          ADD KEY `journal_title` (`journal_title`),
          ADD KEY `journal_ISSN` (`journal_ISSN`)
        ;"
        )

    MySQL.execute!(conn, "ALTER TABLE `author_ref`
          ADD KEY `last_name` (`last_name`),
          ADD KEY `last_first_name` (`last_name`, `first_name`),
          ADD KEY `collective` (`collective`)
        ;"
        )

    MySQL.execute!(conn, "ALTER TABLE `pub_type`
          ADD KEY `name` (`name`)
        ;"
        )

    MySQL.execute!(conn, "ALTER TABLE `abstract_structured`
          ADD KEY `label` (`label`),
          ADD KEY `nlm_category` (`nlm_category`)
        ;"
        )

    MySQL.execute!(conn, "ALTER TABLE `mesh_desc`
          ADD KEY `name` (`name`)
        ;"
        )

    MySQL.execute!(conn, "ALTER TABLE `mesh_qual`
          ADD KEY `name` (`name`)
        ;"
        )

   return nothing
end

"""
    drop_mysql_keys!(conn)

Removes keys/indices from MySQL PubMed tables.
"""
function drop_mysql_keys!(conn::MySQL.Connection)

    res = MySQL.Query(conn, "SHOW INDEX FROM basic WHERE key_name = 'pub_year'", DataFrame)
    size(res)[1] == 0 && return nothing

    MySQL.execute!(conn, "ALTER TABLE `basic`
          DROP KEY `pub_year`,
          DROP KEY `pub_month`,
          DROP KEY `pub_year_month`,
          DROP KEY `journal_title`,
          DROP KEY `journal_ISSN`
        ;"
        )

    MySQL.execute!(conn, "ALTER TABLE `author_ref`
          DROP KEY `last_name`,
          DROP KEY `last_first_name`,
          DROP KEY `collective`
        ;"
        )

    MySQL.execute!(conn, "ALTER TABLE `pub_type`
          DROP KEY `name`
        ;"
        )

    MySQL.execute!(conn, "ALTER TABLE `abstract_structured`
          DROP KEY `label`,
          DROP KEY `nlm_category`
        ;"
        )

    MySQL.execute!(conn, "ALTER TABLE `mesh_desc`
          DROP KEY `name`
        ;"
        )

    MySQL.execute!(conn, "ALTER TABLE `mesh_qual`
          DROP KEY `name`
        ;"
        )

    MySQL.execute!(conn, "ALTER TABLE `mesh_heading`
          DROP KEY `pmid_uids`
        ;"
        )
end

"""
    create_pmid_table!(conn; tablename="article")

Creates a table, using either MySQL of SQLite, to store PMIDs from
Entrez related searches. All tables are empty at this point
"""
function create_pmid_table!(conn)

    # Determine engine
    sql_engine = (typeof(conn) == MySQL.Connection) ? MySQL : SQLite

    #purge related tables
    # sql_engine.execute!(conn, "DROP TABLE IF EXISTS $tablename")

    #Create tables to store
    sql_engine.execute!(conn, "CREATE TABLE IF NOT EXISTS basic (
                                pmid INTEGER NOT NULL PRIMARY KEY
                                );"
                        )
end

"""
    all_pmids(db)

Return all PMIDs stored in the *basic* table of the input database
"""
function all_pmids(conn)
    sql_engine = (typeof(conn)== MySQL.Connection) ? MySQL : SQLite
    query = sql_engine.Query(conn, "SELECT pmid FROM basic;") |> DataFrame
    return query[:pmid]
end

"""
    all_mesh(db)

Return all MeSH stored in the *mesh_desc* table of the input database
"""
function all_mesh(db)
    sel = db_query(db, "SELECT name FROM mesh_desc;")
    return sel[:name]
end

"""
    abstracts_by_year(db, pub_year; local_medline=false)

Return all abstracts of article published in the given year.
If local_medline flag is set to true, it is assumed that db contains *article*
table with only PMIDs and all other info is available in a (same host) medline database
"""
function abstracts_by_year(db, pub_year; local_medline=false, uid_str = "pmid")

    #get all abstracts UNIQUE pairs
    query_code = ""
    if local_medline
    query_code = "SELECT article.$uid_str as $uid_str,
                         medline.abstract_full.abstract_text as abstract_text
                    FROM article
              INNER JOIN medline.basic ON medline.basic.pmid = article.pmid
              INNER JOIN medline.abstract_full ON medline.abstract_full.pmid = article.pmid
                   WHERE medline.basic.pub_year = $pub_year "
    else
        query_code = "SELECT at.$uid_str at $uid_str,
                             at.abstract as abstract_text
                        FROM abstract_text as at
                  INNER JOIN basic ar on at.pmid = ar.pmid
                       WHERE ar.pub_year = $pub_year "
    end

    sel = db_query(db, query_code)
    num_abstracts = size(sel)[1]
    println("Retrieved: ", num_abstracts, " abstracts")
    return sel
end

"""
    abstracts(db; local_medline=false)

Return all abstracts related to PMIDs in the database.
If local_medline flag is set to true, it is assumed that db contains *basic*
table with only PMIDs and all other info is available in a (same host) medline database
"""
function abstracts(db; local_medline=false)

    #get all abstracts UNIQUE pairs
    query_code = ""
    if local_medline
        query_code = "SELECT basic.pmid as pmid,
                             medline.abstract_full.abstract_text as abstract_text
                        FROM basic
                  INNER JOIN medline.abstract_full ON medline.abstract_full.pmid = basic.pmid"
    else
        query_code = "SELECT ar.pmid as pmid,
                             ar.abstract_text as abstract_text
                        FROM abstract_full as ar "
    end

    sel = db_query(db, query_code)
    num_abstracts = size(sel)[1]
    @info "Retrieved abstracts: " num_abstracts
    return sel
end

"""
    get_article_mesh(db, pmid)

Get the all mesh-descriptors associated with a given article
"""
function get_article_mesh(db, pmid::Integer)

    query_string = "SELECT md.name
                      FROM mesh_heading as mh,
                           mesh_desc as md
                     WHERE mh.desc_uid = md.uid
                      AND mh.pmid = $pmid"
    query  = db_query(db, query_string)

    #return data array
    return query[:name]

end

"""
    get_article_mesh_by_concept(db, pmid, umls_concepts...; local_medline)

Get the all mesh-descriptors associated with a given article
## Arguments:
* query_string: "" - assumes full set of results were saved by BioMedQuery directly from XML
"""
function get_article_mesh_by_concept(db, pmid::Integer, umls_concepts...; query_string="")

    concept_set_string = "'" * join(umls_concepts,"','") * "'"

    if query_string == ""
        query_string = "SELECT DISTINCT(md.name), m2u.umls
                        FROM mesh_desc AS md
                        JOIN mesh_heading AS mh ON mh.desc_uid = md.uid
                        JOIN mesh2umls AS m2u ON m2u.mesh = md.name
                        WHERE mh.pmid = $pmid
                        AND m2u.umls IN  ($(concept_set_string))"
    end

    query  = db_query(db, query_string)

    #return data array
    return query[:name]

end

"""
    db_insert!(conn, articles::Dict{String,DataFrame}, csv_path=pwd(), csv_prefix="<current date>_PubMed_"; verbose=false, drop_csvs=false)

Writes dictionary of dataframes to a MySQL database.  Tables must already exist (see PubMed.create_tables!).  CSVs that are created during writing can be saved (default) or removed.
"""
function db_insert!(db::MySQL.Connection, articles::Dict{String,DataFrame}, csv_path::String = pwd(), csv_prefix::String = "$(Date(now()))_PubMed_"; verbose=false, drop_csv=false)

    dfs_to_csv(articles, csv_path, csv_prefix)

    #Insert csv prefix into files_meta talbe
    meta_sql = """INSERT INTO file_meta (file_name,ins_start_time) VALUES ('$csv_prefix',CURRENT_TIMESTAMP)"""
    MySQL.execute!(db, meta_sql)

    for (table, df) in articles
        # check if column names all exist in mysql table
        if !col_match(db, table, df)
            error("Each DataFrame column must match the name of a table column. $table had mismatches.")
        end

        path = joinpath(csv_path, "$(csv_prefix)$(table).csv")
        cols_string = assemble_cols(df)

        # Save article data (MySQL.stream from df)
        ins_sql = """LOAD DATA LOCAL INFILE '$path' INTO TABLE $table CHARACTER SET utf8mb4 FIELDS TERMINATED BY ',' OPTIONALLY ENCLOSED BY '"' IGNORE 1 LINES ($cols_string)"""
        MySQL.execute!(db, ins_sql)
    end

    meta_sql = """UPDATE file_meta SET ins_end_time = CURRENT_TIMESTAMP WHERE file_name = '$csv_prefix'"""
    MySQL.execute!(db, meta_sql)

    if drop_csv
        remove_csvs(articles, csv_path, csv_prefix)
    end

    return nothing

end

"""
    db_insert!(conn, csv_path=pwd(), csv_prefix="<current date>_PubMed_"; verbose=false, drop_csvs=false)

Writes CSVs from PubMed parsing to a MySQL database.  Tables must already exist (see PubMed.create_tables!).  CSVs can optionally be removed after being written to DB.
"""
function db_insert!(db::MySQL.Connection, csv_path::String = pwd(), csv_prefix::String = "$(Date(now()))_PubMed_"; verbose=false, drop_csv=false)
    paths = Vector{String}()

    #Insert csv prefix into files_meta talbe
    meta_sql = """INSERT INTO file_meta (file_name,ins_start_time) VALUES ('$csv_prefix',CURRENT_TIMESTAMP)"""
    MySQL.execute!(db, meta_sql)

    for table in select_all_tables(db)
        # for all non-file_meta tables
        if table != "file_meta"
            path = joinpath(csv_path, "$(csv_prefix)$(table).csv")
            drop_csv && push!(paths,path)

            headers = CSV.read(path, limit=1)

            cols = String.(getfield(headers, :colindex).names)
            if !col_match(db, table, cols)
                error("Each CSV column must match the name of a table column. $table had mismatches.")
            end

            cols_string = join(cols, ",")

            # Save article data (MySQL.stream from df)
            ins_sql = """LOAD DATA LOCAL INFILE '$path' INTO TABLE $table CHARACTER SET utf8mb4 FIELDS TERMINATED BY ',' OPTIONALLY ENCLOSED BY '"' IGNORE 1 LINES ($cols_string)"""
            MySQL.execute!(db, ins_sql)
        end
    end

    meta_sql = """UPDATE file_meta SET ins_end_time = CURRENT_TIMESTAMP WHERE file_name = '$csv_prefix'"""
    MySQL.execute!(db, meta_sql)

    if drop_csv
        remove_csvs(paths)
    end

    return nothing

end

function db_insert!(db::MySQL.Connection, pmid::Int64, articles::Dict{String,DataFrame}, csv_path::String = pwd(), csv_prefix::String = "$(Date(now()))_PubMed_"; verbose=false, drop_csv=false)

    dfs_to_csv(articles, csv_path, csv_prefix)

    for (table, df) in articles
        if occursin(r"$mesh*", table)
            # check if column names all exist in mysql table
            if !col_match(db, table, df)
                error("Each DataFrame column must match the name of a table column. $table had mismatches.")
            end

            path = joinpath(csv_path, "$(csv_prefix)$(table).csv")
            cols_string = assemble_cols(df)

            # Save article data (MySQL.stream from df)
            ins_sql = """LOAD DATA LOCAL INFILE '$path' INTO TABLE $table CHARACTER SET utf8mb4 FIELDS TERMINATED BY ',' OPTIONALLY ENCLOSED BY '"' IGNORE 1 LINES ($cols_string)"""
            MySQL.execute!(db, ins_sql)
        end
    end

    if drop_csv
        remove_csvs(articles, csv_path, csv_prefix)
    end

    return nothing

end

"""
    db_insert!(conn, articles::Dict{String,DataFrame}, csv_path=pwd(), csv_prefix="<current date>_PubMed_"; verbose=false, drop_csvs=false)

Writes dictionary of dataframes to a SQLite database.  Tables must already exist (see PubMed.create_tables!).  CSVs that are created during writing can be saved (default) or removed.
"""
function db_insert!(db::SQLite.DB, articles::Dict{String,DataFrame}, csv_path::String = pwd(), csv_prefix::String = "$(Date(now()))_PubMed_"; verbose=false, drop_csv=false)

    #Insert csv prefix into files_meta talbe
    meta_sql = """INSERT INTO file_meta (file_name,ins_start_time) VALUES ('$csv_prefix',CURRENT_TIMESTAMP)"""
    SQLite.execute!(db, meta_sql)

    for (table, df) in articles

        # check if column names all exist in mysql table
        if !col_match(db, table, df)
            error("Each DataFrame column must match the name of a table column. $table had mismatches.")
        end

        for i = 1:size(df)[1]
            col_dict = Dict{Symbol,Any}()
            for col in getfield(df, :colindex).names
                col_dict[col] = df[i,col]
            end
            insert_row!(db, table, col_dict)
        end

    end

    meta_sql = """UPDATE file_meta SET ins_end_time = CURRENT_TIMESTAMP WHERE file_name = '$csv_prefix'"""
    SQLite.execute!(db, meta_sql)

    return nothing

end


text:
const STAGE_TABLE = UInt8(0)
const STAGE_INIT_NOISY = UInt8(1)
const STAGE_GOOD_NOISY = UInt8(2)
const STAGE_KILLER_1 = UInt8(3)
const STAGE_KILLER_2 = UInt8(4)
const STAGE_COUNTER = UInt8(5)
const STAGE_INIT_QUIET = UInt8(6)
const STAGE_QUIET = UInt8(7)
const STAGE_BAD_NOISY = UInt8(8)
const STAGE_DONE = UInt8(9)

const NORMAL_TYPE = UInt8(1)
const NOISY_TYPE = UInt8(2)

const MVVLVA_VALS = @SVector Int32[100, 450, 450, 675, 1300, 5000]

MoveOrder() = MoveOrder(NORMAL_TYPE, STAGE_TABLE, MoveStack(MAX_MOVES), zeros(Int32, MAX_MOVES),
    0, 0, 0, MOVE_NONE, MOVE_NONE, MOVE_NONE, MOVE_NONE)

function setmargin!(moveorder::MoveOrder, margin::Int)
    moveorder.margin = margin
end

function clear!(moveorder::MoveOrder)
    moveorder.type = NORMAL_TYPE
    moveorder.stage = STAGE_TABLE
    clear!(moveorder.movestack)
    moveorder.noisy_size = 0
    moveorder.quiet_size = 0
end

# find the index of the best scored move in the movestack
function idx_bestmove(moveorder::MoveOrder, idx_start::Int, idx_end::Int)
    best_idx = idx_start
    best_val = moveorder.values[best_idx]
    @inbounds for i in (idx_start + 1):idx_end
        if (cv = moveorder.values[i]) >= best_val
            best_idx = i
            best_val = cv
        end
    end
    return best_idx
end

function popmove!(moveorder::MoveOrder, idx::Int)
    move = moveorder.movestack[idx]
    @inbounds for i in idx:moveorder.movestack.idx
        moveorder.movestack.list[i] = moveorder.movestack.list[i + 1]
        moveorder.values[i] = moveorder.values[i + 1]
    end
    moveorder.movestack.idx -= 1
    if idx <= moveorder.noisy_size
        moveorder.noisy_size -= 1
    elseif moveorder.noisy_size < idx
        moveorder.quiet_size -= 1
    end
    return move
end

# function sortmoves!(moveorder::MoveOrder, idx_start::Int, idx_end::Int)
#     idx = Vector{Int}(undef, max(0, idx_end-idx_start+1))
#     m = view(moveorder.movestack, idx_start:idx_end)
#     v = view(moveorder.values, idx_start:idx_end)
#     sortperm!(idx, v, rev = true)
#     permute!(m, idx)
#     permute!(v, idx)
# end

function MVVLVA!(moveorder::MoveOrder, board::Board)
    @inbounds for i in 1:moveorder.noisy_size
        move = moveorder.movestack[i]
        sqr_from = from(move)
        sqr_to = to(move)
        ptype_from = type(board[sqr_from])
        ptype_to = type(board[sqr_to])

        if ptype_to !== VOID
            moveorder.values[i] = MVVLVA_VALS[ptype_to.val] - ptype_from.val
        end

        _flag = flag(move)
        if _flag <= __QUEEN_CASTLE
            continue
        elseif _flag === __ENPASS
            moveorder.values[i] = MVVLVA_VALS[1] - one(Int32)
        elseif _flag === __QUEEN_PROMO
            moveorder.values[i] += MVVLVA_VALS[5]
        end
    end
    return
end

function init_noisy_moveorder!(thread::Thread, ply::Int, margin::Int)
    board = thread.board
    moveorder = thread.moveorders[ply + 1]
    moveorder.type = NOISY_TYPE
    moveorder.margin = margin
    # Set move heuristics
    moveorder.tt_move = MOVE_NONE
    moveorder.counter = MOVE_NONE
    moveorder.killer1 = MOVE_NONE
    moveorder.killer2 = MOVE_NONE
    # Generate moves
    gen_noisy_moves!(moveorder.movestack, board)
    moveorder.noisy_size = moveorder.movestack.idx
    moveorder.stage = STAGE_INIT_NOISY
    return
end

function init_normal_moveorder!(thread::Thread, tt_move::Move, ply::Int)
    board = thread.board
    moveorder = thread.moveorders[ply + 1]
    moveorder.type = NORMAL_TYPE
    moveorder.margin = 0
    # Set move heuristics
    moveorder.tt_move = tt_move
    if ply > 0
        @inbounds previous_move = thread.movestack[ply]
        @inbounds previous_piece = thread.piecestack[ply]
    else
        previous_move = MOVE_NONE
        previous_piece = VOID
    end
    previous_to = to(previous_move)
    if (previous_move === MOVE_NONE || previous_move === NULL_MOVE)
         moveorder.counter = MOVE_NONE
    else
        @inbounds moveorder.counter = thread.cmtable[(!board.turn).val][previous_piece.val][previous_to]
    end
    @inbounds moveorder.killer1 = thread.killer1s[ply + 1]
    @inbounds moveorder.killer2 = thread.killer2s[ply + 1]
    # generate moves
    gen_noisy_moves!(moveorder.movestack, board)
    moveorder.noisy_size = moveorder.movestack.idx
    gen_quiet_moves!(moveorder.movestack, board)
    moveorder.quiet_size = moveorder.movestack.idx - moveorder.noisy_size
    moveorder.stage = STAGE_TABLE
    return
end

function selectmove!(thread::Thread, ply::Int, skipquiets::Bool)::Move
    @inbounds moveorder = thread.moveorders[ply + 1]
    board = thread.board

    # First pick the transposition table move.
    # However, due to collisions and current TT implementation it may not always be legal.
    # Therefore verify legality, by *first* generating possible moves, and checking it exists.
    if moveorder.stage === STAGE_TABLE
        moveorder.stage = STAGE_INIT_NOISY
        if (moveorder.tt_move !== MOVE_NONE) && (moveorder.tt_move ∈ moveorder.movestack)
            return moveorder.tt_move
        end
    end

    # Score the noisy moves
    # Score the moves using MVV-LVA heuristics.
    if moveorder.stage === STAGE_INIT_NOISY
        MVVLVA!(moveorder, board)
        # sortmoves!(moveorder, 1, moveorder.noisy_size)
        moveorder.stage = STAGE_GOOD_NOISY
    end

    # Pick the noisy moves which pass MVV-LVA, and a static_exchange_evaluator runby.
    # Do not play TT moves twice.
    if moveorder.stage === STAGE_GOOD_NOISY
        while moveorder.noisy_size > 0
            idx = idx_bestmove(moveorder, 1, moveorder.noisy_size)
            if moveorder.values[idx] >= 0
                if !ischeck(board) && (static_exchange_evaluator(board, moveorder.movestack[idx], moveorder.margin) == false)
                    moveorder.values[idx] = -1
                    continue
                end
                move = popmove!(moveorder, idx)
                if move === moveorder.tt_move
                    continue
                else
                    # Avoid playing killer or counter moves twice.
                    (move == moveorder.killer1) && (moveorder.killer1 = MOVE_NONE)
                    (move == moveorder.killer2) && (moveorder.killer2 = MOVE_NONE)
                    (move == moveorder.counter) && (moveorder.counter = MOVE_NONE)
                    return move
                end
            else
                break
            end
        end
        moveorder.stage = STAGE_KILLER_1
    end

    # If noisy type, we can stop here.
    if moveorder.type === NOISY_TYPE
        return MOVE_NONE
    end

    # If skipquiets flag is set, now is the time to skip a few stages.
    if skipquiets
        moveorder.stage = STAGE_BAD_NOISY
    end

    # First killer move stage.
    if moveorder.stage === STAGE_KILLER_1
        if (moveorder.killer1 !== moveorder.tt_move) && (moveorder.killer1 ∈ moveorder.movestack)
            moveorder.stage = STAGE_KILLER_2
            return moveorder.killer1
        else
            moveorder.stage = STAGE_KILLER_2
        end
    end

    # Second killer move stage.
    if moveorder.stage === STAGE_KILLER_2
        if (moveorder.killer2 !== moveorder.tt_move) && (moveorder.killer2 ∈ moveorder.movestack)
            moveorder.stage = STAGE_COUNTER
            return moveorder.killer2
        else
            moveorder.stage = STAGE_COUNTER
        end
    end

    # Counter moves stage.
    if moveorder.stage === STAGE_COUNTER
        if (moveorder.counter !== moveorder.tt_move) && (moveorder.counter !== moveorder.killer1) &&
            (moveorder.counter !== moveorder.killer2) && (moveorder.counter ∈ moveorder.movestack)
            moveorder.stage = STAGE_INIT_QUIET
            return moveorder.counter
        else
            moveorder.stage = STAGE_INIT_QUIET
        end
    end

    # Score all quiet moves.
    # Score using the history heuristics.
    if moveorder.stage === STAGE_INIT_QUIET
        gethistoryscores!(thread, moveorder.movestack, moveorder.values, moveorder.noisy_size + 1, moveorder.movestack.idx, ply)
        # sortmoves!(moveorder, moveorder.noisy_size + 1, moveorder.movestack.idx)
        moveorder.stage = STAGE_QUIET
    end

    # Pick the best quiet moves.
    # Do not play TT move twice.
    if moveorder.stage === STAGE_QUIET
        if moveorder.quiet_size !== 0
            idx = idx_bestmove(moveorder, moveorder.noisy_size + 1, moveorder.movestack.idx)
            move = popmove!(moveorder, idx)
            if (move == moveorder.tt_move) || ((move == moveorder.killer1) ||
                (move == moveorder.killer1) || (move == moveorder.counter))
                return selectmove!(thread, ply, skipquiets)
            else
                return move
            end
        else
            moveorder.stage = STAGE_BAD_NOISY
        end
    end

    # Lastly, if we are left with noisy moves that failed STAGE_GOOD_NOISY, we can pick them now.
    if moveorder.stage === STAGE_BAD_NOISY
        if (moveorder.noisy_size !== 0)
            move = popmove!(moveorder, 1)
            if (move == moveorder.tt_move) || (!skipquiets && ((move == moveorder.killer1) ||
                (move == moveorder.killer2) || (move == moveorder.counter)))
                return selectmove!(thread, ply, skipquiets)
            else
                return move
            end
        else
            moveorder.stage = STAGE_DONE
        end
    end

    # We are done
    if moveorder.stage === STAGE_DONE
        return MOVE_NONE
    end

    # If return type == Nothing, something went wrong.
end

max_stars_repo_name: J-Revell/Astellarn

max_issues_repo_name: J-Revell/Astellarn

max_forks_repo_name: J-Revell/Astellarn


text:
@testset "I343" begin
    if VERSION ≥ v"0.6.0-dev.2746"
        @testset "Type Alias Shadowing" begin
            s = """
            const TT = Int64
            @compat SharedVector{X} = SharedArray{X,1}

            type MyType{X}
                t::X
                MyType{X}(x) where X = new(convert(X, x))
            end
            """
            msgs = lintstr(s)
            @test messageset(msgs) == Set([:I343])
            @test msgs[1].variable == "SharedVector"
        end
    end

    msgs = lintstr("""
    e = 1
    """)
    @test messageset(msgs) == Set([:I343])
    @test contains(msgs[1].message, "with same name as export from Base")

    @test isempty(lintstr("""
    import Base: parent
    immutable MyType end
    parent(_::MyType) = MyType()
    """))
end


text:
# FASTQ Reader
# ============

"""
    FASTQReader(input::IO, quality_encoding=:sanger, fill_ambiguous=nothing)

Create a data reader of the FASTQ file format.

# Arguments
* `input`: data source
* `quality_encoding=:sanger`: encoding of base qualities; see the following table for available values
* `fill_ambiguous=nothing`: fill ambiguous nucleotides with the given nucleotide

| Quality encoding | Symbol           | ASCII offset | Quality range |
|:---------------- |:---------------- |:------------:|:-------------:|
|  Sanger          | `:sanger`        |     +33      |      0-93     |
|  Solexa          | `:solexa`        |     +64      |     -5-62     |
|  Illumina 1.3+   | `:illumina13`    |     +64      |      0-62     |
|  Illumina 1.5+   | `:illumina15`    |     +64      |      2-62     |
|  Illumina 1.8+   | `:illumina18`    |     +33      |      0-93     |
"""
type FASTQReader{S<:Sequence} <: Bio.IO.AbstractReader
    state::Ragel.State
    seqbuf::BufferedOutputStream{BufferedStreams.EmptyStream}
    qualbuf::BufferedOutputStream{BufferedStreams.EmptyStream}
    name2buf::StringField
    desc2buf::StringField
    qualcount::Int
    quality_encoding::QualityEncoding
    fill_ambiguous::Nullable{DNA}

    function FASTQReader(input::BufferedInputStream, quality_encoding, fill_ambiguous)
        return new(Ragel.State(fastqparser_start, input),
                   BufferedOutputStream(), BufferedOutputStream(),
                   StringField(), StringField(), 0, quality_encoding, fill_ambiguous)
    end
end

function Bio.IO.stream(reader::FASTQReader)
    return reader.state.stream
end

function Base.eltype{S}(::Type{FASTQReader{S}})
    return FASTQSeqRecord{S}
end

function FASTQReader(input::IO; quality_encoding=:sanger, fill_ambiguous=nothing)
    return FASTQReader{DNASequence}(BufferedInputStream(input), sym2qualenc(quality_encoding), fill_ambiguous)
end

function (::Type{FASTQReader{S}}){S<:Sequence}(input::IO;
                                               quality_encoding=:sanger,
                                               fill_ambiguous=nothing)
    return FASTQReader{S}(BufferedInputStream(input), sym2qualenc(quality_encoding), fill_ambiguous)
end

# Map the quality encoding name to a corresponding QualityEncoding object.
function sym2qualenc(name::Symbol)
    if name == :sanger
        return SANGER_QUAL_ENCODING
    elseif name == :solexa
        return SOLEXA_QUAL_ENCODING
    elseif name == :illumina13
        return ILLUMINA13_QUAL_ENCODING
    elseif name == :illumina15
        return ILLUMINA15_QUAL_ENCODING
    elseif name == :illumina18
        return ILLUMINA18_QUAL_ENCODING
    end
    throw(ArgumentError("quality encoding ':$(name)' is not supported"))
end


text:
"""
    exchange_leastsq!(v, x::BEDFile, y, k) -> ELSQResults

    Execute the exchange algorithm using a `BEDFile` object `x` for one model size `k`.
"""
function exchange_leastsq!{T <: Float}(
    v        :: ELSQVariables{T},
    x        :: BEDFile{T},
    y        :: SharedVector{T},
    kernfile :: String,
    k        :: Int,
    w        :: PlinkGPUVariables{T};
    pids     :: DenseVector{Int} = procs(x),
    n        :: Int  = length(y),
    p        :: Int  = size(x,2),
    window   :: Int  = k,
    max_iter :: Int  = 100,
    tol      :: T    = convert(T, 1e-6),
    quiet    :: Bool = false
)
    # initial value for previous residual sum of squares
    old_rss = oftype(tol,Inf)

    # obtain top r components of bvec in magnitude
    selectperm!(v.perm, v.b, k, by=abs, rev=true, initialized=true)

    # update estimated response
    update_indices!(v.idx, v.b)
    A_mul_B!(v.xb, x, v.b, v.idx, k, v.mask_n)

    # update residuals
    difference!(v.r, y, v.xb)
    mask!(v.r, v.mask_n, 0, zero(T))

    # save value of RSS before starting algorithm
    rss = sumabs2(v.r) / 2

    # compute inner products of x and residuals
    # this is basically the negative gradient
    At_mul_B!(v.df, x, v.r, v.mask_n, w)

    # outer loop controls number of total iterations for algorithm run on one r
    for iter = 1:(max_iter)

        # output algorithm progress to console
        quiet || println("\titer = ", iter, ", RSS = ", rss)

        # middle loop tests each of top r parameters (by magnitude?)
        for i = abs(k-window+1):k

            # save information for current value of i
            l     = v.perm[i]
            betal = v.b[l]
            decompress_genotypes!(v.tempn, x, l) # tempn now holds X[:,l]

            # if necessary, compute inner products of current predictor against all other predictors
            # store this information in Dict inner
            # for current index, hold dot products in memory for duration of inner loop
            get_inner_product!(v.dotprods, v.tempn, v, x, w, l, pids=pids)

            # subroutine compares current predictor i against all predictors k+1, k+2, ..., p
            # these predictors are candidates for inclusion in set
            # _exlstsq_innerloop! find best new predictor r
            a, b, r, adb = _exlstsq_innerloop!(v, k, i, p, tol)

            # now want to update residuals with current best predictor
            m = update_current_best_predictor!(v, x, betal, adb, r)

            # if necessary, compute inner product of current predictor against all other predictors
            # save in our Dict for future reference
            get_inner_product!(v.tempp, v.tempn2, v, x, w, m, pids=pids)

            # also update df
            axpymbz!(v.df, betal, v.dotprods, adb, v.tempp)

            # now swap best predictor with current predictor
            _swap_predictors!(v, i, r, m, adb)

        end # end middle loop over predictors

        # update residual sum of squares
        rss = sumabs2(v.r) / 2

        # test for descent failure
        # if no descent failure, then test for convergence
        # if not converged, then save RSS and check finiteness
        # if not converged and still finite, then save RSS and continue
        ascent    = rss > old_rss + tol
        converged = abs(old_rss - rss) / abs(old_rss + 1) < tol
        old_rss = rss
        ascent && print_descent_error(iter, old_rss, rss)
        converged && return nothing

        check_finiteness(rss)

    end # end outer iteration loop

    # at this point, maximum iterations reached
    # warn and return
    print_maxiter(max_iter, rss)
    return nothing

end # end exchange_leastsq!

"""
    exlstsq(x::BEDFile, y, models) -> ELSQResults

    Execute the exchange algorithm with a `BEDFile` object `x` over a vector `models` of model sizes.
"""
function exlstsq{T <: Float}(
    x        :: BEDFile{T},
    y        :: SharedVector{T},
    kernfile :: String;
    models   :: DenseVector{Int} = collect(1:min(20,size(x,2))),
    pids     :: DenseVector{Int} = procs(x),
    v        :: ELSQVariables{T} = ELSQVariables(x, y),
    w        :: PlinkGPUVariables{T} = PlinkGPUVariables(v.df, x, y, kernfile, v.mask_n),
    window   :: Int  = max(20, min(maximum(models), size(x,2))),
    max_iter :: Int  = 100,
    tol      :: T    = convert(T, 1e-6),
    quiet    :: Bool = true
)
    # dimensions of problem
    nmodels = length(models)
    n,p = size(x)

    # error checking
    errorcheck(x, y, tol, max_iter, window, p)

    # initialize sparse matrix to return
    betas = spzeros(T, p, nmodels)

    # loop through models
    for i in models
        exchange_leastsq!(v, x, y, kernfile, i, w,  window=min(i,window), max_iter=max_iter, tol=tol, quiet=quiet, n=n, p=p, pids=pids)
        betas[:,i] = sparse(sdata(v.b))
    end

    # return matrix of betas
    return betas
end

"""
    one_fold(x::BEDFile ,y, moels, folds, fold) -> Vector{Float}

    Compute one crossvalidation fold with the exchange algorithm using a `BEDFile` object `x`.
"""
function one_fold{T <: Float}(
    x        :: BEDFile{T},
    y        :: SharedVector{T},
    models   :: DenseVector{Int},
    kernfile :: String,
    folds    :: DenseVector{Int},
    fold     :: Int;
    tol      :: T    = convert(T, 1e-6),
    max_iter :: Int  = 100,
    window   :: Int  = max(20, min(maximum(models), size(x,2))),
    quiet    :: Bool = true
)

    # find testing indices
    # how big is training sample?
    test_idx = folds .== fold
    test_size = sum(test_idx)

    # train_idx is the vector that indexes the TRAINING set
    train_idx = convert(Vector{Int}, !test_idx)
    test_idx  = convert(Vector{Int}, test_idx)

    # preallocate temporary variables
    v = ELSQVariables(x, y, train_idx)
    w = PlinkGPUVariables(v.df, x, y, kernfile, v.mask_n)

    # will return a sparse matrix of betas
    betas = exlstsq(x, y, kernfile, models=models, v=v, w=w, window=window, max_iter=max_iter, tol=tol, quiet=quiet)

    # preallocate vector for output
    mses = zeros(T, length(models))

    # compute the mean out-of-sample error for the TEST set
    for i in eachindex(models)
        # set b
        copy!(v.b, view(betas, :, i))

        # update indices of current b
        update_indices!(v.idx, v.b)

        # compute estimated response with current b
        A_mul_B!(v.xb, x, v.b, v.idx, models[i], test_idx)

        # compute residuals
        difference!(v.r, y, v.xb)
        mask!(v.r, test_idx, 0, zero(T))

        # compute out-of-sample error as squared residual averaged over size of test set
        mses[i] = sumabs2(v.r) / (2*test_size)
    end

    return mses :: Vector{T}
end

"""
    pfold(x, y, path, folds, q) -> Vector

This function is the parallel execution kernel in `cv_exlstsq`. It is not meant to be called outside of `cv_exlstsq`.
It will distribute `q` crossvalidation folds across the processes supplied by the optional argument `pids` and call `one_fold` for each fold.
Each fold will compute a regularization path `1:path`.
`pfold` collects the vectors of MSEs returned by calling `one_fold` for each process, reduces them, and returns their average across all folds.
"""
function pfold(
    T        :: Type,
    xfile    :: String,
    xtfile   :: String,
    x2file   :: String,
    yfile    :: String,
    meanfile :: String,
    precfile :: String,
    models   :: DenseVector{Int},
    kernfile :: String,
    folds    :: DenseVector{Int},
    q        :: Int;
    pids     :: DenseVector{Int} = procs(),
    tol      :: Float = convert(T, 1e-6),
    max_iter :: Int   = 100,
    window   :: Int   = 20,
    quiet    :: Bool  = true,
    header   :: Bool  = false,
)
    # how many CPU processes can pfold use?
    np = length(pids)

    # report on CPU processes
    quiet || println("pfold: np = ", np)
    quiet || println("pids = ", pids)

    # set up function to share state (indices of folds)
    i = 1
    nextidx() = (idx=i; i+=1; idx)

    # preallocate cell array for results
    results = Array{Any}(q)

    # master process will distribute tasks to workers
    # master synchronizes results at end before returning
    @sync begin

        # loop over all workers
        for worker in pids

            # exclude process that launched pfold, unless only one process is available
            if worker != myid() || np == 1

                # asynchronously distribute tasks
                @async begin
                    while true

                        # grab next fold
                        current_fold = nextidx()

                        # if current fold exceeds total number of folds then exit loop
                        current_fold > q && break

                        # report distribution of fold to worker and device
                        quiet || print_with_color(:blue, "Computing fold $current_fold on worker $worker.\n\n")

                        # launch job on worker
                        # worker loads data from file paths and then computes the errors in one fold
                        results[current_fold] = remotecall_fetch(worker) do
                            pids = [worker]
                            x = BEDFile(T, xfile, xtfile, x2file, meanfile, precfile, pids=pids, header=header)
                            y = SharedArray(abspath(yfile), T, (x.geno.n,), pids=pids)
                            one_fold(x, y, models, kernfile, folds, current_fold, max_iter=max_iter, window=window, tol=tol, quiet=quiet)
                        end # end remotecall_fetch()
                    end # end while
                end # end @async
            end # end if
        end # end for
    end # end @sync

    # return reduction (row-wise sum) over results
    return (reduce(+, results[1], results) ./ q) :: Vector{T}
end

"""
    cv_exlstsq(x::BEDFile, y, models, q) -> ELSQCrossvalidationResults

Perofrm `q`-fold crossvalidation for the best model size in the exchange algorithm for a `BEDFile` object `x`, response vector `y`, and model sizes specified in `models`.
"""
function cv_exlstsq(
    T        :: Type,
    xfile    :: String,
    xtfile   :: String,
    x2file   :: String,
    yfile    :: String,
    meanfile :: String,
    precfile :: String,
    kernfile :: String;
    q        :: Int = max(3, min(Sys.CPU_CORES, 5)),
    models   :: DenseVector{Int} = begin
           # find p from the corresponding BIM file, then make path
            bimfile = xfile[1:(endof(xfile)-3)] * "bim"
            p       = countlines(bimfile)
            collect(1:min(20,p))
            end,
    folds    :: DenseVector{Int} = begin
           # find n from the corresponding FAM file, then make folds
            famfile = xfile[1:(endof(xfile)-3)] * "fam"
            n       = countlines(famfile)
            cv_get_folds(n, q)
            end,
    pids     :: DenseVector{Int} = procs(),
    tol      :: Float = convert(T, 1e-6),
    max_iter :: Int   = 100,
    window   :: Int   = 20,
    wg_size  :: Int   = 512,
    quiet    :: Bool  = true,
    header   :: Bool  = false,
)

#    1 <= minimum(models) <= maximum(models) <= size(x,2) || throw(ArgumentError("Model sizes must be positive and cannot exceed number of predictors"))

    mses = pfold(T, xfile, xtfile, x2file, yfile, meanfile, precfile, models, kernfile, folds, q, max_iter=max_iter, quiet=quiet, pids=pids, header=header, window=window)

    # what is the best model size?
    k = convert(Int, floor(mean(models[mses .== minimum(mses)])))

    # print results
    quiet || print_cv_results(mses, models, k)

    # recompute ideal model
    # not worth effort to recompute only one model size with GPU, so use CPU instead
    # also extract names of predictors
    b, bidx, bids = refit_exlstsq(T, xfile, xtfile, x2file, yfile, meanfile, precfile, k, models=models, pids=pids, tol=tol, max_iter=max_iter, quiet=quiet, header=header)

    return ELSQCrossvalidationResults{T}(mses, b, bidx, k, sdata(models), bids)
end

# default type for cv_exlstsq is Float64
cv_exlstsq(xfile::String, xtfile::String, x2file::String, yfile::String, meanfile::String, precfile::String, kernfile::String; q::Int = max(3, min(Sys.CPU_CORES, 5)), models::DenseVector{Int} = begin bimfile = xfile[1:(endof(xfile)-3)] * "bim"; p = countlines(bimfile); collect(1:min(20,p)) end, folds::DenseVector{Int} = begin famfile = xfile[1:(endof(xfile)-3)] * "fam"; n = countlines(famfile); cv_get_folds(n, q) end, pids::DenseVector{Int} = procs(), tol::Float64 = 1e-6, max_iter::Int  = 100, window::Int  = 20, quiet::Bool = true, header::Bool = false,
) = cv_exlstsq(Float64, xfile, xtfile, x2file, yfile, meanfile, precfile, kernfile, folds=folds, q=q, models=models, pids=pids, tol=tol, max_iter=max_iter, window=window, quiet=quiet, header=header)

function pfold(
    T          :: Type,
    xfile      :: String,
    x2file     :: String,
    yfile      :: String,
    models     :: DenseVector{Int},
    kernfile   :: String,
    folds      :: DenseVector{Int},
    q          :: Int;
    pids       :: DenseVector{Int} = procs(),
    max_iter   :: Int  = 100,
    quiet      :: Bool = true,
    header     :: Bool = false
)

    # ensure correct type
    T <: Float || throw(ArgumentError("Argument T must be either Float32 or Float64"))

    # how many CPU processes can pfold use?
    np = length(pids)

    # report on CPU processes
    quiet || println("pfold: np = ", np)
    quiet || println("pids = ", pids)

    # set up function to share state (indices of folds)
    i = 1
    nextidx() = (idx=i; i+=1; idx)

    # preallocate cell array for results
    results = cell(q)

    # master process will distribute tasks to workers
    # master synchronizes results at end before returning
    @sync begin

        # loop over all workers
        for worker in pids

            # exclude process that launched pfold, unless only one process is available
            if worker != myid() || np == 1

                # asynchronously distribute tasks
                @async begin
                    while true

                        # grab next fold
                        current_fold = nextidx()

                        # if current fold exceeds total number of folds then exit loop
                        current_fold > q && break

                        # report distribution of fold to worker
                        quiet || print_with_color(:blue, "Computing fold $current_fold on worker $worker.\n\n")

                        # launch job on worker
                        # worker loads data from file paths and then computes the errors in one fold
                        results[current_fold] = remotecall_fetch(worker) do
                                pids = [worker]
                                x = BEDFile(T, xfile, x2file, pids=pids, header=header)
                                y = SharedArray(abspath(yfile), T, (x.geno.n,), pids=pids)
                                one_fold(x, y, models, kernfile, folds, current_fold, max_iter=max_iter, window=window, quiet=quiet, pids=pids)
                        end # end remotecall_fetch()
                    end # end while
                end # end @async
            end # end if
        end # end for
    end # end @sync

    # return reduction (row-wise sum) over results
    return (reduce(+, results[1], results) ./ q) :: Vector{T}
end

function cv_exlstsq(
    T        :: Type,
    xfile    :: String,
    x2file   :: String,
    yfile    :: String,
    kernfile :: String;
    q        :: Int = max(3, min(Sys.CPU_CORES, 5)),
    models   :: DenseVector{Int} = begin
           # find p from the corresponding BIM file, then make path
            bimfile = xfile[1:(endof(xfile)-3)] * "bim"
            p       = countlines(bimfile)
            collect(1:min(20,p))
            end,
    folds    :: DenseVector{Int} = begin
           # find n from the corresponding FAM file, then make folds
            famfile = xfile[1:(endof(xfile)-3)] * "fam"
            n       = countlines(famfile)
            cv_get_folds(n, q)
            end,
    pids     :: DenseVector{Int} = procs(),
    tol      :: Float = convert(T, 1e-4),
    max_iter :: Int   = 100,
    wg_size  :: Int   = 512,
    quiet    :: Bool  = true,
    header   :: Bool  = false
)
    # enforce type
    T <: Float || throw(ArgumentError("Argument T must be either Float32 or Float64"))

    # how many elements are in the path?
    num_models = length(models)

    # want to compute a path for each fold
    # the folds are computed asynchronously
    # only use the worker processes
    mses = pfold(T, xfile, x2file, yfile, models, kernfile, folds, q, max_iter=max_iter, quiet=quiet, pids=pids, header=header, window=window)

    # what is the best model size?
    k = convert(Int, floor(mean(models[mses .== minimum(mses)])))

    # print results
    !quiet && print_cv_results(mses, models, k)

    # recompute ideal model
    # not worth effort to recompute only one model size with GPU, so use CPU instead
    # also extract predictor names
    b, bidx, bids = refit_exlstsq(T, xfile, x2file, yfile, k, models=models, pids=pids, tol=tol, max_iter=max_iter, quiet=quiet, header=header)

    return ELSQCrossvalidationResults{T}(mses, sdata(models), b, bidx, k)
end

# default type for cv_iht is Float64
cv_exlstsq(xfile::String, x2file::String, yfile::String, kernfile::String; q::Int = max(3, min(Sys.CPU_CORES, 5)), models::DenseVector{Int} = begin bimfile=xfile[1:(endof(xfile)-3)] * "bim"; p=countlines(bimfile); collect(1:min(20,p)) end, folds::DenseVector{Int} = begin famfile=xfile[1:(endof(xfile)-3)] * "fam"; n=countlines(famfile); cv_get_folds(n, q) end, pids::DenseVector{Int}=procs(), tol::Float64=1e-4, max_iter::Int=100, quiet::Bool=true, header::Bool=false) = cv_iht(Float64, xfile, x2file, yfile, kernfile, models=models, folds=folds, q=q, pids=pids, tol=tol, max_iter=max_iter, quiet=quiet, header=header)


text:
using IBMQClient, Yao, YaoBlocksQobj
using Test

@testset "YaoBlocksQobj.jl" begin
    qc = chain(
        3,
        put(1 => X),
        put(2 => Y),
        put(3 => Z),
        put(2 => T),
        put(2 => I2),
        swap(1, 2),
        put(3 => Ry(0.7)),
        control(2, 1 => Y),
        control(3, 2 => Z),
    )
    qc1 = chain(
        3,
        put(1 => H),
        put(3 => Rx(0.7)),
        control(2, 1 => Y),
        control(3, 2 => Z),
        put(1 => Yao.Measure(1)),
    )

    header = Dict("description" => "test circuits")
    exp_header = [Dict("description" => "1"), Dict("description" => "2")]
    circuits = [qc, qc1]

    q = convert_to_qobj(circuits, id = "test", header = header, exp_header = exp_header)
    q1 = convert_to_qobj(circuits)
    q2 = convert_to_qobj(circuits, exp_header = exp_header)

    set = [q, q1, q2]
    for i in set
        experiments = i.experiments
        for i = 1:length(experiments)
            exp = experiments[i]
            @test exp isa YaoBlocksQobj.Schema.Experiment
            @test exp.header == Dict("description" => "$i") || exp.header === nothing
            @test exp.config === nothing
            inst = exp.instructions
            ir = convert_to_qbir(inst)
            @test ir == circuits[i]
        end

        qc_inst = chain(
            1,
            put(1 => YaoBlocksQobj.U1{Float64}(2)),
            put(1 => YaoBlocksQobj.U2{Float64}(1, 0.7)),
            put(1 => YaoBlocksQobj.U3{Float64}(0, 1, 0.7)),
        )

        inst = [
            IBMQClient.Schema.Gate("u1", [0], [2], nothing, nothing),
            IBMQClient.Schema.Gate("u2", [0], [1, 0.7], nothing, nothing),
            IBMQClient.Schema.Gate("u3", [0], [0, 1, 0.7], nothing, nothing),
        ]

        @test qc_inst == convert_to_qbir(inst)
    end

    @testset "qbir misc." begin
        @test YaoBlocksQobj.mat(Number, YaoBlocksQobj.U1(0.7)) ≈ Number[
            1 0
            0 exp(im * 0.7)
        ]

        @test YaoBlocksQobj.mat(Number, YaoBlocksQobj.U2(0.7, 0.5)) ≈ Number[
            1/√2 (-exp(im * 0.5))/√2
            (exp(im * 0.7))/√2 (exp(im * (1.2)))/√2
        ]

        @test YaoBlocksQobj.mat(Number, YaoBlocksQobj.U3(0.7, 0.5, 0.2)) ≈ Number[
            cos(0.7 / 2) -sin(0.7 / 2)*exp(im * 0.2)
            sin(0.7 / 2)*exp(im * 0.5) cos(0.7 / 2)*exp(im * (0.7))
        ]

        @test YaoBlocksQobj.iparams_eltype(YaoBlocksQobj.U1(0.7)) == Float64
        @test YaoBlocksQobj.iparams_eltype(YaoBlocksQobj.U2(0.7, 0.5)) == Float64
        @test YaoBlocksQobj.iparams_eltype(YaoBlocksQobj.U3(0.7, 0.5, 0.2)) == Float64

        @test YaoBlocksQobj.getiparams(YaoBlocksQobj.U1(0.7)) == (0.7)
        @test YaoBlocksQobj.getiparams(YaoBlocksQobj.U2(0.7, 0.5)) == (0.7, 0.5)
        @test YaoBlocksQobj.getiparams(YaoBlocksQobj.U3(0.7, 0.5, 0.2)) == (0.7, 0.5, 0.2)
    end
end