text:
@compat for (func, elty) in [(:CLBlastStrmv, Float32), (:CLBlastDtrmv, Float64),
                     (:CLBlastCtrmv, ComplexF32), (:CLBlastZtrmv, ComplexF64)]
    #TODO: (:CLBlastHtrmv, Float16)

    @eval function $func(layout::CLBlastLayout, triangle::CLBlastTriangle,
                         a_transpose::CLBlastTranspose, diagonal::CLBlastDiagonal,
                         n::Integer,
                         a_buffer::cl.CL_mem, a_offset::Integer, a_ld::Integer,
                         x_buffer::cl.CL_mem, x_offset::Integer, x_inc::Integer,
                         queue::cl.CmdQueue, event::cl.Event)
        err = ccall(
            ($(string(func)), libCLBlast),
            cl.CL_int,
            (Cint, Cint, Cint, Cint, Csize_t, Ptr{Cvoid}, Csize_t, Csize_t,
              Ptr{Cvoid}, Csize_t, Csize_t, Ptr{Cvoid}, Ptr{Cvoid}),
            Cint(layout), Cint(triangle), Cint(a_transpose), Cint(diagonal), n,
              a_buffer, a_offset, a_ld, x_buffer, x_offset, x_inc, Ref(queue), Ref(event)
        )
        if err != cl.CL_SUCCESS
            println(stderr, "Calling function $(string($func)) failed!")
            throw(CLBlastError(err))
        end
        return err
    end

    @eval function trmv!(uplo::Char, trans::Char, diag::Char,
                         A::cl.CLArray{$elty,2}, x::cl.CLArray{$elty};
                         queue::cl.CmdQueue=cl.queue(x))
        # check and convert arguments
        m, n = size(A)
        if m != n
            throw(DimensionMismatch("`A` has dimensions $(size(A)) but must be square."))
        end
        if length(x) != n
            throw(DimensionMismatch("`x` has length $(length(x)) while $n is required."))
        end
        if uplo == 'U'
            triangle = CLBlastTriangleUpper
        elseif uplo == 'L'
            triangle = CLBlastTriangleLower
        else
            throw(ArgumentError("Upper/lower marker `uplo` is $(uplo) but only 'U' and 'L' are allowed."))
        end
        if trans == 'N'
            a_transpose = CLBlastTransposeNo
        elseif trans == 'T' || (trans == 'C' && $elty <: Real)
            a_transpose = CLBlastTransposeYes
        elseif trans == 'C' && $elty <: Complex
            a_transpose = CLBlastTransposeConjugate
        else
            throw(ArgumentError("Transpose marker `trans` is $(uplo) but only 'N', 'T', and 'C' are allowed."))
        end
        if diag == 'N'
            diagonal = CLBlastDiagonalNonUnit
        elseif diag == 'U'
            diagonal = CLBlastDiagonalUnit
        else
            throw(ArgumentError("Diagonal marker `diag` is $(diag) but only 'N' and 'U' are allowed."))
        end
        layout = CLBlastLayoutColMajor

        # output event
        event::cl.Event = cl.Event(C_NULL)

        $func(layout, triangle, a_transpose, diagonal,
              n,
              pointer(A), 0, size(A,1),
              pointer(x), 0, 1,
              queue, event)

        # wait for kernel
        cl.wait(event)

        x
    end

end


text:
export getdevice, setdevice, sync

function getdevice()
  p = CUdevice[0]
  cuCtxGetDevice(p)
  Int(p[1])
end

setdevice(dev::Int) = cuCtxSetCurrent(contexts[dev+1])

sync() = cuCtxSynchronize()

ext: jl


text:
using HypothesisTests, Distributions, Base.Test
using HypothesisTests: default_tail

# One sample test
n = 1000
srand(1984948)

x = rand(Normal(), n)
t = OneSampleADTest(x, Normal())
@test t.A² ≈ 0.2013 atol=0.1^4
@test pvalue(t) ≈ 0.99 atol=0.1^4
@test default_tail(t) == :right

x = rand(DoubleExponential(), n)
t = OneSampleADTest(x, Normal())
@test t.A² ≈ 10.7439 atol=0.1^4
@test pvalue(t) ≈ 0.0 atol=0.1^4

x = rand(Cauchy(), n)
t = OneSampleADTest(x, Normal())
@test pvalue(t) ≈ 0.0 atol=0.1^4

x = rand(LogNormal(), n)
t = OneSampleADTest(x, Normal())
@test pvalue(t) ≈ 0.0 atol=0.1^4

x = rand(Uniform(), n)
t = OneSampleADTest(x, Uniform())
@test pvalue(t) ≈ 0.145 atol=0.1^4

x = rand(Uniform(0, 1.8), n)
t = OneSampleADTest(x, Uniform())
@test pvalue(t) ≈ 0.5112 atol=0.1^4

x = rand(Exponential(), n)
t = OneSampleADTest(x, Exponential())
@test pvalue(t) ≈ 0.9427 atol=0.1^4

# k-sample test
srand(948574875)
samples = Any[
    [38.7, 41.5, 43.8, 44.5, 45.5, 46.0, 47.7, 58.0],
    [39.2, 39.3, 39.7, 41.4, 41.8, 42.9, 43.3, 45.8],
    [34.0, 35.0, 39.0, 40.0, 43.0, 43.0, 44.0, 45.0],
    [34.0, 34.8, 34.8, 35.4, 37.2, 37.8, 41.2, 42.8]
]

t = KSampleADTest(samples...)
@test t.A²k ≈ 8.3926 atol=0.1^4
@test t.σ ≈ 1.2038 atol=0.1^4
@test pvalue(t) ≈ 0.0022 atol=0.1^4
@test default_tail(t) == :right

ts = KSampleADTest(samples..., nsim = 20000);
@test pvalue(ts) ≈ 0.00155 atol=0.1^3

t = KSampleADTest(samples..., modified = false)
@test t.A²k ≈ 8.3559 atol=0.1^4
@test t.σ ≈ 1.2038 atol=0.1^4
@test pvalue(t) ≈ 0.0023 atol=0.1^4

ts = KSampleADTest(samples..., modified = false, nsim = 20000);
@test pvalue(ts) ≈ 0.00150 atol=0.1^3

# more tests
srand(31412455)
samples = Any[rand(Normal(), 50), rand(Normal(0.5), 30)]
t = KSampleADTest(samples...)
@test pvalue(t) < 0.05

samples = Any[rand(Normal(), 50), rand(Normal(), 30), rand(Normal(), 20)]
t = KSampleADTest(samples...)
@test pvalue(t) > 0.05

@test pvalue(OneSampleADTest(vcat(rand(Normal(),500), rand(Beta(2,2),500)), Beta(2,2))) ≈ 0.0 atol=0.1^4


text:
## Various helper functions for working with latitude and longitude
const Location2D = SVector{2, Float64}
const LatLonCoords = NamedTuple{(:lat, :lon), Tuple{Float64,Float64}}
LatLonCoords() = (lat = 0.0, lon = 0.0)

"""
Return the vector form of a latitude-longitude point [lat, lon]
"""
function convert_to_vector(c::LatLonCoords)
    return Location2D(c.lat, c.lon)
end

struct EuclideanLatLongMetric <: Metric
end

function Distances.evaluate(::EuclideanLatLongMetric,
                            coords1::Location2D,
                            coords2::Location2D)
    deglen = 110.25
    x = coords1[1]- coords2[1]
    y = (coords1[2] - coords2[2])*cos(coords2[1])
    return deglen*sqrt(x^2 + y^2)
end

@with_kw struct CityParams
    lat_start::Float64
    lat_end::Float64
    lon_start::Float64
    lon_end::Float64
end

function parse_city_params(param_file::String)

    params_dict = TOML.parsefile(param_file)

    return CityParams(lat_start = params_dict["LATSTART"],
                      lat_end = params_dict["LATEND"],
                      lon_start = params_dict["LONSTART"],
                      lon_end = params_dict["LONEND"])
end

@with_kw struct Drone
    avg_speed::Float64      = 0.00777
    depot_number::Int64
    depot_loc::LatLonCoords
end

@with_kw mutable struct DroneProperties
    at_depot::Bool                                  = true
    current_package::String                         = ""
    interaction_events::Vector{InteractionEvent}    = Vector{InteractionEvent}(undef, 0)
    tree::SearchTree                                = SearchTree()
end

const RoutingAllocation = GenericAllocation{Drone,DroneProperties}

struct Package
    delivery::LatLonCoords
    time_window::Tuple{Float64,Float64}
end

"""
Bookkeeping data structure only useful for plotting.
"""
@with_kw struct CurrDroneSiteLocs
    curr_drone_locs_cols::Vector{Tuple{LatLonCoords,Symbol}}    = Tuple{LatLonCoords,Symbol}[]
    curr_sites_locs_cols::Vector{Tuple{LatLonCoords,Symbol}}    = Tuple{LatLonCoords,Symbol}[]
end

@with_kw mutable struct RoutingSimulator <: MRTAEnvironment
    current_time::Float64
    city_params::CityParams
    new_request_prob::Float64
    time_window_duration::Float64
    halton_nn_tree::BallTree
    estimate_matrix::Matrix{Float64}
    delivery_reward::Float64
    active_packages::Dict{String,Package}   = Dict{String,Package}()
    busy_packages::Dict{String,Package}     = Dict{String,Package}()
    done_packages::Dict{String,Package}     = Dict{String,Package}()
    num_total_packages::Int64               = 0
    num_active_packages::Int64              = 0
    late_packages::Int64                    = 0
    delivered_packages::Int64               = 0
    sum_of_delivery_time::Float64           = 0.0
    time_scale::Float64                     = 60.0 # Dealing in 1-minute chunks
    true_delivery_return::Dict{Tuple{String,String},Tuple{Float64,Float64}} = Dict{Tuple{String,String},Tuple{Float64,Float64}}()
    tt_est_std_scale::Float64               = 3.0
    distance_thresh::Float64                = 5.0
    curr_drone_site_locs::CurrDroneSiteLocs = CurrDroneSiteLocs()
end


text:
path = dirname(@__FILE__)

# Initialize arguments to function
h5 = h5open("$path/reference/kalman_filter_args.h5")
for arg in ["data", "TTT", "RRR", "CCC", "QQ", "ZZ", "DD", "EE", "z0", "P0"]
    eval(parse("$arg = read(h5, \"$arg\")"))
end
close(h5)

h5 = h5open("$path/reference/kalman_filter_out.h5")
for arg in ["pred", "vpred", "filt", "vfilt"]
    eval(parse("$arg = read(h5, \"$arg\")"))
end
close(h5)

# Run smoothers
states = Dict{Symbol, Matrix{Float64}}()
shocks = Dict{Symbol, Matrix{Float64}}()

states[:hamilton], shocks[:hamilton] = hamilton_smoother(data, TTT, RRR, CCC, QQ, ZZ, DD, EE, z0, P0)
states[:koopman], shocks[:koopman] = koopman_smoother(data, TTT, RRR, CCC, QQ, ZZ, DD, EE, z0, P0, pred, vpred)
states[:carter_kohn], shocks[:carter_kohn] = carter_kohn_smoother(data, TTT, RRR, CCC, QQ, ZZ, DD, EE, z0, P0; draw_states = false)
states[:durbin_koopman], shocks[:durbin_koopman] = durbin_koopman_smoother(data, TTT, RRR, CCC, QQ, ZZ, DD, EE, z0, P0; draw_states = false)

# Check that last-period smoothed states equal last-period filtered states
@testset "Ensure last smoothed state = the last filtered state" begin
    for smoother in [:hamilton, :koopman, :carter_kohn, :durbin_koopman]
        @test filt[:, end] ≈ states[smoother][:, end]
    end
end

# Compare to expected output
exp_states, exp_shocks = h5open("$path/reference/smoothers_out.h5", "r") do file
    read(file, "exp_states"), read(file, "exp_shocks")
end

@testset "Check expected output for all smoothers" begin
    for smoother in [:hamilton, :koopman, :carter_kohn, :durbin_koopman]
        @test exp_states ≈ states[smoother]
        @test exp_shocks ≈ shocks[smoother]
    end
end

# Make sure that simulation smoothers run with `draw_states` on
carter_kohn_smoother(data, TTT, RRR, CCC, QQ, ZZ, DD, EE, z0, P0; draw_states = true)
durbin_koopman_smoother(data, TTT, RRR, CCC, QQ, ZZ, DD, EE, z0, P0; draw_states = true)

nothing


text:
include("../../../src/includeall.jl")

m = COTHANK()
Γ0, Γ1, C, Ψ, Π = eqcond(m)
#=
rows_zero_Γ0 = []
for i in 1:size(Γ0, 1)
    if all(Γ0[i, :] .== 0.)
        push!(rows_zero_Γ0, i)
    end
end

cols_zero_Γ0 = []
for i in 1:size(Γ0, 1)
    if all(Γ0[:, i] .== 0.)
        push!(rows_zero_Γ0, i)
    end
end
=#

eqcondmats = matopen("eqcond.mat")
#gensysmats = matopen("gensys_out.mat")

Γ0m = read(eqcondmats, "GAM0")
Γ1m = read(eqcondmats, "GAM1")
Cm  = read(eqcondmats, "C")
Ψm  = read(eqcondmats, "PSI")
Πm  = read(eqcondmats, "PPI")
ssmat = read(eqcondmats,"ssvec")
#TTTm = read(gensysmats, "G1")
#RRRm = read(gensysmats, "impact")


text:
################################################################################
#
#  Enumeration of subfields
#
################################################################################

# This is a simple implementation of the algorithms described in
# "Generating subfields" by van Hoeij, Kluners and Novocin
# (see https://doi.org/10.1016/j.jsc.2012.05.010).

#################################################################################
#
#  Principal subfields
#
################################################################################

# Compute the bases of the principal subfields
function _principal_subfields_basis(K::SimpleNumField)
  f = K.pol
  Kx, x = PolynomialRing(K, "x", cached = false)
  n = degree(K)
  #f in Kx
  #fk = Kx([coeff(f,i) for i in 0:n])
  #determine roots
  rts = roots(f, K)
  ar_lin_fac = elem_type(Kx)[x - root for root in rts]
  fK = change_base_ring(K, f, parent = Kx)
  ##divide by roots
  for lin_fac in ar_lin_fac
    fK = div(fK, lin_fac)#divaxa
  end

  fac = factor(fK)

  for (g, e) in fac
    push!(ar_lin_fac, g)
  end

  factor_ar = ar_lin_fac
  k = base_field(K)
  principal_subfields_ar = dense_matrix_type(elem_type(k))[]

  #compute kernel of (phi - id)
  for fi in factor_ar
    M = zero_matrix(k, n, n * degree(fi))
    im_ar = elem_type(Kx)[(mod(x^l,fi)-gen(K)^l) for l in 0:n-1]
    for j in 1:n
      for l in 0:degree(fi)-1
        t = basis_matrix([coeff(im_ar[j], l)])
        for m in 1:n
          M[j, l * n + m] = t[1, m]
        end
      end
    end
    nu, ker = kernel(M, side = :left)

    # This might be expensive for bigger fields?
    if K isa NumField{fmpq}
      ker_rref = fmpq_mat(lll(saturate(FakeFmpqMat(rref(ker)[2]).num)))
    else
      ker_rref = rref(ker)[2]
    end

    if ker_rref in [ b for b in principal_subfields_ar if nrows(b) == nrows(ker_rref)]
      continue
    else
      push!(principal_subfields_ar, ker_rref)
    end
  end
  return principal_subfields_ar
end

@doc Markdown.doc"""
    principal_subfields(L::SimpleNumField) -> Vector{Tuple{NumField, Map}}

Return the principal subfields of $L$ as pairs consisting of a subfield $k$
and an embedding $k \to L$.
"""
function principal_subfields(K::SimpleNumField)
  ba = _principal_subfields_basis(K)
  elts = Vector{Vector{nf_elem}}(undef, length(ba))
  for i in 1:length(ba)
    baf = FakeFmpqMat(ba[i])
    elts[i] = elem_type(K)[elem_from_mat_row(K, baf.num, j, baf.den) for j in 1:nrows(ba[i])]
  end
  T = typeof(K)
  return Tuple{T, morphism_type(T)}[ subfield(K, elts[i], isbasis = true) for i in 1:length(elts)]
end

# Computes the minpoly of a over M if k(a)=K/M/k
# TODO: Improve this.
function _get_sfpoly(Kx, M)
  K = base_ring(Kx)
  k = base_field(K)
  n = degree(K)
  @assert k === base_ring(M)
  M = transpose(M)
  rk = rank(M)
  if rk == 1
      return Kx([K(coeff(minpoly(gen(K)),i)) for i in 0:degree(minpoly(gen(K)))])
  end
  my = div(n, rk)
  ar_basis = Vector{elem_type(K)}(undef, rk)
  for i in 1:rk
    elem_basis = zero(K)
    for j in 1:n
      elem_basis += M[j,i] * gen(K)^(j-1)
    end
    ar_basis[i] = elem_basis
  end
  null_mat = zero_matrix(k, n, 0)
  for i in 1:length(ar_basis)
    null_mat = hcat(null_mat,representation_matrix(ar_basis[i])'[:,1:my])
  end
  temp = zero_matrix(k,n,1)
  temp[my+1,1] = one(k)
  null_mat = hcat(null_mat,temp)
  mat_poly = nullspace(null_mat)[2]
  ar_coeff = Array{elem_type(K),1}(undef,my)
  for i = 1:my
    indx = 1
    temp_coeff = K(0)
    for j = i:my:n
      temp_coeff += mat_poly[j,1] * ar_basis[indx]
      indx += 1
    end
    ar_coeff[i] = temp_coeff
  end
  #w.l. we cannot assume basis ordered in any way
  push!(ar_coeff,K(mat_poly[n+1,1]))
  if !isone(mat_poly[n+1,1])
      error("LC of minpoly not 1")
  end
  return Kx(ar_coeff)
end

# TODO:
# The algorithm works with subfields represented as vector spaces and it needs
# to compute intersections and test containment
# - Improve this by exploiting that everything should be in rref (use reduce_mod)
# - Maybe also cache the pivots

# Computes the intersection of subfields A,B of K/k, represented as k-VS
# TODO (easy): Get rid of the transpose :)
function _intersect_spaces(A::Hecke.AbstractAlgebra.Generic.MatElem, B::Hecke.AbstractAlgebra.Generic.MatElem)
  A = transpose(A)
  B = transpose(B)
  M = nullspace(hcat(A,-B))[2][1:size(A)[2],:]
  intersect_mat = zero_matrix(base_ring(A), nrows(A), ncols(M))
  for i in 1:ncols(intersect_mat)
    intersect_mat[:,i] = A * M[:,i]
  end
  return transpose(intersect_mat)
end

# Computes the intersection of subfields A = [a1,...,an] of K/k, represented as k-VS
function _intersect_spaces(A::Vector{T}) where T
  if length(A) < 1
    throw(error("Number of spaces must be non-zero"))
  elseif length(A) == 1
    return A[1]
  elseif length(A) == 2
    return _intersect_spaces(A[1],A[2])
  else
    intersection_temp = _intersect_spaces(A[1],A[2])
    for i in 3:length(A)
      intersection_temp = _intersect_spaces(intersection_temp,A[i])
    end
    return intersection_temp
  end
end

# Returns true if A is subspace of B,otherwise false, for A,B k-VS
# TODO (easy): Make this great again :)
function _issubspace(A::MatElem, B::MatElem, proper_subspace::Bool = false)     #or cmpr remark7
  intersectAB = _intersect_spaces(A, B)
  Bol = rank(intersectAB) == rank(A)
  if proper_subspace
    return Bol && (rank(intersectAB) < rank(B))
  else
    return Bol
  end
end

function _generating_subfields(S, len::Int = -1)
  # I don't know (yet) why it is necessary, but we have to do it
  if length(S) == 1
    return S
  end
  ar_2delete = Bool[false for i in 1:length(S)]
  for i in 1:length(S)
    if nrows(S[i]) >= len
      ar_intersection = typeof(S)()
      for j in 1:length(S)
        if _issubspace(S[i],S[j],true)
          push!(ar_intersection, S[j])
        end
      end
      #K always principal_subfield
      if length(ar_intersection) == 0
        ar_2delete[i] = true
      else
        intersection = _intersect_spaces(ar_intersection)
        if _issubspace(intersection, S[i]) && _issubspace(S[i], intersection)
          ar_2delete[i] == true
        end
      end
    else
      ar_2delete[i] = true
    end
  end
  _generating_subfields_ar = S
  indx = length(ar_2delete)
  while indx > 0
    if ar_2delete[indx]
      deleteat!(_generating_subfields_ar, indx)
    end
    indx -= 1
  end
  return _generating_subfields_ar
end

function _all_subfields(K, S::Vector{T}, len::Int = -1) where {T}
    Kx, _  = PolynomialRing(K, "x", cached = false)
    if length(S) == 0
        return S
    end
    #compute associated tuple of K
    e = Int[0 for i in 1:length(S)]
    K_mat = identity_matrix(base_field(K), degree(K))
    sf_ar = elem_type(Kx)[_get_sfpoly(Kx, M) for M in S]
    ListSubfields = T[K_mat]
    _next_subfields(ListSubfields, Kx, S ,K_mat,e,0,sf_ar, len)
    return ListSubfields
end

#Hilfsfunktion for _all_subfields
function _next_subfields(ListSubfields, Kx, S::Vector{T}, L::T, e::Array{Int,1}, s::Int, sf_ar, len::Int) where {T}
  i = s + 1
  while i <= length(S)
    if e[i] == 0
      M = _intersect_spaces(L,S[i])
      #check len
      if rank(M) >= len
        #constr associated tuple of M
        ee = Int[0 for j in 1:length(S)]
        for j in 1:length(S)
          #    if _issubspace(M,S[j])
          if iszero(mod(_get_sfpoly(Kx, M),sf_ar[j]))
            ee[j] = 1
          end
        end
        #check whether to append M
        # This should be checked when creating the e'
        appnd = true
        for j in 1:i-1
          if ee[j] > e[j]
            appnd = false
            break
          end
        end
        if appnd
          push!(ListSubfields,M)
          hi = i
          _next_subfields(ListSubfields, Kx, S,M,ee,hi,sf_ar, len)
        end
      end
    end
    i += 1
  end
  #check length
  if len != -1
    indx = length(ListSubfields)
    while indx > 0
      if rank(ListSubfields[indx]) != len
        deleteat!(ListSubfields,indx)
      end
      indx -= 1
    end
  end
end

################################################################################
#
#  Subfield function
#
################################################################################

@doc Markdown.doc"""
    subfields(L::SimpleNumField) -> Vector{Tuple{NumField, Map}}

Given a simple extension $L/K$, returns all subfields of $L$ containing
$K$ as tuples $(k, \iota)$ consisting of a simple extension $k$
and an embedding $\iota k \to K$.
"""
function subfields(K::SimpleNumField; degree::Int = -1)
  n = Hecke.degree(K) # I want to keep the degree keyword
  k = base_field(K)
  #K = k[x]/f
  # K no generating subfield

  # TODO (medium)
  # I don't know why we have to do this.
  # This needs to be fixed properly
  T = typeof(K)
  if n == 1
    return Tuple{T, morphism_type(T)}[(K, id_hom(K))]
  end

  if degree == n
    return Tuple{T, morphism_type(T)}[(K, id_hom(K))]
  end

  if isprime(n)
    res = Tuple{T, morphism_type(T)}[]
    if degree == n
      push!(res, (K, id_hom(K)))
    elseif degree == 1
      kt, t = PolynomialRing(k, "t", cached = false)
      k_as_field = number_field(t-1, check = false, cached = false)[1]
      push!(res, (k_as_field, hom(k_as_field, K, one(K))))
    elseif degree == -1
      kt, t = PolynomialRing(k, "t", cached = false)
      k_as_field = number_field(t-1, check = false, cached = false)[1]
      push!(res, (K, id_hom(K)))
      push!(res, (k_as_field, hom(k_as_field, K, one(K))))
    end
    return res
  end
  princ_subfields = _principal_subfields_basis(K)
  gg = _generating_subfields(princ_subfields)
  sf_asmat_ar = _all_subfields(K, gg, degree)
  #compute embedding
  Res = Vector{Tuple{typeof(K), morphism_type(K)}}()
  #get minimal polynomial of primitive elem k(pe) = M, over k
  for sf_mat in sf_asmat_ar
    #interpret column vectors as field elems
    #sf_mat = transpose(sf_mat)
    #sf_mat_f = FakeFmpqMat(sf_mat)#
    sf_mat_f = sf_mat
    basis_ar = Array{elem_type(K),1}(undef, nrows(sf_mat_f))
    for i in 1:nrows(sf_mat_f)
      if K isa NumField{fmpq}
        _t = FakeFmpqMat(sf_mat_f)
        basis_ar[i] = elem_from_mat_row(K, _t.num, i, _t.den)
      else
        basis_ar[i] = elem_from_mat_row(K, sf_mat_f, i)
      end

    end
    push!(Res, subfield(K, basis_ar, isbasis = true))
  end
  return Res
end

# TODO: Write a dedicated function for the normal case and use the subgroup functions
function subfields_normal(K::SimpleNumField, classes::Bool = false)
  G, mG = automorphism_group(K)
  subs = subgroups(G, conjugacy_classes = classes)
  res = Vector{Tuple{typeof(K), morphism_type(typeof(K))}}()
  for (i, (H, mH)) in enumerate(subs)
    auts = morphism_type(typeof(K))[ mG(mH(h)) for h in H ]
    push!(res, fixed_field(K, auts))
  end
  return res
end


text:
module MaxHelpingHandCustomizableCrumblePlatform

using ..Ahorn, Maple

@mapdef Entity "MaxHelpingHand/CustomizableCrumblePlatform" CustomizableCrumblePlatform(x::Integer, y::Integer, width::Integer=Maple.defaultBlockWidth,
    texture::String="default", oneUse::Bool=false, respawnDelay::Number=2.0, grouped::Bool=false, minCrumbleDurationOnTop::Number=0.2,
    maxCrumbleDurationOnTop::Number=0.6, crumbleDurationOnSide::Number=1.0, outlineTexture::String="objects/crumbleBlock/outline", onlyEmitSoundForPlayer::Bool=false)

const placements = Ahorn.PlacementDict(
    "Crumble Blocks ($(uppercasefirst(texture)), Customizable)\n(max480's Helping Hand)" => Ahorn.EntityPlacement(
        CustomizableCrumblePlatform,
        "rectangle",
        Dict{String, Any}(
            "texture" => texture
        )
    ) for texture in Maple.crumble_block_textures
)

Ahorn.editingOptions(entity::CustomizableCrumblePlatform) = Dict{String, Any}(
    "texture" => Maple.crumble_block_textures
)

Ahorn.editingOrder(entity::CustomizableCrumblePlatform) = String["x", "y", "width", "crumbleDurationOnSide", "minCrumbleDurationOnTop", "maxCrumbleDurationOnTop"]

Ahorn.minimumSize(entity::CustomizableCrumblePlatform) = 8, 0
Ahorn.resizable(entity::CustomizableCrumblePlatform) = true, false

function Ahorn.selection(entity::CustomizableCrumblePlatform)
    x, y = Ahorn.position(entity)
    width = Int(get(entity.data, "width", 8))

    return Ahorn.Rectangle(x, y, width, 8)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::CustomizableCrumblePlatform, room::Maple.Room)
    texture = get(entity.data, "texture", "default")
    texture = "objects/crumbleBlock/$texture"

    # Values need to be system specific integer
    x = Int(get(entity.data, "x", 0))
    y = Int(get(entity.data, "y", 0))

    width = Int(get(entity.data, "width", 8))
    tilesWidth = div(width, 8)

    Ahorn.Cairo.save(ctx)

    Ahorn.rectangle(ctx, 0, 0, width, 8)
    Ahorn.clip(ctx)

    for i in 0:ceil(Int, tilesWidth / 4)
        Ahorn.drawImage(ctx, texture, 32 * i, 0, 0, 0, 32, 8)
    end

    Ahorn.restore(ctx)
end

end

max_stars_repo_name: swoolcock/MaxHelpingHand

max_issues_repo_name: swoolcock/MaxHelpingHand

max_forks_repo_name: swoolcock/MaxHelpingHand


text:
# Chris Rackaukas' macro generator
"""
Chris Rackaukas' macro generator

e.g.
```julia
@def mymac begin
   a = 2
end

@mymac
```

Used e.g. to define long lists of imports that are used in several
submodules.
"""
macro def(name, definition)
   return quote
      macro $(esc(name))()
         esc($(Expr(:quote, definition)))
      end
   end
end


text:
export  addorderby!,
        clearlimit!,
        setlimit!,
        clearoffset!,
        setoffset!

struct CompoundQuery <: SQLQuery
    larg::String
    rarg::String
    cop::Symbol
    isall::Bool
    orderby::Vector{String}
    asc::Vector{Bool}
    nulllast::Vector{Bool}
    limit::Int
    offset::Int

    CompoundQuery( larg::AbstractString, rarg::AbstractString, cop::Symbol, isall::Bool ) = new( larg, rarg, cop, isall, Var[], Bool[], Bool[], -1, 0 )
end

function Base.show( io::IO, cq::CompoundQuery )
    print( io, "(", cq.larg, ") ", uppercase(cq.cop |> string), " " )
    cq.isall && print( io, "ALL " )
    print( io, "(", cq.rarg, ")" )

    if !isempty(cq.orderby)
        print( io, "\n    ORDER BY " )
        ostr = map( eachindex(cq.orderby) ) do ii
            print( io, cq.orderby[ii], ' ', cq.asc[ii] ? "ASC" : "DESC", " NULLS ", cq.nulllast[ii] ? "LAST" : "FIRST" )
        end
        print( io, join( ostr, ", " ) )
    end

    if (cq.limit >= 0) || (cq.offset > 0)
        print( io, "\n   " )
        cq.limit < 0 || print( io, " LIMIT ", cq.limit )
        cq.offset > 0 && print( io, " OFFSET ", cq.offset )
    end
end

for cop in [:union, :intersect, :except]
    Core.eval( @__MODULE__, "export query$cop" |> Meta.parse )
    Core.eval( @__MODULE__, """\"\"\"
    ```
    query$cop(
        larg::Union{Select, CompoundQuery},
        rarg::Union{Select, CompoundQuery},
        isall::Bool=false )
    ```
    This function creates a `CompoundQuery` object consisting of the SQL operator `$(lowercase( cop |> string ))` applied to `larg` and `rarg` as left and right arguments respectively. The argument `isall` determines whether the flag `ALL` will be set, that is, whether duplicate rows will be kept or not.

    Note that this function does not check of the results of the two queries are compatible.
    \"\"\"
    query$cop( larg::Union{Select, CompoundQuery}, rarg::Union{Select, CompoundQuery}, isall::Bool=false ) = CompoundQuery( larg |> string, rarg |> string, :$cop, isall )""" |> Meta.parse )
end

"""
```
addorderby!(
    column::SVar;
    asc::Bool=true,
    nullorder::Symbol=:default )
```
This function permits the user to add a command that adds the column `column` to the `ORDER BY` clause of a `Select` or `CompoundQuery` struct. The argument `asc` controls the sorting order of the column, ascending or descending, and the argument `nullorder` controls where to place `NULL` entries in the column. This argument can take only the values `:first`, `:last`, or `:default`.

The output of this function is a function that takes an argument of type `Select` or `CompoundQuery` and returns the updated argument, such that Julia's function chaining syntax can be used.

Example:
```
sq |> addorderby!( "var2" ) |> addorderby!( "var1", asc=false )
```
sets the command to order the result of `sq` first by `var2` in ascending order, then by `var2` in descending order.
"""
function addorderby!( column::SVar; asc::Bool=true, nullorder::Symbol=:default )
    sq::SQLQuery -> begin
        if nullorder ∉ [:first, :last, :default]
            @warn "Unknown nullorder, must be one of :first, :last, or :default"
            return sq
        end

        push!( sq.orderby, Var(column) )
        push!( sq.asc, asc )
        push!( sq.nulllast, nullorder === :default ? asc : nullorder === :first )
        sq
    end
end

"""
```
clearlimit!( sq::Union{Select, CompoundQuery} )
```
This function clears the `LIMIT` clause of the query `sq` and returns the updated query.
"""
function clearlimit!( sq::SQLQuery )
    sq.limit = -1
    sq
end

"""
```
setlimit!( limit::Integer )
```
This function permits the user to add a command that sets the `LIMIT` clause of a `Select` or `CompoundQuery` struct to `limit`.

The output of this function is a function that takes an argument of type `Select` or `CompoundQuery` and returns the updated argument, such that Julia's function chaining syntax can be used.

Example:
```
sq |> setlimit!( 100 )
```
sets the `LIMIT` clause of `sq` to 100.
"""
setlimit!( limit::Integer ) = sq::SQLQuery -> begin
    limit < 0 && return clearlimit!(sq)
    sq.limit = limit |> Int
    sq
end

"""
```
clearoffset!( sq::Union{Select, CompoundQuery} )
```
This function clears the `OFFSET` clause of the query `sq` and returns the updated query.
"""
function clearoffset!( sq::SQLQuery )
    sq.offset = 0
    sq
end

"""
```
setoffset!( offset::Integer )
```
This function permits the user to add a command that sets the `OFFSET` clause of a `Select` or `CompoundQuery` struct to `offset`.

The output of this function is a function that takes an argument of type `Select` or `CompoundQuery` and returns the updated argument, such that Julia's function chaining syntax can be used.

Example:
```
sq |> setoffset!( 10 )
```
sets the `OFFSET` clause of `sq` to 10.
"""
setoffset!( offset::Integer ) = sq::SQLQuery -> begin
    offset > 0 || return clearoffset!(sq)
    sq.offset = offset |> Int
    sq
end

"""
```
compoundquery( conn::Conn, df::Bool=true )
```
This function permits the user to execute the SQL equivalent of a `CompoundQuery` struct on the database that `conn` connects to. The flag `df` determines the format of the output (`DataFrame` or `NamedTuple`) if the statement will return an output.

The output of this function is a function that takes a `Struct` as argument and returns the result of the executed SQL statement.

Example:
```
sq |> compoundquery(conn)
```
"""
compoundquery( conn::Conn, df::Bool=true ) = cq::CompoundQuery -> begin
    res = execute( conn, string( cq, ";" ) )
    res |> makeunique! |> (df ? DataFrame : columntable)
end


text:
using DataInterpolations, Test

t = [1,2,3]
u = [0, 1, 0]

for di in [LinearInterpolation, QuadraticInterpolation, ConstantInterpolation]
    li = di(copy(u), copy(t))
    append!(li, u, t)
    @test li == di(vcat(u, u), vcat(t, t))

    li = di(copy(u), copy(t))
    push!(li, 1, 4)
    @test li == di(vcat(u, 1), vcat(t, 4))
end


text:
"""
The ColorSchemes module provides simple access to color schemes (arrays of colors).

Use `get(cscheme, n)` with n varying from 0 to 1 to find a color at that point
in the color scheme `cscheme`, where 0 is the leftmost color, and 1 the rightmost.

Access the scheme's colors directly using `cscheme.colors`.

See also `getinverse()`.
"""
module ColorSchemes

import Base.get, Base.reverse

using Colors, ColorTypes, FixedPointNumbers, StaticArrays

export ColorScheme,
       get,
       getinverse,
       colorschemes,
       loadcolorscheme,
       findcolorscheme,
       ColorSchemeCategory

"""
    ColorScheme(colors, category, notes)

Create a ColorScheme from an array of colors. You can also name it, assign a
category to it, and add notes.

```julia
myscheme = ColorScheme([Colors.RGB(0.0, 0.0, 0.0), Colors.RGB(1.0, 1.0, 1.0)],
    "custom", "twotone, black and white")
```

"""
struct ColorScheme{V <: AbstractVector{<:Colorant},S1 <: AbstractString,S2 <: AbstractString}
    colors::V
    category::S1
    notes::S2
end

ColorScheme(colors::V, category::S1="", notes::S2="") where {V,S1,S2} =
    ColorScheme{V,S1,S2}(colors, category, notes)

# displaying swatches
Base.show(io::IO, m::MIME"image/svg+xml", cscheme::ColorScheme) = show(io, m, cscheme.colors)

# for making the catalogue pages in the docs
struct ColorSchemeCategory
    name :: String
end

Base.getindex(cs::ColorScheme, i::AbstractFloat) = get(cs, i)
Base.getindex(cs::ColorScheme, i::AbstractVector{<: AbstractFloat}) = get(cs, i)

"""
    loadcolorscheme(vname, colors, cat="", notes="")

Define a ColorScheme from a symbol and an array of colors,
and add it to the `colorschemes` dictionary. Optionally
speify a category and some notes.

### Example

This is an example (from `ColorsSchemes/data/flags.jl`) of how the supplied colorschemes are
loaded into the ColorSchemes dictionary.

```
loadcolorscheme(:flag_nu, [
        RGB(0.7843137254901961, 0.06274509803921569, 0.1803921568627451),
        RGB(0.00392156862745098, 0.12941176470588237, 0.4117647058823529),
        RGB(1.0, 1.0, 1.0),
        RGB(0.996078431372549, 0.8666666666666667, 0.0),
    ], "flags", "The flag of Niue")
```
"""
function loadcolorscheme(vname, colors, cat="", notes="")
    haskey(colorschemes, vname) && println("$vname overwritten")
    colorschemes[vname] = ColorSchemes.ColorScheme(colors, cat, notes)
    return colorschemes[vname]
end

"""
    colorschemes

An exported dictionary of pre-defined colorschemes:

```
colorschemes[:summer] |> show
   ColorScheme(
      ColorTypes.RGB{Float64}[
          RGB{Float64}(0.0,0.5,0.4), RGB{Float64}(0.01,0.505,0.4), RGB{Float64}(0.02,0.51,0.4), RGB{Float64}(0.03,0.515,0.4),
          ...
```

To choose a random ColorScheme:

```
scheme = rand(keys(colorschemes))
```

"""
const colorschemes = Dict{Symbol,ColorScheme}()

function loadallschemes()
    # load the installed schemes
    datadir = joinpath(dirname(@__DIR__), "data")
    include(joinpath(datadir, "allcolorschemes.jl"))
    include(joinpath(datadir, "colorbrewerschemes.jl"))
    include(joinpath(datadir, "matplotlib.jl"))
    include(joinpath(datadir, "cmocean.jl"))
    include(joinpath(datadir, "sampledcolorschemes.jl"))
    include(joinpath(datadir, "gnu.jl"))
    include(joinpath(datadir, "colorcetdata.jl"))
    include(joinpath(datadir, "scicolor.jl"))
    include(joinpath(datadir, "seaborn.jl"))
    include(joinpath(datadir, "tableau.jl"))
    include(joinpath(datadir, "cvd.jl"))
    include(joinpath(datadir, "flags.jl"))

    # create them as constants...
    for key in keys(colorschemes)
        @eval const $key = colorschemes[$(QuoteNode(key))]
    end
end

loadallschemes()

"""
    findcolorscheme(str)

Find all color schemes matching `str`. `str` is interpreted as a regular expression (case-insensitive).

To read the notes of built-in colorscheme `cscheme`:

```
colorschemes[:cscheme].notes
```
"""
function findcolorscheme(str)
    println("\ncolorschemes containing \"$str\"\n")
    counter = 0
    for (k, v) in colorschemes
        if occursin(Regex(str, "i"), string(k))
            printstyled("$(rpad(k, 20))\n", bold=true)
            counter += 1
        elseif occursin(Regex(str, "i"), string(v.category))
            printstyled("$(rpad(k, 20))", bold=true)
            println(" (category) $(v.category)")
            counter += 1
        end
        if occursin(Regex(str, "i"), string(v.notes))
            printstyled("$(rpad(k, 20))", bold=true)
            l = min(30, length(v.notes))
            println(" (notes) $(v.notes[1:l])...")
            counter += 1
        end
    end
    counter > 0 ?
        println("\n\nfound $counter result$(counter > 1 ? "s" : "") for \"$str\"") :
        println("\n\nnothing found for \"$str\"")
    return nothing
end

# Interfaces
Base.length(cscheme::ColorScheme) = length(cscheme.colors)
Base.getindex(cscheme::ColorScheme, I) = cscheme.colors[I]
Base.size(cscheme::ColorScheme) = length(cscheme.colors)
Base.IndexStyle(::Type{<:ColorScheme}) = IndexLinear()

# iterable
function Base.iterate(cscheme::ColorScheme)
    if length(cscheme.colors) == 0
        return nothing
    end
    return (first(cscheme.colors), 1)
end

function Base.iterate(cscheme::ColorScheme, state)
    if state >= length(cscheme.colors)
        return nothing
    else
        return (cscheme.colors[state + 1], state + 1)
    end
end

Base.lastindex(cscheme::ColorScheme) = lastindex(cscheme.colors)

# utility lerping function to convert a value between oldmin/oldmax to equivalent value between newmin/newmax
remap(value, oldmin, oldmax, newmin, newmax) =
    ((value .- oldmin) ./ (oldmax .- oldmin)) .* (newmax .- newmin) .+ newmin

const AllowedInput = Union{<:Real,AbstractArray{<:Real}}

defaultrange(x::Real) = zero(x), oneunit(x)
defaultrange(x::AbstractArray) = zero(eltype(x)), oneunit(eltype(x))

"""
    get(cscheme::ColorScheme, inData :: Array{Number, 2}, rangescale=:clamp)
    get(cscheme::ColorScheme, inData :: Array{Number, 2}, rangescale=(minVal, maxVal))

Return an RGB array of colors generated by applying the
colorscheme to the 2D input data.

If `rangescale` is `:clamp` the colorscheme is applied to
values between 0.0-1.0, and values outside this range get
clamped to the ends of the colorscheme.

Else, if `rangescale` is `:extrema`, the colorscheme is
applied to the range `minimum(indata)..maximum(indata)`.

TODO: this function expects the colorscheme to consist of
RGB [0.0-1.0] values. It should work with more colortypes.

# Examples

```
img = get(colorschemes[:leonardo], rand(10,10)) # displays in Juno Plots window, but
save("testoutput.png", img)      # you'll need FileIO or similar to do this

img2 = get(colorschemes[:leonardo], 10.0 * rand(10, 10), :extrema)
img3 = get(colorschemes[:leonardo], 10.0 * rand(10, 10), (1.0, 9.0))

# Also works with PerceptualColourMaps
using PerceptualColourMaps # warning, installs PyPlot, PyCall, LaTeXStrings
img4 = get(PerceptualColourMaps.cmap("R1"), rand(10,10))
```
"""
function get(cscheme::ColorScheme, x::AllowedInput, rangemode::Symbol)
    rangescale = if rangemode == :clamp
        defaultrange(x)
    elseif rangemode == :extrema
        extrema(x)
    else
        error("rangescale ($rangemode) not supported, should be :clamp, :extrema or tuple (minVal, maxVal)")
    end
    get(cscheme, x, rangescale)
end
function get(cscheme::ColorScheme, x::AllowedInput, rangescale::NTuple{2,<:Real}=defaultrange(x))
    x isa AbstractRange && (x = collect(x))
    # check for empty range
    if iszero(first(rangescale) - last(rangescale))
        rangescale = zero(first(rangescale)), oneunit(last(rangescale))
    end
    x = clamp.(x, rangescale...)
    before_fp = remap(x, rangescale..., 1, length(cscheme))
    before = round.(Int, before_fp, RoundDown)
    after =  min.(before .+ 1, length(cscheme))
    #  blend between the two colors adjacent to the point
    cpt = before_fp - before
    return weighted_color_mean.(1 .- cpt, cscheme.colors[before], cscheme.colors[after])
end
# Boolean just takes the start and end values of the scheme. Also avoid using a branch.
function get(cscheme::ColorScheme, x::Union{Bool,AbstractArray{Bool}})
    i = x .* (length(cscheme) .- 1) .+ 1
    return getindex.(Ref(cscheme.colors), i)
end

"""
    get(cs::ColorScheme, g::Color{T,1} where T<:Union{Bool, AbstractFloat, FixedPoint})

Return the color in `cs` that corresponds to the gray value `g`.
"""
function get(cs::ColorScheme, g::Color{T,1} where T <: Union{Bool,AbstractFloat,FixedPoint})
    get(cs, ColorTypes.gray(g)) # don't confuse with the 'grays' colorscheme
end

"""
    getinverse(cscheme::ColorScheme, c, range=(0.0, 1.0))

Computes where the provided Color `c` would fit in `cscheme`.

This is the inverse of `get()` — it returns the value `x` in the provided `range`
for which `get(scheme, x)` would most closely match the provided Color `c`.

# Examples

The first example asks: "where in the `leonardo` colorscheme will I find the nearest
color to red?":

```
julia> getinverse(colorschemes[:leonardo], RGB(1, 0, 0))
0.6248997995654847

julia> getinverse(ColorScheme([Colors.RGB(0,0,0), Colors.RGB(1,1,1)]),  Colors.RGB(.5,.5,.5))
0.5432555858022048

julia> cs = ColorScheme(range(Colors.RGB(0,0,0), stop=Colors.RGB(1,1,1), length=5))

julia> getinverse(cs, cs[3])
0.5
```
"""
function getinverse(cscheme::ColorScheme, c, rangescale::Tuple{Number,Number}=(0.0, 1.0))
    # TODO better error handling please
    length(cscheme) <= 1 && throw(ErrorException("ColorScheme of length $(length(cscheme)) is not long enough"))
    cdiffs = map(c_i -> colordiff(promote(c, c_i)...), cscheme.colors)
    closest = argmin(cdiffs)
    left = right = 0
    if closest == 1 ; left = closest; right = closest + 1;
    elseif closest == length(cscheme) ; left = closest - 1; right = closest;
    else
        next_closest = cdiffs[closest - 1] < cdiffs[closest + 1] ? closest - 1 : closest + 1
        left = min(closest, next_closest)
        right = max(closest, next_closest)
    end
    v = left
    if cdiffs[left] != cdiffs[right] # prevents divide by zero
        v += ( cdiffs[left] / (cdiffs[left] + cdiffs[right]))
    end
    return remap(v, 1, length(cscheme), rangescale...)
end

"""
    reverse(cscheme)

Make a new ColorScheme with the same colors as `cscheme` but in reverse order.
"""
Base.reverse(cscheme::ColorScheme) =
    ColorScheme(reverse(cscheme.colors), cscheme.category, cscheme.notes)

end


text:
include("test_pycall.jl")
include("test_zygote.jl")
include("test_chainrules.jl")


text:
#===============================================================================

    mc_peak_analysis.jl

    melting curve peak analysis

===============================================================================#

import DataStructures.OrderedDict
import DataFrames.DataFrame
import StatsBase: rle
import Dierckx: Spline1D, derivative
import Memento: debug, info

#==============================================================================#

"Finds peak in the melting curve data, specifically the negative of the derivative
of fluorescence with respect to temperature, for each well and channel."
function mc_peak_analysis(
    output_type                 ::Type{P},
    i                           ::McInput,
    df                          ::DataFrame,
) where {P <: McPeakOutput}

    ## functions to parse input data >>

    "Filter out data points that are separated by narrow temperature intervals."
    filter_too_close(df ::DataFrame) =
        df[df[:temperature] |> mold(Float_T) |> temperature_intervals |> not_too_close, :]

    "Calculate temperature intervals."
    temperature_intervals(temps_orig ::Vector{<: Float_T}) =
        vcat(diff(temps_orig), Inf_T)

    "Criterion for narrow temperature intervals."
    not_too_close(temp_intervals ::Vector{<: Float_T}) =
        temp_intervals .> i.temp_too_close_frac * median(temp_intervals)

    jitter_temps!(df ::DataFrame) =
        (mutate_dups(df[:temperature], i.jitter_constant), df[:fluorescence] |> mold(Float_T))

    ## functions used to calculate `span_smooth` >>

    "Choose the value for the smoothing parameter `span_smooth`."
    choose_span_smooth() =
        if i.auto_span_smooth
            info(logger, "automatic selection of `span_smooth`...")
            calc_span_smooth(fu_rle())
        else
            info(logger, "no automatic selection, use span_smooth_default " *
                string(i.span_smooth_default) * " as `span_smooth`")
            i.span_smooth_default
        end ## if

    "Calculate the smoothing parameter `span_smooth`."
    @inline function calc_span_smooth(fu ::Tuple{Vector{Bool},Vector{<: Int_T}})

        larger_span(span_smooth_product ::Real) =
            if span_smooth_product > i.span_smooth_default
                info(logger, "`span_smooth` was selected as $span_smooth_product")
                return span_smooth_product
            else
                info(logger, "`span_smooth_product` $span_smooth_product < " *
                    "`span_smooth_default`, use `span_smooth_default` " *
                    string(i.span_smooth_default))
                return i.span_smooth_default
            end ## if

        "Calculate `span_smooth_product` as the longest temperature span
        in which fluorescence rises with temperature."
        span_smooth_product() =
            i.span_smooth_factor * max_fu_temp_span(fu_idc()) / i.temp_span

        max_fu_temp_span(fu ::AbstractVector{Int_T}) =
            maximum(temps[fu_idc[2:2:end]] .- temps[fu[1:2:end]])

        fu_idc() = (cumsum(fu[1][1] ? vcat(0,fu[2]) : fu[2]) .+ 1)[1:2*sum(fu[1])]

        ## << end of function definitions nested within calc_span_smooth()

        if fu[1] == [false]
            info(logger, "no fluorescence increase with temperature increase was detected: " *
                "using `span_smooth_default` $(i.span_smooth_default)")
            return i.span_smooth_default
        else
            info(logger, "fluorescence increase with temperature increase was detected")
            return larger_span(span_smooth_product())
        end ## if
    end ## calc_span_smooth()

    ## `fu_rle` - fluorescence up run length encoding
    "Find the region(s) where there is a positive gradient such that fluorescence
    increases with temperature."
    fu_rle() = calc_fu_rle(temps_in_span(max_fluo_decrease()))

    is_increasing(x ::AbstractVector) =
        diff(x) .> zero(x)

    calc_fu_rle(idc ::AbstractVector{Int_T}) =
        rle(is_increasing(fluos[idc]))

    "Find the region of length 2 * `temperature_bandwidth` showing the steepest
    fluorescence decrease between start and end."
    max_fluo_decrease() =
        indmax(fluo_decrease(temps_in_span(i)) for i in 1:len_raw)

    fluo_decrease(sel_idc_int ::AbstractVector{Int_T}) =
        fluos[sel_idc_int[1]] - fluos[sel_idc_int[end]]

    "Find nearby data points in a vector."
    find_in_span(
        X           ::AbstractVector,
        i           ::Int_T,
        half_span   ::Real
    ) =
        find(X) do x
            X[i] - half_span <= x <= X[i] + half_span
        end

    temps_in_span(i ::Int_T) =
        find_in_span(temps, i, i.temperature_bandwidth)

    ## smoothing functions >>

    "Truncate elements in a vector to the length of the shortest element."
    shorten(x) =
        x |> mold(index(x |> mold(length) |> minimum |> from(1)))

    "Model fluorescence as a cubic spline function of temperature using the Dierckx
    Fortran library. The default parameter `s = 0.0` interpolates without smoothing."
    spline_model(x) =
        Spline1D(shorten(x)...; k = 3)

    "Smooth raw fluorescence values."
    smooth_raw_fluo() =
        (temps, supsmu(temps, fluos, span_smooth / i.denser_factor))

    "Model fluorescence as a cubic spline function of temperature, then recalculates
    the fluorescence and its derivative with respect to temperature using a denser
    sequence of interpolated temperatures."
    smoothed_data_matrix() =
        hcat(
            ## collates processed, interpolated data into matrix
            ## note: memory intensive
            denser_temps,
            fluo_spl_blsub(spl(denser_temps)),
            negderiv_smu(-derivative(spl, denser_temps)))

    "Calculate baseline-subtracted spline-smoothed fluorescence data, and optionally
    smooth the output of the spline function."
    fluo_spl_blsub(fluo_spl ::AbstractVector{<: Float_T}) =
        ## assumes constant baseline == minimum fluorescence value
        sweep(minimum)(-)(
            ## optionally, smooth the output of the spline function
            i.smooth_fluo_spline ?
                supsmu(denser_temps, fluo_spl, span_smooth) :
                fluo_spl)

    "Calculate the negative derivative at interpolated temperatures and smooth
    the output using `supsmu`."
    negderiv_smu(negderiv ::AbstractVector{<: Float_T}) =
        supsmu(denser_temps, negderiv, span_smooth)

    ## note: DataArray format doesn't work for `derivative` by "Dierckx"
    "Create a dense sequence of interpolated temperature values."
    interpolated_temperatures() =
        Array(colon(
            min_temp,
            temp_span / (len_raw * i.denser_factor - 1),
            max_temp))

    ## peak finding functions >>

    "Find the indices in a vector where the value at the index equals the summary
    value of a sliding window centered at the index. This can be used to find local
    summits and nadirs (peaks and troughs) in a 1-dimensional data series."
    function find_local(
        summary_func    ::Function,
        vals            ::AbstractVector,
        half_width      ::Int_T,
    )
        vals_in_window(i ::Int_T) = vals_padded[i : i + half_width * 2]
        match_summary_val(v ::AbstractVector) = summary_func(v) == v[half_width + 1]
        #
        padding = -summary_func(-vals) |> furnish(half_width)
        vals_padded = [padding; vals; padding]
        vals |> length |> from(1) |> collect |>
            mold(match_summary_val ∘ vals_in_window) |> find
    end

    # half_peak_span_temperature = 0.5 * i.peak_span_temperature
    half_peak_window() =
        0.3 * (i.peak_span_temperature / temp_span) * len_denser |>
            roundoff(0) |> Int_T

    "Find the indices of summits and nadirs the smoothed -df/dT data series."
    find_summits_and_nadirs() =
        [maximum, minimum] |> mold(f -> find_local(f, negderiv_smu, half_peak_window()))

    summits_and_nadirs(sn_idc) =
        OrderedDict(
            zip([:summit, :nadir],
                sn_idc |> mold(idc -> smoothed_data[idc, :])))

    ## peak constructors used in find_peaks()
    make_peak(::Void) = nothing
    function make_peak(element ::PeakIndicesElement)
        normalize_range(x ::AbstractArray) =
            sweep(minimum)(-)(x) |> sweep(maximum)(/)

        nd_s=normalize_range(negderiv_smu)[element.summit_idx]
        nd_rn=normalize_range(negderiv_smu)[element.right_nadir_idx]
        nd_ln=normalize_range(negderiv_smu)[element.left_nadir_idx]
        t_s=normalize_range(denser_temps)[element.summit_idx]
        t_rn=normalize_range(denser_temps)[element.right_nadir_idx]
        t_ln=normalize_range(denser_temps)[element.left_nadir_idx]

        equal_dist=minimum([element.summit_idx-element.left_nadir_idx,element.right_nadir_idx-element.summit_idx])
        equal_dist<0?equal_dist=5:nothing
        nd_rn_=normalize_range(negderiv_smu)[element.summit_idx+equal_dist]
        nd_ln_=normalize_range(negderiv_smu)[element.summit_idx-equal_dist]
        t_rn_=normalize_range(denser_temps)[element.summit_idx+equal_dist]
        t_ln_=normalize_range(denser_temps)[element.summit_idx-equal_dist]

        a_=[(nd_s-nd_ln_)/(t_s-t_ln_),(nd_s-nd_rn_)/(t_rn_-t_s)]
        return Peak(
            element.summit_idx,                 ## summit_idx
            denser_temps[element.summit_idx],   ## temperature at peak
            minimum([nd_s-nd_rn,nd_s-nd_ln])>0.1 && maximum(a_)/minimum(a_)<1.9 ? peak_bounds(element) |> calc_area : 0)  ## area
    end
    "Iterate over summits and nadirs to create a vector of `Peak`s."
    function find_peaks(
        summit_idc      ::AbstractVector,
        nadir_idc       ::AbstractVector
    )
        peak_finder =
            PeakIndices(
                negderiv_smu[summit_idc],
                summit_idc,
                nadir_idc)
        Vector{Peak}(peak_finder |> collect |> mold(make_peak) |> sift(thing))
    end ## find_peaks()

    "Find the shoulders of a peak."
    @inline function peak_bounds(element ::PeakIndicesElement)
        left_nadir_idx, summit_idx, right_nadir_idx =
            element |> fieldnames |> mold(curry(getfield)(element))
        # nadir_vec = [left_nadir_idx, right_nadir_idx]
        # low_nadir_idx, high_nadir_idx = map(
        #     func -> nadir_vec[func(negderiv_smu[nadir_vec])],
        #     [indmin, indmax])
        (low_nadir_idx, high_nadir_idx) =
            (negderiv_smu[left_nadir_idx] < negderiv_smu[right_nadir_idx]) ?
                (left_nadir_idx, right_nadir_idx) :
                (right_nadir_idx, left_nadir_idx)
        hn_ns = negderiv_smu[high_nadir_idx]
        #
        ## find the nearest location to `summit_idx`
        ## on the `low_nadir_idx` side of the peak
        ## where the line `negderiv = negderiv_smu[high_nadir_idx]`
        ## is crossed by the `negderiv_smu` curve
        idx_step = -1 + 2 * (summit_idx < low_nadir_idx)
        about2cross_idx = low_nadir_idx - idx_step ## peak slightly narrower than using actual crossing point
        idx = summit_idx
        negderiv_smu1 = negderiv_smu[summit_idx]
        limit_idx = low_nadir_idx - 2 * idx_step
        while (idx != limit_idx)
            negderiv_smu0 = negderiv_smu1
            idx1 = idx + idx_step
            negderiv_smu1 = negderiv_smu[idx1]
            if negderiv_smu1 <= hn_ns <= negderiv_smu0
                about2cross_idx = idx
                break
            end ## if
            idx = idx1
        end ## while
        return (high_nadir_idx, about2cross_idx)
    end ## peak_bounds()

    ## == integrated -df/dT peak area elevated from x-axis minus
    ## trapezium-shaped baseline area elevated from x-axis
    "Calculate the area of the peak above its baseline."
    @inline function calc_area(peak_bounds_idc ::Tuple{Int_T, Int_T})
        #
        area_func(temp_lo ::Real, temp_hi ::Real) =
            -sum(negderiv_smu[[peak_bounds_idc...]]) * 0.5 * (temp_hi - temp_lo) -
                (spl(temp_hi) - spl(temp_lo)) ## 'baseline' of peak
        #
        ordered_tuple(x, y) = (x < y) ? (x, y) : (y, x)
        #
        ## << end of function definitions nested in peak_area()

        # temp_lo, temp_hi = map(
        #     func -> func(denser_temps[peak_bound_idc]),
        #     [minimum, maximum])
        return area_func(ordered_tuple(denser_temps[[peak_bounds_idc...]]...)...)
    end ## calc_area()

    ## count cross points
    @inline function count_cross_points()
        ## vectorized version
        # num_cross_points = sum(map(1:(len_denser-1)) do i
        #     (negderiv_smu[i] - negderiv_midrange) * (negderiv_smu[i+1] - negderiv_midrange) <= 0
        # end)
        #
        ## devectorized version
        num_cross_points = 0
        negderiv_smu_centred1 = sign(negderiv_smu[1] - negderiv_midrange)
        local i = 0
        while (i < len_denser)
            i += 1
            negderiv_smu_centred0 = negderiv_smu_centred1
            negderiv_smu_centred1 = sign(negderiv_smu[i] - negderiv_midrange)
            if negderiv_smu_centred0 != negderiv_smu_centred1
                num_cross_points += 1
            end
        end ## while
        return num_cross_points
    end ## count_cross_points()

    @inline function count_cross_points_dist()
        normalize_range(x ::AbstractArray) =
            sweep(minimum)(-)(x) |> sweep(maximum)(/)
        num_cross_points = 0
        cross_loc=zeros(0)
        negderiv_smu_centred1 = sign(negderiv_smu[1] - negderiv_midrange)
        local ii = 0
        while (ii < len_denser)
            ii += 1
            negderiv_smu_centred0 = negderiv_smu_centred1
            negderiv_smu_centred1 = sign(negderiv_smu[ii] - negderiv_midrange)
            if negderiv_smu_centred0 != negderiv_smu_centred1
                num_cross_points += 1
                push!(cross_loc,ii/len_denser)
            end
        end ## while
        if num_cross_points>min(i.max_num_cross_points, len_raw * i.noise_factor) && mean(cross_loc)>0.4 && mean(cross_loc)<0.6
            println(cross_loc)
            return true
        else
            return false
        end
    end ## count_cross_points()
    split_vector_and_return_larger_quantile(
        x                   ::AbstractVector,
        len                 ::Int_T,      ## == length(x)
        idx                 ::Int_T,      ## 1 <= idx <= len
        q                   ::Float_T     ## 0 <= p <= 1
    ) = (1:idx, idx:len) |> mold(range -> quantile(x[range], q)) |> maximum

    ## note: largest_peak_idx calculated incorrectly in original code
    "Filter out small or randomly-occurring peaks."
    function peak_filter()
        ## function: normalize values to a range from 0 to 1
        normalize_range(x ::AbstractArray) =
            sweep(minimum)(-)(x) |> sweep(maximum)(/)
        #
        largest_peak = first(area_order)
        if  (split_vector_and_return_larger_quantile(
                normalize_range(negderiv_smu), ## originally normalize_range(negderiv), but why ???
                len_denser,
                peaks_raw[largest_peak].idx, ## index of peak with largest area
                i.norm_negderiv_quantile) ## larger_normd_qtv_of_two_sides
                    > i.max_norm_negderiv) ||
            (count_cross_points() > min(i.max_num_cross_points, len_raw * i.noise_factor)*1.9) ||
            linreg(temps, fluos)[2] >= 0.0 || normalize_range(negderiv_smu)[1]>0.6 || count_cross_points_dist() ||
            abs(maximum(negderiv_smu)-minimum(negderiv_smu))<3500

            return []
        end ## if
        ## else
        #
        ## Disabled because it caused false suppression of Tm peak reporting for
        ## `db_name_ = "20160309_chaipcr"; exp_id_ = 7`, well A2.
        ## Needs to be enabled because it is the only criterion that
        ## suppress false reporting of not-real peaks in
        ## `db_name_ = "20161103_chaipcr_ip152"; exp_id_ = 45`, well A1
        # temperature_max_negderiv = denser_temps[max_negderiv_smu[2]]
        # largest_peak_idx = summit_idc[indmax(areas_raw)] ## peak with largest area
        # top1_from_max = abs(denser_temps[largest_peak_idx] - temperature_max_negderiv)
        # if (top1_from_max > top1_from_max_ub)
        #     return []
        # end
        #
        ## when all the peaks are noise, no peak is expected
        ## to have a substantially larger area than the other peaks.
        ## Provided `i.max_num_peaks + 1 < num_peaks`, this is observed as
        ## `areas_raw[area_order[i.max_num_peaks + 1]] >= area_lb`
        ## implying `length(filtered_idc_topNp1) > i.max_num_peaks`
        ## If `num_peaks >= i.max_num_peaks` there is no problem.
        min_area = areas_raw[largest_peak] * i.min_normalized_area
        largest_idc = area_order[min(i.max_num_peaks + 1, num_peaks) |> from(1)]
        filtered_idc = largest_idc |> sift(idx -> areas_raw[idx] >= min_area)
        filtered_idc = filtered_idc |> sift(idx_ -> normalize_range(negderiv_smu)[peaks_raw[idx_].idx] > 0.3)
        return filtered_idc
    end ## peak_filter()

    ## << end of function definitions nested in mc_peak_analysis()

    debug(logger, "at mc_peak_analysis()")
    #
    ## filter out data points separated by narrow temperature intervals
    ## return temperatures as array, assuming no missing values
    (temps, fluos)    = df |> filter_too_close |> jitter_temps!
    len_raw           = length(temps)
    #
    ## negative derivative by central finite differencing (cfd)
    ## only used if the data array is too short to find peaks
	if (len_raw == 0)
        return output_type == McPeakLongOutput ?
			McPeakLongOutput() :
            McPeakShortOutput()
	end

    if (len_raw <= 3)
        slope = finite_diff(temps, fluos; nu = 1)

        return output_type == McPeakLongOutput ?
			McPeakLongOutput(
                observed_data       = i.reporting(hcat(temps, fluos, -slope))) :
            McPeakShortOutput(
                observed_data       = i.reporting(hcat(temps, fluos, -slope)))
    end ## if
    #
    ## else
    min_temp          = minimum(temps)
    max_temp          = maximum(temps)
    temp_span         = max_temp - min_temp
    #
    ## choose smoothing parameter
    span_smooth       = choose_span_smooth()
    #
    ## fit a cubic spline model in order to interpolate data points
    ## then smooth data and calculate slope at denser sequence of temperatures
    spl               = spline_model(smooth_raw_fluo())
    denser_temps      = interpolated_temperatures()
    smoothed_data     = smoothed_data_matrix()
    negderiv_smu      = smoothed_data[:, 3]
    max_negderiv_smu  = findmax(negderiv_smu)
    negderiv_midrange =
        mean([quantile(negderiv_smu, i.negderiv_range_low_quantile), max_negderiv_smu[1]])
    #
    ## extract data at observed temperatures
    len_denser        = length(denser_temps)
    observed_data     = smoothed_data[1:i.denser_factor:len_denser, :]
    #
    ## find peak and trough locations
    ## `sn` - summits and nadirs
    sn_idc            = find_summits_and_nadirs()
    sn_dict           = summits_and_nadirs(sn_idc)
    #
    peaks_raw         = find_peaks(sn_idc...)
    num_peaks         = length(peaks_raw)
    #
    ## return smoothed data if no peaks

    if num_peaks == 0
        return output_type == McPeakLongOutput ?
			McPeakLongOutput(
                observed_data       = i.reporting(observed_data),
                smoothed_data       = smoothed_data,
                negderiv_midrange   = negderiv_midrange,
                extremes            = sn_dict) :
            McPeakShortOutput(
                observed_data       = i.reporting(observed_data))
    end ## if no peaks
    #
    ## keep only the biggest peak(s)
    areas_raw = peaks_raw |> their(:area)
    area_order = sortperm(areas_raw, rev = true)
    peaks_filtered = peaks_raw[peak_filter()]
    sort!(peaks_filtered,by=x->negderiv_smu[x.idx],rev=true)
    length(peaks_filtered)>i.max_num_peaks?peaks_filtered=peaks_filtered[1:i.max_num_peaks]:nothing
    return output_type == McPeakLongOutput ?
        McPeakLongOutput(
            i.reporting(observed_data),
            i.reporting(peaks_filtered),
            smoothed_data,
            negderiv_midrange,
            sn_dict,
            i.reporting(peaks_raw[area_order]),
            length(peaks_filtered) > 0) :
        McPeakShortOutput(
            i.reporting(observed_data),
            i.reporting(peaks_filtered))

end ## mc_peak_analysis()

#==============================================================================#

## function called by mc_peak_analysis() >>

"Jitter duplicated temperature values so that all the elements become unique."
function mutate_dups(
    vec_2mut        ::AbstractVector,
    jitter_constant ::Float_T,
)
    debug(logger, "at mutate_dups()")
    vec_len       = vec_2mut |> length
    vec_uniq      = vec_2mut |> unique |> sort
    vec_uniq_len  = vec_uniq |> length
    ## return if no ties
    (vec_len == vec_uniq_len) && return vec_2mut
    ## find ties
    order_to   = sortperm(vec_2mut)
    order_back = sortperm(order_to)
    vec_sorted = (vec_2mut + 0.0)[order_to]
    dups = find(diff(vec_sorted) .== 0.0) .+ 1
    ## calculate value of jitter
    add1 = jitter * median(diff(vec_uniq))
    ## break ties
    accumulator1 = 0
    accumulator2 = 0
    local i
    for i in eachindex(dups)
        multiway_tie  = i > 1 && dups[i] - dups[i - 1] == 1
        accumulator1 += multiway_tie
        accumulator2  = multiway_tie ? accumulator2 : accumulator1
        rank = accumulator1 - accumulator2 + 1
        vec_sorted[dups[i]] += add1 * rank
    end
    return vec_sorted[order_back]
end ## mutate_dups()

max_stars_repo_name: MakerButt/chaipcr

max_issues_repo_name: MakerButt/chaipcr

max_forks_repo_name: MakerButt/chaipcr


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

module TestFunctions

using Test
using MathOptInterface
const MOI = MathOptInterface

"""
    test_isbits()

Test isbit-ness of VariableIndex and terms.

It is important that these struct remain isbits as otherwise, the performance of
the corresponding function will be deteriored. These tests explicit this design
choice and makes sure that it remains the case.
"""
function test_isbits()
    x = @inferred MOI.VariableIndex(1)
    @test isbits(x)
    at = @inferred MOI.ScalarAffineTerm(1.0, x)
    @test isbits(at)
    @test isbits(@inferred MOI.VectorAffineTerm(1, at))
    qt = @inferred MOI.ScalarQuadraticTerm(1.0, x, x)
    @test isbits(qt)
    @test isbits(@inferred MOI.VectorQuadraticTerm(1, qt))
end

function test_functions_broadcast()
    x = MOI.VariableIndex(1)
    y = MOI.VariableIndex(2)
    z = MOI.VariableIndex(3)
    function sum_indices(sv1::MOI.VariableIndex, sv2::MOI.VariableIndex)
        return sv1.value + sv2.value
    end
    @test sum_indices.(x, [y, z]) == [3, 4]
end

function test_functions_copy_VectorOfVariables()
    x = MOI.VariableIndex(1)
    y = MOI.VariableIndex(2)
    z = MOI.VariableIndex(3)
    f = MOI.VectorOfVariables([x, y])
    f_copy = copy(f)
    f_copy.variables[2] = z
    @test f.variables[2] == y
end

function test_functions_convert_VariableIndex()
    model = MOI.Utilities.Model{Float64}()
    x = MOI.add_variable(model)
    f_vov = convert(MOI.VectorOfVariables, x)
    @test f_vov ≈ MOI.VectorOfVariables([x])
    f_vaf = convert(MOI.VectorAffineFunction{Float64}, x)
    @test f_vaf ≈ MOI.VectorAffineFunction(
        [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, x))],
        [0.0],
    )
    f_vqf = convert(MOI.VectorQuadraticFunction{Float64}, x)
    @test f_vqf ≈ MOI.VectorQuadraticFunction(
        MOI.VectorQuadraticTerm{Float64}[],
        [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, x))],
        [0.0],
    )
end

function test_functions_convert_ScalarAffineFunction()
    model = MOI.Utilities.Model{Float64}()
    x = MOI.add_variable(model)
    f = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(2.0, x)], 1.0)
    @test_throws(MethodError, convert(MOI.VectorOfVariables, f))
    f_vaf = convert(MOI.VectorAffineFunction{Float64}, f)
    @test f_vaf ≈ MOI.VectorAffineFunction(
        [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(2.0, x))],
        [1.0],
    )
    f_vqf = convert(MOI.VectorQuadraticFunction{Float64}, f)
    @test f_vqf ≈ MOI.VectorQuadraticFunction(
        MOI.VectorQuadraticTerm{Float64}[],
        [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(2.0, x))],
        [1.0],
    )
end

function test_functions_convert_ScalarQuadraticFunction()
    model = MOI.Utilities.Model{Float64}()
    x = MOI.add_variable(model)
    f = MOI.ScalarQuadraticFunction(
        [MOI.ScalarQuadraticTerm(3.0, x, x)],
        [MOI.ScalarAffineTerm(2.0, x)],
        1.0,
    )
    @test_throws(MethodError, convert(MOI.VectorOfVariables, f))
    @test_throws(MethodError, convert(MOI.VectorAffineFunction{Float64}, f))
    f_vqf = convert(MOI.VectorQuadraticFunction{Float64}, f)
    @test f_vqf ≈ MOI.VectorQuadraticFunction(
        MOI.VectorQuadraticTerm{Float64}[MOI.VectorQuadraticTerm(
            1,
            MOI.ScalarQuadraticTerm(3.0, x, x),
        )],
        [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(2.0, x))],
        [1.0],
    )
end

function test_isapprox_VectorOfVariables()
    x = MOI.VariableIndex(1)
    y = MOI.VariableIndex(2)
    @test MOI.VectorOfVariables([x, y]) == MOI.VectorOfVariables([x, y])
    @test MOI.VectorOfVariables([y, x]) != MOI.VectorOfVariables([x, y])
    @test MOI.VectorOfVariables([x, x]) != MOI.VectorOfVariables([x])
    @test MOI.VectorOfVariables([x]) != MOI.VectorOfVariables([y])
    return
end

function test_isapprox_ScalarAffineFunction()
    x = MOI.VariableIndex(1)
    y = MOI.VariableIndex(2)
    z = MOI.VariableIndex(3)
    @test MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1, 1], [x, z]), 1) ≈
          MOI.ScalarAffineFunction(
        MOI.ScalarAffineTerm.([1.0, 1e-7, 1.0], [x, y, z]),
        1.0,
    ) atol = 1e-6
    @test MOI.ScalarAffineFunction(
        MOI.ScalarAffineTerm.([1.0, 1e-7], [x, y]),
        1.0,
    ) ≈ MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1, x)], 1) atol = 1e-6
    f = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2, 4], [x, y]), 6)
    g = deepcopy(f)
    @test g ≈ f
    @test f ≈ g
    f.terms[2] = MOI.ScalarAffineTerm(3, y)
    @test !(g ≈ f)
    @test !(f ≈ g)
    return
end

function test_isapprox_VectorAffineFunction()
    x = MOI.VariableIndex(1)
    y = MOI.VariableIndex(2)
    f = MOI.VectorAffineFunction(
        MOI.VectorAffineTerm.(
            [1, 1, 2],
            MOI.ScalarAffineTerm.([2, 4, 3], [x, y, y]),
        ),
        [6, 8],
    )
    g = deepcopy(f)
    @test f ≈ g
    f.terms[3] = MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(9, y))
    @test !(f ≈ g)
    @test !(g ≈ f)
    push!(f.terms, MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(-6, y)))
    @test f ≈ g
    @test g ≈ f
    return
end

function test_isapprox_ScalarQuadraticFunction()
    x = MOI.VariableIndex(1)
    y = MOI.VariableIndex(2)
    z = MOI.VariableIndex(3)
    f = MOI.ScalarQuadraticFunction(
        MOI.ScalarQuadraticTerm.([1, 2, 3], [x, y, x], [x, y, y]),
        [MOI.ScalarAffineTerm(3, x)],
        8,
    )
    g = deepcopy(f)
    @test f ≈ g
    push!(f.affine_terms, MOI.ScalarAffineTerm(2, y))
    @test !(f ≈ g)
    g = deepcopy(f)
    push!(f.quadratic_terms, MOI.ScalarQuadraticTerm(2, y, x))
    @test !(f ≈ g)
    @test !(g ≈ f)
    push!(f.quadratic_terms, MOI.ScalarQuadraticTerm(-2, y, x))
    @test f ≈ g
    @test g ≈ f
    return
end

function test_isapprox_VectorQuadraticFunction()
    x = MOI.VariableIndex(1)
    y = MOI.VariableIndex(2)
    z = MOI.VariableIndex(3)
    f = MOI.VectorQuadraticFunction(
        MOI.VectorQuadraticTerm.(
            [1, 1, 2],
            MOI.ScalarQuadraticTerm.([1, 2, 3], [x, y, x], [x, y, y]),
        ),
        MOI.VectorAffineTerm.(
            [1, 2, 1],
            MOI.ScalarAffineTerm.([3, 1, 1], [x, x, y]),
        ),
        [10, 11, 12],
    )
    g = deepcopy(f)
    @test f ≈ g
    @test g ≈ f
    f.affine_terms[1] = MOI.VectorAffineTerm(3, MOI.ScalarAffineTerm(4, x))
    @test !(f ≈ g)
    @test !(g ≈ f)
    push!(g.affine_terms, MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(-3, x)))
    push!(g.affine_terms, MOI.VectorAffineTerm(3, MOI.ScalarAffineTerm(4, x)))
    @test f ≈ g
    @test g ≈ f
    f.quadratic_terms[1] =
        MOI.VectorQuadraticTerm(3, MOI.ScalarQuadraticTerm(1, x, x))
    @test !(f ≈ g)
    @test !(g ≈ f)
    return
end

function test_isapprox_issue_1483()
    x = MOI.ScalarQuadraticFunction(
        MOI.ScalarQuadraticTerm{Float16}[],
        MOI.ScalarAffineTerm{Float16}[],
        Float16(0.0),
    )
    y = MOI.ScalarQuadraticFunction(
        MOI.ScalarQuadraticTerm{Float16}[],
        MOI.ScalarAffineTerm{Float16}[MOI.ScalarAffineTerm(
            Float16(1.0),
            MOI.VariableIndex(1234),
        )],
        Float16(0.0),
    )
    z = MOI.ScalarQuadraticFunction(
        MOI.ScalarQuadraticTerm{Float16}[],
        MOI.ScalarAffineTerm{Float16}[],
        Float16(0.0),
    )
    @test !(x ≈ y)
    @test !(y ≈ x)
    @test x ≈ z
    @test z ≈ x
    @test !(y ≈ z)
    @test !(z ≈ y)
    return
end

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$name", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
end

end

TestFunctions.runtests()

max_stars_repo_name: JuliaOpt/MathOptInterface

max_issues_repo_name: JuliaOpt/MathOptInterface

max_forks_repo_name: JuliaOpt/MathOptInterface


text:
using FastGroupBy,InternedStrings
import FastGroupBy.mask16bit
const N = 100_000_000; const K = 100
srand(1);
@time sample_space = [string(rand(Char.(97:97+25), rand(1:16))...) for k in 1:N÷K];
@time svec = rand(sample_space, N);
@time radixsort!(UInt192,svec)
issorted(svec)

using FastGroupBy
import FastGroupBy.mask16bit
const N = 100_000_000; const K = 100
srand(1);
@time sample_space = [string(rand(Char.(97:97+25), rand(1:24))...) for k in 1:N÷K];
@time svec = rand(sample_space, N);
@time radixsort!(UInt192,svec) #88
issorted(svec)

const N = 100_000_000; const K = 100
srand(1);
@time sample_space = [string(rand(Char.(97:97+25), rand(1:24))...) for k in 1:N÷K];
@time svec = rand(sample_space, N);
@time sort!(svec) #78 seconds
issorted(svec)

@time radixsort!(svec)
issorted(svec)

T = UInt192
i = 1
@time vs = FastGroupBy.load_bits.(T, svec, Int(i-1)*sizeof(T));
@time FastGroupBy.sorttwo_lsd16!(vs1, svec);

srand(1);
@time sample_space = [string(rand(Char.(97:97+25), rand(1:24))...) for k in 1:N÷K];
@time svec = rand(sample_space, N);
@time vs1 = FastGroupBy.load_bits_fast.(T, svec);
@time FastGroupBy.sorttwo_lsd!(vs1, svec)

using BenchmarkTools
T = UInt128
i = 1
srand(1);
@time sample_space = [string(rand(Char.(97:97+25), rand(1:16))...) for k in 1:N÷K];
@time svec = rand(sample_space, N);
function fn1(svec, T)
    @time vs1 = FastGroupBy.load_bits_fast.(T, svec); # 11
    @time FastGroupBy.sorttwo_lsd16!(vs1, svec); # 21
end;
@time fn1(svec, T); #34 #37
issorted(svec)

srand(1);
@time sample_space = [string(rand(Char.(97:97+25), rand(1:16))...) for k in 1:N÷K];
@time svec = rand(sample_space, N);
function fn2(svec, T)
    @time vs1 = FastGroupBy.load_bits_fast.(T, svec); # 11
    @time FastGroupBy.sorttwo_lsd!(vs1, svec); # 21
end
@time fn2(svec, T);
issorted(svec, T) #31 #44

srand(1);
@time sample_space = [string(rand(Char.(97:97+25), rand(1:16))...) for k in 1:N÷K];
@time svec = rand(sample_space, N);
function fn3(svec, T)
    @time vs1 = FastGroupBy.load_bits_fast.(T, svec); # 11
    @time FastGroupBy.sorttwo!(ntoh.(vs1), svec); #29
end
@time fn3(svec, T); #40 #30
issorted(svec)

srand(1);
@time sample_space = [string(rand(Char.(97:97+25), rand(1:16))...) for k in 1:N÷K];
@time svec = rand(sample_space, N);
function fn4(svec, T)
    @time vs1 = FastGroupBy.load_bits_fast_ntoh.(T, svec); #2
    @time FastGroupBy.sorttwo!(vs1, svec) #18
end
@time fn4(svec, T); #47 #40
issorted(svec)

@time sorttwo!(vs, svec)
@time issorted(svec)

@time radixsort!(UInt192,svec)
issorted(svec)

T = UInt
i = 1
@time vs = FastGroupBy.load_bits.(T, svec, Int(i-1)*sizeof(T));
@time vs1 = load_bits.(T, svec, Int(i-1)*sizeof(T));

@time radixsort!(UInt128,svec)
issorted(svec)
index = svec;


text:
module  TestUtils

using PopGenCore
using Test

x = @nancycats

@testset "Utilities" begin
    @testset "GeneralUtils.jl" begin
        @test typeof(copy(x)) == PopData
        @test size(x) == (;samples = 237, loci = 9)
        @test sort((2,1,3)) == (1,2,3)
        @test convertcoord("-41 31.52") == -41.5253f0
        @test convertcoord.(["-41:31.52", "25 11:54S"]) == Float32[-41.5253, -25.1983]
        @test length(loci(dropmonomorphic(x))) == 9
        @test length(loci(dropmultiallelic(x))) == 0
    end

    @testset "GenotypeUtils.jl" begin
        @test allelecount([(1,2), (1,4)]) == 3
        @test allelecount(x.genodata.genotype[1:10]) == 4
        @test alleles([(1,2), (1,4)]) == [1,2,1,4]
        @test alleles([(1,2), (1,4)], false) == [1,2,1,4]
        @test alleles([(1,2), (1,4), missing], false) == [1,2,1,4]
        @test all(alleles([(1,2), (1,4), missing], true) .=== [1,2,1,4,missing])
        @test uniquealleles([(1,2), (1,4), missing]) == [1,2,4]
        @test size(locidataframe(x)) == (9,237)
        @test size(locimatrix(x)) == (237,9)
        @test length(phasedmatrix(x)) == 2
        @test size.(phasedmatrix(x)) == [(237,9),(237,9)]
        @test eltype(phasedmatrix(x)) == Matrix{Union{Missing, Int16}}
        @test eltype(phasedmatrix(x)[1]) == Union{Missing, Int16}
    end

    @testset "ioUtils.jl" begin
        @test findploidy(x.genodata.genotype[1:20]) == 2
        @test findploidy([(1,2,3), (1,2,2)]) == 3
        @test phase("128114", Int16, 3) == (114,128)
        @test phase("128114", Int16, 2) == (12,14,81)
        @test phase("128", Int16, 3) == (128,)
        @test eltype(phase("128", Int16, 3)) == Int16
        @test eltype(phase("110", Int8, 3)) == Int8
        @test phase(128114, Int16, 3) == (114,128)
        @test phase(128114, Int16, 2) == (12,14,81)
        @test phase(128, Int16, 3) == (128,)
        @test eltype(phase(128, Int16, 3)) == Int16
        @test eltype(phase(110, Int8, 3)) == Int8
        @test phase(missing, Int8, 124) === missing
        @test unphase((1,2,3,4), digits = 3) == "001002003004"
        @test unphase(missing, digits = 2, ploidy = 2, miss = -9) == "-9"
        @test unphase(missing, digits = 2, ploidy = 2, miss = 0) == "0000"
    end

    @testset "MathUtils.jl" begin
        @test countnonzeros([1,2,0,0,3,4,5]) == 5
        @test reciprocal(0) == 0.0
        @test reciprocal(2) == 1/2
        @test reciprocalsum([1,2,5,0,0,4]) == (1/1 + 1/2 + 1/5 + 1/4)
    end

    @testset "MissingUtils.jl" begin
        @test nonmissing(fill(missing, 5)) == 0
        @test nonmissing([1,2,3,missing, missing]) == 3
        @test nonmissing(x, "fca8") == 217
        @test nonmissings([1,2,missing, missing, 3], [3,missing, missing, 4, 5]) == [1,5]
    end
end
end


text:
using Automa
using Automa.RegExp: @re_str
const re = Automa.RegExp

ab = re"ab*"
ab.when = :cond
c = re"c"
pattern = re.cat(ab, c)

write("preconditions.dot", Automa.machine2dot(Automa.compile(pattern)))
run(`dot -Tpng -o src/figure/preconditions.png preconditions.dot`)


text:
struct NewtonResult{T, N}
    converged::Bool
    iterations::Int
    solution::SVector{N, T}
    residual::T
end

"""
    newton(F, x0, maxiterations, ftol)

Use a trust-region newton iteration[^1] to compute a root of the fiber `F` using maximal `iterations` iterations.
Convergence is declared if the infinity norm is less than `ftol`.

[^1]: Chapter 11.2 in Wright, Stephen J., and Jorge Nocedal. "Numerical optimization." Springer Science 35.67-68 (1999): 7.
"""
function newton(
    F::AbstractFiber{T, N, N},
    x0::SVector{N, T},
    iterations::Integer,
    ftol::Float64,
    abort=always_continue) where {N, T}

    factor = 1.0

    x = x0     # Current point

    r = evaluate(F, x)
    J = jacobian(F, x, true)

    it = 0
    norm_residual = maxnorm(r)
    converged = norm_residual < ftol

    Δ = factor * norm(x)
    if Δ == 0.0
        Δ = factor
    end

    η = 1e-4
    force_exit = false
    while !(converged || force_exit) && it < iterations
        it += 1
        # Compute proposed iteration step
        p, is_ok = dogleg(r, J, Δ)
        if !is_ok
            break
        end
        x += p

        r_new = evaluate(F, x)

        # Ratio of actual to predicted reduction (equation 11.45 in N&W)
        r_predict = r + J * p
        norm_r2 = sum(abs2,r)
        ρ = (norm_r2 - sum(abs2,r_new)) / (norm_r2 - sum(abs2,r_predict))

        if ρ > η
            # Successful iteration
            r = r_new
            J = jacobian(F, x, true)

            norm_residual = maxnorm(r)
            converged = norm_residual < ftol

            if abort(x)
                force_exit = true
            end
        else
            x -= p
            converged = false
        end

        # Update size of trust region
        if ρ < 0.1
            Δ *= 0.5
        elseif ρ >= 0.9
            Δ = 2 * norm(p)
        elseif ρ >= 0.5
            Δ = max(Δ, 2 * norm(p))
        end
    end

    NewtonResult(converged, it,  x, norm_residual)
end

always_continue(x) = false

"""
    dogleg(r, d, J, Δ::Real)
See Procedure 11.6 in "Numerical optimization."[^1] from Wright and Nocedal

[^1]: Wright, Stephen J., and Jorge Nocedal. "Numerical optimization." Springer Science 35.67-68 (1999): 7.
"""
function dogleg(r::SVector{N, T}, J::SMatrix{N, N, T}, Δ::Real) where {N,T}
    p_i = r
    # try
    p_i, is_ok = mysolve(J, r)
    if !is_ok
        return (p_i, false)
    end
    # catch
    #     # use Moore-Penrose inverse instead
    #     println("pinv")
    #     p_i = -SMatrix{N, N}(pinv(Matrix(J))) * r
    # end

    # Test if Gauss-Newton step is within the region
    if norm(p_i) <= Δ
        return p_i, true   # accepts equation 4.13 from N&W for this step
    else
        g = transposed_mul(J, r)

        # compute Cauchy point
        norm_g2 = sum(abs2, g)
        p_c = - norm_g2 / sum(abs2, J*g) .* g

        if norm(p_c) >= Δ
            # Cauchy point is out of the region, take the largest step along
            # gradient direction
            return g .* (-Δ / sqrt(norm_g2)), true
        else
            p_diff = p_i - p_c
            # Compute the optimal point on dogleg path
            b = 2 * dot(p_c, p_diff)
            a = dot(p_diff, p_diff)
            τ = (-b+sqrt(b^2 - 4a*(dot(p_c, p_c) - Δ^2)))/(2a)
            return p_c .+ τ .* p_diff, true
        end
    end
end

# We have a variant of Newton's method for underdetermined system
function newton(
    F::AbstractFiber{T, M, N},
    x0::SVector{N, T},
    iterations::Integer,
    ftol::Float64, abort=always_continue) where {M, N, T}

    x = x0     # Current point

    r = evaluate(F, x)
    norm_residual = maxnorm(r)
    J = jacobian(F, x, true)
    it = 0
    converged = norm_residual < ftol
    force_exit = false
    while !(converged || force_exit) && it < iterations
        it += 1
        # Compute proposed iteration step
        p, is_ok = mysolve(J, r)
        if !is_ok
            break
        end
        x -= p
        if abort(x)
            force_exit = true
        end

        r = evaluate(F, x)
        norm_residual = maxnorm(r)
        converged = norm_residual < ftol
    end

    NewtonResult(converged, it,  x, norm_residual)
end

function mysolve(M::SMatrix{2, 2, T}, u::SVector{2, T}) where T
    @inbounds minus_m3 = -M[3]
    @inbounds denom = muladd(M[1], M[4], minus_m3 * M[2])
    if abs(denom) < eps(T) * 10
        return u, false
    end
    invdenom = -inv(denom)
    @inbounds res = SVector(muladd(M[4], u[1], minus_m3 * u[2]) * invdenom, muladd(M[1], u[2], -M[2]*u[1]) * invdenom)
    return res, true
end

function mysolve(a::SMatrix{3, 3, T}, b::SVector{3, T}) where T
    d = det(a)
    if abs(d) < eps(T) * 10
        return b, false
        # throw(Base.LinAlg.SingularException(1))
    end
    invd = -inv(d)
    @inbounds out = SVector{3, T}(
        ((a[2,2]*a[3,3] - a[2,3]*a[3,2])*b[1] +
            (a[1,3]*a[3,2] - a[1,2]*a[3,3])*b[2] +
            (a[1,2]*a[2,3] - a[1,3]*a[2,2])*b[3]) * invd,
        ((a[2,3]*a[3,1] - a[2,1]*a[3,3])*b[1] +
            (a[1,1]*a[3,3] - a[1,3]*a[3,1])*b[2] +
            (a[1,3]*a[2,1] - a[1,1]*a[2,3])*b[3]) * invd,
        ((a[2,1]*a[3,2] - a[2,2]*a[3,1])*b[1] +
            (a[1,2]*a[3,1] - a[1,1]*a[3,2])*b[2] +
            (a[1,1]*a[2,2] - a[1,2]*a[2,1])*b[3]) * invd)
    return out, true
end

function mysolve(A::SMatrix{M, N, T}, u::SVector{M, T}) where {T, M, N}
    Q, R = qr(A')
    z = forward_substitution(R', u)
    Q * z, true
end

function forward_substitution(R::SMatrix{1, 1, T}, u::SVector{1, T}) where T
    SVector(u[1] / R[1, 1])
end
function forward_substitution(R::SMatrix{2, 2, T}, u::SVector{2, T}) where T
    x1 = u[1] / R[1, 1]
    x2 = (u[2] - R[2, 1] * x1) / R[2, 2]
    SVector(x1, x2)
end

function forward_substitution(R::SMatrix{3, 3, T}, u::SVector{3, T}) where T
    x1 = u[1] / R[1, 1]
    x2 = (u[2] - R[2, 1] * x1) / R[2, 2]
    x3 = (u[3] - R[3, 1] * x1 - R[3, 2] * x2) / R[3, 3]
    SVector(x1, x2, x3)
end

"""
    transposed_mul(J, r)

Compute `J'*r`.
"""
function transposed_mul(J::SMatrix{2, 2, T}, r::SVector{2, T}) where T
    @inbounds a = muladd(J[1, 1], r[1], J[2, 1] * r[2])
    @inbounds b = muladd(J[1, 2], r[1], J[2, 2] * r[2])
    SVector{2, T}(a, b)
end

function transposed_mul(J::SMatrix{3, 3, T}, r::SVector{3, T}) where T
    @inbounds a = J[1, 1] * r[1] + J[2, 1] * r[2] + J[3, 1] * r[3]
    @inbounds b = J[1, 2] * r[1] + J[2, 2] * r[2] + J[3, 2] * r[3]
    @inbounds c = J[1, 3] * r[1] + J[2, 3] * r[2] + J[3, 3] * r[3]
    SVector{3, T}(a, b, c)
end

function transposed_mul(J::SMatrix{2, 3, T}, r::SVector{2, T}) where T
    @inbounds a = J[1, 1] * r[1] + J[2, 1] * r[2]
    @inbounds b = J[1, 2] * r[1] + J[2, 2] * r[2]
    @inbounds c = J[1, 3] * r[1] + J[2, 3] * r[2]
    SVector{3, T}(a, b, c)
end

@inline maxnorm(x::SVector{2}) = max(abs(x[1]), abs(x[2]))
@inline maxnorm(x::SVector{3}) = max(abs(x[1]), abs(x[2]), abs(x[3]))


text:
const B = 1
const KB = 1024
const MB = 1024KB
const GB = 1024MB
const TB = 1024GB

SCHEDULER_MAX_CPU = nthreads() > 1 ? nthreads()-1 : Sys.CPU_THREADS
SCHEDULER_MAX_MEM = round(Int, Sys.total_memory() * 0.9)
SCHEDULER_UPDATE_SECOND = 0.6

const JOB_QUEUE = Vector{Job}()
JOB_QUEUE_LOCK = SpinLock()
const JOB_QUEUE_OK = Vector{Job}()  # jobs not queuing
JOB_QUEUE_MAX_LENGTH = 10000

SCHEDULER_BACKUP_FILE = ""

const QUEUING = :queuing
const RUNNING = :running
const DONE = :done
const FAILED = :failed
const CANCELLED = :cancelled

const PAST = :past # super set of DONE, FAILED, CANCELLED

function release_lock()
    global JOB_QUEUE_LOCK
    unlock(JOB_QUEUE_LOCK)
end

function wait_for_lock()
    global JOB_QUEUE_LOCK
    while !trylock(JOB_QUEUE_LOCK)
        sleep(0.05)
    end
end

"""
    submit!(job::Job)

Submit the job. If `job.create_time == 0000-01-01T00:00:00 (default)`, it will change to the time of submission.
"""
function submit!(job::Job)
    global JOB_QUEUE
    global QUEUING

    if scheduler_status(verbose=false) == :not_running
        @error "Scheduler is not running. Please start scheduler by using scheduler_start()"
        return job
    end

    # cannot run a backup task
    if isnothing(job.task)
        if job.state in (RUNNING, QUEUING)
            job.state = CANCELLED
        end
        @error "Cannot submit a job recovered from backup!" job
        return job
    end

    # check task state
    if istaskstarted(job.task) || job.state !== QUEUING
        @error "Cannot submit running/done/failed/canceled job!" job
        return job
    end

    wait_for_lock()
    try
        # check duplicate submission (queuing)
        for existing_job in JOB_QUEUE
            if existing_job === job
                error("Cannot re-submit the same job! Job ID: $(job.id)")
            end
        end

        # job.state = QUEUING
        if job.create_time == DateTime(0)
            job.create_time = now()
        end

        push!(JOB_QUEUE, job)
    catch e
        rethrow(e)
    finally
        release_lock()
    end

    return job
end

"""
    unsafe_run!(job::Job) :: Bool

Jump the queue and run `job` immediately, no matter what other jobs are running or waiting. If successful initiating to run, return `true`, else `false`.
"""
function unsafe_run!(job::Job) :: Bool
    global QUEUING
    global RUNNING
    global FAILED
    global DONE
    global CANCELLED

    # cannot run a backup task
    if isnothing(job.task)
        @error "Cannot run a job recovered from backup!" job
        if job.state in (RUNNING, QUEUING)
            job.state = CANCELLED
        end
        return false
    end

    try
        # set_non_sticky!(job)
        # schedule(job.task)
        schedule_thread(job)
        job.start_time = now()
        job.state = RUNNING
        true
    catch e
        # check status when fail
        if istaskfailed(job.task)
            if job.state !== CANCELLED
                job.state = FAILED
                @error "A job has failed: $(job.id)" exception=job.task.result
            end
            false
        elseif istaskdone(job.task)
            job.state = DONE
            false
        elseif istaskstarted(job.task)
            job.state = RUNNING
            false
        else
            rethrow(e)
            false
        end
    end
end

"""
    cancel!(job::Job)

Cancel `job`, stop queuing or running.
"""
function cancel!(job::Job)
    wait_for_lock()
    res = try
        unsafe_cancel!(job)
    catch e
        rethrow(e)
    finally
        release_lock()
    end
    res
end

"""
    unsafe_cancel!(job::Job)

Caution: it is unsafe and should only be called within lock. Do not call from other module.
"""
function unsafe_cancel!(job::Job)
    global CANCELLED
    global RUNNING

    # cannot cancel a backup task (can never be run)
    if isnothing(job.task)
        if job.state in (RUNNING, QUEUING)
            job.state = CANCELLED
        end
        return job.state
    end

    # no need to cancel finished ones
    if istaskfailed(job.task)
        if job.state !== CANCELLED
            job.state = FAILED
            # job.task.result isa Exception, notify errors
            @error "A job has failed: $(job.id)" exception=job.task.result
        end
        free_thread(job)
        return job.state
    elseif istaskdone(job.task)
        if job.state !== DONE
            job.state = DONE
        end
        free_thread(job)
        return job.state
    end

    try
        schedule(job.task, InterruptException(), error=true)
        job.stop_time = now()
        job.state = CANCELLED
    catch e
        unsafe_update_state!(job)
        if job.state === RUNNING
            @error "unsafe_cancel!(job): cannot cancel a job." job
        end
        job.state
    finally
        free_thread(job)
        return job.state
    end
end

### Scheduler main

function scheduler()
    global SCHEDULER_UPDATE_SECOND
    while true
        update_queue!()
        sleep(SCHEDULER_UPDATE_SECOND)
    end
end

function update_queue!()
    global SCHEDULER_MAX_CPU
    # step 0: cancel jobs if run time > wall time
    cancel_jobs_reaching_wall_time!()
    # step 1: update running jobs' states
    update_state!()
    # step 2: migrate finished jobs to JOB_QUEUE_OK from JOB_QUEUE
    migrate_finished_jobs!()
    # step 3: compute current CPU/MEM usage
    used_ncpu, used_mem = current_usage()
    ncpu_available = SCHEDULER_MAX_CPU - used_ncpu
    mem_available = SCHEDULER_MAX_MEM - used_mem
    # step 4: sort by priority
    update_queue_priority!()
    # step 5: run queuing jobs
    run_queuing_jobs(ncpu_available, mem_available)
end

function cancel_jobs_reaching_wall_time!()
    global JOB_QUEUE
    global RUNNING
    wait_for_lock()
    try
        for job in JOB_QUEUE
            if job.state === RUNNING
                elapsed_time =  now() - job.start_time
                remaining_time = Millisecond(job.wall_time) - elapsed_time
                if remaining_time.value < 0
                    unsafe_cancel!(job)  # cannot call cancel! because lock is on
                end
            end
        end
    catch e
        rethrow(e)
    finally
        release_lock()
    end
end

"""
    unsafe_update_state!(job::Job)

Update the state of `job` from `job.task` when `job.state === :running`.

Caution: it is unsafe and should only be called within lock.
"""
function unsafe_update_state!(job::Job)
    global RUNNING
    global DONE
    global FAILED
    if job.state === RUNNING
        task_state = job.task.state
        if task_state === DONE
            job.stop_time = now()
            free_thread(job)
            job.state = DONE
        elseif task_state === FAILED
            job.stop_time = now()
            free_thread(job)
            @error "A job has failed: $(job.id): $(job.name)" exception=job.task.result
            job.state = FAILED
        end
    else
        job.state
    end
end

"""
    update_state!()

Update the state of each `job` in JOB_QUEUE from `job.task` when `job.state === :running`.
"""
function update_state!()
    global JOB_QUEUE
    wait_for_lock()
    try
        foreach(unsafe_update_state!, JOB_QUEUE)
    catch e
        rethrow(e)
    finally
        release_lock()
    end
    return
end

function migrate_finished_jobs!()
    global JOB_QUEUE
    global JOB_QUEUE_OK
    global JOB_QUEUE_MAX_LENGTH
    wait_for_lock()
    try
        finished_indices = map(j -> !(j.state === QUEUING || j.state === RUNNING), JOB_QUEUE)
        append!(JOB_QUEUE_OK, JOB_QUEUE[finished_indices])
        deleteat!(JOB_QUEUE, finished_indices)
    catch e
        rethrow(e)
    finally
        release_lock()
    end

    # delete JOB_QUEUE_OK if too many
    n_delete = length(JOB_QUEUE_OK) - JOB_QUEUE_MAX_LENGTH
    n_delete > 0 && deleteat!(JOB_QUEUE_OK, 1:n_delete)
    return
end

function current_usage()
    global JOB_QUEUE
    global RUNNING
    cpu_usage = 0
    mem_usage = 0
    wait_for_lock()
    try
        for job in JOB_QUEUE
            if job.state === RUNNING
                cpu_usage += job.ncpu
                mem_usage += job.mem
            end
        end
    catch e
        rethrow(e)
    finally
        release_lock()
    end
    return cpu_usage, mem_usage
end

function update_queue_priority!()
    global JOB_QUEUE
    wait_for_lock()
    try
        sort!(JOB_QUEUE, by=get_priority)
    catch e
        rethrow(e)
    finally
        release_lock()
    end
    return
end
get_priority(job::Job) = job.priority

function run_queuing_jobs(ncpu_available::Int, mem_available::Int)
    global JOB_QUEUE
    global QUEUING
    wait_for_lock()
    try
        for job in JOB_QUEUE
            ncpu_available < 1 && break
            mem_available  < 1 && break
            if job.state === QUEUING && job.schedule_time < now() && job.ncpu <= ncpu_available && job.mem <= mem_available && is_dependency_ok(job)
                if unsafe_run!(job)
                    ncpu_available -= job.ncpu
                    mem_available  -= job.mem
                end
            end
        end
    catch e
        rethrow(e)
    finally
        release_lock()
    end
    return ncpu_available, mem_available
end

"""
Caution: run it within lock only.
"""
function is_dependency_ok(job::Job)
    for dep in job.dependency
        state = dep.first
        dep_id = dep.second
        dep_job = job_query_by_id_no_lock(dep_id)

        if isnothing(dep_job)
            return false
        end

        state == dep_job.state && continue

        state == PAST &&
            dep_job.state in (FAILED, CANCELLED, DONE) && continue

        if dep_job.state in (FAILED, CANCELLED)
            # change job state to cancelled
            @warn "Cancel job ($(job.id)) because one of its dependency ($(dep_job.id)) is failed or cancelled."
            job.state = CANCELLED
            return false
        elseif dep_job.state == DONE
            # change job state to cancelled
            @warn "Cancel job ($(job.id)) because one of its dependency ($(dep_job.id)) is done, but the required state is $(state)."
            job.state = CANCELLED
            return false
        end

        return false
    end
    return true
end


text:
import Pkg
Pkg.activate(".")

using GBIF
using SimpleSDMLayers
using StatsBase
using Plots
using ArchGDAL
using Images
using StatsPlots.PlotMeasures
using ProgressMeter

include("bioclim.jl")
include("landcover.jl")

# Load raccoon emoji
img = load("images/raccoon_freepik_cropped.png")

## Get occurrences & climatic data

# Get the GBIF taxon
raccoon = GBIF.taxon("Procyon lotor", rank=:SPECIES, strict=true)

# Get the initial occurrences
raccoon_occ = GBIF.occurrences(raccoon, "limit" => 100)

# Get the rest of the occurrences
while length(raccoon_occ) <= 2000
    try
        GBIF.occurrences!(raccoon_occ)
    catch#
        break
    end
end

# Get the temperature and precipitation
temperature   = worldclim(1)[(bottom = -60.0)]
precipitation = worldclim(12)[(bottom = -60.0)]
# Get landcover data
tree  = landcover(1)[(bottom = -60.0)]
urban = landcover(2)[(bottom = -60.0)]
water = landcover(3)[(bottom = -60.0)]
# Combine variables
vars = [temperature, precipitation, tree, urban, water]

# Map climate variables
temp_map = heatmap(temperature, c = :inferno, xlab= "Longitude", ylab = "Latitude",
                   colorbar_title = "Temperature (° C)")
prec_map = heatmap(precipitation, c = :blues, xlab= "Longitude", ylab = "Latitude",
                   colorbar_title = "Precipitation (mm)")

# Map raccoon occurrences
# Option 1: Regular markers & raccoon emoji on the side
occ_map  = heatmap(temperature, c = :lightgrey, xlab= "Longitude", ylab = "Latitude",
                   colorbar = :none) |> x ->
    scatter!(x, longitudes(raccoon_occ), latitudes(raccoon_occ),
         lab = "Raccoons",
         legend = :outerbottomright,
         foreground_color_legend = nothing,
         size = (600, 400)
    ) |> x ->
    plot!(x, img,
        yflip = true,
        inset = bbox(0.4, -0.05, 150px, 75px, :center),
        subplot = 2,
        grid = false, axis = false,
        bg_inside = nothing
    )
# Option 2: Raccoon emojis as markers (max 200 for now)
function emojimap(layer, lons, lats, n)
    lons = lons[1:n]
    lats = lats[1:n]
    emoji_map = heatmap(layer, c = :lightgrey, # xlab= "Longitude", ylab = "Latitude",
        colorbar = :none, size = (360,150).*2,
        ticks = false, margin = -1.9mm, box = :none)
    @time @showprogress for (lon, lat, i) in zip(lons, lats, 1:n)
        plot!(emoji_map, img, yflip = true, grid = false, axis = false, bg_inside = nothing,
            inset = (1, bbox(((lon + 180)/360)w - 15px, ((lat+60)/150)h - 8px, 30px, 15px, :bottom, :left)),
            subplot = i+1
        )
    end
    return emoji_map
end
occ_map2 = emojimap(temperature, longitudes(raccoon_occ), latitudes(raccoon_occ), 200)

# Option 3: Unique sites
using DataFrames

df = DataFrame(raccoon_occ)
select!(df, :longitude, :latitude)
filter!(x -> !ismissing(x.latitude) || !ismissing(x.longitude), df)
df.gridlon = [SimpleSDMLayers._match_longitude(temperature, lon) for lon in df.longitude]
df.gridlat = [SimpleSDMLayers._match_latitude(temperature, lat) for lat in df.latitude]
uniquedf = unique(df, [:gridlon, :gridlat])

occ_map2 = emojimap(temperature, uniquedf.longitude, uniquedf.latitude, 125)

## Bioclim model

# Get prediction for each variable
vars_predictions = [bioclim(v, raccoon_occ) for v in vars]

# Get minimum prediction per site
sdm_raccoon = reduce(min, vars_predictions)

# Set value to NaN if prediction is zero
replace!(x -> iszero(x) ? NaN : x, sdm_raccoon.grid)

# Filter predictions with threshold
threshold = quantile(filter(!isnan, sdm_raccoon.grid), 0.05)
replace!(x -> x <= threshold ? NaN : x, sdm_raccoon.grid)

# Group predictions in categories
lim1 = 0.08
lim2 = 0.25
replace!(x -> x >= lim2 ? 1.0 : x, sdm_raccoon.grid)
replace!(x -> lim1 < x < lim2 ? 0.5 : x, sdm_raccoon.grid)
replace!(x -> x <= lim1 ? 0.0 : x, sdm_raccoon.grid)

# Custom colorpicking function
colorpick(cg::ColorGradient, n::Int) = RGB[cg[i] for i in LinRange(0, 1, n)]

# Map predictions
begin
    pred_map = heatmap(temperature, c = :lightgrey, colorbar = :none,
                       xlab= "Longitude", ylab = "Latitude",
                       framestyle = :box, size = (600, 300))
    heatmap!(pred_map, sdm_raccoon, c = :viridis, clim = (0,1), colorbar_title = "Suitability for raccoons")
    scatter!(pred_map, [NaN NaN NaN NaN],
                 c = [colorpick(cgrad(:viridis), 3)... :lightgrey],
                 labels = ["Low" "Medium" "High" "Not suitable"],
                 legend = :outerbottomright, legendtitle = "Suitability",
                 foreground_color_legend = nothing,
                 legendtitlefontsize = 10)
    plot!(pred_map, img,
            yflip = true,
            inset = bbox(0.36, -0.23, 150px, 75px, :center),
            subplot = 2,
            grid = false, axis = false,
            bg_inside = nothing
        )
end

## Export figures
# Occurrence figure
savefig(plot(occ_map2, dpi = 600), joinpath("fig", "occurrences_emojis.png"))
savefig(plot(occ_map2, dpi = 600), joinpath("fig", "occurrences_emojis.svg"))
# Predictions figure
savefig(plot(pred_map, dpi = 600), joinpath("fig", "predictions.png"))
savefig(plot(pred_map, dpi = 600), joinpath("fig", "predictions.svg"))
savefig(plot(pred_map, dpi = 1200), joinpath("fig", "predictions_1200.png"))
# Extra figures
savefig(plot(temp_map, dpi = 150), joinpath("fig", "temperature.png"))
savefig(plot(prec_map, dpi = 150), joinpath("fig", "precipitation.png"))
savefig(plot(occ_map,  dpi = 150), joinpath("fig", "occurrences.png"))


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule libass_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("libass")
JLLWrappers.@generate_main_file("libass", UUID("0ac62f75-1d6f-5e53-bd7c-93b484bb37c0"))
end  # module libass_jll