text:
module test_bukdu_html5_form

using Test # @test
using Bukdu # Changeset
import Bukdu.HTML5.Form: @tags
import Bukdu.HTML5.Form: form_for, label_for
import Bukdu.HTML5.Form: text_input, text_area, radio_button, checkbox
import Bukdu.HTML5.Form: submit

@tags input button form label textarea

struct User
end

f = Changeset(User, (name="Alex",))

@test text_input(f, :name, "Chen")                  ==  input[:id => "user_name", :name => "user_name", :type => "text", :value => "Chen"]

@test submit("Submit")                              ==  button[:type => "submit"]("Submit")

@test form_for((f) -> nothing, f, "/post")          ==  form[:action => "/post", :method => "POST"]

@test label_for(radio_button(f, :choice, "Choice")) ==  [input[:id => "user_choice_Choice", :name => "user_choice", :type => "radio", :value => "Choice"],
                                                         label[:for => "user_choice_Choice"]]

@test checkbox(f, :famous, true)                    ==  [input[:name => "user_famous", :type => "hidden", :value => "false"],
                                                         input[:checked => "checked", :id => "user_famous", :name => "user_famous", :type => "checkbox", :value => "true"]]

@test checkbox(f, :famous, false)                   ==  [input[:name => "user_famous", :type => "hidden", :value => "false"],
                                                         input[:id => "user_famous", :name => "user_famous", :type => "checkbox", :value => "true"]]

@test text_area(f, :intro, nothing)                  ==  textarea[:id => "user_intro", :name => "user_intro"]("")

@tags a
@test string(a[:href => "/"]()) == """"""

end  # module test_bukdu_html5_form

max_stars_repo_name: wookay/TestBukdu

max_issues_repo_name: wookay/TestBukdu

max_forks_repo_name: wookay/TestBukdu


text:
# ---
# title: 219. Contains Duplicate II
# id: problem219
# author: zhwang
# date: 2022-01-19
# difficulty: Easy
# categories: Array, Hash Table
# link: <https://leetcode.com/problems/contains-duplicate-ii/description/>
# hidden: true
# ---
#
# Given an array of integers and an integer _k_ , find out whether there are two
# distinct indices _i_ and _j_ in the array such that **nums[i] = nums[j]** and
# the **absolute** difference between _i_ and _j_ is at most _k_.
#
# **Example 1:**
#
#
#
#     Input: nums = [1,2,3,1], k = 3
#     Output: true
#
#
# **Example 2:**
#
#
#
#     Input: nums = [1,0,1,1], k = 1
#     Output: true
#
#
# **Example 3:**
#
#
#
#     Input: nums = [1,2,3,1,2,3], k = 2
#     Output: false
#
#
#
## @lc code=start
using LeetCode

function contains_nearby_duplicate(nums::Vector{Int}, k::Int)::Bool
    pos = Dict{Int,Int}()
    for (i, num) in enumerate(nums)
        haskey(pos, num) && i - pos[num] <= k && return true
        pos[num] = i
    end
    return false
end

## @lc code=end

ext: jl


text:
export tridiag_real
using LinearAlgebra

"""
    tridiag_real(A::AbstractMatrix, d::Integer)

The symmetric real or Hermitian input matrix `A` is transformed to a real
`SymTridiagonal` matrix.
For general matrices only `d` superdiagonals are processed.
Works efficiently for `BandedMatrices`.
"""
function tridiag_real(A::AbstractMatrix, d::Integer)
    m, n = size(A)
    m == n || throw(DimensionMismatch("matrix is not square: dimensions are $((m, n))"))
    0 < d <= n || throw(ArgumentError("number of diagonals $d not in 1:$n"))

    B = copybands(A, d)
    td = make_tridiag!(B)
    dv = [real(B[i,1]) for i in 1:n]
    ev = d > 1 ? [abs(B[i,2]) for i = 2:n] : zeros(real(eltype(A)), n-1)
    SymTridiagonal(dv, ev)
end

# B[1:q,1:n] has the bands of a (2q-1)-banded hermitian nxn-matrix.
# superdiagonal i valid in [i+1:n].
function make_tridiag!(B)
    n, d = size(B)
    0 < d <= n || throw(ArgumentError("number of diagonals $d not in 1:$n"))

    @inbounds for bm = d-1:-1:2
        for k = 1:n-bm
            kp = k
            apiv = B[bm+kp,bm+1]
            iszero(apiv) && continue
            for i = bm+k-1:bm:n-1
                b = B[i, i-kp+1]
                c, s, r, α = givens2(b, apiv)
                u, v = B[i,1], B[i+1,1]
                upx = (u + v) / 2
                B[i,1] = (u - v) / 2
                B[i,i-kp+1] = r
                for j = kp+1:i
                    u = B[i,i-j+1]
                    v = B[i+1,i-j+2] * α'
                    B[i,i-j+1], B[i+1,i-j+2] = u * c + v * s, -u * s + v * c
                end
                B[i+1,1] = -(B[i,1] * α)'
                ip = i + bm
                for j = i+1:min(ip, n)
                    u = B[j,j-i+1]
                    v = B[j,j-i] * α
                    B[j,j-i+1], B[j,j-i] = u * c + v * s, -u * s + v * c
                end
                w = real(B[i+1,1])
                B[i,1] = upx - w
                B[i+1,1] = upx + w
                if ip < n
                    v = B[ip+1,ip-i+1] * α
                    apiv, B[ip+1,ip-i+1] = v * s, v * c
                end
                kp = i
            end
        end
    end
    B
end

function givens2(a::T, b::T) where T <: AbstractFloat
    r = _hypot(a, b)
    s = b / r
    c = a / r
    c, s, r, one(T)
end
@inline function givens2(a::T, b::T) where T <: Complex
    aa = _abs(a)
    ba = _abs(b)
    ra = _hypot(aa, ba)
    s = ba / ra
    c = aa / ra
    ua = a / aa
    α = ua' * (b / ba)
    r = ua * ra
    c, s, r, α
end

# speedy versions of abs and hypot (in the common case)
function _abs(x::Complex)
    s = sqrt(abs2(real(x)) + abs2(imag(x)))
    _isclean(s) ? s : abs(x)
end
_abs(x) = abs(x)
function _hypot(x, y)
    s = sqrt(abs2(x) + abs2(y))
    _isclean(s) ? s : hypot(x, y)
end

function _isclean(a::T) where T<:AbstractFloat
    isfinite(a) &&
    a >= sqrt(2 * floatmin(T))
end

copybands(A, d) = copybands!(zeros(eltype(A), size(A,1), d), A)

function copybands!(B::AbstractMatrix, A::AbstractMatrix)
    n, d = size(B)
    for i = 1:d
        for j = i:n
            B[j,i] = A[j-i+1,j]
        end
    end
    B
end


text:
module spin

export SpinBasis, sigmax, sigmay, sigmaz, sigmap, sigmam, spinup, spindown

import Base: ==

using Compat
using ..bases, ..states, ..operators, ..operators_sparse

"""
    SpinBasis(n)

Basis for spin-n particles.

The basis can be created for arbitrary spinnumbers by using a rational number,
e.g. `SpinBasis(3//2)`. The Pauli operators are defined for all possible
spin numbers.
"""
mutable struct SpinBasis <: Basis
    shape::Vector{Int}
    spinnumber::Rational{Int}
    function SpinBasis(spinnumber::Rational{Int})
        n = numerator(spinnumber)
        d = denominator(spinnumber)
        @assert d==2 || d==1
        @assert n > 0
        N = numerator(spinnumber*2 + 1)
        new([N], spinnumber)
    end
end
SpinBasis(spinnumber::Int) = SpinBasis(convert(Rational{Int}, spinnumber))

==(b1::SpinBasis, b2::SpinBasis) = b1.spinnumber==b2.spinnumber

"""
    sigmax(b::SpinBasis)

Pauli ``σ_x`` operator for the given Spin basis.
"""
function sigmax(b::SpinBasis)
    N = length(b)
    diag = Complex128[complex(sqrt(real((b.spinnumber + 1)*2*a - a*(a+1)))) for a=1:N-1]
    data = spdiagm(diag, 1, N, N) + spdiagm(diag, -1, N, N)
    SparseOperator(b, data)
end

"""
    sigmay(b::SpinBasis)

Pauli ``σ_y`` operator for the given Spin basis.
"""
function sigmay(b::SpinBasis)
    N = length(b)
    diag = Complex128[1im*complex(sqrt(real((b.spinnumber + 1)*2*a - a*(a+1)))) for a=1:N-1]
    data = spdiagm(diag, -1, N, N) - spdiagm(diag, 1, N, N)
    SparseOperator(b, data)
end

"""
    sigmaz(b::SpinBasis)

Pauli ``σ_z`` operator for the given Spin basis.
"""
function sigmaz(b::SpinBasis)
    N = length(b)
    diag = Complex128[complex(2*m) for m=b.spinnumber:-1:-b.spinnumber]
    data = spdiagm(diag, 0, N, N)
    SparseOperator(b, data)
end

"""
    sigmap(b::SpinBasis)

Raising operator ``σ_+`` for the given Spin basis.
"""
function sigmap(b::SpinBasis)
    N = length(b)
    S = (b.spinnumber + 1)*b.spinnumber
    diag = Complex128[complex(sqrt(float(S - m*(m+1)))) for m=b.spinnumber-1:-1:-b.spinnumber]
    data = spdiagm(diag, 1, N, N)
    SparseOperator(b, data)
end

"""
    sigmam(b::SpinBasis)

Lowering operator ``σ_-`` for the given Spin basis.
"""
function sigmam(b::SpinBasis)
    N = length(b)
    S = (b.spinnumber + 1)*b.spinnumber
    diag = [complex(sqrt(float(S - m*(m-1)))) for m=b.spinnumber:-1:-b.spinnumber+1]
    data = spdiagm(diag, -1, N, N)
    SparseOperator(b, data)
end

"""
    spinup(b::SpinBasis)

Spin up state for the given Spin basis.
"""
spinup(b::SpinBasis) = basisstate(b, 1)

"""
    spindown(b::SpinBasis)

Spin down state for the given Spin basis.
"""
spindown(b::SpinBasis) = basisstate(b, b.shape[1])

end #module


text:
# --------------------------------------------------------------------------
# ACE1.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# Licensed under ASL - see ASL.md for terms and conditions.
# --------------------------------------------------------------------------

# ----------------------------------------------------

export PolyPairPot

struct PolyPairPot{T,TJ,NZ} <: PairPotential
   coeffs::Vector{T}
   basis::PolyPairBasis{TJ, NZ}
end

@pot PolyPairPot

PolyPairPot(pB::PolyPairBasis, coeffs::Vector) = PolyPairPot(coeffs, pB)

JuLIP.MLIPs.combine(pB::PolyPairBasis, coeffs::AbstractVector) =
            PolyPairPot(identity.(collect(coeffs)), pB)

JuLIP.cutoff(V::PolyPairPot) = cutoff(V.basis)

z2i(V::PolyPairPot, z::AtomicNumber) = z2i(V.basis, z)
i2z(V::PolyPairPot, iz::Integer) = i2z(V.basis, iz)
numz(V::PolyPairPot) = numz(V.basis)

==(V1::PolyPairPot, V2::PolyPairPot) =
            ( (V1.basis == V2.basis) && (V1.coeffs == V2.coeffs) )

write_dict(V::PolyPairPot{T}) where {T} = Dict(
      "__id__" => "ACE1_PolyPairPot",
      "T" => write_dict(T),
      "coeffs" => V.coeffs,
      "basis" => write_dict(V.basis)
      )

read_dict(::Val{:ACE1_PolyPairPot}, D::Dict, T = read_dict(D["T"])) =
      PolyPairPot(read_dict(D["basis"]), T.(D["coeffs"]))

alloc_temp(V::PolyPairPot{T}, N::Integer) where {T} =
      ( R = zeros(JVec{T}, N),
        Z = zeros(AtomicNumber, N),
        alloc_temp(V.basis)... )

alloc_temp_d(V::PolyPairPot{T}, N::Integer) where {T} =
      ( dV = zeros(JVec{T}, N),
         R = zeros(JVec{T}, N),
         Z = zeros(AtomicNumber, N),
         alloc_temp_d(V.basis)... )

function _dot_zij(V, B, z, z0)
   i0 = _Bidx0(V.basis, z, z0)  # cf. pair_basis.jl
   return sum( V.coeffs[i0 + n] * B[n]  for n = 1:length(V.basis, z0) )
end

evaluate!(tmp, V::PolyPairPot, r::Number, z, z0) =
      _dot_zij(V, evaluate!(tmp.J, tmp.tmp_J, V.basis.J, r, z, z0), z, z0)

evaluate_d!(tmp, V::PolyPairPot, r::Number, z, z0) =
      _dot_zij(V, evaluate_d!(tmp.J, tmp.dJ, tmp.tmpd_J, V.basis.J, r, z, z0), z, z0)

function evaluate!(tmp, V::PolyPairPot, r::Number)
   @assert numz(V) == 1
   z = i2z(V, 1)
   return evaluate!(tmp, V, r, z, z)
end

function evaluate_d!(tmp, V::PolyPairPot, r::Number)
   @assert numz(V) == 1
   z = i2z(V, 1)
   return evaluate_d!(tmp, V, r, z, z)
end

evaluate(V::PolyPairPot, r::Number, args...) = evaluate!(alloc_temp(V, 1), V, r, args...)
evaluate_d(V::PolyPairPot, r::Number, args...) = evaluate_d!(alloc_temp_d(V, 1), V, r, args...)


text:
using JuMP
using Ipopt
using CSV, DataFrames, DataStructures, MathOptFormat
using Statistics,Distributions, StatsFuns, StatsBase
using PyPlot

function critical_fractile(dist_type)
    #closed form of solution for newsvendor for normal distribution
    # dist_type - str, "norm" or "exp" for normal or exponential distribution, respectively

    b = 1 #underage cost
    h = 1 #overage cost

    if dist_type == "norm"
        mean = 50
        var = 50

        beta = b/(b+h)

        b_inv = norminvcdf(beta) #inverse of the cdf for normal distribution

        return mean + b_inv * var

    elseif dist_type == "exp"
        # plugged and chugged with cdf of exponential function to get warm start value (solver kept converging to point of infeasibility)
        cdf_func(x) = cdf(Exponential(50),x)

        # print(cdf_func(34.65))

        model3 = Model(with_optimizer(Ipopt.Optimizer,print_level=0))

        @variable(model3, x_var, start=34) #use warm start to avoid converge to infeasible point (see above)

        JuMP.register(model3, :cdf_func, 1, cdf_func; autodiff = true)  #to use cdf function in constraint

        @NLconstraint(model3,con1, cdf_func(x_var) >= b/(b+h))

        @objective(model3,Min,x_var)
        optimize!(model3)

        return value(x_var)

    else
        return nothing
    end
end

function scarf(df)
    #closed form scarf solution
    # df - dataframe, data of observed demand

    b = 1 #underage cost
    h = 1 #overage cost

    mu_hat = mean(df[:,"Column1"]) #sample mean
    sigma_sq_hat = var(df[:,"Column1"]) #sample variance

    return mu_hat + sigma_sq_hat/2 * (sqrt(b / h) - sqrt(h / b))
end

function SAA(df)
    # compute SAA solution (empirical risk minimization)
    # df - dataframe, data of observed demand

    b = 1 #underage cost
    h = 1 #overage cost

    model = Model(with_optimizer(Ipopt.Optimizer,print_level=0))

    data_size = length(df[:,"Column1"])

    @variable(model, x_var) #stocking quantity
    @variable(model, y_vars[1:data_size]) #variable to linearize max(d_i - x, x - d_i) in objective function

    @constraint(model,con1[i = 1:data_size], y_vars[i] >= x_var - df[:,"Column1"][i])
    @constraint(model, con2[i = 1:data_size], y_vars[i] >= df[:,"Column1"][i] - x_var )

    @NLobjective(model,Min, 1 / data_size * sum(y_vars[i] for i = 1:data_size))

    optimize!(model)

    return value(x_var)
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   SGP4 orbit propagator model.
#
#   This is a independent implementation of the algorithm presented in [1].
#   Notice that the readability of the code was the major concern about the
#   implementation here. Algorithms with better performance can be found at
#   Vallado's repository in:
#
#       https://celestrak.com/software/vallado-sw.asp
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] Hoots, F. R., Roehrich, R. L (1980). Models for Propagation of NORAD
#       Elements Set. Spacetrack Report No. 3.
#
#   [2] Vallado, D. A., Crawford, P., Hujsak, R., Kelso, T. S (2006). Revisiting
#       Spacetrack Report #3: Rev1. AIAA.
#
#   [3] SGP4 Source code of STRF: https://github.com/cbassa/strf
#       The SGP4 C code available on STRF was converted by Paul. S. Crawford and
#       Andrew R. Brooks.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export sgp4_gc_wgs72, sgp4_gc_wgs84
export sgp4_init, sgp4!
export dsinit, dsper!, dssec!

################################################################################
#                                  Overloads
################################################################################

# Copy for SGP4_Structure.
Base.copy(m::SGP4_Structure) = SGP4_Structure([ getfield(m, k) for k = 1:length(fieldnames(m)) ]...)

# Deepcopy for SGP4_Structure.
Base.deepcopy(m::SGP4_Structure) = SGP4_Structure([ deepcopy(getfield(m, k)) for k = 1:length(fieldnames(m)) ]...)

################################################################################
#                                  Constants
################################################################################

# WGS-84 / EGM-08 Gravitational constants.
const sgp4_gc_wgs84 = SGP4_GravCte{Float64}(
        R0/1000,
        60.0/sqrt(6378.137^3/398600.5),
         0.00108262998905,
        -0.00000253215306,
        -0.00000161098761
       )

# WGS-72 Gravitational constants.
const sgp4_gc_wgs72 = SGP4_GravCte{Float64}(
        6378.135,
        60.0/sqrt(6378.135^3/398600.8),
         0.001082616,
        -0.00000253881,
        -0.00000165597
       )

# Constants realted to the Deep Space perturbations.
# ==================================================

const STEP          = 720.0
const MAX_INTEGRATE = STEP * 10000
const ZNS           = 1.19459E-5
const C1SS          = 2.9864797e-6
const ZES           = 0.01675
const ZNL           = 1.5835218e-4
const ZEL           = 0.05490
const C1L           = 4.7968065e-7
const ZSINIS        = 0.39785416
const ZCOSIS        = 0.91744867
const ZCOSGS        = 0.1945905
const ZSINGS        = -0.98088458
const Q22           = 1.7891679e-6
const Q31           = 2.1460748e-6
const Q33           = 2.2123015e-7
const G22           = 5.7686396
const G32           = 0.95240898
const G44           = 1.8014998
const G52           = 1.0508330
const G54           = 4.4108898
const ROOT22        = 1.7891679e-6
const ROOT32        = 3.7393792e-7
const ROOT44        = 7.3636953e-9
const ROOT52        = 1.1428639e-7
const ROOT54        = 2.1765803e-9
const THDT          = 4.37526908801129966e-3

################################################################################
#                                  Functions
################################################################################

"""
    sgp4_init(spg4_gc::SGP4_GravCte{T}, epoch::Number, n_0::Number, e_0::Number, i_0::Number, Ω_0::Number, ω_0::Number, M_0::Number, bstar::Number) where T
    sgp4_init(tle::TLE, sgp4_gc::SGP4_GravCte = sgp4_gc_wgs84)

Initialize the data structure of SGP4 algorithm.

# Args

* `spg4_gc`: SPG4 gravitational constants (see `SGP4_GravCte`).
* `epoch`: Epoch of the orbital elements [Julian Day].
* `n_0`: SGP type "mean" mean motion at epoch [rad/min].
* `e_0`: "Mean" eccentricity at epoch.
* `i_0`: "Mean" inclination at epoch [rad].
* `Ω_0`: "Mean" longitude of the ascending node at epoch [rad].
* `ω_0`: "Mean" argument of perigee at epoch [rad].
* `M_0`: "Mean" mean anomaly at epoch [rad].
* `bstar`: Drag parameter (B*).
* `tle`: TLE to initialize the SPG4 (see `TLE`).

# Returns

The structure `SGP4_Structure` with the initialized parameters.

"""
function sgp4_init(tle::TLE, sgp4_gc::SGP4_GravCte = sgp4_gc_wgs84)
    return sgp4_init(sgp4_gc, tle.epoch, tle.n*2*pi/(24*60), tle.e,
                     tle.i*pi/180, tle.Ω*pi/180, tle.ω*pi/180, tle.M*pi/180,
                     tle.bstar)
end

function sgp4_init(sgp4_gc::SGP4_GravCte{T},
                   epoch::Number,
                   n_0::Number,
                   e_0::Number,
                   i_0::Number,
                   Ω_0::Number,
                   ω_0::Number,
                   M_0::Number,
                   bstar::Number) where T

    # Unpack the gravitational constants to improve code readability.
    @unpack_SGP4_GravCte sgp4_gc

    # Constants
    # =========
    #
    # Note: [er] = Earth radii.

    # Distance units / Earth radii.
    AE = T(1)

    k_2  = +1/2*J2*AE^2
    k_4  = -3/8*J4*AE^4
    A_30 = -J3*AE^3

    # Kilometers / Earth radii.
    XKMPER = R0

    # Parameters for the SGP4 density function.
    s   =  78/XKMPER + 1
    q_0 = 120/XKMPER + 1

    # (q_0-s)^4 [er]^4
    QOMS2T = (q_0-s)^4

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

    # Auxiliary variables to improve the performance.
    # ===============================================

    e_0² = e_0^2

    sin_i_0, θ = sincos(i_0)
    θ²         = θ^2
    θ³         = θ^3
    θ⁴         = θ^4

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

    # Recover the original mean motion (nll_0) and semi-major axis (all_0) from
    # the input elements.

    aux = (3*θ²-1)/(1-e_0²)^(3/2)

    a_1 = (XKE/n_0)^(2/3)
    δ_1 = 3/2*k_2/a_1^2*aux
    a_0 = a_1*(1 - 1/3*δ_1 - δ_1^2 - 134/81*δ_1^3)
    δ_0 = 3/2*k_2/a_0^2*aux

    nll_0 = n_0/(1 + δ_0)

    # Vallado's implementation of SGP4 [2] compute the semi-major axis
    # considering the new angular velocity, which is called `no_unkozai`. In the
    # original SGP4 technical report [1], the semi-major axis was computed
    # considering:
    #
    #   all_0 = a_0/(1 - δ_0)
    #
    all_0 = (XKE/nll_0)^(2/3)

    # Initialization
    # ==============

    # Compute the orbit perigee [ER].
    perigee = (all_0*(1-e_0)-AE)*XKMPER

    # For perigee below 156 km, the values of S and QOMS2T are altered.
    if perigee < 156
        if perigee < 98
            s = 20/XKMPER + AE
        # Perigee between 98km and 156km.
        else
            s = all_0*(1-e_0) - s + AE
        end

        QOMS2T = (q_0 - s)^4
    end

    # Compute SGP4 constants.
    ξ    = 1/(all_0-s)
    β_0  = sqrt(1-e_0^2)
    η    = all_0*e_0*ξ

    # Vallado's implementation of SGP4 [2] considers the absolute value of
    # (1-η^2) here and in the C2 and C4 computation. Notice that, if (1-η^2) <
    # 0, then aux1 cannot be computed. The original SGP4 technical report [1]
    # does not mention anything about this.

    aux0 = abs(1-η^2)
    aux1 = aux0^(-7/2)
    aux2 = ξ^4*all_0*β_0^2*aux1

    C2 = QOMS2T*ξ^4*nll_0*aux1*
         ( all_0*( 1 + (3/2)η^2 + 4e_0*η + e_0*η^3) +
           3/2*(k_2*ξ)/aux0*(-1/2 + (3/2)θ²)*(8 + 24η^2 + 3η^4) )

    C1 = bstar*C2

    C3 = (e_0 > 1e-4) ? QOMS2T*ξ^5*A_30*nll_0*AE*sin_i_0/(k_2*e_0) : T(0)

    C4 = 2nll_0*QOMS2T*aux2*
         ( 2η*(1+e_0*η) + (1/2)e_0 + (1/2)η^3 -
           2k_2*ξ/(all_0*aux0)*
                (3(1-3θ²)*(1 + (3/2)η^2 - 2e_0*η - (1/2)e_0*η^3) +
                3/4*(1-θ²)*(2η^2 - e_0*η - e_0*η^3)*cos(2ω_0)))

    C5 = 2QOMS2T*aux2*(1 + (11/4)η*(η+e_0) + e_0*η^3)

    D2 = 4all_0*ξ*C1^2

    D3 = 4/3*all_0*ξ^2*( 17all_0 +   s)*C1^3

    # Vallado's implementation of SGP4 [2] uses all_0^2, instead of only all_0
    # that is seen in the original SGP4 Technical Report [1].
    D4 = 2/3*all_0^2*ξ^3*(221all_0 + 31s)*C1^4

    # Compute the time-derivative of some orbital elements.
    dotM = ( 1 + 3k_2*(-1 + 3θ²)/(2all_0^2*β_0^3) +
             3k_2^2*(13 - 78θ² + 137θ⁴)/(16all_0^4*β_0^7) )*nll_0

    dotω = ( -3k_2*(1-5θ²)/(2all_0^2*β_0^4) +
              3k_2^2*(7 - 114θ² + 395θ⁴)/(16all_0^4*β_0^8) +
              5k_4*(3 - 36θ² + 49θ⁴)/(4all_0^4*β_0^8) )*nll_0

    dotΩ1 = -3k_2*θ/(all_0^2*β_0^4)*nll_0

    dotΩ  = dotΩ1 + ( 3k_2^2*(4θ - 19θ³)/(2all_0^4*β_0^8) +
                        5k_4*θ*(3 - 7θ²)/(2all_0^4*β_0^8) )*nll_0

    # The current orbital parameters are obtained from the TLE.
    a_k = all_0
    e_k = e_0
    i_k = i_0
    Ω_k = Ω_0
    ω_k = ω_0
    M_k = M_0
    n_k = nll_0

    sgp4_ds::SGP4_DeepSpace{T} = SGP4_DeepSpace{T}()

    # If the orbit period is higher than 225 min., then we must consider the
    # deep space perturbations. This is indicated by selecting the algorithm
    # `:sdp4`.
    if 2π / n_0 >= 225.0
        algorithm = :sdp4

        # Initialize the values for the SDP4 (deep space) algorithm.
        sgp4_ds = dsinit(epoch,
                         nll_0, all_0, e_0, i_0, Ω_0, ω_0, M_0,
                         dotM, dotω, dotΩ)
    else
        # For perigee lower than 220 km, the equations are truncated to a linear
        # variation in `sqrt(a)` and quadratic variation in mean anomaly. Also,
        # the C5 term, the δω term, and the δM term are dropped. This is
        # indicated by selecting the algorithm `:sgp4_lowper`. Otherwise, if
        # perigee is higher or equal 220 km and the orbit period is lower than
        # 225 min., then we use the normal SGP4 algorithm by selecting `:sgp4`.
        algorithm = (perigee/AE >= (220 + (AE-1)*XKMPER)) ? :sgp4 : :sgp4_lowper
    end

    # Create the output structure with the data.
    SGP4_Structure{T}(

        epoch, n_0, e_0, i_0, Ω_0, ω_0, M_0, bstar, 0, a_k, e_k, i_k, Ω_k, ω_k,
        M_k, n_k, all_0, nll_0, AE, QOMS2T, β_0, ξ, η, sin_i_0, θ, θ², A_30,
        k_2, k_4, C1, C3, C4, C5, D2, D3, D4, dotM, dotω, dotΩ1, dotΩ,
        algorithm, sgp4_gc, sgp4_ds

       )
end

"""
    sgp4!(sgp4d::SGP4_Structure{T}, t::Number) where T

Propagate the orbit defined in `sgp4d` (see `SGP4_Structure`) until the time `t`
[min]. Notice that the values in `sgp4d` will be modified.

# Returns

* The position vector represented in TEME frame at time `t` [km].
* The velocity vector represented in TEME frame at time `t` [km/s].

"""
function sgp4!(sgp4d::SGP4_Structure{T}, t::Number) where T
    # Unpack variables.
    @unpack_SGP4_Structure sgp4d
    @unpack_SGP4_GravCte   sgp4_gc

    # After unpacking sgp4d, we have two sets of orbit elements:
    #
    #   (n_0, e_0, i_0, Ω_0, ω_0, M_0),
    #
    # and
    #
    #   (n_k, e_k, i_k, Ω_k, ω_k, M_k).
    #
    # The first are those initial elements from the orbit defined in `sgp4_init`
    # function. The second are the current elements. During this functions, the
    # second set is updated by adding the many effects considered in SGP4.

    # Time elapsed since epoch.
    #
    # We convert to `T` to avoid numerical problems with very big numbers as
    # pointed out in:
    #
    #   https://github.com/JuliaLang/julia/issues/27355
    Δt = T(t)

    # Initialization of the current elements with the values of the epoch.
    n_k = nll_0
    a_k = all_0
    e_k = e_0
    i_k = i_0
    Ω_k = Ω_0
    ω_k = ω_0
    M_k = M_0

    # Auxiliary variables to improve code performance.
    sin_i_k = sin_i_0

    # Secular effects of atmospheric drag and gravitation.
    # ====================================================

    M_k = M_0 + dotM*Δt
    Ω_k = Ω_0 + dotΩ*Δt - 21/2*(nll_0*k_2*θ)/(all_0^2*β_0^2)*C1*Δt^2
    ω_k = ω_0 + dotω*Δt

    # Check if we need to use SDP4 (deep space) algorithm.
    if algorithm == :sdp4
        # Compute the elements perturbed by the secular effects.
        (n_k, e_k, i_k, Ω_k, ω_k, M_k) = dssec!(sgp4_ds,
                                                nll_0,
                                                e_0,
                                                i_0,
                                                ω_0,
                                                Ω_k,
                                                ω_k,
                                                M_k,
                                                dotω,
                                                Δt)

        a_k  = (XKE / n_k)^(2/3) * ( 1 - C1*Δt )^2
        e_k += -bstar*C4*Δt
        M_k += nll_0*(1.5C1*Δt^2)

    # Check if perigee is above 220 km.
    elseif algorithm == :sgp4

        sin_M_0, cos_M_0 = sincos(M_0)

        δω  = bstar*C3*cos(ω_0)*Δt

        δM  = (e_0 > 1e-4) ?
            -2/3*QOMS2T*bstar*ξ^4*AE/(e_0*η)*( (1 + η*cos(M_k))^3 -
                                               (1 + η*cos_M_0 )^3 ) : T(0)

        M_k += +δω + δM

        ω_k += -δω - δM

        e_k = e_0 - bstar*C4*Δt - bstar*C5*(sin(M_k) - sin_M_0)

        a_k = all_0*( @evalpoly(Δt, 1, -C1, -D2, -D3, -D4) )^2

        IL  = M_k + ω_k + Ω_k + nll_0*
              ( @evalpoly(Δt, 0, 0, 3/2*C1, +(D2 + 2C1^2),
                          +(1/4)*(3D3 + 12C1*D2 + 10C1^3),
                          +(1/5)*(3D4 + 12C1*D3 + 6D2^2 + 30C1^2*D2 + 15C1^4) ) )

    elseif algorithm == :sgp4_lowper
        # If so, then
        #     1. Drop all terms after C1 in `a` and `IL`.
        #     2. Drop all terms involving C5.
        #     3. Drop δω.
        #     4. Drop δM.
        e_k = e_0 - bstar*C4*Δt

        a_k = all_0*(1 - C1*Δt)^2

        IL  = M_k + ω_k + Ω_k + nll_0*(3/2)*C1*Δt^2
    else
        error("Unknown algorithm :$algorithm. Possible values are :sgp4, :sgp4_lowper, :sdp4.")
    end

    # TODO: Vallado's implementation [2] apply this normalization to the mean
    # anomaly. It is necessary to verify the reason for that.
    M_k_aux = M_k + ω_k + Ω_k
    Ω_k     = rem(Ω_k, 2π)
    ω_k     = rem(ω_k, 2π)
    M_k_aux = rem(M_k_aux, 2π)
    M_k     = rem(M_k_aux - ω_k - Ω_k, 2π)

    # Lunar-Solar Periodics for Deep Space Orbits
    # ===========================================

    # This is only necessary if we are using SDP4 algorithm.
    if algorithm == :sdp4
        # Compute the elements perturbed by the Lunar-Solar periodics.
        (e_k, i_k, Ω_k, ω_k, M_k) = dsper!(sgp4_ds,
                                           e_k,
                                           i_k,
                                           Ω_k,
                                           ω_k,
                                           M_k,
                                           Δt)

        IL = M_k + ω_k + Ω_k

        # Make sure that the inclination is always positive.
        if i_k < 0
            i_k = -i_k
            Ω_k += π
            ω_k -= π
        end

        # The inclination was changed, hence some auxiliary variables must be
        # recomputed.
        sin_i_k, θ = sincos(i_k)
        θ²         = θ^2
    end

    # Vallado's code does not let the eccentricity to be smaller than 1e-6.
    #
    # TODO: Verify why this is necessary. I did not find any reason for that.
    (e_k < 1e-6) && (e_k = T(1e-6))

    β = sqrt(1-e_k^2)

    # Compute the angular velocity [rad/min].
    n_k = XKE/a_k^(3/2)

    # Long-period periodic terms.
    # ===========================

    sin_ω_k, cos_ω_k = sincos(ω_k)

    a_xN = e_k*cos_ω_k

    # TODO: Vallado's implementation of SGP4 uses another equation here.
    # However, both produces the same result. Verify which one is better.
    #
    a_yNL = A_30*sin_i_k/(4*k_2*a_k*β^2)
    a_yN = e_k*sin_ω_k + a_yNL

    IL_L =  (1/2)a_yNL*a_xN*(3 + 5θ)/(1 + θ)

    IL_T = IL + IL_L

    # Solve Kepler's equation for (E + ω).
    # ====================================

    U = rem(IL_T - Ω_k, 2π)

    E_ω = U

    # Define the following variables that will be modified inside the loop so
    # that we can use them after the loop.
    sin_E_ω = T(0)
    cos_E_ω = T(0)

    for k = 1:10
        sin_E_ω, cos_E_ω = sincos(E_ω)

        ΔE_ω = (U - a_yN*cos_E_ω + a_xN*sin_E_ω - E_ω)/
               (1 - a_yN*sin_E_ω - a_xN*cos_E_ω)

        # Vallado proposes to limit the maximum increment.
        (abs(ΔE_ω) >= 0.95) && (ΔE_ω = sign(ΔE_ω)*0.95)

        E_ω += ΔE_ω

        # If the increment is less than a threshold, break the loop.
        #
        # Vallado proposes a threshold of 10^-12 instead of 10^-6.
        (abs(ΔE_ω) < 1e-12) && (break)
    end

    # Short-term periodic terms.
    # ==========================

    # Auxiliary variables.
    #
    # Note: the sine and cosine of E+ω was already computed in the previous
    # loop.

    e_cos_E = a_xN*cos_E_ω + a_yN*sin_E_ω
    e_sin_E = a_xN*sin_E_ω - a_yN*cos_E_ω

    e_L     = sqrt(a_xN^2 + a_yN^2)

    p_L     = a_k*(1-e_L^2)

    r       = a_k*(1-e_cos_E)

    dot_r   = XKE*sqrt(a_k)*e_sin_E/r

    r_dot_f = XKE*sqrt(p_L)/r

    aux     = e_sin_E/(1 + sqrt(1 - e_L^2))
    cos_u   = a_k/r*( cos_E_ω - a_xN + a_yN*aux )
    sin_u   = a_k/r*( sin_E_ω - a_yN - a_xN*aux )
    cos_2u  = 1 - 2*sin_u^2
    sin_2u  = 2cos_u*sin_u
    u       = atan(sin_u, cos_u)

    # Short-term periodic terms.

    Δr       = +k_2/(2p_L)*(1 - θ²)*cos_2u

    Δu       = -k_2/(4p_L^2)*(7θ² - 1)*sin_2u

    ΔΩ       = +3k_2*θ/(2p_L^2)*sin_2u

    Δi       = +3k_2*θ/(2p_L^2)*sin_i_k*cos_2u

    Δdot_r   = -k_2*n_k/p_L*(1 - θ²)*sin_2u

    Δr_dot_f = +k_2*n_k/p_L*( (1 - θ²)*cos_2u - 3/2*(1 - 3θ²) )

    # The short-term periodics are added to give the osculating quantities.

    r_k       = r*( 1 - (3/2)k_2*sqrt(1 - e_L^2)/p_L^2*(3θ² - 1) ) + Δr

    u_k       = u + Δu

    Ω_k       = Ω_k + ΔΩ

    i_k       = i_k + Δi

    dot_r_k   = dot_r + Δdot_r

    r_dot_f_k = r_dot_f + Δr_dot_f

    # Orientation vectors.
    sin_Ω_k, cos_Ω_k = sincos(Ω_k)
    sin_i_k, cos_i_k = sincos(i_k)
    sin_u_k, cos_u_k = sincos(u_k)

    M = SVector{3}(-sin_Ω_k*cos_i_k,
                   +cos_Ω_k*cos_i_k,
                            sin_i_k )

    N = SVector{3}(+cos_Ω_k,
                   +sin_Ω_k,
                   +T(0)   )

    Uv = M*sin_u_k + N*cos_u_k
    Vv = M*cos_u_k - N*sin_u_k

    r_TEME = r_k*Uv*R0
    v_TEME = (dot_r_k*Uv + r_dot_f_k*Vv)*R0/60

    # Update the variables.
    sgp4d.Δt  = Δt
    sgp4d.a_k = a_k
    sgp4d.e_k = e_k
    sgp4d.i_k = i_k
    sgp4d.Ω_k = Ω_k
    sgp4d.ω_k = ω_k
    sgp4d.M_k = M_k
    sgp4d.n_k = n_k

    (r_TEME, v_TEME)
end

################################################################################
#                             Deep Space Functions
################################################################################

"""
    dsinit(epoch::T, nll_0::T, all_0::T, e_0::T, i_0::T, Ω_0::T, ω_0::T, M_0::T, dotM::T, dotω::T, dotΩ::T) where T<:Number

Initialize the deep space structure. This function performs the initial
computations and save the values at an instance of the structure
`SGP4_DeepSpace`. Those will be used when calling the functions `dsper!` and
`dpsec!`.

# Args

* `epoch`: Epoch of the initial orbit [Julian Day].
* `nll_0`: Initial mean motion [rad/min].
* `all_0`: Initial semi-major axis [ER].
* `e_0`: Initial eccentricity.
* `i_0`: Initial inclination [rad].
* `Ω_0`: Initial right ascencion of the ascending node [rad].
* `ω_0`: Initial argument of perigee [rad].
* `M_0`: Initial mean motion [rad].
* `dotM`: Time-derivative of the mean motion [rad/min].
* `dotω`: Time-derivative of the argument of perigee [rad/min].
* `dotΩ`: Time-derivative of the RAAN [rad/min].

# Returns

An instance of the structure `SGP4_DeepSpace` with the initalized values.

"""
function dsinit(epoch::T,
                nll_0::T,
                all_0::T,
                e_0::T,
                i_0::T,
                Ω_0::T,
                ω_0::T,
                M_0::T,
                dotM::T,
                dotω::T,
                dotΩ::T) where T<:Number

    sgp4_ds::SGP4_DeepSpace{T} = SGP4_DeepSpace{T}()
    @unpack_SGP4_DeepSpace sgp4_ds

    #                         Auxiliary Variables
    # ==========================================================================
    e_0²        = e_0^2
    e_0³        = e_0 * e_0²
    sqrt_1_e_0² = sqrt((1 - e_0²))
    inv_all_0   = 1/all_0
    inv_nll_0   = 1/nll_0
    se          = T(0)
    si          = T(0)
    sl          = T(0)
    sgh         = T(0)
    shdq        = T(0)

    sin_i_0, cos_i_0 = sincos(i_0)
    sin_Ω_0, cos_Ω_0 = sincos(Ω_0)
    sin_ω_0, cos_ω_0 = sincos(ω_0)

    sin_i_0²    = sin_i_0^2
    cos_i_0²    = cos_i_0^2
    xpidot      = dotω + dotΩ

    #                        Initial Configuration
    # ==========================================================================

    # Drop terms if inclination is smaller than 3 deg.
    ishq = (i_0 >= 3π/180) ? true : false

    # Do not let `sin_i_0` be 0.
    (abs(sin_i_0) < 1e-12) && (sin_i_0 = sign(sin_i_0)*1e-12)

    # Compute the Greenwhich Mean Sidereal Time at epoch.
    gmst = JDtoGMST(epoch)

    #                      Initialize Lunar Solar Terms
    # ==========================================================================

    # `day` is the number of days since Jan 0, 1900 at 12h.
    day    = epoch - (datetime2julian(DateTime(1900,1,1,12,0,0))-1)

    xnodce     = mod(4.5236020 - 9.2422029e-4day, 2π)
    stem, ctem = sincos(xnodce)

    zcosil = 0.91375164 - 0.03568096ctem
    zsinil = sqrt(1 - zcosil^2)

    zsinhl = 0.089683511stem/zsinil
    zcoshl = sqrt(1 - zsinhl^2)

    gam    = 5.8351514 + 0.0019443680day

    zx     = 0.39785416stem/zsinil
    zy     = zcoshl*ctem + 0.91744867zsinhl*stem
    zx     = atan(zx,zy)
    zx     = gam + zx - xnodce

    zsingl, zcosgl = sincos(zx)

    zmol   = mod(4.7199672 + 0.22997150day - gam, 2π)
    zmos   = mod(6.2565837 + 0.017201977day,      2π)

    #                            Do Solar Terms
    # ==========================================================================
    zcosg = ZCOSGS
    zsing = ZSINGS
    zcosi = ZCOSIS
    zsini = ZSINIS
    zcosh = cos_Ω_0
    zsinh = sin_Ω_0
    cc    = C1SS
    zn    = ZNS
    ze    = ZES
    zmo   = zmos

    for ls = 0:1
        a1  = +zcosg*zcosh + zsing*zcosi*zsinh
        a3  = -zsing*zcosh + zcosg*zcosi*zsinh
        a7  = -zcosg*zsinh + zsing*zcosi*zcosh
        a8  = +zsing*zsini
        a9  = +zsing*zsinh + zcosg*zcosi*zcosh
        a10 = +zcosg*zsini
        a2  = +cos_i_0*a7  + sin_i_0*a8
        a4  = +cos_i_0*a9  + sin_i_0*a10
        a5  = -sin_i_0*a7  + cos_i_0*a8
        a6  = -sin_i_0*a9  + cos_i_0*a10

        x1 = +a1*cos_ω_0 + a2*sin_ω_0
        x2 = +a3*cos_ω_0 + a4*sin_ω_0
        x3 = -a1*sin_ω_0 + a2*cos_ω_0
        x4 = -a3*sin_ω_0 + a4*cos_ω_0
        x5 = +a5*sin_ω_0
        x6 = +a6*sin_ω_0
        x7 = +a5*cos_ω_0
        x8 = +a6*cos_ω_0

        z31 = 12x1^2    - 3x3^2
        z32 = 24x1 * x2 - 6x3 * x4
        z33 = 12x2^2    - 3x4^2
        z1  = 3(   a1^2 + a2^2   ) + z31 * e_0²
        z2  = 6(a1 * a3 + a2 * a4) + z32 * e_0²
        z3  = 3(   a3^2 + a4^2   ) + z33 * e_0²
        z11 = -6a1 * a5 + e_0² * (-24x1 * x7 - 6x3 * x5)
        z12 = -6(a1 * a6 + a3 * a5) + e_0² * ( -24(x2 * x7 + x1 * x8) -
                                                 6(x3 * x6 + x4 * x5) )
        z13 = -6a3 * a6 + e_0² * (-24x2 * x8 - 6x4 * x6)
        z21 = +6a2 * a5 + e_0² * (+24x1 * x5 - 6x3 * x7)
        z22 = +6(a4 * a5 + a2 * a6) + e_0² * (24(x2 * x5 + x1 * x6) -
                                               6(x4 * x7 + x3 * x8) )
        z23 = +6a4 * a6 + e_0² * (24x2 * x6 - 6x4 * x8)
        z1  = +2z1 + (1 - e_0²) * z31
        z2  = +2z2 + (1 - e_0²) * z32
        z3  = +2z3 + (1 - e_0²) * z33
        s3  = +cc * inv_nll_0
        s2  = -0.5s3 / sqrt_1_e_0²
        s4  = +s3 * sqrt_1_e_0²
        s1  = -15e_0 * s4
        s5  = +x1 * x3 + x2 * x4
        s6  = +x2 * x3 + x1 * x4
        s7  = +x2 * x4 - x1 * x3
        se  = +s1 * zn * s5
        si  = +s2 * zn * (z11 + z13)
        sl  = -zn * s3 * (z1 + z3 - 14 - 6e_0²)
        sgh = +s4 * zn * (z31 + z33 - 6)

        shdq = zero(T)

        if ishq
            sh   = -zn * s2 * (z21 + z23);
			shdq = sh / sin_i_0;
        end

        ee2  =  +2s1 * s6
        e3   =  +2s1 * s7
        xi2  =  +2s2 * z12
        xi3  =  +2s2 * (z13 - z11)
        xl2  =  -2s3 * z2
        xl3  =  -2s3 * (z3 - z1)
        xl4  =  -2s3 * (-21 - 9e_0²) * ze
        xgh2 =  +2s4 * z32
        xgh3 =  +2s4 * (z33 - z31)
        xgh4 = -18s4 * ze
        xh2  =  -2s2 * z22
        xh3  =  -2s2 * (z23 - z21)

        (ls == 1) && break

        #                        Do Lunar Terms
        # ======================================================================
        sse   = se
        ssi   = si
        ssl   = sl
        ssh   = shdq
        ssg   = sgh - cos_i_0 * ssh
        se2   = ee2
        si2   = xi2
        sl2   = xl2
        sgh2  = xgh2
        sh2   = xh2
        se3   = e3
        si3   = xi3
        sl3   = xl3
        sgh3  = xgh3
        sh3   = xh3
        sl4   = xl4
        sgh4  = xgh4
        zcosg = zcosgl
        zsing = zsingl
        zcosi = zcosil
        zsini = zsinil
        zcosh = cos_Ω_0 * zcoshl + sin_Ω_0 * zsinhl
        zsinh = sin_Ω_0 * zcoshl - cos_Ω_0 * zsinhl
        zn    = ZNL
        cc    = C1L
        ze    = ZEL
        zmo   = zmol
    end

    sse += se
    ssi += si
    ssl += sl
    ssg += sgh - cos_i_0 * shdq
    ssh += shdq

    if (nll_0 < 0.0052359877) && (nll_0 > 0.0034906585)
        #        24h Synchronous Resonance Terms Initialization
        # ======================================================================

        iresfl = true;
        isynfl = true;

        g200    = e_0² * (0.8125e_0² - 2.5) + 1
        g310    = 2e_0² + 1
        g300    = e_0² * (6.60937e_0² - 6) + 1
        f220    = 0.75(cos_i_0 + 1)^2
        f311    = 0.9375(3.0cos_i_0 + 1)*sin_i_0^2 - 0.75(cos_i_0 + 1)
        f330    = 1.875 * (cos_i_0 + 1)^3
        del1    = 3(nll_0^2 * inv_all_0^2)
        del2    = 2del1 * f220 * g200 * Q22
        del3    = 3del1 * f330 * g300 * Q33 * inv_all_0
        del1    =  del1 * f311 * g310 * Q31 * inv_all_0
        fasx2   = 0.13130908
        fasx4   = 2.8843198
        fasx6   = 0.37448087
        xlamo   = mod(M_0 + Ω_0 + ω_0 - gmst, 2pi)
        bfact   = dotM + xpidot - THDT + ssl + ssg + ssh

    elseif (nll_0 >= 0.00826) && (nll_0 <= 0.00924) && (e_0 >= 0.5)
        #   Geopotential Resonance Initialization for 12 Hour Orbits
        # ======================================================================

        iresfl = true
        isynfl = false

        g201 = -0.306 - 0.44(e_0 - 0.64)

        if e_0 <= 0.65
            g211 = @evalpoly(e_0, + 3.6160, - 13.2470, + 16.29000)
            g310 = @evalpoly(e_0, -19.3020, +117.3900, -228.4190, +156.5910)
            g322 = @evalpoly(e_0, -18.9068, +109.7927, -214.6334, +146.5816)
            g410 = @evalpoly(e_0, -41.1220, +242.6940, -471.0940, +313.9530)
            g422 = @evalpoly(e_0, -146.407, +841.8800, -1629.014, +1083.435)
            g520 = @evalpoly(e_0, -532.114, +3017.977, -5740.032, +3708.276)
        else
            g211 = @evalpoly(e_0, -  72.099, +331.8190, -508.7380, +266.7240)
            g310 = @evalpoly(e_0, - 346.844, +1582.851, -2415.925, +1246.113)
            g322 = @evalpoly(e_0, - 342.585, +1554.908, -2366.899, +1215.972)
            g410 = @evalpoly(e_0, -1052.797, +4758.686, -7193.992, +3651.957)
            g422 = @evalpoly(e_0, -3581.690, +16178.11, -24462.77, +12422.52)

            if e_0 <= 0.715
                g520 = @evalpoly(e_0, +1464.74, -4664.75, +3763.64)
            else
                g520 = @evalpoly(e_0, -5149.66, +29936.92, -54087.36, +31324.56)
            end
        end

        if e_0 < 0.7
            g533 = @evalpoly(e_0, -919.22770, +4988.6100, -9064.7700, +5542.210)
            g521 = @evalpoly(e_0, -822.71072, +4568.6173, -8491.4146, +5337.524)
            g532 = @evalpoly(e_0, -853.66600, +4690.2500, -8624.7700, +5341.400)
        else
            g533 = @evalpoly(e_0, -37995.780, +161616.52, -229838.20, +109377.94)
            g521 = @evalpoly(e_0, -51752.104, +218913.95, -309468.16, +146349.42)
            g532 = @evalpoly(e_0, -40023.880, +170470.89, -242699.48, +115605.82)
        end

        f220 = +0.75(1 + 2cos_i_0 + cos_i_0²)
        f221 = +1.5sin_i_0²
        f321 = +1.875sin_i_0 * (1 - 2cos_i_0 - 3cos_i_0²)
        f322 = -1.875sin_i_0 * (1 + 2cos_i_0 - 3cos_i_0²)
        f441 = +35sin_i_0² * f220
        f442 = +39.375sin_i_0²^2
        f522 = +9.84375sin_i_0 * (sin_i_0² * (+1 - 2cos_i_0 - 5cos_i_0²) +
                                   0.33333333(-2 + 4cos_i_0 + 6cos_i_0²) )
        f523 = sin_i_0 * (4.92187512sin_i_0² * (-2 - 4cos_i_0 + 10cos_i_0²) +
                                     6.56250012(+1 + 2cos_i_0 -  3cos_i_0²) )
        f542 = 29.53125sin_i_0 * (+2 - 8cos_i_0 + cos_i_0² *
                                    (-12 + 8cos_i_0 + 10cos_i_0²))
        f543 = 29.53125sin_i_0 * (-2 - 8cos_i_0 + cos_i_0² *
                                    (+12 + 8cos_i_0 - 10cos_i_0²))

        temp1   = 3(nll_0 * inv_all_0)^2
        temp0   = temp1 * ROOT22
        d2201   = temp0 * f220 * g201
        d2211   = temp0 * f221 * g211
        temp1  *= inv_all_0
        temp0   = temp1 * ROOT32
        d3210   = temp0 * f321 * g310
        d3222   = temp0 * f322 * g322
        temp1  *= inv_all_0
        temp0   = 2temp1 * ROOT44
        d4410   = temp0 * f441 * g410
        d4422   = temp0 * f442 * g422
        temp1  *= inv_all_0
        temp0   = temp1 * ROOT52
        d5220   = temp0 * f522 * g520
        d5232   = temp0 * f523 * g532
        temp0   = 2temp1 * ROOT54
        d5421   = temp0 * f542 * g521
        d5433   = temp0 * f543 * g533
        xlamo   = mod(M_0 + 2Ω_0 - 2gmst, 2pi)
        bfact   = dotM + 2dotΩ - 2THDT + ssl + 2ssh
    else
        #                     Non Resonant Orbits
        # ======================================================================

        iresfl = false
        isynfl = false
    end

	if iresfl
        #                   Initialize the Integrator
        # ======================================================================
		xfact = bfact - nll_0
		xli   = xlamo
		atime = 0.0

        # TODO: Check if this variable can be removed from SGP4_DeepSpace.
        xni   = nll_0

        # Compute the "dot" terms.
        # ========================

        if (isynfl)
            sin_1, cos_1 = sincos(  (xli - fasx2) )
            sin_2, cos_2 = sincos( 2(xli - fasx4) )
            sin_3, cos_3 = sincos( 3(xli - fasx6) )

            xndot = del1 * sin_1 +  del2 * sin_2 +  del3 * sin_3
            xnddt = del1 * cos_1 + 2del2 * cos_2 + 3del3 * cos_3
        else
            ω     = ω_0 + dotω * atime

            sin_1,  cos_1  = sincos(2ω + xli  - G22)
            sin_2,  cos_2  = sincos(   + xli  - G22)
            sin_3,  cos_3  = sincos(+ω + xli  - G32)
            sin_4,  cos_4  = sincos(-ω + xli  - G32)
            sin_5,  cos_5  = sincos(+ω + xli  - G52)
            sin_6,  cos_6  = sincos(-ω + xli  - G52)
            sin_7,  cos_7  = sincos(2ω + 2xli - G44)
            sin_8,  cos_8  = sincos(     2xli - G44)
            sin_9,  cos_9  = sincos(+ω + 2xli - G54)
            sin_10, cos_10 = sincos(-ω + 2xli - G54)

            xndot = d2201 * sin_1 + d2211 * sin_2 + d3210 * sin_3 +
                    d3222 * sin_4 + d5220 * sin_5 + d5232 * sin_6 +
                    d4410 * sin_7 + d4422 * sin_8 + d5421 * sin_9 +
                    d5433 * sin_10

            xnddt =   d2201 * cos_1 + d2211 * cos_2 + d3210 * cos_3 +
                      d3222 * cos_4 + d5220 * cos_5 + d5232 * cos_6 +
                    2(d4410 * cos_7 + d4422 * cos_8 + d5421 * cos_9 +
                      d5433 * cos_10)
        end

        xldot  = xni + xfact
        xnddt *= xldot
    end

	# Set up for original mode (LS terms at epoch non-zero).
	pgh0 = ph0 = pe0 = pinc0 = pl0 = 0.0

    @pack! sgp4_ds = atime, xli, xni, xnq, xfact, ssl, ssg, ssh, sse, ssi,
                     xlamo, omegaq, omgdt, gmst, del1, del2, del3, fasx2, fasx4,
                     fasx6, d2201, d2211, d3210, d3222, d4410, d4422, d5220,
                     d5232, d5421, d5433, xnddt, xndot, xldot, zmos, se2, se3,
                     si2, si3, sl2, sl3, sl4, sgh2, sgh3, sgh4, sh2, sh3, zmol,
                     ee2, e3, xi2, xi3, xl2, xl3, xl4, xgh2, xgh3, xgh4, xh2,
                     xh3, pe, pinc, pgh, ph, pl, pgh0, ph0, pe0, pinc0, pl0,
                     isynfl, iresfl, ilsz

    return sgp4_ds
end

"""
    dssec!(sgp4_ds::SGP4_DeepSpace{T}, nll_0::T, e_0::T, i_0::T, ω_0::T, Ω_k::T, ω_k::T, M_k::T, dotω::T, Δt::Number) where T<:Number

Compute the secular effects.

Notice that the values in the structure `sgp4_ds` **will be modified**.

# Args

* `sgp4_ds`: Deep space structure (see `SGP4_DeepSpace`).
* `nll_0`: Initial mean motion [rad/min].
* `e_0`: Initial eccentricity.
* `i_0`: Initial inclination [rad].
* `ω_0`: Initial argument of perigee [rad].
* `Ω_k`: Current right ascension of the ascending node [rad].
* `ω_k`: Current argument of perigee [rad].
* `M_k`: Current mean anomaly [rad].
* `dotω`: Time-derivative of the argument of perigee [rad/min].
* `Δt`: Time interval since the epoch [min].

# Returns

The following elements perturbed by the secular effects:

* Mean motion [rad/min].
* Eccentricity.
* Inclination [rad].
* Right ascension of the ascending node [rad].
* Argument of perigee [rad].
* Mean anomaly [rad].

"""
function dssec!(sgp4_ds::SGP4_DeepSpace{T},
                nll_0::T,
                e_0::T,
                i_0::T,
                ω_0::T,
                Ω_k::T,
                ω_k::T,
                M_k::T,
                dotω::T,
                Δt::Number) where T<:Number
    # Unpack variables.
    @unpack_SGP4_DeepSpace sgp4_ds

    M_sec = M_k + ssl * Δt
    e_sec = e_0 + sse * Δt
    i_sec = i_0 + ssi * Δt
    Ω_sec = Ω_k + ssh * Δt
    ω_sec = ω_k + ssg * Δt

    # TODO: Verify what this variable means. This is found in `dspace.m` of
    # Vallado's implementation [2].
    θ = mod(gmst + THDT * Δt, 2π)

    # If the orbit is not resonant, then nothing more should be computed.
    (!iresfl) && return (nll_0, e_sec, i_sec, Ω_sec, ω_sec, M_sec)

    #   Update Resonances using Numerical (Euler-Maclaurin) Integration
    # ==========================================================================

    # Epoch restart
    # =============

    # This verification is different between Vallado's [2] and [3]. We will use
    # [2] since it seems more recent.
    if  (atime == 0) || (Δt * atime <= 0) || (abs(Δt) < abs(atime))
        atime = T(0)
        xni   = nll_0
        xli   = xlamo
    end

    # Integration
    # ===========

    ft = Δt - atime

    # In [3], the integration process is performed only if `ft` is larger than
    # `STEP`. However, Vallado's implementation [2] does not verify this and the
    # integration is performed every time. This behavior was chose because it
    # seems that [3] is a more recent version of the algorithm.

    # Check integration direction.
    delt = (Δt >= atime) ? STEP : -STEP

    # Perform the integration with step `delt` until the difference between
    # the time `Δt` and `atime` is less then `STEP`.
    while true
        # Compute the dot terms.
        if (isynfl)

            sin_1, cos_1 = sincos(  (xli - fasx2) )
            sin_2, cos_2 = sincos( 2(xli - fasx4) )
            sin_3, cos_3 = sincos( 3(xli - fasx6) )

            xndot = del1 * sin_1 +  del2 * sin_2 +  del3 * sin_3
            xnddt = del1 * cos_1 + 2del2 * cos_2 + 3del3 * cos_3
        else
            ω     = ω_0 + dotω * atime

            sin_1,  cos_1  = sincos(2ω + xli  - G22)
            sin_2,  cos_2  = sincos(   + xli  - G22)
            sin_3,  cos_3  = sincos(+ω + xli  - G32)
            sin_4,  cos_4  = sincos(-ω + xli  - G32)
            sin_5,  cos_5  = sincos(+ω + xli  - G52)
            sin_6,  cos_6  = sincos(-ω + xli  - G52)
            sin_7,  cos_7  = sincos(2ω + 2xli - G44)
            sin_8,  cos_8  = sincos(     2xli - G44)
            sin_9,  cos_9  = sincos(+ω + 2xli - G54)
            sin_10, cos_10 = sincos(-ω + 2xli - G54)

            xndot = d2201 * sin_1 + d2211 * sin_2 + d3210 * sin_3 +
                    d3222 * sin_4 + d5220 * sin_5 + d5232 * sin_6 +
                    d4410 * sin_7 + d4422 * sin_8 + d5421 * sin_9 +
                    d5433 * sin_10

            xnddt =   d2201 * cos_1 + d2211 * cos_2 + d3210 * cos_3 +
                      d3222 * cos_4 + d5220 * cos_5 + d5232 * cos_6 +
                    2(d4410 * cos_7 + d4422 * cos_8 + d5421 * cos_9 +
                      d5433 * cos_10)
        end

        xldot  = xni + xfact
        xnddt *= xldot

        ft = Δt - atime
        (abs(ft) < STEP) && break

        # In Vallado's implementation [2], this is in the final of the loop
        # instead of at the beginning.
        xli   += delt * (xldot + delt * xndot / 2)
        xni   += delt * (xndot + delt * xnddt / 2)
        atime += delt
    end

    xl    = xli + ft * (xldot + ft * xndot / 2)
    n_sec = xni + ft * (xndot + ft * xnddt / 2)
    M_sec = !isynfl ? xl - 2Ω_sec + 2θ : xl - Ω_sec - ω_sec + θ

    @pack! sgp4_ds = atime, xni, xli, xnddt, xndot, xldot

    (n_sec, e_sec, i_sec, Ω_sec, ω_sec, M_sec)
end

"""
    dsper!(sgp4_ds::SGP4_DeepSpace{T}, e_k::T, i_k::T, Ω_k::T, ω_k::T, M_k::T, Δt:Number) where T<:Number

Compute the effects caused by Lunar-Solar periodics.

Notice that the values in the structure `sgp4_ds` **will be modified**.

# Args

* `sgp4_ds`: Deep space structure (see `SGP4_DeepSpace`).
* `e_k`: Current eccentricity.
* `i_k`: Current inclination [rad].
* `Ω_k`: Current right ascension of the ascending node [rad].
* `ω_k`: Current argument of perigee [rad].
* `M_k`: Current mean anomaly [rad].
* `Δt`: Time interval since the epoch [min].

# Returns

The following elements perturbed by lunar-solar periodics.

* Eccentricity.
* Inclination [rad].
* Right ascension of the ascending node [rad].
* Argument of perigee [rad].
* Mean anomaly [rad].

"""
function dsper!(sgp4_ds::SGP4_DeepSpace{T},
                e_k::T,
                i_k::T,
                Ω_k::T,
                ω_k::T,
                M_k::T,
                Δt::Number) where T<:Number
    # Unpack variables.
    @unpack_SGP4_DeepSpace sgp4_ds

    #                          Update Solar Terms
    # ==========================================================================

	zm    = zmos +  ZNS * Δt
	zf    = zm   + 2ZES * sin(zm)

    sinzf, coszf = sincos(zf)

	f2    = +sinzf * sinzf / 2 - 0.25
	f3    = -sinzf * coszf / 2
	ses   = se2 * f2 + se3 * f3
	sis   = si2 * f2 + si3 * f3
	sls   = sl2 * f2 + sl3 * f3 + sl4 * sinzf

	sghs  = sgh2 * f2 + sgh3 * f3 + sgh4 * sinzf
	shs   = sh2  * f2 + sh3  * f3

    #                          Update Lunar Terms
    # ==========================================================================

	zm    = zmol +  ZNL * Δt
	zf    = zm   + 2ZEL * sin(zm)

    sinzf, coszf = sincos(zf)

	f2    = +sinzf * sinzf / 2 - 0.25
	f3    = -sinzf * coszf / 2
	sel   = ee2 * f2 + e3 * f3
	sil   = xi2 * f2 + xi3 * f3
	sll   = xl2 * f2 + xl3 * f3 + xl4 * sinzf

	sghl  = xgh2 * f2 + xgh3 * f3 + xgh4 * sinzf
	shl   = xh2  * f2 + xh3  * f3

    #                         Save computed values
    # ==========================================================================

	pgh  = sghs + sghl
	ph   = shs  + shl
	pe   = ses  + sel
	pinc = sis  + sil
	pl   = sls  + sll

    # Update inclination and eccentricity.
    e_per = e_k + pe
    i_per = i_k + pinc

    sinis, cosis = sincos(i_per)

    # The original algorithm considered the original inclination to select the
    # Lyddane Lunar-Solar perturbations algorithm. However, Vallado's
    # implementation [2] test the perturbed inclination to select this. It is
    # mentioned that this is the behavior selected in GSFC source code.
    if i_per >= 0.2
        tmp_ph = ph / sinis;
		ω_per  = ω_k + pgh - cosis * tmp_ph;
		Ω_per  = Ω_k + tmp_ph;
		M_per  = M_k + pl;
    else
        sinok   = sin(Ω_k)
        cosok   = cos(Ω_k)

        #                       |----------    dalf     ----------|
		alfdp   = sinis * sinok + ph * cosok + pinc * cosis * sinok
        #                       |----------    dbet     ----------|
		betdp   = sinis * cosok - ph * sinok + pinc * cosis * cosok

        # For the following computation, in which `Ω_per` is used without a
        # trigonometric function, it is advisable to make sure that it stays in
        # the interval [0, 2π].
        Ω_per   = mod(Ω_k, 2π)

        #                                   |----------    dls    ----------|
		xls     = M_k + ω_k + cosis * Ω_per + pl + pgh - pinc * Ω_per * sinis
        Ω_aux   = Ω_per
        Ω_per   = mod(atan(alfdp, betdp), 2π)

        if abs(Ω_aux - Ω_per) > π
            Ω_per = (Ω_per < Ω_aux) ? Ω_per + 2π : Ω_per - 2π
        end

		M_per   = M_k + pl;
		ω_per   = xls - M_per - cosis * Ω_per
    end

    @pack! sgp4_ds = pgh, ph, pe, pinc, pl

    (e_per, i_per, Ω_per, ω_per, M_per)
end


text:
import JSON

function readJSON(json_file="./config.json")
    open(json_file) do f
        r = JSON.parse(f)
        return r
    end
end

# using JSON

# function readJSON(fpath)
#   jcont = JSON.parsefile(fpath)
#   for (a,b) in jcont
#     println(jcont[a],b)
#   end
#   println(jcont["input"]["art_sig1"])
# end
# readJSON("./config.json")

max_stars_repo_name: av-maslov/Recurrency

max_issues_repo_name: av-maslov/Recurrency

max_forks_repo_name: av-maslov/Recurrency


text:
#possible issues
#1. intermediate yy needed; MeshArrays errors otherwise
#2. lack of exch_z, which would allow using corner points
#4. Float32 / Float64 dispatch may be an issue in e.g.
#   `pi/2 .-Γ.YC*pi/180`

#    XC=γ.read(γ.write(Γ.XC),MeshArray(γ,Float64))
#    YC=γ.read(γ.write(Γ.YC),MeshArray(γ,Float64))

using MeshArrays
import Meshes, MeshViz
import GLMakie as Mkie
#import CairoMakie as Mkie

col=[:blue,:greenyellow,:magenta,:yellow2,:tomato,:black]

"""
    plot_as_sphere(Γ)

Plot mesh (Γ.XC,Γ.YC) and depth (Γ.Depth) on the surface of the sphere in 3D.
"""
function plot_as_sphere(Γ; title="", az=-0.2π, el=0.2π)
    yy=pi/2 .-Γ.YC*pi/180
    x=sin.(yy)*cos.(Γ.XC*pi/180)
    y=sin.(yy)*sin.(Γ.XC*pi/180)
    z=cos.(yy)
    d=Γ.Depth
    crng=(minimum(d),maximum(d))
    minimum(d)==maximum(d) ? crng=(0.8*crng[1],1.1*crng[1]) : nothing

    fig = Mkie.Figure(resolution = (900,900), backgroundcolor = :grey95)
    ax = Mkie.Axis3(fig, aspect = :data, viewmode = :fitzoom, #perspectiveness = 0.5,
        azimuth = az, elevation = el)
    for i=1:length(z.fSize)
        Mkie.surface!(ax, x[i], y[i], z[i], color = d[i], colorrange = crng, colormap = :grays)
        Mkie.wireframe!(x[i],y[i],z[i], overdraw = false, linewidth = 0.25,color=col[i])
    end

    Mkie.hidedecorations!(ax)
    Mkie.hidespines!(ax)
    ax.title=title
    fig[1,1] = ax
    fig
end

"""
    plot_as_plane(Γ)

Plot mesh (Γ.XC,Γ.YC) and depth (Γ.Depth) on a 2D plane.
"""
function plot_as_plane(Γ; title="")
    x=Float64.(Γ.XC)
    y=Float64.(Γ.YC)
	z=0*y
    d=Float64.(Γ.Depth)
    crng=(minimum(d),maximum(d))
    minimum(d)==maximum(d) ? crng=(0.8*crng[1],1.1*crng[1]) : nothing

    fig = Mkie.Figure(resolution = (900,900), backgroundcolor = :grey95)

    ax = Mkie.Axis(fig[1,1])
    for i=1:length(x.fSize)
        Mkie.surface!(ax, x[i], y[i], z[i], color = d[i], colorrange = crng, colormap = :grays)
		#contourf!(x[i], y[i], d[i], colorrange = (0.0, 5e3), colormap = :grays)
		Mkie.wireframe!(x[i],y[i], z[i], overdraw = false, linewidth = 0.25,color=col[i])
    end
    #hidedecorations!(ax)
    #hidespines!(ax)
    ax.title=title
    fig[1,1] = ax
    fig
end

"""
    plot_as_scatter(Γ; title="")

Plot mesh (Γ.XC,Γ.YC) as scatter of points.
"""
function plot_as_scatter(Γ; title="")
    x=Float64.(Γ.XC)
    y=Float64.(Γ.YC)

    fig = Figure(resolution = (900,900), backgroundcolor = :grey95)

    ax = Axis(fig[1,1])
    for i=1:length(x.fSize)
        lines!(ax, x[i][:], y[i][:], color = col[i])
        scatter!(ax, x[i][:], y[i][:], color = col[i], marker = 'o')
    end
    #hidedecorations!(ax)
    #hidespines!(ax)
    ax.title=title
    fig[1,1] = ax
    fig
end

#Γ=GridLoad(GridSpec("LatLonCap",MeshArrays.GRID_LLC90))
function plot_as_plane_b()
    Mkie.set_theme!(Mkie.theme_light())
    fig = Mkie.Figure(resolution = (900,600),markersize=0.1)
    ax = Mkie.Axis(fig[1,1],xlabel="longitude",ylabel="latitude")

    col=[:red,:green,:blue,:magenta,:cyan,:black]

    for f in 1:length(Γ.XC)
        tmp=[Meshes.Point2(Γ.XC[f][i],Γ.YC[f][i]) for i in eachindex(Γ.XC[f])]
        MeshViz.viz!(ax,tmp,elementcolor=col[f],markersize=2.0)
    end

    fig,ax
end

#Γ=GridLoad(GridSpec("LatLonCap",MeshArrays.GRID_LLC90))
function plot_as_sphere_b()
    Mkie.set_theme!(Mkie.theme_light())
    fig = Mkie.Figure(resolution = (900,600),markersize=0.1)
    ax = Mkie.Axis3(fig[1,1])

    col=[:red,:green,:blue,:magenta,:cyan,:black]

    yy=pi/2 .-Γ.YC*pi/180
    x=sin.(yy)*cos.(Γ.XC*pi/180)
    y=sin.(yy)*sin.(Γ.XC*pi/180)
    z=cos.(yy)
    #d=Γ.Depth

    #sphr=sample(Meshes.Sphere((0.0,0.0,0.0), 0.99), Meshes.RegularSampling(360, 180))
    #sphr=Meshes.Sphere((0.0,0.0,0.0), 0.99)
    #MeshViz.viz!(ax,sphr)

    sphere = Mkie.Sphere(Mkie.Point3f(0.0), 0.99)
    Mkie.mesh!(ax,sphere, color=:white)

    for f in 1:length(Γ.XC)
        tmp=[Meshes.Point3(x[f][i],y[f][i],z[f][i]) for i in eachindex(Γ.XC[f])]
        MeshViz.viz!(ax,tmp,color=col[f],pointsize=1000.0)
    end

    Mkie.hidedecorations!(ax)
    Mkie.hidespines!(ax)
    Mkie.ylims!(-0.7,0.7)
    Mkie.xlims!(-0.7,0.7)
    Mkie.zlims!(-0.7,0.7)

    fig,ax
end


text:
using Manifolds
using Manifolds: +ₚ, *ₚ, offsetframe, dimincrease, dim, codim
using Test

using LinearAlgebra
using Random
using StaticArrays
using Suppressor

# Randomized equality checking. Also tests call and rand.
function ≅(m1::AbstractManifold{D}, m2::AbstractManifold{D}, n=1000) where D
    params = rand(D, n)
    vals1 = hcat(m1(zeros(D)...), mapslices(x -> m1(x...), params, dims = 1))
    vals2 = hcat(m2(zeros(D)...), mapslices(x -> m2(x...), params, dims = 1))
    if vals1 ≈ vals2
        true
    else
        bad = findall([vals1[:, i] ≉ vals2[:, i] for i in axes(vals1, 1)])
        println("Inconsistencies found in $(length(bad)) parameters. Examples:")
        foreach(Iterators.take(bad, 10)) do b
            println("  ", params[:, b], ", ", vals1[:, b], ", ", vals2[:, b])
        end
        false
    end
end
function ≅(m1::AbstractManifold, m2::AbstractManifold, n=1000)
    println("Mismatched dimensions.")
    false
end

⟂(v1::AbstractVector, v2::AbstractVector) = iszero(v1 ⋅ v2)

include("base.jl")
include("modifiers.jl")
#include("products.jl")


text:
import Pkg

function print_welcome_message(a::AbstractDict = ENV)::Nothing
    predictmd_version::VersionNumber = version()
    predictmd_pkgdir::String = package_directory()
    predictmdapi_version::VersionNumber = api_version()
    predictmdapi_pkgdir::String = api_package_directory()
    if !_suppress_welcome_message(a)
        @info("This is PredictMD, version $(predictmd_version)")
        @info("For help, please visit https://predictmd.net")
        @debug("PredictMD package directory: $(predictmd_pkgdir)")
        @debug("PredictMDAPI version: $(predictmdapi_version)")
        @debug("PredictMDAPI package directory: $(predictmdapi_pkgdir)")
    end
    return nothing
end

function _suppress_welcome_message(a::AbstractDict = ENV)::Bool
    return get(a, "SUPPRESS_PREDICTMD_WELCOME_MESSAGE", "false") == "true"
end


text:
module InventoryManagement

using DataFrames
using Reexport
using Random
using IntervalSets
using Distributions

@reexport using Graphs
@reexport using MetaGraphs

include("environment.jl")
include("utils.jl")
include("actions.jl")
include("policy.jl")
include("spaces.jl")

export SupplyChainEnv, reset!, is_terminated, action_space, show_action
export reorder_policy, simulate_policy!, service_measures

end


text:
"The perfect sidekick to your scientific inquiries"
module DrWatson
import Pkg, LibGit2

# Pure Julia implementation
include("project_setup.jl")
include("naming.jl")
include("saving_tools.jl")

# Functionality that saves/loads
import FileIO
wsave = FileIO.save
wload = FileIO.load
include("saving_files.jl")

# Functionality that requires Dataframes and other heavy dependencies:
using Requires
function __init__()
    @require DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" begin
        include("result_collection.jl")
    end
end

function greet()
    println(
    """
    DrWatson is currently in beta. Help us make it better by opening
    issues on GitHub or submitting feature requests!

    Have fun with your new project!
    """
    )
end

end


text:
# const PointCollection = StructVector{TimedPoint}
# pointcollection(x::Missing, t₀) = StructVector{TimedPoint}(undef, 0)
# pointcollection(x, t₀) = StructVector(TimedPoint(Space(i[1], i[2]), i[3] - t₀) for i in eachrow(x.data))

struct Track
    coords::Vector{Space}
    t::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}
    tp::Int
    rawcoords
end

function findturningpoint(::Missing, spl, tl, _)
    tp = gettpknot(spl)
    i = findfirst(≥(tp), tl)
    isnothing(i) && return length(tl)
    i
end

function findturningpoint(t, _, tl, t₀)
    i = findfirst(≥(t - t₀), tl)
    return isnothing(i) ? length(tl) : i
end

function Track(xyt, tp)
    sort!(xyt, by = row -> row.t)
    xyt.t .-= xyt.t[1]
    unique!(x -> x.xy, xyt)
    Δt = mean(trim(diff(xyt.t), prop = 0.1))
    t, xy = filterdance(xyt.xy, Δt)
    spl = ParametricSpline(t, reduce(hcat, xy); s = 10, k = 1) # s = 500
    tl = range(0.0, step = Δt, stop = t[end])
    xyl = Space.(spl.(tl))
    i = findturningpoint(tp, spl, tl, xyt[1].t)
    Track(xyl, tl, i, xyt)
end

# homing(t::Track) = t.coords[1:t.tp]
# searching(t::Track) = t.coords[t.tp:end]
# center_of_search(t::Track) = mean(searching(t))
# turning_point(t::Track) = t.coords[t.tp]
Base.getproperty(x::Track, k::Symbol) = k === :homing ? x.coords[1:x.tp] :
                                         k === :searching ?  x.coords[x.tp:end] :
                                         k === :center_of_search ? mean(x.coords[x.tp:end]) :
                                         k === :turning_point ? x.coords[x.tp] :
                                         getfield(x, k)

const ignorefirst = 10 # cm
const bigturn = π/3 # 60°

_getv(spl, k) = SVector{2, Float64}(derivative(spl, k))

function gettpindex(spl, ks)
    tp = ks[1]
    vlast = _getv(spl, ks[1])
    for k in Iterators.drop(ks, 1)
        v = _getv(spl, k)
        Δ = angle(vlast, v)
        tp = k
        Δ > bigturn && break
        vlast = v
    end
    return tp
end

function gettpknot(spl)
    ks = Dierckx.get_knots(spl)
    filter!(k -> norm(spl(k) - spl(0)) > ignorefirst, ks)
    isempty(ks) && return 0.0
    tp2 = gettpindex(spl, ks)
    # return tp2
    tp1 = copy(tp2)
    for k in ks
        k == tp2 && break
        tp1 = k
    end
    tp1 += 0.1
    if tp1 < tp2
        main = _getv(spl, tp1)
        for t in tp2:-0.3:tp1
            v = _getv(spl, t)
            Δ = angle(main, v)
            Δ < bigturn && return t
        end
    end
    return tp2
end

function filterdance(xy, Δt)
    xy2 = [xy[1]]
    t = [0.0]
    for p in xy
        if norm(p - xy2[end]) > 4
            push!(xy2, p)
            push!(t, t[end] + Δt)
        else
            t[end] += Δt
        end
    end
    t .-= t[1]
    return t, xy2
end


text:
using KeyedFrames
using DataFrames
using Test

@testset "KeyedFrames" begin
    df1 = DataFrame(; a=1:10, b=2:11, c=3:12)
    df2 = DataFrame(; a=1:5, d=4:8)
    df3 = DataFrame(; a=[4, 2, 1], e=[2, 5, 2], f=1:3)

    @testset "constructor" begin
        kf1 = KeyedFrame(df1, [:a, :b])
        @test KeyedFrames.frame(kf1) === df1
        @test keys(kf1) == [:a, :b]

        kf2 = KeyedFrame(df2, :a)
        @test KeyedFrames.frame(kf2) === df2
        @test keys(kf2) == [:a]

        kf3 = KeyedFrame(df3, ["e", "a"])
        @test KeyedFrames.frame(kf3) === df3
        @test keys(kf3) == [:e, :a]

        @test_throws ArgumentError KeyedFrame(df1, [:a, :b, :d])
        @test_throws ArgumentError KeyedFrame(df1, :d)

        @test keys(KeyedFrame(df1, [:a, :a, :b, :a])) == [:a, :b]
    end

    kf1 = KeyedFrame(df1, [:a, :b])
    kf2 = KeyedFrame(df2, :a)
    kf3 = KeyedFrame(df3, ["e", "a"])

    @testset "equality" begin
        cp = deepcopy(kf1)
        cpd = deepcopy(df1)

        @test kf1 == kf1
        @test isequal(kf1, kf1)
        @test hash(kf1) == hash(kf1)

        @test kf1 == cp
        @test isequal(kf1, cp)
        @test hash(kf1) == hash(cp)

        @test kf1 == KeyedFrame(cpd, [:a, :b])
        @test isequal(kf1, KeyedFrame(cpd, [:a, :b]))
        @test hash(kf1) == hash(KeyedFrame(cpd, [:a, :b]))

        @test kf1 == KeyedFrame(cpd, [:b, :a])
        @test !isequal(kf1, KeyedFrame(cpd, [:b, :a]))
        @test hash(kf1) != hash(KeyedFrame(cpd, [:b, :a]))

        @test kf1 == df1
        @test df1 == kf1
        @test !isequal(kf1, df1)
        @test !isequal(df1, kf1)
        @test hash(kf1) != hash(df1)

        @test kf1 != KeyedFrame(cpd, [:a, :b, :c])
        @test !isequal(kf1, KeyedFrame(cpd, [:a, :b, :c]))
        @test hash(kf1) != hash(KeyedFrame(cpd, [:a, :b, :c]))

        @test kf2 != KeyedFrame(df3, :a)
        @test !isequal(kf2, KeyedFrame(df3, :a))
        @test hash(kf2) != hash(KeyedFrame(df3, :a))
    end

    @testset "copy" begin
        @test kf1 == deepcopy(kf1)
        @test kf1 !== deepcopy(kf1)
        @test DataFrame(kf1) !== DataFrame(deepcopy(kf1))
        @test keys(kf1) !== keys(deepcopy(kf1))
    end

    @testset "convert" begin
        for (kf, df) in [(kf1, df1), (kf2, df2), (kf3, df3)]
            @test convert(DataFrame, kf) == df
            @test DataFrame(kf) == df
            @test SubDataFrame(kf, 1:3, :) == SubDataFrame(df, 1:3, :)
        end
    end

    @testset "size" begin
        for (kf, df) in [(kf1, df1), (kf2, df2), (kf3, df3)]
            @test nrow(kf) == nrow(df)
            @test ncol(kf) == ncol(df)
            @test size(kf) == size(df)
        end
    end

    @testset "names/index/key" begin
        for (kf, df) in [(kf1, df1), (kf2, df2), (kf3, df3)]
            @test KeyedFrames.names(kf) == DataFrames.names(df)
            @test KeyedFrames.index(kf) == DataFrames.index(df)
            @test keys(kf) == keys(kf)
        end
    end

    @testset "getindex" begin
        @test isequal(kf1[:], kf1)
        @test isequal(kf1[:, :], kf1)

        @test isequal(kf1[1, :], KeyedFrame(DataFrame(; a=1, b=2, c=3), [:a, :b]))
        @test isequal(kf1[8:10, :], KeyedFrame(DataFrame(; a=8:10, b=9:11, c=10:12), [:a, :b]))

        @test isequal(kf1[1:2], KeyedFrame(DataFrame(; a=1:10, b=2:11), [:a, :b]))
        @test isequal(kf1[[:a, :b]], KeyedFrame(DataFrame(; a=1:10, b=2:11), [:a, :b]))

        @test isequal(kf1[1, [:a, :b]], KeyedFrame(DataFrame(; a=1, b=2), [:a, :b]))
        @test isequal(kf1[8:10, 1:2], KeyedFrame(DataFrame(; a=8:10, b=9:11), [:a, :b]))

        @test isequal(kf1[:, 1:2], KeyedFrame(DataFrame(; a=1:10, b=2:11), [:a, :b]))
        @test isequal(kf1[:, [:a, :b]], KeyedFrame(DataFrame(; a=1:10, b=2:11), [:a, :b]))

        # When :a column disappears it is removed from the key
        @test isequal(kf1[2:3], KeyedFrame(DataFrame(; b=2:11, c=3:12), :b))
        @test isequal(kf1[[:b, :c]], KeyedFrame(DataFrame(; b=2:11, c=3:12), :b))

        @test isequal(kf1[:, 2:3], KeyedFrame(DataFrame(; b=2:11, c=3:12), :b))
        @test isequal(kf1[:, [:b, :c]], KeyedFrame(DataFrame(; b=2:11, c=3:12), :b))

        # Returns a column/value instead of a KeyedFrame
        @test isequal(kf1[1, :a], 1)
        @test isequal(kf1[8:10, 1], [8, 9, 10])

        @test isequal(kf1[:, :b], collect(2:11))
        @test isequal(kf1[:, 2], collect(2:11))
    end

    @testset "setindex!" begin
        cp = deepcopy(kf1)
        cp[:b] = collect(11:20)     # Need to collect on a setindex! with DataFrames.
        @test isequal(cp, KeyedFrame(DataFrame(; a=1:10, b=11:20, c=3:12), [:a, :b]))

        cp = deepcopy(kf1)
        cp[:b] = 3
        @test isequal(cp, KeyedFrame(DataFrame(; a=1:10, b=3, c=3:12), [:a, :b]))

        cp = deepcopy(kf1)
        cp[2] = 3
        @test isequal(cp, KeyedFrame(DataFrame(; a=1:10, b=3, c=3:12), [:a, :b]))

        cp = deepcopy(kf1)
        cp[[:b, :c]] = 3
        @test isequal(cp, KeyedFrame(DataFrame(; a=1:10, b=3, c=3), [:a, :b]))

        cp = deepcopy(kf1)
        cp[2:3] = 3
        @test isequal(cp, KeyedFrame(DataFrame(; a=1:10, b=3, c=3), [:a, :b]))

        cp = deepcopy(kf2)
        cp[1, :d] = 10
        @test isequal(cp, KeyedFrame(DataFrame(; a=1:5, d=[10, 5, 6, 7, 8]), :a))

        cp = deepcopy(kf2)
        cp[1, 1:2] = 10
        @test isequal(cp, KeyedFrame(DataFrame(;a=[10, 2, 3, 4, 5],d=[10, 5, 6, 7, 8]), :a))
    end

    @testset "first/last" begin
        # Don't assume that n will always equal 6
        @test first(kf1) isa KeyedFrame
        @test isequal(first(kf1, 1), kf1[1, :])
        @test isequal(first(kf1, 3), kf1[1:3, :])
        @test isequal(first(kf1, 6), kf1[1:6, :])

        @test last(kf1) isa KeyedFrame
        @test isequal(last(kf1, 1), kf1[end, :])
        @test isequal(last(kf1, 3), kf1[end - 2:end, :])
        @test isequal(last(kf1, 6), kf1[end - 5:end, :])
    end

    @testset "sort" begin
        # No columns specified (sort by key)
        expected = KeyedFrame(DataFrame(; a=[1, 4, 2], e=[2, 2, 5], f=[3, 1, 2]), [:e, :a])
        reversed = KeyedFrame(DataFrame(; a=[2, 4, 1], e=[5, 2, 2], f=[2, 1, 3]), [:e, :a])

        @test sort(kf3) == expected
        @test sort(kf3; rev=true) == reversed

        cp = deepcopy(kf3)
        @test sort!(cp) == expected
        @test cp == expected

        cp = deepcopy(kf3)
        @test sort!(cp; rev=true) == reversed
        @test cp == reversed

        @test issorted(kf1)
        @test !issorted(kf3)
        @test issorted(expected)
        @test !issorted(reversed)
        @test issorted(reversed; rev=true)

        # Columns specified
        expected = KeyedFrame(DataFrame(; a=[1, 2, 4], e=[2, 5, 2], f=[3, 2, 1]), [:e, :a])
        reversed = KeyedFrame(DataFrame(; a=[4, 2, 1], e=[2, 5, 2], f=[1, 2, 3]), [:e, :a])

        @test sort(kf3, :a) == expected
        @test sort(kf3, :a; rev=true) == reversed

        cp = deepcopy(kf3)
        @test sort!(cp, :a) == expected
        @test cp == expected

        cp = deepcopy(kf3)
        @test sort!(cp, :a; rev=true) == reversed
        @test cp == reversed

        @test !issorted(expected)
        @test issorted(expected, :a)
        @test !issorted(reversed)
        @test !issorted(reversed, :a)
        @test issorted(reversed, :a; rev=true)

        # Test return type of `sort!`
        @test isa(sort!(deepcopy(kf3), :a), KeyedFrame)
    end

    @testset "push!" begin
        cp = deepcopy(kf2)

        push!(cp, [6, 9])
        @test cp == KeyedFrame(DataFrame(; a=1:6, d=4:9), :a)

        # Test return type of `push!`
        @test isa(push!(deepcopy(kf2), [6, 9]), KeyedFrame)
    end

    @testset "append!" begin
        cp = deepcopy(kf2)
        append!(cp, DataFrame(; a=6:8, d=9:11))
        @test cp == KeyedFrame(DataFrame(; a=1:8, d=4:11), :a)

        # With append! we discard the key from the second KeyedFrame, which I think is fine.
        cp = deepcopy(kf2)
        append!(cp, KeyedFrame(DataFrame(; a=6:8, d=9:11), [:a, :d]))
        @test cp == KeyedFrame(DataFrame(; a=1:8, d=4:11), :a)

        # Test return type of `append!`
        @test isa(
            append!(deepcopy(kf2), KeyedFrame(DataFrame(; a=6:8, d=9:11), [:a, :d])),
            KeyedFrame
        )
    end

    @testset "deleterows!" begin
        cp = deepcopy(kf1)
        deleterows!(cp, 1)
        @test cp == KeyedFrame(DataFrame(; a=2:10, b=3:11, c=4:12), [:a, :b])

        cp = deepcopy(kf1)
        deleterows!(cp, 1:4)
        @test cp == KeyedFrame(DataFrame(; a=5:10, b=6:11, c=7:12), [:a, :b])

        cp = deepcopy(kf1)
        deleterows!(cp, [1, 10])
        @test cp == KeyedFrame(DataFrame(; a=2:9, b=3:10, c=4:11), [:a, :b])

        # Test return type of `deleterows!`
        @test isa(deleterows!(deepcopy(kf1), 1), KeyedFrame)
    end

    @testset "deletecols!" begin
        for ind in (:b, 2, [:b], [2])
            cp = deepcopy(kf1)
            deletecols!(cp, ind)
            @test cp == KeyedFrame(DataFrame(; a=1:10, c=3:12), [:a])
        end
        for ind in ([:a, :c], [1, 3])
            cp = deepcopy(kf1)
            deletecols!(cp, ind)
            @test cp == KeyedFrame(DataFrame(; b=2:11), [:b])
        end
        for ind in ([:a, :b], [1, 2])
            cp = deepcopy(kf1)
            deletecols!(cp, ind)
            @test cp == KeyedFrame(DataFrame(; c=3:12), Symbol[])
        end
        for ind in (:d, 4, [:a, :d], [1, 4])
            cp = deepcopy(kf1)
            @test_throws Exception deletecols!(cp, ind)
        end

        # Test return type of `deletecols!`
        @test isa(deletecols!(deepcopy(kf1), :b), KeyedFrame)
    end

    @testset "rename" begin
        initial = copy(kf1)
        expected = KeyedFrame(DataFrame(; new_a=1:10, b=2:11, new_c=3:12), [:new_a, :b])

        @test rename(initial, :a => :new_a, :c => :new_c) == expected
        @test initial == kf1
        @test rename(initial, [:a => :new_a, :c => :new_c]) == expected
        @test initial == kf1
        @test rename(initial, Dict(:a => :new_a, :c => :new_c)) == expected
        @test initial == kf1
        @test rename(x -> x == :b ? x : Symbol("new_$x"), initial) == expected
        @test initial == kf1

        @test rename!(initial, :a => :new_a, :c => :new_c) == expected
        @test initial == expected
        initial = copy(kf1)
        @test rename!(initial, [:a => :new_a, :c => :new_c]) == expected
        @test initial == expected
        initial = copy(kf1)
        @test rename!(initial, Dict(:a => :new_a, :c => :new_c)) == expected
        @test initial == expected
        initial = copy(kf1)
        @test rename!(x -> x == :b ? x : Symbol("new_$x"), initial) == expected
        @test initial == expected
    end

    @testset "unique" begin
        kf4 = KeyedFrame(DataFrame(; a=[1, 2, 3, 1, 2], b=[1, 2, 3, 4, 2], c=1:5), [:a, :b])

        @test nonunique(kf4) == [false, false, false, false, true]
        @test nonunique(kf4, :a) == [false, false, false, true, true]

        # Use default columns (key)
        expected = KeyedFrame(DataFrame(; a=[1, 2, 3, 1], b=1:4, c=1:4), [:a, :b])
        @test isequal(unique(kf4), expected)
        cp = deepcopy(kf4)
        unique!(cp)
        @test isequal(cp, expected)

        # Specify columns
        expected = KeyedFrame(DataFrame(; a=1:3, b=1:3, c=1:3), [:a, :b])
        @test isequal(unique(kf4, :a), expected)
        cp = deepcopy(kf4)
        unique!(cp, :a)
        @test isequal(cp, expected)

        # Test return type of `unique!`
        @test isa(unique!(deepcopy(kf1)), KeyedFrame)
        @test isa(unique!(deepcopy(kf1), :a), KeyedFrame)
    end

    @testset "join" begin
        expected = KeyedFrame(DataFrame(; a=1:5, b=2:6, c=3:7, d=4:8), [:a, :b])
        @test isequal(join(kf1, kf2), expected)
        @test isequal(join(kf1, df2), expected)               # Join a KeyedFrame and a DF
        @test isequal(join(df1, kf2), DataFrame(expected))    # Join a DF and a KeyedFrame

        @test isequal(join(kf1, kf2; kind=:right), expected)
        @test isequal(join(kf1, df2; kind=:right), expected)
        @test isequal(join(df1, kf2; kind=:right), DataFrame(expected))

        expected = KeyedFrame(DataFrame(; a=1:5, d=4:8, b=2:6, c=3:7), [:a, :b])
        @test isequal(join(kf2, kf1; kind=:left), expected)
        expected = KeyedFrame(DataFrame(; a=1:5, d=4:8, b=2:6, c=3:7), [:a])
        @test isequal(join(kf2, df1; kind=:left), expected)
        @test isequal(join(df2, kf1; kind=:left), DataFrame(expected))

        expected = KeyedFrame(
            DataFrame(; a=1:10, b=2:11, c=3:12, d=[4:8; fill(missing, 5)]), [:a, :b]
        )
        @test isequal(join(kf1, kf2; kind=:outer), expected)
        @test isequal(join(kf1, df2; kind=:outer), expected)
        @test isequal(join(df1, kf2; kind=:outer), DataFrame(expected))

        @test isequal(join(kf1, kf2; kind=:left), expected)
        @test isequal(join(kf1, df2; kind=:left), expected)
        @test isequal(join(df1, kf2; kind=:left), DataFrame(expected))

        expected = KeyedFrame(
            DataFrame(; a=1:10, d=[4:8; fill(missing, 5)], b=2:11, c=3:12), [:a, :b]
        )
        @test isequal(join(kf2, kf1; kind=:right), expected)
        expected = KeyedFrame(
            DataFrame(; a=1:10, d=[4:8; fill(missing, 5)], b=2:11, c=3:12), :a
        )
        @test isequal(join(kf2, df1; kind=:right), expected)
        expected = DataFrame(; a=1:10, d=[4:8; fill(missing, 5)], b=2:11, c=3:12)
        @test isequal(join(df2, kf1; kind=:right), expected)

        expected = KeyedFrame(DataFrame(; a=1:5, b=2:6, c=3:7), [:a, :b])
        @test isequal(join(kf1, kf2; kind=:semi), expected)
        @test isequal(join(kf1, df2; kind=:semi), expected)
        @test isequal(join(df1, kf2; kind=:semi), DataFrame(expected))

        expected = KeyedFrame(DataFrame(; a=1:5, d=4:8), :a)
        @test isequal(join(kf2, kf1; kind=:semi), expected)
        @test isequal(join(kf2, df1; kind=:semi), expected)
        @test isequal(join(df2, kf1; kind=:semi), DataFrame(expected))

        expected = KeyedFrame(DataFrame(; a=6:10, b=7:11, c=8:12), [:a, :b])
        @test isequal(join(kf1, kf2; kind=:anti), expected)
        @test isequal(join(kf1, df2; kind=:anti), expected)
        @test isequal(join(df1, kf2; kind=:anti), DataFrame(expected))

        expected = KeyedFrame(DataFrame(; a=[], d=[]), :a)
        @test isequal(join(kf2, kf1; kind=:anti), expected)
        @test isequal(join(kf2, df1; kind=:anti), expected)
        @test isequal(join(df2, kf1; kind=:anti), DataFrame(expected))

        expected = KeyedFrame(
            DataFrame(; a=[1, 2, 4], d=[4, 5, 7], e=[2, 5, 2], f=[3, 2, 1]), [:a, :e]
        )
        @test isequal(join(kf2, kf3), expected)

        expected = KeyedFrame(
            DataFrame(;
                a=[1, 2, 3, 4, 5, 4, 1],
                d=[4, 5, 6, 7, 8, 2, 2],
                f=[missing, 2, missing, missing, missing, 1, 3],
            ),
            :a,     # Key :e disappears, because it's renamed :d by the join
        )
        @test isequal(join(kf2, kf3; on=[:a => :a, :d => :e], kind=:outer), expected)
    end

    @testset "permutecols!" begin
        cp = deepcopy(kf1)
        permutecols!(cp, [1, 3, 2])

        @test isequal(cp, KeyedFrame(DataFrame(; a=1:10, c=3:12, b=2:11), [:a, :b]))
        permutecols!(cp, [2, 3, 1])
        @test isequal(cp, KeyedFrame(DataFrame(; c=3:12, b=2:11, a=1:10), [:a, :b]))

        @test_throws Exception permutecols!(cp, [1, 3])
        @test_throws Exception permutecols!(cp, [1, 2, 3, 4])

        # Test return type of `permutecols!`
        @test isa(permutecols!(deepcopy(kf1), [1, 2, 3]), KeyedFrame)
    end
end


text:
include("./chroma.jl")
include("./melspectrogram.jl")
include("./mfcc.jl")

using CSV, DataFrames, MP3, DSP, Statistics, Random, ProgressMeter

function shuffle_rows(tagsfile)

	annotations = CSV.File(tagsfile)
	new_order = shuffle(1:size(annotations)[1])
	annotations = annotations[new_order]
	CSV.write(tagsfile, annotations)

end

function preprocess_tags(tagsfile, pathcolumn, excludedcolumns; partitions=nothing)

	annotations = DataFrame(CSV.File(tagsfile))

	# Sort the tags by frequency...
	tag_frequencies = [
		column => sum(annotations[:, column])
		for column in names(annotations) if column ∉ [pathcolumn, excludedcolumns...]
	]
	sort!(tag_frequencies, lt = (pair1, pair2) -> pair1.second > pair2.second)

	#  ...and get the 50 most frequent tags
	select!(annotations, All(pathcolumn, [pair.first for pair in tag_frequencies[1:50]]))

	datasetsize, = size(annotations)
	filepaths = Array(select(annotations, pathcolumn))[:]
	tags = permutedims(Array(select!(annotations, Not(pathcolumn))))
	tagnames = names(annotations)

	# If a partition tuple is specified, divide the dataset in the given proportions
	if partitions isa Tuple
		partitionsizes = [datasetsize * partition ÷ sum(partitions)
			for partition in partitions]

		partitionends = cumsum(partitionsizes)
		partitionstarts = [0; partitionends[1:end-1]] .+ 1

		return (map(partitionstarts, partitionends) do pstart, pend
			return (filepaths[pstart:pend], tags[:, pstart:pend])
		end, tagnames)
	end

	return ((filepaths, tags), tagnames)

end

function preprocess_audio(feature_type, tagsfile, pathcolumn, samplespath, outputpath)

	annotations = [row[Symbol(pathcolumn)] for row in CSV.File(tagsfile, select=[pathcolumn])]
	progress = Progress(length(annotations))

	for track in annotations
		# Pre-processes the audio...
		audio = load(joinpath(samplespath, track))
		feature = feature_list[feature_type](audio.data[:], audio.samplerate)
		normalized_feature = normalize_feature(feature)

		# ...and write the result to file
		(trackdir, ) = splitdir(track)
		mkpath(joinpath(outputpath, trackdir))
		write(joinpath(outputpath, track[1:end-4]), normalized_feature)
		next!(progress)
	end

end

function normalize_feature(feature)
	featuremean = mean(feature[:])
	return Float32.(feature .- featuremean ./ std(feature, mean=featuremean) .+ eps(Float32))
end

feature_list = Dict(
	"chroma" => (signal, sr) ->
		chroma(signal, sr, hanning, Int(sr ÷ 20), 0, 0.1, 0.1),
	"melspectrogram" => (signal, sr) ->
		melspectrogram(signal, sr, hanning, Int(sr ÷ 20), 0, 128),
	"mfcc" => (signal, sr) ->
		mfcc(signal, sr, hanning, Int(sr ÷ 20), 0, 26),
)

max_stars_repo_name: gabriel-ss/music-autotagger

max_issues_repo_name: gabriel-ss/music-autotagger

max_forks_repo_name: gabriel-ss/music-autotagger


text:
for N in [Float64, Rational{Int}, Float32]
    c = N[1, 1]
    r = N[1, 1]
    H = Hyperrectangle(c, r)
    M = N[2 1; 1 -1]
    R = RotatedHyperrectangle(M, H)

    # dimension
    @test dim(R) == 2
    # center
    @test center(R) == N[3, 0]

    # generators
    gens = collect(generators(R))
    @test ngens(R) == 2
    @test length(gens) == 2 &&
        (N[2, 1] ∈ gens || N[-2, -1] ∈ gens) &&
        (N[1, -1] ∈ gens || N[-1, 1] ∈ gens)
    @test genmat(R) ∈ [N[2 1; 1 -1], N[-2 1; -1 -1], N[2 -1; 1 1], N[-2 -1; -1 1]]

    # support vector and support function
    d = N[1, 0]
    @test σ(d, R) == N[6, 0]
    @test ρ(d, R) == N(6)
    d = N[-1, 0]
    @test σ(d, R) == N[0, 0]
    @test ρ(d, R) == N(0)
    d = N[0, 1]
    @test σ(d, R) == N[4, 2]
    @test ρ(d, R) == N(2)
    d = N[0, -1]
    @test σ(d, R) == N[2, -2]
    @test ρ(d, R) == N(2)

    # boundedness
    @test isbounded(R)

    # is_polyhedral
    @test is_polyhedral(R)

    # isempty
    @test !isempty(R)

    # universality
    answer, w = isuniversal(R, true)
    @test !isuniversal(R) && !answer && w ∉ R

    # membership
    @test N[4, 1] ∈ R
    @test N[6, 1] ∉ R

    # an_element function
    @test an_element(R) ∈ R

    # linear map
    R2 = linear_map(M, R)
    @test R2 == RotatedHyperrectangle(M^2, H)
    R2 = linear_map(N[1 1;], R)
    @test R2 == RotatedHyperrectangle(N[3 0;], H)

    # list of vertices
    @test ispermutation(vertices_list(R),
                        [N[6, 0], N[2, -2], N[0, 0], N[4, 2]])

    if test_suite_polyhedra
        # list of constraints (only the length, see #2626)
        @test length(constraints_list(R)) == 4
    end
end


text:
using DataFrames
using MySQL

#------------------ BioMedQuery -------------------
    @testset "Testing Entrez" begin
    #testset "globals"
    narticles =10
    email=""
    ids = []
    efetch_dict = Dict()
    verbose = false
    articles = []

    try
        email = ENV["NCBI_EMAIL"]
    catch
        println("Entrez tests require the following enviroment variables:")
        println("NCBI_EMAIL")
    end

    @testset "Testing ESearch" begin
        println("-----------------------------------------")
        println("       Testing ESearch     ")

        search_term="obstructive sleep apnea[MeSH Major Topic]"
        search_dic = Dict("db"=>"pubmed","term" => search_term,
        "retstart" => 0, "retmax"=>narticles, "tool" =>"BioJulia",
        "email" => email)

        esearch_response = BioMedQuery.Entrez.esearch(search_dic)

        #convert xml to dictionary
        esearch_dict = BioMedQuery.Entrez.eparse(esearch_response)

        #examine how many ids were returned
        @test haskey(esearch_dict, "IdList")

        for id_node in esearch_dict["IdList"][1]["Id"]
            push!(ids, id_node)
        end

        @test length(ids)==narticles

    end

    @testset "Testing EFetch"    begin
        println("-----------------------------------------")
        println("       Testing EFetch      ")

        fetch_dic = Dict("db"=>"pubmed","tool" =>"BioJulia", "email" => email,
                         "retmode" => "xml", "rettype"=>"null")
        efetch_response = BioMedQuery.Entrez.efetch(fetch_dic, ids)
        efetch_dict = BioMedQuery.Entrez.eparse(efetch_response)

        @test haskey(efetch_dict, "PubmedArticle")

        articles = efetch_dict["PubmedArticle"]

        #articles should be an array of lenght narticles
        @test isa(articles, Array{Any, 1})
        @test length(articles) == narticles

        if true
            XMLconvert.xmlASCII2file(efetch_response, "./efetch.xml")
        end
    end

    @testset "Testing Citations" begin
        println("-----------------------------------------")
        println("       Testing Citations      ")
        #parse 1st article
        # art = BioMedQuery.Entrez.MedlineArticle(articles[1])
        # println(art)
        config = Dict(:type => "endnote", :output_file => "./citations_temp.endnote", :overwrite=>true)
        BioMedQuery.Entrez.save_article_citations(efetch_dict, config, verbose)

        #test that citations are the same as the ones already stored
        ref_lines=[]
        open("./pubmed_test.endnote") do f
           ref_lines = readlines(f)
        end

        temp_lines = []
        open("./citations_temp.endnote") do f
           temp_lines = readlines(f)
        end

        @test length(ref_lines)==length(temp_lines)

        for li=1:length(ref_lines)
            if ref_lines[li] != temp_lines[li]
                println("Line number: $li not equal")
                println(ref_lines[li])
                println("-----------------")
                println(temp_lines[li])
                @test 1==2
            end
        end
        rm("./citations_temp.endnote")
    end

    # save the results of an entrez fetch to a sqlite database
    @testset "Testing SQLite Saving" begin
        println("-----------------------------------------")
        println("       Testing SQLite Saving")

        db_path = "./test_db.db"

        config = Dict(:db_path=> db_path, :overwrite=>true)
        db = BioMedQuery.Entrez.save_efetch_sqlite(efetch_dict, config, verbose)

        #query the article table and make sure the count is correct
        all_pmids = BioMedQuery.Entrez.DB.all_pmids(db)
        @test length(all_pmids) == narticles

        #check we can get the MESH descriptor for an article
        mesh = BioMedQuery.Entrez.DB.get_article_mesh(db, all_pmids[1])
        @test length(mesh) > 0

        #check that reminder of tables are not empty
        tables = ["author", "author2article", "mesh_descriptor",
        "mesh_qualifier", "mesh_heading"]

        for t in tables
            query_str = "SELECT count(*) FROM "*t*";"
            q = BioMedQuery.DBUtils.db_query(db, query_str)
            count = get(q[1][1])
            @test count > 0
        end

        # remove temp files
        if isfile(db_path)
            rm(db_path)
        end
    end

    @testset "Testing MySQL Saving" begin
        println("-----------------------------------------")
        println("       Testing MySQL Saving")

        dbname = "entrez_test"
        config = Dict(:host=>"localhost", :dbname=>dbname, :username=>"root",
        :pswd=>"", :overwrite=>true)
        @time db = BioMedQuery.Entrez.save_efetch_mysql(efetch_dict, config, verbose)

        #query the article table and make sure the count is correct
        all_pmids = BioMedQuery.Entrez.DB.all_pmids(db)
        @test length(all_pmids) == narticles

        #check we can get the MESH descriptor for an article
        mesh = BioMedQuery.Entrez.DB.get_article_mesh(db, all_pmids[1])
        println(mesh)
        @test length(mesh) > 0

        #check that reminder of tables are not empty
        tables = ["author", "author2article", "mesh_descriptor",
        "mesh_qualifier", "mesh_heading"]

        for t in tables
            query_str = "SELECT count(*) FROM "*t*";"
            q = BioMedQuery.DBUtils.db_query(db, query_str)
            count = q[1][1]
            @test count > 0
        end

        #clean-up
        db_query(db, "DROP DATABASE IF EXISTS $dbname;")

    end

    @testset "Testing ELink" begin
        println("-----------------------------------------")
        println("       Testing ELink")

        pmid = "19304878"
        elink_dict = Dict("dbfrom" =>"pubmed", "id" => pmid,
                          "linkname" => "pubmed_pubmed", "email"=>email)
        elink_response = BioMedQuery.Entrez.elink(elink_dict)

        elink_response_dict = BioMedQuery.Entrez.eparse(elink_response)

        @test haskey( elink_response_dict, "LinkSet")
    end

    @testset "Testing ESummary" begin
        println("-----------------------------------------")
        println("       Testing ESummary")

        pmid = "30367"
        esummary_dict = Dict("db" =>"pubmed", "id" => pmid, "email"=>email)
        esummary_response = BioMedQuery.Entrez.esummary(esummary_dict)

        esummary_response_dict = BioMedQuery.Entrez.eparse(esummary_response)

        @test haskey( esummary_response_dict, "DocSum")
    end

    println("------------End Test Entrez--------------")
    println("-----------------------------------------")

end


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

JLLWrappers.@generate_main_file_header("Telescope")
JLLWrappers.@generate_main_file("Telescope", UUID("ae456e04-aaf2-55ce-80c6-09193969412e"))
end  # module Telescope_jll


text:
function bayes_optimization(n, f, fdata, lb, ub, x, minf, params)
    cf = cfunction(f, Cdouble, (Cuint, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Void}))
    ret = ccall(
        (:bayes_optimization, libbayesopt),
        Cint,
        (Cuint, Ptr{Void}, Ptr{Void},
         Ptr{Cdouble}, Ptr{Cdouble},
         Ptr{Cdouble}, Ptr{Cdouble},
         bopt_params),
        n, cf, fdata, lb, ub, x, minf, params)
    return ret
end


text:
function jump_nlp_01(optimizer; mode = BilevelJuMP.ProductMode(1e-8), config = Config())

    atol = config.atol

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)
    BilevelJuMP.set_copy_names(model)

    @variable(Upper(model), x >= 0, start = 2)
    @variable(Lower(model), y >= 0, start = 0)

    @objective(Upper(model), Max, x -3y)

    @NLconstraint(Upper(model), x^2 <= 4)

    @objective(Lower(model), Min, y)

    @constraint(Lower(model), y >= 0)

    optimize!(model)

    primal_status(model)

    @test termination_status(model) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED]

    @test objective_value(model) ≈ 2 atol=atol

    @test value(x) ≈  2 atol=atol
    @test value(y) ≈  0 atol=atol

end

function jump_nlp_02(optimizer; mode = BilevelJuMP.ProductMode(1e-8), config = Config())

    atol = config.atol

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)
    BilevelJuMP.set_copy_names(model)

    @variable(Upper(model), x >= 0, start = 2)
    @variable(Lower(model), y >= 0, start = 0)

    @objective(Upper(model), Max, x -3y)

    @NLconstraint(Upper(model), x^2 <= 4)
    @NLconstraint(Upper(model), x^2 + y + 2 >= 4)
    @NLconstraint(Upper(model), x^3 + 18*y >= 1)
    @NLconstraint(Upper(model), -1 <= x*y <= 1)

    @objective(Lower(model), Min, y)

    @constraint(Lower(model), y >= 0)

    optimize!(model)

    primal_status(model)

    @test termination_status(model) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED]

    @test objective_value(model) ≈ 2 atol=atol

    @test value(x) ≈  2 atol=atol
    @test value(y) ≈  0 atol=atol

end

function jump_nlp_03(optimizer; mode = BilevelJuMP.ProductMode(1e-8), config = Config())

    atol = config.atol

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)
    BilevelJuMP.set_copy_names(model)

    @variable(Upper(model), x >= 0, start = 2)
    @variable(Lower(model), y >= 0, start = 0)
    @variable(LowerOnly(model), z >= 0, start = 0)

    @objective(Upper(model), Max, x)
    # TODO NL obj
    @NLobjective(Upper(model), Max, x^2 -3y)

    @test_throws ErrorException @NLconstraint(Lower(model), x^2 + y + 2 >= 4)
    @NLconstraint(Upper(model), x^2 <= 4)
    @test_throws ErrorException @NLconstraint(Lower(model), x^3 + 18*y >= 1)
    @NLconstraint(Upper(model), -1 <= x*y <= 1)
    @test_throws ErrorException @NLconstraint(Upper(model), -1 <= z*y <= 1)

    @test_throws ErrorException @NLobjective(Lower(model), Min, y^3)
    @objective(Lower(model), Min, y^2 + z)

    @constraint(Lower(model), y >= 0)

    optimize!(model)

    primal_status(model)

    @test termination_status(model) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED]

    @test objective_value(model) ≈ 4 atol=atol

    @test value(x) ≈  2 atol=atol
    @test value(y) ≈  0 atol=atol

end

function jump_nlp_04(optimizer; mode = BilevelJuMP.ProductMode(1e-8), config = Config())

    atol = config.atol

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)
    BilevelJuMP.set_copy_names(model)

    @variable(Upper(model), x >= 0, start = 2)
    @variable(Lower(model), y >= 0, start = 0)

    @test_throws ErrorException @NLparameter(model, rhs == 4.0)
    @test_throws ErrorException @NLparameter(Lower(model), rhs == 4.0)
    @NLparameter(Upper(model), rhs == 4.0)

    @objective(Upper(model), Max, x -3y)

    @NLconstraint(Upper(model), x^2 <= rhs)

    @objective(Lower(model), Min, y)

    @constraint(Lower(model), y >= 0)

    optimize!(model)

    primal_status(model)

    @test termination_status(model) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED]

    @test objective_value(model) ≈ 2 atol=atol

    @test value(x) ≈  2 atol=atol
    @test value(y) ≈  0 atol=atol

    #=
        change param
    =#

    JuMP.set_value(rhs, 9.0)

    optimize!(model)

    primal_status(model)

    @test termination_status(model) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED]

    @test objective_value(model) ≈ 3 atol=atol

    @test value(x) ≈  3 atol=atol
    @test value(y) ≈  0 atol=atol

end


text:
abstract RandomProcessCarrier
abstract QuasiRandomCarrier
abstract IterativeMapCarrier <: QuasiRandomCarrier
abstract FSPWLCarrier <: IterativeMapCarrier # Fully-stretching piece-wise linear maps

typealias Carrier Union(Distribution, RandomProcessCarrier, QuasiRandomCarrier)