text:
"""
    LinearMixedModel

Linear mixed-effects model representation

## Fields

* `formula`: the formula for the model
* `reterms`: a `Vector{AbstractReMat{T}}` of random-effects terms.
* `Xymat`: horizontal concatenation of a full-rank fixed-effects model matrix `X` and response `y` as an `FeMat{T}`
* `feterm`: the fixed-effects model matrix as an `FeTerm{T}`
* `sqrtwts`: vector of square roots of the case weights.  Can be empty.
* `parmap` : Vector{NTuple{3,Int}} of (block, row, column) mapping of θ to λ
* `dims` : NamedTuple{(:n, :p, :nretrms),NTuple{3,Int}} of dimensions.  `p` is the rank of `X`, which may be smaller than `size(X, 2)`.
* `A`: a `Vector{AbstractMatrix}` containing the row-major packed lower triangle of `hcat(Z,X,y)'hcat(Z,X,y)`
* `L`: the blocked lower Cholesky factor of `Λ'AΛ+I` in the same Vector representation as `A`
* `optsum`: an [`OptSummary`](@ref) object

## Properties

* `θ` or `theta`: the covariance parameter vector used to form λ
* `β` or `beta`: the fixed-effects coefficient vector
* `λ` or `lambda`: a vector of lower triangular matrices repeated on the diagonal blocks of `Λ`
* `σ` or `sigma`: current value of the standard deviation of the per-observation noise
* `b`: random effects on the original scale, as a vector of matrices
* `u`: random effects on the orthogonal scale, as a vector of matrices
* `lowerbd`: lower bounds on the elements of θ
* `X`: the fixed-effects model matrix
* `y`: the response vector
"""
struct LinearMixedModel{T<:AbstractFloat} <: MixedModel{T}
    formula::FormulaTerm
    reterms::Vector{AbstractReMat{T}}
    Xymat::FeMat{T}
    feterm::FeTerm{T}
    sqrtwts::Vector{T}
    parmap::Vector{NTuple{3,Int}}
    dims::NamedTuple{(:n, :p, :nretrms),NTuple{3,Int}}
    A::Vector{AbstractMatrix{T}}            # cross-product blocks
    L::Vector{AbstractMatrix{T}}
    optsum::OptSummary{T}
end

function LinearMixedModel(
    f::FormulaTerm, tbl; contrasts=Dict{Symbol,Any}(), wts=[], σ=nothing
)
    return LinearMixedModel(f::FormulaTerm, Tables.columntable(tbl); contrasts, wts, σ)
end

function LinearMixedModel(
    f::FormulaTerm, tbl::Tables.ColumnTable; contrasts=Dict{Symbol,Any}(), wts=[], σ=nothing
)
    # TODO: perform missing_omit() after apply_schema() when improved
    # missing support is in a StatsModels release
    tbl, _ = StatsModels.missing_omit(tbl, f)
    sch = try
        schema(f, tbl, contrasts)
    catch e
        if isa(e, OutOfMemoryError)
            @warn "Random effects grouping variables with many levels can cause out-of-memory errors.  Try manually specifying `Grouping()` contrasts for those variables."
        end
        rethrow(e)
    end
    form = apply_schema(f, sch, LinearMixedModel)
    # tbl, _ = StatsModels.missing_omit(tbl, form)

    y, Xs = modelcols(form, tbl)

    return LinearMixedModel(y, Xs, form, wts, σ)
end

"""
    LinearMixedModel(y, Xs, form, wts=[], σ=nothing)

Private constructor for a LinearMixedModel.

To construct a model, you only need the response (`y`), already assembled
model matrices (`Xs`), schematized formula (`form`) and weights (`wts`).
Everything else in the structure can be derived from these quantities.

!!! note
    This method is internal and experimental and so may change or disappear in
    a future release without being considered a breaking change.
"""
function LinearMixedModel(
    y::AbstractArray,
    Xs::Tuple, # can't be more specific here without stressing the compiler
    form::FormulaTerm,
    wts=[],
    σ=nothing,
)
    T = promote_type(Float64, float(eltype(y)))  # ensure eltype of model matrices is at least Float64

    reterms = AbstractReMat{T}[]
    feterms = FeTerm{T}[]
    for (i, x) in enumerate(Xs)
        if isa(x, AbstractReMat{T})
            push!(reterms, x)
        elseif isa(x, ReMat) # this can occur in weird situation where x is a ReMat{U}
            # avoid keeping a second copy if unweighted
            z = convert(Matrix{T}, x.z)
            wtz = x.z === x.wtz ? z : convert(Matrix{T}, x.wtz)
            S = size(z, 1)
            x = ReMat{T,S}(
                x.trm,
                x.refs,
                x.levels,
                x.cnames,
                z,
                wtz,
                convert(LowerTriangular{Float64,Matrix{Float64}}, x.λ),
                x.inds,
                convert(SparseMatrixCSC{T,Int32}, x.adjA),
                convert(Matrix{T}, x.scratch),
            )
            push!(reterms, x)
        else
            cnames = coefnames(form.rhs[i])
            push!(feterms, FeTerm(x, isa(cnames, String) ? [cnames] : collect(cnames)))
        end
    end
    return LinearMixedModel(convert(Array{T}, y), only(feterms), reterms, form, wts, σ)
end

"""
    LinearMixedModel(y, feterm, reterms, form, wts=[], σ=nothing)

Private constructor for a `LinearMixedModel` given already assembled fixed and random effects.

To construct a model, you only need a vector of `FeMat`s (the fixed-effects
model matrix and response), a vector of `AbstractReMat` (the random-effects
model matrices), the formula and the weights. Everything else in the structure
can be derived from these quantities.

!!! note
    This method is internal and experimental and so may change or disappear in
    a future release without being considered a breaking change.
"""
function LinearMixedModel(
    y::AbstractArray,
    feterm::FeTerm{T},
    reterms::AbstractVector{<:AbstractReMat{T}},
    form::FormulaTerm,
    wts=[],
    σ=nothing,
) where {T}
    # detect and combine RE terms with the same grouping var
    if length(reterms) > 1
        reterms = amalgamate(reterms)
    end

    sort!(reterms; by=nranef, rev=true)
    Xy = FeMat(feterm, vec(y))
    sqrtwts = sqrt.(convert(Vector{T}, wts))
    reweight!.(reterms, Ref(sqrtwts))
    reweight!(Xy, sqrtwts)
    A, L = createAL(reterms, Xy)
    lbd = foldl(vcat, lowerbd(c) for c in reterms)
    θ = foldl(vcat, getθ(c) for c in reterms)
    optsum = OptSummary(θ, lbd, :LN_BOBYQA; ftol_rel=T(1.0e-12), ftol_abs=T(1.0e-8))
    optsum.sigma = isnothing(σ) ? nothing : T(σ)
    fill!(optsum.xtol_abs, 1.0e-10)
    return LinearMixedModel(
        form,
        reterms,
        Xy,
        feterm,
        sqrtwts,
        mkparmap(reterms),
        (n=length(y), p=feterm.rank, nretrms=length(reterms)),
        A,
        L,
        optsum,
    )
end

function fit(
    ::Type{LinearMixedModel},
    f::FormulaTerm,
    tbl;
    wts=[],
    contrasts=Dict{Symbol,Any}(),
    progress::Bool=true,
    REML::Bool=false,
    σ=nothing,
    thin=typemax(Int),
)
    return fit(
        LinearMixedModel,
        f,
        Tables.columntable(tbl);
        wts,
        contrasts,
        progress,
        REML,
        σ,
        thin,
    )
end

function fit(
    ::Type{LinearMixedModel},
    f::FormulaTerm,
    tbl::Tables.ColumnTable;
    wts=[],
    contrasts=Dict{Symbol,Any}(),
    progress=true,
    REML=false,
    σ=nothing,
    thin=typemax(Int),
)
    return fit!(LinearMixedModel(f, tbl; contrasts, wts, σ); progress, REML, thin)
end

function _offseterr()
    return throw(
        ArgumentError(
            "Offsets are not supported in linear models. You can simply shift the response by the offset.",
        ),
    )
end

function fit(
    ::Type{MixedModel},
    f::FormulaTerm,
    tbl;
    wts=[],
    contrasts=Dict{Symbol,Any}(),
    progress::Bool=true,
    REML::Bool=false,
    offset=[],
    σ=nothing,
    thin=typemax(Int),
)
    return if !isempty(offset)
        _offseterr()
    else
        fit(LinearMixedModel, f, tbl; wts, contrasts, progress, REML, σ, thin)
    end
end

function fit(
    ::Type{MixedModel},
    f::FormulaTerm,
    tbl,
    d::Normal,
    l::IdentityLink;
    wts=[],
    contrasts=Dict{Symbol,Any}(),
    progress::Bool=true,
    REML::Bool=false,
    offset=[],
    σ=nothing,
    thin=typemax(Int),
    fast=nothing,
    nAGQ=nothing,
)
    return if !isempty(offset)
        _offseterr()
    else
        if !isnothing(fast) || !isnothing(nAGQ)
            @warn "fast and nAGQ arguments are ignored when fitting a LinearMixedModel"
        end
        fit(LinearMixedModel, f, tbl; wts, contrasts, progress, REML, σ, thin)
    end
end

function StatsBase.coef(m::LinearMixedModel{T}) where {T}
    piv = m.feterm.piv
    return invpermute!(fixef!(similar(piv, T), m), piv)
end

βs(m::LinearMixedModel) = NamedTuple{(Symbol.(coefnames(m))...,)}(coef(m))

function StatsBase.coefnames(m::LinearMixedModel)
    Xtrm = m.feterm
    return invpermute!(copy(Xtrm.cnames), Xtrm.piv)
end

function StatsBase.coeftable(m::LinearMixedModel)
    co = coef(m)
    se = stderror!(similar(co), m)
    z = co ./ se
    pvalue = ccdf.(Chisq(1), abs2.(z))
    names = coefnames(m)

    return CoefTable(
        hcat(co, se, z, pvalue),
        ["Coef.", "Std. Error", "z", "Pr(>|z|)"],
        names,
        4, # pvalcol
        3, # teststatcol
    )
end

"""
    condVar(m::LinearMixedModel)

Return the conditional variances matrices of the random effects.

The random effects are returned by `ranef` as a vector of length `k`,
where `k` is the number of random effects terms.  The `i`th element
is a matrix of size `vᵢ × ℓᵢ`  where `vᵢ` is the size of the
vector-valued random effects for each of the `ℓᵢ` levels of the grouping
factor.  Technically those values are the modes of the conditional
distribution of the random effects given the observed data.

This function returns an array of `k` three dimensional arrays,
where the `i`th array is of size `vᵢ × vᵢ × ℓᵢ`.  These are the
diagonal blocks from the conditional variance-covariance matrix,

    s² Λ(Λ'Z'ZΛ + I)⁻¹Λ'
"""
function condVar(m::LinearMixedModel{T}) where {T}
    return [condVar(m, fnm) for fnm in fnames(m)]
end

function condVar(m::LinearMixedModel{T}, fname) where {T}
    Lblk = LowerTriangular(densify(sparseL(m; fname=fname)))
    blk = findfirst(isequal(fname), fnames(m))
    λt = Array(m.λ[blk]') .* sdest(m)
    vsz = size(λt, 2)
    ℓ = length(m.reterms[blk].levels)
    val = Array{T}(undef, (vsz, vsz, ℓ))
    scratch = Matrix{T}(undef, (size(Lblk, 1), vsz))
    for b in 1:ℓ
        fill!(scratch, zero(T))
        copyto!(view(scratch, (b - 1) * vsz .+ (1:vsz), :), λt)
        ldiv!(Lblk, scratch)
        mul!(view(val, :, :, b), scratch', scratch)
    end
    return val
end

function _cvtbl(arr::Array{T,3}, trm) where {T}
    return merge(
        NamedTuple{(fname(trm),)}((trm.levels,)),
        columntable([
            NamedTuple{(:σ, :ρ)}(sdcorr(view(arr, :, :, i))) for i in axes(arr, 3)
        ]),
    )
end

"""
    condVartables(m::LinearMixedModel)

Return the conditional covariance matrices of the random effects as a `NamedTuple` of columntables
"""
function condVartables(m::MixedModel{T}) where {T}
    return NamedTuple{fnames(m)}((map(_cvtbl, condVar(m), m.reterms)...,))
end

function _pushALblock!(A, L, blk)
    push!(L, blk)
    return push!(A, deepcopy(isa(blk, BlockedSparse) ? blk.cscmat : blk))
end

function createAL(reterms::Vector{AbstractReMat{T}}, Xy::FeMat{T}) where {T}
    k = length(reterms)
    vlen = kchoose2(k + 1)
    A = sizehint!(AbstractMatrix{T}[], vlen)
    L = sizehint!(AbstractMatrix{T}[], vlen)
    for i in eachindex(reterms)
        for j in 1:i
            _pushALblock!(A, L, densify(reterms[i]' * reterms[j]))
        end
    end
    for j in eachindex(reterms)   # can't fold this into the previous loop b/c block order
        _pushALblock!(A, L, densify(Xy' * reterms[j]))
    end
    _pushALblock!(A, L, densify(Xy'Xy))
    for i in 2:k      # check for fill-in due to non-nested grouping factors
        ci = reterms[i]
        for j in 1:(i - 1)
            cj = reterms[j]
            if !isnested(cj, ci)
                for l in i:k
                    ind = block(l, i)
                    L[ind] = Matrix(L[ind])
                end
                break
            end
        end
    end
    return A, L
end

StatsBase.deviance(m::LinearMixedModel) = objective(m)

GLM.dispersion(m::LinearMixedModel, sqr::Bool=false) = sqr ? varest(m) : sdest(m)

GLM.dispersion_parameter(m::LinearMixedModel) = true

"""
    feL(m::LinearMixedModel)

Return the lower Cholesky factor for the fixed-effects parameters, as an `LowerTriangular`
`p × p` matrix.
"""
function feL(m::LinearMixedModel)
    XyL = m.L[end]
    k = size(XyL, 1)
    inds = Base.OneTo(k - 1)
    return LowerTriangular(view(XyL, inds, inds))
end

"""
    fit!(m::LinearMixedModel; progress::Bool=true, REML::Bool=false,
                              σ::Union{Real, Nothing}=nothing,
                              thin::Int=typemax(Int))

Optimize the objective of a `LinearMixedModel`.  When `progress` is `true` a
`ProgressMeter.ProgressUnknown` display is shown during the optimization of the
objective, if the optimization takes more than one second or so.

At every `thin`th iteration  is recorded in `fitlog`, optimization progress is
saved in `m.optsum.fitlog`.
"""
function fit!(
    m::LinearMixedModel{T};
    progress::Bool=true,
    REML::Bool=false,
    σ::Union{Real,Nothing}=nothing,
    thin::Int=typemax(Int),
) where {T}
    optsum = m.optsum
    # this doesn't matter for LMM, but it does for GLMM, so let's be consistent
    if optsum.feval > 0
        throw(ArgumentError("This model has already been fitted. Use refit!() instead."))
    end
    opt = Opt(optsum)
    optsum.REML = REML
    prog = ProgressUnknown("Minimizing"; showspeed=true)
    # start from zero for the initial call to obj before optimization
    iter = 0
    fitlog = optsum.fitlog
    function obj(x, g)
        isempty(g) || throw(ArgumentError("g should be empty for this objective"))
        val = objective(updateL!(setθ!(m, x)))
        iszero(rem(iter, thin)) && push!(fitlog, (copy(x), val))
        progress && ProgressMeter.next!(prog; showvalues=[(:objective, val)])
        iter += 1
        return val
    end
    NLopt.min_objective!(opt, obj)
    try
        optsum.finitial = obj(optsum.initial, T[])
    catch PosDefException
        if all(==(first(m.y)), m.y)
            throw(
                ArgumentError("The response is constant and thus model fitting has failed")
            )
        else
            rethrow()
        end
    end
    empty!(fitlog)
    push!(fitlog, (copy(optsum.initial), optsum.finitial))
    fmin, xmin, ret = NLopt.optimize!(opt, copyto!(optsum.final, optsum.initial))
    ProgressMeter.finish!(prog)
    ## check if small non-negative parameter values can be set to zero
    xmin_ = copy(xmin)
    lb = optsum.lowerbd
    for i in eachindex(xmin_)
        if iszero(lb[i]) && zero(T) < xmin_[i] < T(0.001)
            xmin_[i] = zero(T)
        end
    end
    loglength = length(fitlog)
    if xmin ≠ xmin_
        if (zeroobj = obj(xmin_, T[])) ≤ (fmin + 1.e-5)
            fmin = zeroobj
            copyto!(xmin, xmin_)
        elseif length(fitlog) > loglength
            # remove unused extra log entry
            pop!(fitlog)
        end
    end
    ## ensure that the parameter values saved in m are xmin
    updateL!(setθ!(m, xmin))

    optsum.feval = opt.numevals
    optsum.final = xmin
    optsum.fmin = fmin
    optsum.returnvalue = ret
    _check_nlopt_return(ret)
    return m
end

function fitted!(v::AbstractArray{T}, m::LinearMixedModel{T}) where {T}
    ## FIXME: Create and use `effects(m) -> β, b` w/o calculating β twice
    Xtrm = m.feterm
    vv = mul!(vec(v), Xtrm.x, fixef!(similar(Xtrm.piv, T), m))
    for (rt, bb) in zip(m.reterms, ranef(m))
        mul!(vv, rt, bb, one(T), one(T))
    end
    return v
end

StatsBase.fitted(m::LinearMixedModel{T}) where {T} = fitted!(Vector{T}(undef, nobs(m)), m)

"""
    fixef!(v::Vector{T}, m::MixedModel{T})

Overwrite `v` with the pivoted fixed-effects coefficients of model `m`

For full-rank models the length of `v` must be the rank of `X`.  For rank-deficient models
the length of `v` can be the rank of `X` or the number of columns of `X`.  In the latter
case the calculated coefficients are padded with -0.0 out to the number of columns.
"""
function fixef!(v::AbstractVector{Tv}, m::LinearMixedModel{T}) where {Tv,T}
    Xtrm = m.feterm
    fill!(v, -zero(Tv))
    XyL = m.L[end]
    L = feL(m)
    k = size(XyL, 1)
    r = size(L, 1)
    for j in 1:r
        v[j] = XyL[k, j]
    end
    ldiv!(feL(m)', length(v) == r ? v : view(v, 1:r))
    return v
end

"""
    fixef(m::MixedModel)

Return the fixed-effects parameter vector estimate of `m`.

In the rank-deficient case the truncated parameter vector, of length `rank(m)` is returned.
This is unlike `coef` which always returns a vector whose length matches the number of
columns in `X`.
"""
fixef(m::LinearMixedModel{T}) where {T} = fixef!(Vector{T}(undef, m.feterm.rank), m)

"""
    fixefnames(m::MixedModel)

Return a (permuted and truncated in the rank-deficient case) vector of coefficient names.
"""
function fixefnames(m::LinearMixedModel{T}) where {T}
    Xtrm = m.feterm
    return Xtrm.cnames[1:(Xtrm.rank)]
end

"""
    fnames(m::MixedModel)

Return the names of the grouping factors for the random-effects terms.
"""
fnames(m::MixedModel) = (fname.(m.reterms)...,)

"""
    getθ(m::LinearMixedModel)

Return the current covariance parameter vector.
"""
getθ(m::LinearMixedModel{T}) where {T} = getθ!(Vector{T}(undef, length(m.parmap)), m)

function getθ!(v::AbstractVector{Tv}, m::LinearMixedModel{T}) where {Tv,T}
    pmap = m.parmap
    if length(v) ≠ length(pmap)
        throw(
            DimensionMismatch(
                "length(v) = $(length(v)) ≠ length(m.parmap) = $(length(pmap))"
            ),
        )
    end
    reind = 1
    λ = first(m.reterms).λ
    for (k, tp) in enumerate(pmap)
        tp1 = first(tp)
        if reind ≠ tp1
            reind = tp1
            λ = m.reterms[tp1].λ
        end
        v[k] = λ[tp[2], tp[3]]
    end
    return v
end

function Base.getproperty(m::LinearMixedModel{T}, s::Symbol) where {T}
    if s == :θ || s == :theta
        getθ(m)
    elseif s == :β || s == :beta
        coef(m)
    elseif s == :βs || s == :betas
        βs(m)
    elseif s == :λ || s == :lambda
        getproperty.(m.reterms, :λ)
    elseif s == :σ || s == :sigma
        sdest(m)
    elseif s == :σs || s == :sigmas
        σs(m)
    elseif s == :σρs || s == :sigmarhos
        σρs(m)
    elseif s == :b
        ranef(m)
    elseif s == :objective
        objective(m)
    elseif s == :corr
        vcov(m; corr=true)
    elseif s == :vcov
        vcov(m; corr=false)
    elseif s == :PCA
        NamedTuple{fnames(m)}(PCA.(m.reterms))
    elseif s == :pvalues
        ccdf.(Chisq(1), abs2.(coef(m) ./ stderror(m)))
    elseif s == :stderror
        stderror(m)
    elseif s == :u
        ranef(m; uscale=true)
    elseif s == :lowerbd
        m.optsum.lowerbd
    elseif s == :X
        modelmatrix(m)
    elseif s == :y
        let xy = m.Xymat.xy
            view(xy, :, size(xy, 2))
        end
    elseif s == :rePCA
        rePCA(m)
    else
        getfield(m, s)
    end
end

StatsBase.islinear(m::LinearMixedModel) = true

"""
    _3blockL(::LinearMixedModel)

returns L in 3-block form:
- a Diagonal or UniformBlockDiagonal block
- a dense rectangular block
- and a dense lowertriangular block
"""
function _3blockL(m::LinearMixedModel{T}) where {T}
    L = m.L
    reterms = m.reterms
    isone(length(reterms)) &&
        return first(L), L[block(2, 1)], LowerTriangular(L[block(2, 2)])
    rows = sum(k -> size(L[kp1choose2(k + 1)], 1), axes(reterms, 1))
    cols = size(first(L), 2)
    B2 = Matrix{T}(undef, (rows, cols))
    B3 = Matrix{T}(undef, (rows, rows))
    rowoffset = 0
    for i in 1 .+ axes(reterms, 1)
        Li1 = L[block(i, 1)]
        rows = rowoffset .+ axes(Li1, 1)
        copyto!(view(B2, rows, :), Li1)
        coloffset = 0
        for j in 2:i
            Lij = L[block(i, j)]
            copyto!(view(B3, rows, coloffset .+ axes(Lij, 2)), Lij)
            coloffset += size(Lij, 2)
        end
        rowoffset += size(Li1, 1)
    end
    return first(L), B2, LowerTriangular(B3)
end

# use dispatch to distinguish Diagonal and UniformBlockDiagonal in first(L)
_ldivB1!(B1::Diagonal{T}, rhs::AbstractVector{T}, ind) where {T} = rhs ./= B1.diag[ind]
function _ldivB1!(B1::UniformBlockDiagonal{T}, rhs::AbstractVector{T}, ind) where {T}
    return ldiv!(LowerTriangular(view(B1.data, :, :, ind)), rhs)
end

"""
    leverage(::LinearMixedModel)

Return the diagonal of the hat matrix of the model.

For a linear model, the sum of the leverage values is the degrees of freedom
for the model in the sense that this sum is the dimension of the span of columns
of the model matrix.  With a bit of hand waving a similar argument could be made
for linear mixed-effects models. The hat matrix is of the form ``[ZΛ X][L L']⁻¹[ZΛ X]'``.
"""
function StatsBase.leverage(m::LinearMixedModel{T}) where {T}
    # To obtain the diagonal elements solve L⁻¹[ZΛ X]'eⱼ
    # where eⱼ is the j'th basis vector in Rⁿ and evaluate the squared length of the solution.
    # The fact that the [1,1] block of L is always UniformBlockDiagonal
    # or Diagonal makes it easy to obtain the first chunk of the solution.
    B1, B2, B3 = _3blockL(m)
    reterms = m.reterms
    re1 = first(reterms)
    re1z = re1.z
    r1sz = size(re1z, 1)
    re1λ = re1.λ
    re1refs = re1.refs
    Xy = m.Xymat
    rhs1 = zeros(T, size(re1z, 1))   # for the first block only the nonzeros are stored
    rhs2 = zeros(T, size(B2, 1))
    value = similar(m.y)
    for i in eachindex(value)
        re1ind = re1refs[i]
        _ldivB1!(B1, mul!(rhs1, adjoint(re1λ), view(re1z, :, i)), re1ind)
        off = (re1ind - 1) * r1sz
        fill!(rhs2, 0)
        rhsoffset = 0
        for j in 2:length(reterms)
            trm = reterms[j]
            z = trm.z
            stride = size(z, 1)
            mul!(
                view(rhs2, (rhsoffset + (trm.refs[i] - 1) * stride) .+ Base.OneTo(stride)),
                adjoint(trm.λ),
                view(z, :, i),
            )
            rhsoffset += length(trm.levels) * stride
        end
        copyto!(view(rhs2, rhsoffset .+ Base.OneTo(size(Xy, 2))), view(Xy, i, :))
        ldiv!(B3, mul!(rhs2, view(B2, :, off .+ Base.OneTo(r1sz)), rhs1, 1, -1))
        rhs2[end] = 0
        value[i] = sum(abs2, rhs1) + sum(abs2, rhs2)
    end
    return value
end

function StatsBase.loglikelihood(m::LinearMixedModel)
    if m.optsum.REML
        throw(ArgumentError("loglikelihood not available for models fit by REML"))
    end
    return -objective(m) / 2
end

lowerbd(m::LinearMixedModel) = m.optsum.lowerbd

function mkparmap(reterms::Vector{AbstractReMat{T}}) where {T}
    parmap = NTuple{3,Int}[]
    for (k, trm) in enumerate(reterms)
        n = LinearAlgebra.checksquare(trm.λ)
        for ind in trm.inds
            d, r = divrem(ind - 1, n)
            push!(parmap, (k, r + 1, d + 1))
        end
    end
    return parmap
end

nθ(m::LinearMixedModel) = length(m.parmap)

"""
    objective(m::LinearMixedModel)

Return negative twice the log-likelihood of model `m`
"""
function objective(m::LinearMixedModel{T}) where {T}
    wts = m.sqrtwts
    denomdf = T(ssqdenom(m))
    σ = m.optsum.sigma
    val = if isnothing(σ)
        logdet(m) + denomdf * (one(T) + log2π + log(pwrss(m) / denomdf))
    else
        denomdf * (log2π + 2 * log(σ)) + logdet(m) + pwrss(m) / σ^2
    end
    return isempty(wts) ? val : val - T(2.0) * sum(log, wts)
end

function Base.propertynames(m::LinearMixedModel, private::Bool=false)
    return (
        fieldnames(LinearMixedModel)...,
        :θ,
        :theta,
        :β,
        :beta,
        :βs,
        :betas,
        :λ,
        :lambda,
        :stderror,
        :σ,
        :sigma,
        :σs,
        :sigmas,
        :σρs,
        :sigmarhos,
        :b,
        :u,
        :lowerbd,
        :X,
        :y,
        :corr,
        :vcov,
        :PCA,
        :rePCA,
        :objective,
        :pvalues,
    )
end

"""
    pwrss(m::LinearMixedModel)

The penalized, weighted residual sum-of-squares.
"""
pwrss(m::LinearMixedModel) = abs2(last(last(m.L)))

"""
    ranef!(v::Vector{Matrix{T}}, m::MixedModel{T}, β, uscale::Bool) where {T}

Overwrite `v` with the conditional modes of the random effects for `m`.

If `uscale` is `true` the random effects are on the spherical (i.e. `u`) scale, otherwise
on the original scale
"""
function ranef!(
    v::Vector, m::LinearMixedModel{T}, β::AbstractArray{T}, uscale::Bool
) where {T}
    (k = length(v)) == length(m.reterms) || throw(DimensionMismatch(""))
    L = m.L
    lind = length(L)
    for j in k:-1:1
        lind -= 1
        Ljkp1 = L[lind]
        vj = v[j]
        length(vj) == size(Ljkp1, 2) || throw(DimensionMismatch(""))
        pp1 = size(Ljkp1, 1)
        copyto!(vj, view(Ljkp1, pp1, :))
        mul!(vec(vj), view(Ljkp1, 1:(pp1 - 1), :)', β, -one(T), one(T))
    end
    for i in k:-1:1
        Lii = L[kp1choose2(i)]
        vi = vec(v[i])
        ldiv!(adjoint(isa(Lii, Diagonal) ? Lii : LowerTriangular(Lii)), vi)
        for j in 1:(i - 1)
            mul!(vec(v[j]), L[block(i, j)]', vi, -one(T), one(T))
        end
    end
    if !uscale
        for (t, vv) in zip(m.reterms, v)
            lmul!(t.λ, vv)
        end
    end
    return v
end

ranef!(v::Vector, m::LinearMixedModel, uscale::Bool) = ranef!(v, m, fixef(m), uscale)

"""
    ranef(m::LinearMixedModel; uscale=false)

Return, as a `Vector{Matrix{T}}`, the conditional modes of the random effects in model `m`.

If `uscale` is `true` the random effects are on the spherical (i.e. `u`) scale, otherwise on
the original scale.

For a named variant, see [`raneftables`](@ref).
"""
function ranef(m::LinearMixedModel{T}; uscale=false) where {T}
    reterms = m.reterms
    v = [Matrix{T}(undef, size(t.z, 1), nlevs(t)) for t in reterms]
    return ranef!(v, m, uscale)
end

LinearAlgebra.rank(m::LinearMixedModel) = m.feterm.rank

"""
    rePCA(m::LinearMixedModel; corr::Bool=true)

Return a named tuple of the normalized cumulative variance of a principal components
analysis of the random effects covariance matrices or correlation
matrices when `corr` is `true`.

The normalized cumulative variance is the proportion of the variance for the first
principal component, the first two principal components, etc.  The last element is
always 1.0 representing the complete proportion of the variance.
"""
function rePCA(m::LinearMixedModel; corr::Bool=true)
    pca = PCA.(m.reterms, corr=corr)
    return NamedTuple{fnames(m)}(getproperty.(pca, :cumvar))
end

"""
    PCA(m::LinearMixedModel; corr::Bool=true)

Return a named tuple of the principal components analysis of the random effects
covariance matrices or correlation matrices when `corr` is `true`.
"""

function PCA(m::LinearMixedModel; corr::Bool=true)
    return NamedTuple{fnames(m)}(PCA.(m.reterms, corr=corr))
end

"""
    reevaluateAend!(m::LinearMixedModel)

Reevaluate the last column of `m.A` from `m.Xymat`.  This function should be called
after updating the response.
"""
function reevaluateAend!(m::LinearMixedModel)
    A = m.A
    reterms = m.reterms
    nre = length(reterms)
    trmn = reweight!(m.Xymat, m.sqrtwts)
    ind = kp1choose2(nre)
    for trm in m.reterms
        ind += 1
        mul!(A[ind], trmn', trm)
    end
    mul!(A[end], trmn', trmn)
    return m
end

"""
    refit!(m::LinearMixedModel[, y::Vector]; REML=m.optsum.REML, kwargs...)

Refit the model `m` after installing response `y`.

If `y` is omitted the current response vector is used.
`kwargs` are the same as [`fit!`](@ref).
"""
function refit!(m::LinearMixedModel; REML=m.optsum.REML, kwargs...)
    return fit!(unfit!(m); REML=REML, kwargs...)
end

function refit!(m::LinearMixedModel, y; kwargs...)
    resp = m.y
    length(y) == length(resp) || throw(DimensionMismatch(""))
    copyto!(resp, y)
    return refit!(m; kwargs...)
end

function reweight!(m::LinearMixedModel, weights)
    sqrtwts = map!(sqrt, m.sqrtwts, weights)
    reweight!.(m.reterms, Ref(sqrtwts))
    reweight!(m.Xymat, sqrtwts)
    updateA!(m)
    return updateL!(m)
end

"""
    sdest(m::LinearMixedModel)

Return the estimate of σ, the standard deviation of the per-observation noise.
"""
sdest(m::LinearMixedModel) = something(m.optsum.sigma, √varest(m))

"""
    setθ!(m::LinearMixedModel, v)

Install `v` as the θ parameters in `m`.
"""
function setθ!(m::LinearMixedModel{T}, θ::AbstractVector) where {T}
    parmap, reterms = m.parmap, m.reterms
    length(θ) == length(parmap) || throw(DimensionMismatch())
    reind = 1
    λ = first(reterms).λ
    for (tv, tr) in zip(θ, parmap)
        tr1 = first(tr)
        if reind ≠ tr1
            reind = tr1
            λ = reterms[tr1].λ
        end
        λ[tr[2], tr[3]] = tv
    end
    return m
end

function Base.setproperty!(m::LinearMixedModel, s::Symbol, y)
    return s == :θ ? setθ!(m, y) : setfield!(m, s, y)
end

function Base.show(io::IO, ::MIME"text/plain", m::LinearMixedModel)
    if m.optsum.feval < 0
        @warn("Model has not been fit")
        return nothing
    end
    n, p, q, k = size(m)
    REML = m.optsum.REML
    println(io, "Linear mixed model fit by ", REML ? "REML" : "maximum likelihood")
    println(io, " ", m.formula)
    oo = objective(m)
    if REML
        println(io, " REML criterion at convergence: ", oo)
    else
        nums = Ryu.writefixed.([-oo / 2, oo, aic(m), aicc(m), bic(m)], 4)
        fieldwd = max(maximum(textwidth.(nums)) + 1, 11)
        for label in ["  logLik", "-2 logLik", "AIC", "AICc", "BIC"]
            print(io, rpad(lpad(label, (fieldwd + textwidth(label)) >> 1), fieldwd))
        end
        println(io)
        print.(Ref(io), lpad.(nums, fieldwd))
        println(io)
    end
    println(io)

    show(io, VarCorr(m))

    print(io, " Number of obs: $n; levels of grouping factors: ")
    join(io, nlevs.(m.reterms), ", ")
    println(io)
    println(io, "\n  Fixed-effects parameters:")
    return show(io, coeftable(m))
end

Base.show(io::IO, m::LinearMixedModel) = Base.show(io, MIME"text/plain"(), m)

"""
    _coord(A::AbstractMatrix)

Return the positions and values of the nonzeros in `A` as a
`NamedTuple{(:i, :j, :v), Tuple{Vector{Int32}, Vector{Int32}, Vector{Float64}}}`
"""
function _coord(A::Diagonal)
    return (i=Int32.(axes(A, 1)), j=Int32.(axes(A, 2)), v=A.diag)
end

function _coord(A::UniformBlockDiagonal)
    dat = A.data
    r, c, k = size(dat)
    blk = repeat(r .* (0:(k - 1)); inner=r * c)
    return (
        i=Int32.(repeat(1:r; outer=c * k) .+ blk),
        j=Int32.(repeat(1:c; inner=r, outer=k) .+ blk),
        v=vec(dat),
    )
end

function _coord(A::SparseMatrixCSC{T,Int32}) where {T}
    rv = rowvals(A)
    cv = similar(rv)
    for j in axes(A, 2), k in nzrange(A, j)
        cv[k] = j
    end
    return (i=rv, j=cv, v=nonzeros(A))
end

_coord(A::BlockedSparse) = _coord(A.cscmat)

function _coord(A::Matrix)
    m, n = size(A)
    return (
        i=Int32.(repeat(axes(A, 1); outer=n)),
        j=Int32.(repeat(axes(A, 2); inner=m)),
        v=vec(A),
    )
end

"""
    sparseL(m::LinearMixedModel; fname::Symbol=first(fnames(m)), full::Bool=false)

Return the lower Cholesky factor `L` as a `SparseMatrix{T,Int32}`.

`full` indicates whether the parts of `L` associated with the fixed-effects and response
are to be included.

`fname` specifies the first grouping factor to include. Blocks to the left of the block corresponding
 to `fname` are dropped. The default is the first, i.e., leftmost block and hence all blocks.
"""
function sparseL(
    m::LinearMixedModel{T}; fname::Symbol=first(fnames(m)), full::Bool=false
) where {T}
    L, reterms = m.L, m.reterms
    sblk = findfirst(isequal(fname), fnames(m))
    if isnothing(sblk)
        throw(ArgumentError("fname = $fname is not the name of a grouping factor"))
    end
    blks = sblk:(length(reterms) + full)
    rowoffset, coloffset = 0, 0
    val = (i=Int32[], j=Int32[], v=T[])
    for i in blks, j in first(blks):i
        Lblk = L[block(i, j)]
        cblk = _coord(Lblk)
        append!(val.i, cblk.i .+ Int32(rowoffset))
        append!(val.j, cblk.j .+ Int32(coloffset))
        append!(val.v, cblk.v)
        if i == j
            coloffset = 0
            rowoffset += size(Lblk, 1)
        else
            coloffset += size(Lblk, 2)
        end
    end
    return dropzeros!(tril!(sparse(val...)))
end

"""
    ssqdenom(m::LinearMixedModel)

Return the denominator for penalized sums-of-squares.

For MLE, this value is the number of observations. For REML, this
value is the number of observations minus the rank of the fixed-effects matrix.
The difference is analagous to the use of n or n-1 in the denominator when
calculating the variance.
"""
function ssqdenom(m::LinearMixedModel)::Int
    n = m.dims.n
    return m.optsum.REML ? n - m.dims.p : n
end

"""
    std(m::MixedModel)

Return the estimated standard deviations of the random effects as a `Vector{Vector{T}}`.

FIXME: This uses an old convention of isfinite(sdest(m)).  Probably drop in favor of m.σs
"""
function Statistics.std(m::LinearMixedModel)
    rl = rowlengths.(m.reterms)
    s = sdest(m)
    return isfinite(s) ? rmul!(push!(rl, [1.0]), s) : rl
end

"""
    stderror!(v::AbstractVector, m::LinearMixedModel)

Overwrite `v` with the standard errors of the fixed-effects coefficients in `m`

The length of `v` should be the total number of coefficients (i.e. `length(coef(m))`).
When the model matrix is rank-deficient the coefficients forced to `-0.0` have an
undefined (i.e. `NaN`) standard error.
"""
function stderror!(v::AbstractVector{Tv}, m::LinearMixedModel{T}) where {Tv,T}
    L = feL(m)
    scr = Vector{T}(undef, size(L, 2))
    s = sdest(m)
    fill!(v, zero(Tv) / zero(Tv))  # initialize to appropriate NaN for rank-deficient case
    for i in eachindex(scr)
        fill!(scr, false)
        scr[i] = true
        v[i] = s * norm(ldiv!(L, scr))
    end
    invpermute!(v, m.feterm.piv)
    return v
end

function StatsBase.stderror(m::LinearMixedModel{T}) where {T}
    return stderror!(similar(m.feterm.piv, T), m)
end

"""
    updateA!(m::LinearMixedModel)

Update the cross-product array, `m.A`, from `m.reterms` and `m.Xymat`

This is usually done after a reweight! operation.
"""
function updateA!(m::LinearMixedModel)
    reterms = m.reterms
    k = length(reterms)
    A = m.A
    ind = 1
    for (i, trmi) in enumerate(reterms)
        for j in 1:i
            mul!(A[ind], trmi', reterms[j])
            ind += 1
        end
    end
    Xymattr = adjoint(m.Xymat)
    for trm in reterms
        mul!(A[ind], Xymattr, trm)
        ind += 1
    end
    mul!(A[end], Xymattr, m.Xymat)
    return m
end

"""
    unfit!(model::MixedModel)

Mark a model as unfitted.
"""
function unfit!(model::LinearMixedModel{T}) where {T}
    model.optsum.feval = -1
    model.optsum.initial_step = T[]
    reevaluateAend!(model)

    return model
end

"""
    updateL!(m::LinearMixedModel)

Update the blocked lower Cholesky factor, `m.L`, from `m.A` and `m.reterms` (used for λ only)

This is the crucial step in evaluating the objective, given a new parameter value.
"""
function updateL!(m::LinearMixedModel{T}) where {T}
    A, L, reterms = m.A, m.L, m.reterms
    k = length(reterms)
    for (l, a) in zip(L, A) # copy each element of A to corresponding element of L
        copyto!(l, a)       # For some reason copyto!.(L, A) allocates a lot of memory
    end
    for j in eachindex(reterms) # pre- and post-multiply by Λ, add I to diagonal
        cj = reterms[j]
        scaleinflate!(L[kp1choose2(j)], cj)
        for i in (j + 1):(k + 1)     # postmultiply column by Λ
            rmulΛ!(L[block(i, j)], cj)
        end
        for jj in 1:(j - 1)        # premultiply row by Λ'
            lmulΛ!(cj', L[block(j, jj)])
        end
    end
    for j in 1:(k + 1)             # blocked Cholesky
        Ljj = L[kp1choose2(j)]
        for jj in 1:(j - 1)
            rankUpdate!(Hermitian(Ljj, :L), L[block(j, jj)], -one(T), one(T))
        end
        cholUnblocked!(Ljj, Val{:L})
        LjjT = isa(Ljj, Diagonal) ? Ljj : LowerTriangular(Ljj)
        for i in (j + 1):(k + 1)
            Lij = L[block(i, j)]
            for jj in 1:(j - 1)
                mul!(Lij, L[block(i, jj)], L[block(j, jj)]', -one(T), one(T))
            end
            rdiv!(Lij, LjjT')
        end
    end
    return m
end

"""
    varest(m::LinearMixedModel)

Returns the estimate of σ², the variance of the conditional distribution of Y given B.
"""
function varest(m::LinearMixedModel)
    return isnothing(m.optsum.sigma) ? pwrss(m) / ssqdenom(m) : m.optsum.sigma
end

function StatsBase.weights(m::LinearMixedModel)
    rtwts = m.sqrtwts
    return isempty(rtwts) ? ones(eltype(rtwts), nobs(m)) : abs2.(rtwts)
end


text:
struct HessianCache{T,fdtype,inplace}
    xpp::T
    xpm::T
    xmp::T
    xmm::T
end

function HessianCache(xpp,xpm,xmp,xmm,
                      fdtype::Type{T1}=Val{:hcentral},
                      inplace::Type{Val{T2}} = x isa StaticArray ? Val{false} : Val{true}) where {T1,T2}
    HessianCache{typeof(xpp),fdtype,inplace}(xpp,xpm,xmp,xmm)
end

function HessianCache(x,fdtype::Type{T1}=Val{:hcentral},
                        inplace::Type{Val{T2}} = x isa StaticArray ? Val{false} : Val{true}) where {T1,T2}
    HessianCache{typeof(x),fdtype,inplace}(copy(x),copy(x),copy(x),copy(x))
end

function finite_difference_hessian(f, x::AbstractArray{<:Number},
    fdtype     :: Type{T1}=Val{:hcentral},
    inplace    :: Type{Val{T2}} = x isa StaticArray ? Val{false} : Val{true};
    relstep=default_relstep(fdtype, eltype(x)),
    absstep=relstep) where {T1,T2}

    cache = HessianCache(x, fdtype, inplace)
    finite_difference_hessian(f, x, cache; relstep=relstep, absstep=absstep)
end

function finite_difference_hessian(
    f,x,
    cache::HessianCache{T,fdtype,inplace};
    relstep=default_relstep(fdtype, eltype(x)),
    absstep=relstep) where {T,fdtype,inplace}
    _H = false .* x .* x'
    _H isa SMatrix ? H = MArray(_H) : H = _H
    finite_difference_hessian!(H, f, x, cache; relstep=relstep, absstep=absstep)
    Symmetric(_H isa SMatrix ? SArray(H) : H)
end

function finite_difference_hessian!(H::AbstractMatrix,f,
    x::AbstractArray{<:Number},
    fdtype     :: Type{T1}=Val{:hcentral},
    inplace    :: Type{Val{T2}} = x isa StaticArray ? Val{false} : Val{true};
    relstep=default_relstep(fdtype, eltype(x)),
    absstep=relstep) where {T1,T2}

    cache = HessianCache(x,fdtype,inplace)
    finite_difference_hessian!(H, f, x, cache; relstep=relstep, absstep=absstep)
end

function finite_difference_hessian!(H,f,x,
                                    cache::HessianCache{T,fdtype,inplace};
                                    relstep = default_relstep(fdtype, eltype(x)),
                                    absstep = relstep) where {T,fdtype,inplace}

    @assert fdtype == Val{:hcentral}
    n = length(x)
    xpp, xpm, xmp, xmm = cache.xpp, cache.xpm, cache.xmp, cache.xmm
    fx = f(x)

    if inplace === Val{true}
        _xpp, _xpm, _xmp, _xmm = xpp, xpm, xmp, xmm
    end

    for i = 1:n
        xi = ArrayInterface.allowed_getindex(x,i)
        epsilon = compute_epsilon(Val{:hcentral}, xi, relstep, absstep)

        if inplace === Val{true}
            ArrayInterface.allowed_setindex!(xpp,xi + epsilon,i)
            ArrayInterface.allowed_setindex!(xmm,xi - epsilon,i)
        else
            _xpp = Base.setindex(xpp,xi + epsilon, i)
            _xmm = Base.setindex(xmm,xi - epsilon, i)
        end

        ArrayInterface.allowed_setindex!(H,(f(_xpp) - 2*fx + f(_xmm)) / epsilon^2,i,i)
        epsiloni = compute_epsilon(Val{:central}, xi, relstep, absstep)
        xp = xi + epsiloni
        xm = xi - epsiloni

        if inplace === Val{true}
            ArrayInterface.allowed_setindex!(xpp,xp,i)
            ArrayInterface.allowed_setindex!(xpm,xp,i)
            ArrayInterface.allowed_setindex!(xmp,xm,i)
            ArrayInterface.allowed_setindex!(xmm,xm,i)
        else
            _xpp = Base.setindex(xpp,xp,i)
            _xpm = Base.setindex(xpm,xp,i)
            _xmp = Base.setindex(xmp,xm,i)
            _xmm = Base.setindex(xmm,xm,i)
        end

        for j = i+1:n
            xj = ArrayInterface.allowed_getindex(x,j)
            epsilonj = compute_epsilon(Val{:central}, xj, relstep, absstep)
            xp = xj + epsilonj
            xm = xj - epsilonj

            if inplace === Val{true}
                ArrayInterface.allowed_setindex!(xpp,xp,j)
                ArrayInterface.allowed_setindex!(xpm,xm,j)
                ArrayInterface.allowed_setindex!(xmp,xp,j)
                ArrayInterface.allowed_setindex!(xmm,xm,j)
            else
                _xpp = Base.setindex(_xpp,xp,j)
                _xpm = Base.setindex(_xpm,xm,j)
                _xmp = Base.setindex(_xmp,xp,j)
                _xmm = Base.setindex(_xmm,xm,j)
            end

            ArrayInterface.allowed_setindex!(H,(f(_xpp) - f(_xpm) - f(_xmp) + f(_xmm))/(4*epsiloni*epsilonj),i,j)

            if inplace === Val{true}
                ArrayInterface.allowed_setindex!(xpp,xj,j)
                ArrayInterface.allowed_setindex!(xpm,xj,j)
                ArrayInterface.allowed_setindex!(xmp,xj,j)
                ArrayInterface.allowed_setindex!(xmm,xj,j)
            else
                _xpp = Base.setindex(_xpp,xj,j)
                _xpm = Base.setindex(_xpm,xj,j)
                _xmp = Base.setindex(_xmp,xj,j)
                _xmm = Base.setindex(_xmm,xj,j)
            end
        end

        if inplace === Val{true}
            ArrayInterface.allowed_setindex!(xpp,xi,i)
            ArrayInterface.allowed_setindex!(xpm,xi,i)
            ArrayInterface.allowed_setindex!(xmp,xi,i)
            ArrayInterface.allowed_setindex!(xmm,xi,i)
        end
    end
    LinearAlgebra.copytri!(H,'U')
end

ext: jl


text:
## file formats ##

module DataFmt

importall Base
import Base: _default_delims

export countlines, readdlm, readcsv, writedlm, writecsv

const invalid_dlm = char(0xfffffffe)
const offs_chunk_size = 5000

countlines(nameorfile) = countlines(nameorfile, '\n')
function countlines(filename::String, eol::Char)
    open(filename) do io
        countlines(io, eol)
    end
end
function countlines(io::IO, eol::Char)
    if !isascii(eol)
        error("only ASCII line terminators are supported")
    end
    a = Array(Uint8, 8192)
    nl = 0
    preceded_by_eol = true
    while !eof(io)
        nb = readbytes!(io, a)
        for i=1:nb
            if a[i] == eol
                preceded_by_eol = true
            elseif preceded_by_eol
                preceded_by_eol = false
                nl+=1
            end
        end
    end
    nl
end

readdlm(input, T::Type; opts...) = readdlm(input, invalid_dlm, T, '\n'; opts...)
readdlm(input, dlm::Char, T::Type; opts...) = readdlm(input, dlm, T, '\n'; opts...)

readdlm(input; opts...) = readdlm(input, invalid_dlm, '\n'; opts...)
readdlm(input, dlm::Char; opts...) = readdlm(input, dlm, '\n'; opts...)

readdlm(input, dlm::Char, eol::Char; opts...) = readdlm_auto(input, dlm, Float64, eol, true; opts...)
readdlm(input, dlm::Char, T::Type, eol::Char; opts...) = readdlm_auto(input, dlm, T, eol, false; opts...)

function readdlm_auto(input, dlm::Char, T::Type, eol::Char, auto::Bool; opts...)
    optsd = val_opts(opts)
    use_mmap = get(optsd, :use_mmap, @windows ? false : true)
    isa(input, String) && (fsz = filesize(input); input = use_mmap && (fsz > 0) && fsz < typemax(Int) ? as_mmap(input,fsz) : readall(input))
    sinp = isa(input, Vector{Uint8}) ? ccall(:jl_array_to_string, ByteString, (Array{Uint8,1},), input) :
           isa(input, IO) ? readall(input) :
           input
    readdlm_string(sinp, dlm, T, eol, auto, optsd)
end

function as_mmap(fname::String, fsz::Int64)
    open(fname) do io
        mmap_array(Uint8, (int(fsz),), io)
    end
end

function ascii_if_possible(sbuff::String)
    isa(sbuff, ASCIIString) && return sbuff

    asci = true
    d = sbuff.data
    for idx in 1:length(d)
        (d[idx] < 0x80) ? continue : (asci = false; break)
    end
    asci ? ASCIIString(sbuff.data) : sbuff
end

#
# Handlers act on events generated by the parser.
# Parser calls store_cell on the handler to pass events.
#
# DLMOffsets: Keep offsets (when result dimensions are not known)
# DLMStore: Store values directly into a result store (when result dimensions are known)
abstract DLMHandler

type DLMOffsets <: DLMHandler
    oarr::Vector{Vector{Int}}
    offidx::Int
    thresh::Int
    bufflen::Int

    function DLMOffsets(sbuff::String)
        offsets = Array(Array{Int,1}, 1)
        offsets[1] = Array(Int, offs_chunk_size)
        thresh = int(Base.Sys.total_memory() / sizeof(Int) / 5)
        new(offsets, 1, thresh, length(sbuff.data))
    end
end

function store_cell(dlmoffsets::DLMOffsets, row::Int, col::Int, quoted::Bool, startpos::Int, endpos::Int)
    offidx = dlmoffsets.offidx
    (offidx == 0) && return     # offset collection stopped to avoid choking on memory

    oarr = dlmoffsets.oarr
    offsets = oarr[end]
    if length(offsets) < offidx
        offlen = offs_chunk_size * length(oarr)
        if (offlen + offs_chunk_size) > dlmoffsets.thresh
            est_tot = int(offlen * dlmoffsets.bufflen / endpos)
            if (est_tot - offlen) > offs_chunk_size    # allow another chunk
                # abandon offset collection
                dlmoffsets.oarr = Vector{Int}[]
                dlmoffsets.offidx = 0
                return
            end
        end
        offsets = Array(Int, offs_chunk_size)
        push!(oarr, offsets)
        offidx = 1
    end
    offsets[offidx] = row
    offsets[offidx+1] = col
    offsets[offidx+2] = int(quoted)
    offsets[offidx+3] = startpos
    offsets[offidx+4] = endpos
    dlmoffsets.offidx = offidx + 5
    nothing
end

function result(dlmoffsets::DLMOffsets)
    trimsz = (dlmoffsets.offidx-1)%offs_chunk_size
    ((trimsz > 0) || (dlmoffsets.offidx == 1)) && resize!(dlmoffsets.oarr[end], trimsz)
    dlmoffsets.oarr
end

type DLMStore{T,S<:String} <: DLMHandler
    hdr::Array{String, 2}
    data::Array{T, 2}

    nrows::Int
    ncols::Int
    lastrow::Int
    lastcol::Int
    hdr_offset::Int
    sbuff::S
    auto::Bool
    eol::Char
    tmp64::Array{Float64,1}
end

function DLMStore{T,S<:String}(::Type{T}, dims::NTuple{2,Integer}, has_header::Bool, sbuff::S, auto::Bool, eol::Char)
    (nrows,ncols) = dims
    ((nrows == 0) || (ncols == 0)) && error("Empty input")
    ((nrows < 0) || (ncols < 0)) && error("Invalid dimensions")
    hdr_offset = has_header ? 1 : 0
    DLMStore{T,S}(fill(SubString(sbuff,1,0), 1, ncols), Array(T, nrows-hdr_offset, ncols), nrows, ncols, 1, 0, hdr_offset, sbuff, auto, eol, Array(Float64,1))
end

function store_cell{T,S<:String}(dlmstore::DLMStore{T,S}, row::Int, col::Int, quoted::Bool, startpos::Int, endpos::Int)
    drow = row - dlmstore.hdr_offset

    ncols = dlmstore.ncols
    lastcol = dlmstore.lastcol
    lastrow = dlmstore.lastrow
    cells::Array{T,2} = dlmstore.data
    sbuff::S = dlmstore.sbuff
    tmp64 = dlmstore.tmp64

    endpos = prevind(sbuff, nextind(sbuff,endpos))
    (endpos > 0) && ('\n' == dlmstore.eol) && ('\r' == sbuff[endpos]) && (endpos = prevind(sbuff, endpos))
    sval = quoted ? SubString(sbuff, startpos+1, endpos-1) : SubString(sbuff, startpos, endpos)

    if drow > 0
        # fill missing elements
        while ((drow - lastrow) > 1) || ((drow > lastrow) && (lastcol < ncols))
            (lastcol == ncols) && (lastcol = 0)
            for cidx in (lastcol+1):ncols
                if (T <: String) || (T == Any)
                    cells[lastrow,cidx] = SubString(sbuff, 1, 0)
                elseif ((T <: Number) || (T <: Char)) && dlmstore.auto
                    throw(TypeError(:store_cell, "", Any, T))
                else
                    error("missing value at row $lastrow column $cidx")
                end
            end
            lastcol = ncols
            lastrow += 1
        end

        # fill data
        if quoted && ('"' in sval)
            fail = colval(replace(sval, r"\"\"", "\""), cells, drow, col, tmp64)
        else
            fail = colval(sval, cells, drow, col, tmp64)
        end
        if fail
            ((T <: Number) && dlmstore.auto) ? throw(TypeError(:store_cell, "", Any, T)) : error("file entry \"$(sval)\" cannot be converted to $T")
        end

        dlmstore.lastrow = drow
        dlmstore.lastcol = col
    else
        # fill header
        colval((quoted && ('"' in sval)) ? replace(sval, r"\"\"", "\"") : sval, dlmstore.hdr, 1, col, tmp64)
    end

    nothing
end

function result{T}(dlmstore::DLMStore{T})
    nrows = dlmstore.nrows - dlmstore.hdr_offset
    ncols = dlmstore.ncols
    lastcol = dlmstore.lastcol
    lastrow = dlmstore.lastrow
    cells = dlmstore.data
    sbuff = dlmstore.sbuff

    if (lastcol < ncols) || (lastrow < nrows)
        while lastrow <= nrows
            (lastcol == ncols) && (lastcol = 0; lastrow += 1)
            for cidx in (lastcol+1):ncols
                if (T <: String) || (T == Any)
                    cells[lastrow,cidx] = SubString(sbuff, 1, 0)
                elseif ((T <: Number) || (T <: Char)) && dlmstore.auto
                    throw(TypeError(:store_cell, "", Any, T))
                else
                    error("missing value at row $lastrow column $cidx")
                end
            end
            lastcol = ncols
            (lastrow == nrows) && break
        end
        dlmstore.lastrow = lastrow
        dlmstore.lastcol = ncols
    end
    (dlmstore.hdr_offset > 0) ? (dlmstore.data, dlmstore.hdr) : dlmstore.data
end

function readdlm_string(sbuff::String, dlm::Char, T::Type, eol::Char, auto::Bool, optsd::Dict)
    ign_empty = (dlm == invalid_dlm)
    quotes = get(optsd, :quotes, true)
    comments = get(optsd, :comments, true)
    comment_char = get(optsd, :comment_char, '#')
    dims = get(optsd, :dims, nothing)
    has_header = get(optsd, :has_header, false)

    offset_handler = (dims == nothing) ? DLMOffsets(sbuff) : DLMStore(T, dims, has_header, sbuff, auto, eol)

    for retry in 1:2
        try
            dims = dlm_parse(sbuff, eol, dlm, '"', comment_char, ign_empty, quotes, comments, offset_handler)
            break
        catch ex
            if isa(ex, TypeError) && (ex.func == :store_cell)
                T = ex.expected
            elseif get(optsd, :ignore_invalid_chars, false)
                sbuff = ascii_if_possible(convert(typeof(sbuff), sbuff.data, ""))
            else
                rethrow(ex)
            end
            offset_handler = (dims == nothing) ? DLMOffsets(sbuff) : DLMStore(T, dims, has_header, sbuff, auto, eol)
        end
    end

    isa(offset_handler, DLMStore) && (return result(offset_handler))

    offsets = result(offset_handler)
    !isempty(offsets) && (return dlm_fill(T, offsets, dims, has_header, sbuff, auto, eol))

    optsd[:dims] = dims
    return readdlm_string(sbuff, dlm, T, eol, auto, optsd)
end

const valid_opts = [:has_header, :ignore_invalid_chars, :use_mmap, :quotes, :comments, :dims, :comment_char]
const valid_opt_types = [Bool, Bool, Bool, Bool, Bool, NTuple{2,Integer}, Char]
function val_opts(opts)
    d = Dict{Symbol,Union(Bool,NTuple{2,Integer},Char)}()
    for opt in opts
        !in(opt[1], valid_opts) && error("unknown option $(opt[1])")
        opt_typ = valid_opt_types[findfirst(valid_opts, opt[1])]
        !isa(opt[2], opt_typ) && error("$(opt[1]) should be of type $opt_typ")
        d[opt[1]] = opt[2]
    end
    d
end

function dlm_fill(T::DataType, offarr::Vector{Vector{Int}}, dims::NTuple{2,Integer}, has_header::Bool, sbuff::String, auto::Bool, eol::Char)
    idx = 1
    offidx = 1
    offsets = offarr[1]
    try
        dh = DLMStore(T, dims, has_header, sbuff, auto, eol)
        while idx <= length(offsets)
            row = offsets[idx]
            col = offsets[idx+1]
            quoted = bool(offsets[idx+2])
            startpos = offsets[idx+3]
            endpos = offsets[idx+4]

            ((idx += 5) > offs_chunk_size) && (offidx < length(offarr)) && (idx = 1; offsets = offarr[offidx += 1])

            store_cell(dh, row, col, quoted, startpos, endpos)
        end
        return result(dh)
    catch ex
        isa(ex, TypeError) && (ex.func == :store_cell) && (return dlm_fill(ex.expected, offarr, dims, has_header, sbuff, auto, eol))
        rethrow(ex)
    end
end

colval{T<:Number, S<:String}(sval::S, cells::Array{T,2}, row::Int, col::Int, tmp64::Array{Float64,1}) = (float64_isvalid(sval, tmp64) ? ((cells[row,col] = tmp64[1]); false) : true)
colval{T<:String, S<:String}(sval::S, cells::Array{T,2}, row::Int, col::Int, tmp64::Array{Float64,1}) = ((cells[row,col] = sval); false)
colval{S<:String}(sval::S, cells::Array{Any,2}, row::Int, col::Int, tmp64::Array{Float64,1}) = ((cells[row,col] = float64_isvalid(sval, tmp64) ? tmp64[1] : sval); false)
colval{T<:Char, S<:String}(sval::S, cells::Array{T,2}, row::Int, col::Int, tmp64::Array{Float64,1}) = ((length(sval) == 1) ? ((cells[row,col] = next(sval,1)[1]); false) : true)
colval{S<:String}(sval::S, cells::Array, row::Int, col::Int, tmp64::Array{Float64,1}) = true

dlm_parse(s::ASCIIString, eol::Char, dlm::Char, qchar::Char, cchar::Char, ign_adj_dlm::Bool, allow_quote::Bool, allow_comments::Bool, dh::DLMHandler) = dlm_parse(s.data, uint8(eol), uint8(dlm), uint8(qchar), uint8(cchar), ign_adj_dlm, allow_quote, allow_comments, dh)
function dlm_parse{T,D}(dbuff::T, eol::D, dlm::D, qchar::D, cchar::D, ign_adj_dlm::Bool, allow_quote::Bool, allow_comments::Bool, dh::DLMHandler)
    all_ascii = (D <: Uint8) || (isascii(eol) && isascii(dlm) && (!allow_quote || isascii(qchar)) && (!allow_comments || isascii(cchar)))
    (T <: UTF8String) && all_ascii && (return dlm_parse(dbuff.data, uint8(eol), uint8(dlm), uint8(qchar), uint8(cchar), ign_adj_dlm, allow_quote, allow_comments, dh))
    ncols = nrows = col = 0
    is_default_dlm = (dlm == convert(D, invalid_dlm))
    error_str = ""
    state = 0   # 0: begin field, 1: quoted field, 2: unquoted field, 3: second quote (could either be end of field or escape character), 4: comment
    is_eol = is_dlm = is_cr = is_quote = is_comment = expct_col = false
    idx = 1
    try
        slen = sizeof(dbuff)
        col_start_idx = 1
        was_cr = false
        while idx <= slen
            val,idx = next(dbuff, idx)
            if (is_eol = (val == eol))
                is_dlm = is_comment = is_cr = is_quote = false
            elseif (is_dlm = (is_default_dlm ? in(val, _default_delims) : (val == dlm)))
                is_comment = is_cr = is_quote = false
            elseif (is_quote = (val == qchar))
                is_comment = is_cr = false
            elseif (is_comment = (val == cchar))
                is_cr = false
            else
                is_cr = (eol == '\n') && (val == '\r')
            end

            if 2 == state   # unquoted field
                if is_dlm
                    state = 0
                    col += 1
                    store_cell(dh, nrows+1, col, false, col_start_idx, idx-2)
                    col_start_idx = idx
                    !ign_adj_dlm && (expct_col = true)
                elseif is_eol
                    nrows += 1
                    col += 1
                    store_cell(dh, nrows, col, false, col_start_idx, idx - (was_cr ? 3 : 2))
                    col_start_idx = idx
                    ncols = max(ncols, col)
                    col = 0
                    state = 0
                elseif (is_comment && allow_comments)
                    nrows += 1
                    col += 1
                    store_cell(dh, nrows, col, false, col_start_idx, idx - 2)
                    ncols = max(ncols, col)
                    col = 0
                    state = 4
                end
            elseif 1 == state   # quoted field
                is_quote && (state = 3)
            elseif 4 == state   # comment line
                if is_eol
                    col_start_idx = idx
                    state = 0
                end
            elseif 0 == state   # begin field
                if is_quote
                    state = (allow_quote && !was_cr) ? 1 : 2
                    expct_col = false
                elseif is_dlm
                    if !ign_adj_dlm
                        expct_col = true
                        col += 1
                        store_cell(dh, nrows+1, col, false, col_start_idx, idx-2)
                    end
                    col_start_idx = idx
                elseif is_eol
                    nrows += 1
                    if expct_col
                        col += 1
                        store_cell(dh, nrows, col, false, col_start_idx, idx - (was_cr ? 3 : 2))
                    end
                    col_start_idx = idx
                    ncols = max(ncols, col)
                    col = 0
                    expct_col = false
                elseif is_comment && allow_comments
                    if col > 0
                        nrows += 1
                        if expct_col
                            col += 1
                            store_cell(dh, nrows, col, false, col_start_idx, idx - 2)
                        end
                        ncols = max(ncols, col)
                        col = 0
                    end
                    expct_col = false
                    state = 4
                elseif !is_cr
                    state = 2
                    expct_col = false
                end
            elseif 3 == state   # second quote
                if is_quote && !was_cr
                    state = 1
                elseif is_dlm && !was_cr
                    state = 0
                    col += 1
                    store_cell(dh, nrows+1, col, true, col_start_idx, idx-2)
                    col_start_idx = idx
                    !ign_adj_dlm && (expct_col = true)
                elseif is_eol
                    nrows += 1
                    col += 1
                    store_cell(dh, nrows, col, true, col_start_idx, idx - (was_cr ? 3 : 2))
                    col_start_idx = idx
                    ncols = max(ncols, col)
                    col = 0
                    state = 0
                elseif is_comment && allow_comments && !was_cr
                    nrows += 1
                    col += 1
                    store_cell(dh, nrows, col, true, col_start_idx, idx - 2)
                    ncols = max(ncols, col)
                    col = 0
                    state = 4
                elseif (is_cr && was_cr) || !is_cr
                    error_str = escape_string("unexpected character '$(char(val))' after quoted field at row $(nrows+1) column $(col+1)")
                    break
                end
            end
            was_cr = is_cr
        end

        if isempty(error_str)
            if 1 == state       # quoted field
                error_str = "truncated column at row $(nrows+1) column $(col+1)"
            elseif (2 == state) || (3 == state) || ((0 == state) && is_dlm)   # unquoted field, second quote, or begin field with last character as delimiter
                col += 1
                nrows += 1
                store_cell(dh, nrows, col, (3 == state), col_start_idx, idx-1)
                ncols = max(ncols, col)
            end
        end
    catch ex
        if isa(ex, TypeError) && (ex.func == :store_cell)
            rethrow(ex)
        else
            error("at row $(nrows+1), column $col : $ex)")
        end
    end
    !isempty(error_str) && error(error_str)

    return (nrows, ncols)
end

readcsv(io; opts...)          = readdlm(io, ','; opts...)
readcsv(io, T::Type; opts...) = readdlm(io, ',', T; opts...)

# todo: keyword argument for # of digits to print
writedlm_cell(io::IO, elt::FloatingPoint, dlm, quotes) = print_shortest(io, elt)
function writedlm_cell{T}(io::IO, elt::String, dlm::T, quotes::Bool)
    if quotes && !isempty(elt) && (('"' in elt) || ('\n' in elt) || ((T <: Char) ? (dlm in elt) : contains(elt, dlm)))
        print(io, '"', replace(elt, r"\"", "\"\""), '"')
    else
        print(io, elt)
    end
end
writedlm_cell(io::IO, elt, dlm, quotes) = print(io, elt)
function writedlm(io::IO, a::AbstractVecOrMat, dlm; opts...)
    optsd = val_opts(opts)
    quotes = get(optsd, :quotes, true)
    pb = PipeBuffer()
    nr = size(a,1)
    nc = size(a,2)
    for i = 1:nr
        for j = 1:nc
            writedlm_cell(pb, a[i,j], dlm, quotes)
            j == nc ? write(pb,'\n') : print(pb,dlm)
        end
        (nb_available(pb) > (16*1024)) && write(io, takebuf_array(pb))
    end
    write(io, takebuf_array(pb))
    nothing
end

writedlm{T}(io::IO, a::AbstractArray{T,0}, dlm; opts...) = writedlm(io, reshape(a,1), dlm; opts...)

function writedlm(io::IO, a::AbstractArray, dlm; opts...)
    tail = size(a)[3:end]
    function print_slice(idxs...)
        writedlm(io, sub(a, 1:size(a,1), 1:size(a,2), idxs...), dlm; opts...)
        if idxs != tail
            print(io, "\n")
        end
    end
    cartesianmap(print_slice, tail)
end

function writedlm(io::IO, itr, dlm; opts...)
    optsd = val_opts(opts)
    quotes = get(optsd, :quotes, true)
    pb = PipeBuffer()
    for row in itr
        state = start(row)
        while !done(row, state)
            (x, state) = next(row, state)
            writedlm_cell(pb, x, dlm, quotes)
            done(row, state) ? write(pb,'\n') : print(pb,dlm)
        end
        (nb_available(pb) > (16*1024)) && write(io, takebuf_array(pb))
    end
    write(io, takebuf_array(pb))
    nothing
end

function writedlm(fname::String, a, dlm; opts...)
    open(fname, "w") do io
        writedlm(io, a, dlm; opts...)
    end
end

writedlm(io, a; opts...) = writedlm(io, a, '\t'; opts...)
writecsv(io, a; opts...) = writedlm(io, a, ','; opts...)

writemime(io::IO, ::MIME"text/csv", a::AbstractVecOrMat) = writedlm(io, a, ',')
writemime(io::IO, ::MIME"text/tab-separated-values", a::AbstractVecOrMat) = writedlm(io, a, '\t')

end # module DataFmt

max_stars_repo_name: porterjamesj/julia

max_issues_repo_name: porterjamesj/julia

max_forks_repo_name: porterjamesj/julia


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

include("replace.jl")


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Functions related to the conversion from DCM to quaternion.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export angle_to_quat

"""
    angle_to_quat(θ₁::T1[, θ₂::T2[, θ₃::T3]], rot_seq::Symbol = :ZYX) where {T1<:Number, T2<:Number, T3<:Number}
    angle_to_quat(eulerang::EulerAngles)

Create a quaternion that perform a set of rotations (`θ₁`, `θ₂`, `θ₃`) about the
coordinate axes specified in `rot_seq`.

The input values of the origin Euler angles can also be passed inside the
structure `Θ` (see [`EulerAngles`](@ref)).

The rotation sequence is defined by a `Symbol` specifing the rotation axes. The
possible values depends on the number of rotations as follows:

- **1 rotation** (`θ₁`): `:X`, `:Y`, or `:Z`.
- **2 rotations** (`θ₁`, `θ₂`): `:XY`, `:XZ`, `:YX`, `:YZ`, `:ZX`, or `:ZY`.
- **3 rotations** (`θ₁`, `θ₂`, `θ₃`): `:XYX`, `XYZ`, `:XZX`, `:XZY`, `:YXY`,
    `:YXZ`, `:YZX`, `:YZY`, `:ZXY`, `:ZXZ`, `:ZYX`, or `:ZYZ`

!!! note
    The type of the new quaternion will be obtained by promiting `T1`, `T2`, and
    `T3`.

# Remarks

This function assigns `q = q1 * q2 * q3` in which `qi` is the quaternion related
with the *i*-th rotation, `i Є [1,2,3]`. If the *i*-th rotation is not specified,
then `qi = Quaternion(I)`.

# Example

```jldoctest
julia> angle_to_quat(pi / 2, :X)
Quaternion{Float64}:
  + 0.707107 + 0.707107⋅i + 0.0⋅j + 0.0⋅k

julia> angle_to_quat(pi / 5, pi / 7, :YZ)
Quaternion{Float64}:
  + 0.927212 + 0.0687628⋅i + 0.301269⋅j + 0.21163⋅k

julia> angle_to_quat(pi / 5, pi / 7, 0, :YZX)
Quaternion{Float64}:
  + 0.927212 + 0.0687628⋅i + 0.301269⋅j + 0.21163⋅k

julia> angle_to_quat(pi / 2, pi / 3, pi / 4, :ZYX)
Quaternion{Float64}:
  + 0.701057 - 0.092296⋅i + 0.560986⋅j + 0.430459⋅k
```
"""
function angle_to_quat(θ::Number, rot_seq::Symbol)
    # Compute the sines and cosines of half angle.
    s, c = sincos(θ / 2)

    # Make sure that the real part is always positive.
    if c < 0
        c = -c
        s = -s
    end

    if rot_seq == :X
        return Quaternion(c, s, 0, 0)
    elseif rot_seq == :Y
        return Quaternion(c, 0, s, 0)
    elseif rot_seq == :Z
        return Quaternion(c, 0, 0, s)
    else
        throw(ArgumentError("rot_seq must be :X, :Y, or :Z"))
    end
end

function angle_to_quat(
    θ₁::T1,
    θ₂::T2,
    rot_seq::Symbol
) where {T1<:Number, T2<:Number}
    T = promote_type(T1, T2)

    # Compute the sines and cosines of half angle.
    s₁, c₁ = sincos(T(θ₁) / 2)
    s₂, c₂ = sincos(T(θ₂) / 2)

    # When we have two rotations, the `q0` component is always the same.
    q0 = c₁ * c₂

    s = (q0 < 0) ? -1 : +1

    if rot_seq == :ZY
        return Quaternion(
            s * q0,
            s * (-s₁ * s₂),
            s * ( c₁ * s₂),
            s * ( s₁ * c₂)
        )
    elseif rot_seq == :XY
        return Quaternion(
            s * q0,
            s * (s₁ * c₂),
            s * (c₁ * s₂),
            s * (s₁ * s₂)
        )
    elseif rot_seq == :XZ
        return Quaternion(
            s * q0,
            s * ( s₁ * c₂),
            s * (-s₁ * s₂),
            s * ( c₁ * s₂)
        )
    elseif rot_seq == :YX
        return Quaternion(
            s * q0,
            s * ( c₁ * s₂),
            s * ( s₁ * c₂),
            s * (-s₁ * s₂)
        )
    elseif rot_seq == :YZ
        return Quaternion(
            s * q0,
            s * (s₁ * s₂),
            s * (s₁ * c₂),
            s * (c₁ * s₂)
        )
    elseif rot_seq == :ZX
        return Quaternion(
            s * q0,
            s * (c₁ * s₂),
            s * (s₁ * s₂),
            s * (s₁ * c₂)
        )
    elseif rot_seq == :ZY
        return Quaternion(
            s * q0,
            s * (-s₁ * s₂),
            s * ( c₁ * s₂),
            s * ( s₁ * c₂)
        )
    else
        throw(ArgumentError("The rotation sequence :$rot_seq is not valid."))
    end
end

function angle_to_quat(
    θ₁::T1,
    θ₂::T2,
    θ₃::T3,
    rot_seq::Symbol = :ZYX
) where {T1<:Number, T2<:Number, T3<:Number}
    T = promote_type(T1, T2, T3)

    # Compute the sines and cosines of half angle.
    s₁, c₁ = sincos(T(θ₁) / 2)
    s₂, c₂ = sincos(T(θ₂) / 2)
    s₃, c₃ = sincos(T(θ₃) / 2)

    if rot_seq == :ZYX
        q0 = c₁ * c₂ * c₃ + s₁ * s₂ * s₃

        s = (q0 < 0) ? -1 : +1

        return Quaternion(
            s * q0,
            s * (c₁ * c₂ * s₃ - s₁ * s₂ * c₃),
            s * (c₁ * s₂ * c₃ + s₁ * c₂ * s₃),
            s * (s₁ * c₂ * c₃ - c₁ * s₂ * s₃)
        )
    elseif rot_seq == :XYX
        q0 = c₁ * c₂ * c₃ - s₁ * c₂ * s₃

        s = (q0 < 0) ? -1 : +1

        return Quaternion(
            s * q0,
            s * (c₁ * c₂ * s₃ + s₁ * c₂ * c₃),
            s * (c₁ * s₂ * c₃ + s₁ * s₂ * s₃),
            s * (s₁ * s₂ * c₃ - c₁ * s₂ * s₃)
        )
    elseif rot_seq == :XYZ
        q0 = c₁ * c₂ * c₃ - s₁ * s₂ * s₃

        s = (q0 < 0) ? -1 : +1

        return Quaternion(
            s * q0,
            s * (s₁ * c₂ * c₃ + c₁ * s₂ * s₃),
            s * (c₁ * s₂ * c₃ - s₁ * c₂ * s₃),
            s * (c₁ * c₂ * s₃ + s₁ * s₂ * c₃)
        )
    elseif rot_seq == :XZX
        q0 = c₁ * c₂ * c₃ - s₁ * c₂ * s₃

        s = (q0 < 0) ? -1 : +1

        return Quaternion(
            s * q0,
            s * (c₁ * c₂ * s₃ + s₁ * c₂ * c₃),
            s * (c₁ * s₂ * s₃ - s₁ * s₂ * c₃),
            s * (c₁ * s₂ * c₃ + s₁ * s₂ * s₃)
        )
    elseif rot_seq == :XZY
        q0 = c₁ * c₂ * c₃ + s₁ * s₂ * s₃

        s = (q0 < 0) ? -1 : +1

        return Quaternion(
            s * q0,
            s * (s₁ * c₂ * c₃ - c₁ * s₂ * s₃),
            s * (c₁ * c₂ * s₃ - s₁ * s₂ * c₃),
            s * (c₁ * s₂ * c₃ + s₁ * c₂ * s₃)
        )
    elseif rot_seq == :YXY
        q0 = c₁ * c₂ * c₃ - s₁ * c₂ * s₃

        s = (q0 < 0) ? -1 : +1

        return Quaternion(
            s * q0,
            s * (c₁ * s₂ * c₃ + s₁ * s₂ * s₃),
            s * (c₁ * c₂ * s₃ + s₁ * c₂ * c₃),
            s * (c₁ * s₂ * s₃ - s₁ * s₂ * c₃)
        )
    elseif rot_seq == :YXZ
        q0 = c₁ * c₂ * c₃ + s₁ * s₂ * s₃

        s = (q0 < 0) ? -1 : +1

        return Quaternion(
            s * q0,
            s * (c₁ * s₂ * c₃ + s₁ * c₂ * s₃),
            s * (s₁ * c₂ * c₃ - c₁ * s₂ * s₃),
            s * (c₁ * c₂ * s₃ - s₁ * s₂ * c₃)
        )
    elseif rot_seq == :YZX
        q0 = c₁ * c₂ * c₃ - s₁ * s₂ * s₃

        s = (q0 < 0) ? -1 : +1

        return Quaternion(
            s * q0,
            s * (c₁ * c₂ * s₃ + s₁ * s₂ * c₃),
            s * (s₁ * c₂ * c₃ + c₁ * s₂ * s₃),
            s * (c₁ * s₂ * c₃ - s₁ * c₂ * s₃)
        )
    elseif rot_seq == :YZY
        q0 = c₁ * c₂ * c₃ - s₁ * c₂ * s₃

        s = (q0 < 0) ? -1 : +1

        return Quaternion(
            s * q0,
            s * (s₁ * s₂ * c₃ - c₁ * s₂ * s₃),
            s * (c₁ * c₂ * s₃ + s₁ * c₂ * c₃),
            s * (c₁ * s₂ * c₃ + s₁ * s₂ * s₃)
        )
    elseif rot_seq == :ZXY
        q0 = c₁ * c₂ * c₃ - s₁ * s₂ * s₃

        s = (q0 < 0) ? -1 : +1

        return Quaternion(
            s * q0,
            s * (c₁ * s₂ * c₃ - s₁ * c₂ * s₃),
            s * (c₁ * c₂ * s₃ + s₁ * s₂ * c₃),
            s * (s₁ * c₂ * c₃ + c₁ * s₂ * s₃)
        )
    elseif rot_seq == :ZXZ
        q0 = c₁ * c₂ * c₃ - s₁ * c₂ * s₃

        s = (q0 < 0) ? -1 : +1

        return Quaternion(
            s * q0,
            s * (c₁ * s₂ * c₃ + s₁ * s₂ * s₃),
            s * (s₁ * s₂ * c₃ - c₁ * s₂ * s₃),
            s * (c₁ * c₂ * s₃ + s₁ * c₂ * c₃)
        )
    elseif rot_seq == :ZYZ
        q0 = c₁ * c₂ * c₃ - s₁ * c₂ * s₃

        s = (q0 < 0) ? -1 : +1

        return Quaternion(
            s * q0,
            s * (c₁ * s₂ * s₃ - s₁ * s₂ * c₃),
            s * (c₁ * s₂ * c₃ + s₁ * s₂ * s₃),
            s * (c₁ * c₂ * s₃ + s₁ * c₂ * c₃)
        )
    else
        throw(ArgumentError("The rotation sequence :$rot_seq is not valid."))
    end
end

angle_to_quat(Θ::EulerAngles) = angle_to_quat(Θ.a1, Θ.a2, Θ.a3, Θ.rot_seq)


text:
using Test
using DiffEqBase
using Sundials

@testset "Species Selection" begin
    outcore = readinput("../src/testing/minimal.rms")
    outedge = readinput("../src/testing/minimal_edge.rms")
    corespcs = outcore["phase"]["Species"]
    edgeishspcs = outedge["phase"]["Species"]
    corerxns = outcore["phase"]["Reactions"]
    edgeishrxns = outedge["phase"]["Reactions"]

    coreedgespcs = copy(corespcs)
    index = length(coreedgespcs)
    for spc in edgeishspcs
        if nothing === findfirst(x->x.name==spc.name,coreedgespcs)
            newspc = Species(;name=spc.name,index=index+1,inchi=spc.inchi,smiles=spc.smiles,
                adjlist=spc.adjlist,thermo=spc.thermo,atomnums=spc.atomnums,diffusion=spc.diffusion,
                radius=spc.radius,radicalelectrons=spc.radicalelectrons,molecularweight=spc.molecularweight)
            index += 1
            push!(coreedgespcs,newspc)
        end
    end
    coreedgespcsnames = getfield.(coreedgespcs,:name)

    coreedgerxns = copy(corerxns)
    index = length(coreedgerxns)
    for rxn in edgeishrxns
        out = findfirst(x->getrxnstr(x)==getrxnstr(rxn),coreedgerxns)
        if nothing === out
            inds = findall(x->getrxnstr(x)==getrxnstr(rxn),edgeishrxns)
            for ind in inds
                rxnout = edgeishrxns[ind]
                reactants = [coreedgespcs[findfirst(x->spc.name==x,coreedgespcsnames)] for spc in rxnout.reactants]
                products = [coreedgespcs[findfirst(x->spc.name==x,coreedgespcsnames)] for spc in rxnout.products]
                reactantinds = [findfirst(x->spc.name==x,coreedgespcsnames) for spc in rxnout.reactants]
                productinds = [findfirst(x->spc.name==x,coreedgespcsnames) for spc in rxnout.products]
                newrxn = ElementaryReaction(;index=index+1,reactants=reactants,reactantinds=reactantinds,products=products,
                    productinds=productinds,kinetics=rxnout.kinetics,electronchange=rxnout.electronchange,
                    radicalchange=rxnout.radicalchange,reversible=rxnout.reversible,pairs=rxnout.pairs)
                push!(coreedgerxns,newrxn)
            end
        end
    end

    coregas = IdealGas(corespcs,corerxns);
    coreedgegas = IdealGas(coreedgespcs,coreedgerxns);

    initialconds = Dict(["T"=>1350.0,"P"=>1.0e5,"ethane"=>1.0]);
    spc = coregas.species[5] #ethane
    termination = [TerminationConversion(spc,0.9),TerminationTime(1e6)];
    coredomain,y0,corep = ConstantTPDomain(phase=coregas,initialconds=initialconds);
    react = Reactor(coredomain,y0,(0.0,1e6);p=corep);
    coreedgedomain,coreedgey0,coreedgep = ConstantTPDomain(phase=coreedgegas,initialconds=initialconds);

    (terminated,resurrected,invalidobjects,unimolecularthreshold,bimolecularthreshold,
    trimolecularthreshold,maxedgespeciesrateratios) = selectobjects(react,coreedgedomain,[],coredomain,
        [],corep,coreedgep,0.03,0.03,false,true,5,0.005,1.0,1.0,true,termination,1.0e8)

    @test terminated == false
    @test invalidobjects[1].name == "[CH2]CCC"
    @test unimolecularthreshold[5] == true
    @test maxedgespeciesrateratios[5] ≈ 0.008333897672002308 rtol=1e-3

end;


text:
using Gadfly
using MLLibrary
using Colors
path =Pkg.dir("MLLibrary","resources","datasets")
week3=joinpath("Week3")
# resources = #[["res"=>"regression_test.csv","path"=>joinpath(path,week3),"proxy"=>[]]]
resources = [["res"=>"regression-0.05.csv","path"=>joinpath(path,week3),"proxy"=>[]],["res"=>"regression-A.csv","path"=>joinpath(path,week3),"proxy"=>[]],["res"=>"regression-B.csv","path"=>joinpath(path,week3),"proxy"=>[]],["res"=>"regression-C.csv","path"=>joinpath(path,week3),"proxy"=>[]],["res"=>"abalone.data","path"=>joinpath(path), "proxy"=>[Dict("column"=>1,"proxy"=>["M" "F" "I"])]]]
# [["res"=>"abalone.data","path"=>joinpath(path), "proxy"=>[Dict("column"=>1,"proxy"=>["M" "F" "I"])]]]

errors = Dict()
for res in resources
  seed=123
  training, test = MLLibrary.DataSampling.ThresholdSampling.getSets(joinpath(res["path"],res["res"]),",",.1,seed)

  proxy = res["proxy"]

  training=[MLLibrary.DataManipulation.DataTransform.addProxyColumns(training[:,1:end-1],proxy) training[:,end]]

  test=[MLLibrary.DataManipulation.DataTransform.addProxyColumns(test[:,1:end-1],proxy) test[:,end]]

# println(training)
# println(test)

  training = MLLibrary.DataManipulation.DataTransform.convertToFloat(training)
  test= MLLibrary.DataManipulation.DataTransform.convertToFloat(test)
  println("Data is ready")
  println(typeof(training))

  betacol = MLLibrary.LearningAlgorithms.LinearRegression.linearRegression(training[:,1:end-1],training[:,end],MLLibrary.LinearAlgebra.HouseHolders.householdDecompose)

  println("finished dcomposing")
  test_y = test[:,end]
  test = test[:,1:end-1]

  observed = test*betacol

  #println((observed),(test_y))

  rms = round(MLLibrary.Utils.ErrorAnalysis.calculateRootMeanStandardError(observed,test_y),2)
  errors[res["res"]]=rms
end

println(errors)

p = plot(

  Guide.ylabel("Root Mean Square Error"),
  Guide.xlabel("Data Set"),
  Guide.title("Root Mean Square Errors"),
  x=collect(keys(errors)),
  y=collect(values(errors)),
  Geom.bar,
  Geom.label(position=:above),
  label=map(string,(values(errors)))
  #color="DataSet"
)
# index = 1
# color_keys = []
# colors = distinguishable_colors(length(errors))
#
# for (k,v) in errors
#   a = layer(
#     x=k,
#     y=v,
#     Geom.histogram,
#     Theme(default_color=(colors[index]))
#     #Guide.colorkey(k)
#   )
#   push!(color_keys,k)
#
#   index+=1
#   push!(p,a)
# end
#
# push!(p,  Guide.manual_color_key("Legend", color_keys, colors))

draw(PNG(joinpath(Pkg.dir("MLLibrary"),"resources","paResources","Week3","RMSEErrors.png"), 800px, 600px), p)

max_stars_repo_name: ajermaky/CustomMLLibrary

max_issues_repo_name: ajermaky/CustomMLLibrary

max_forks_repo_name: ajermaky/CustomMLLibrary


text:
using Documenter, TermoQuimica

DocMeta.setdocmeta!(TermoQuimica, :DocTestSetup, :(using TermoQuimica); recursive=true)

makedocs(;
    modules = [TermoQuimica],
    authors  = "Emilio Alvizo Velázquez",
    repo="https://github.com/EmilioAlvizo/TermoQuimica.jl/blob/{commit}{path}#{line}",
    sitename="TermoQuimica.jl",
    format = Documenter.HTML(;
        prettyurls = get(ENV, "CI", "false") == "true",
        canonical="https://EmilioAlvizo.github.io/TermoQuimica.jl",
        assets = String[]),
    pages = [
        "Inicio" => "index.md",
        "Background" => "background.md",
        "Guia de usuario" => "userguide.md",
        "API" => "api.md",
        "Servidor" => "servidor.md",
    ])

deploydocs(
    repo="github.com/EmilioAlvizo/TermoQuimica.jl",
)

max_stars_repo_name: EmilioAlvizo/TermoQuimica

max_issues_repo_name: EmilioAlvizo/TermoQuimica

max_forks_repo_name: EmilioAlvizo/TermoQuimica


text:
"""

# struct Observations

## Fields

- datas : list of data vector one for each observations
- value\\_at\\_datas
For regressions problems for example length(datas) is number of observations.
    and length(datas[1]) is 1+dimension of observations data beccause of interception terms.
"""
struct Observations
    datas :: Vector{Vector{Float64}}
    value_at_data :: Vector{Float64}
    #
    function Observations(datas::Vector{Vector{Float64}}, values::Vector{Float64})
        new(datas, values)
    end
end


text:
const snapsvgjs = joinpath(@__DIR__, "..", "templates", "snap.svg-min.js")
const viewerjs = joinpath(@__DIR__, "viewer.js")

function escape_script(js::AbstractString)
    return replace(js, "]]" => "] ]")
end

function svgheader(io::IO, fig_id::AbstractString; width=1200, height=706, font="Verdana")

    y_msg = height - 17
    print(io, """<?xml version="1.0" standalone="no"?>
<!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN" "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd">
<svg version="1.1" width="$(width)" height="$(height)" viewBox="0 0 $(width) $(height)" xmlns="http://www.w3.org/2000/svg" >
<defs >
    <linearGradient id="background" y1="0" y2="1" x1="0" x2="0" >
        <stop stop-color="#eeeeee" offset="5%" />
        <stop stop-color="#eeeeb0" offset="95%" />
    </linearGradient>
    <clipPath id="$fig_id-image-frame">
      <rect id="$fig_id-clip-rect" x="0" y="0" width="$(width)" height="$(height)" />
    </clipPath>
</defs>
<style type="text/css">
    rect[rx]:hover { stroke:black; stroke-width:1; }
    text:hover { stroke:black; stroke-width:1; stroke-opacity:0.35; }
</style>
<g id="$fig_id-frame" clip-path="url(#$fig_id-image-frame)">
<rect class="pvbackground" x="0.0" y="0" width="$(width).0" height="$(height).0" fill="url(#background)"  />
<text class="pvbackground" text-anchor="middle" x="600" y="24" font-size="17" font-family="$(font)" fill="rgb(0,0,0)"  >Profile results</text>
<text text-anchor="left" x="10" y="$y_msg" font-size="12" font-family="$(font)" fill="rgb(0,0,0)"  >Function:</text>
<text text-anchor="" x="70" y="$y_msg" font-size="12" font-family="$(font)" fill="rgb(0,0,0)" id="$fig_id-details" > </text>
<g id="$fig_id-viewport" transform="scale(1)">
""")
end

function svgfinish(io::IO, fig_id)
    print(io, """
        </g></g>
        <script><![CDATA[$(escape_script(read(snapsvgjs, String)))]]></script>
        <script><![CDATA[
            $(escape_script(read(viewerjs, String)))
            (function (glob, factory) {
                if (typeof require === "function" && typeof define === "function" && define.amd) {
                    require(["Snap.svg", "ProfileSVG"], function (Snap, ProfileSVG) {
                        factory(Snap, ProfileSVG);
                    });
              } else {
                  factory(glob.Snap, glob.ProfileSVG);
              }
            })(window, function (Snap, ProfileSVG) {
                var svg = Snap.select('svg').node;
                var pt = svg.createSVGPoint();

                var fig = {};

                fig.viewport = Snap.select('#$fig_id-viewport');
                fig.frame = Snap.select('#$fig_id-frame');

                fig.viewport_cx = fig.viewport.getBBox().cx;

                fig.rects = Snap.selectAll('#$fig_id-viewport rect');
                fig.texts = Snap.selectAll('#$fig_id-viewport text');
                fig.clip = Snap.select('#$fig_id-clip-rect');
                fig.clip_width  = fig.clip.getBBox().w;
                fig.clip_middle = fig.clip.getBBox().cy;
                fig.details = document.getElementById("$fig_id-details").firstChild;

                fig.scale = 1.0;
                fig.shift = fig.viewport_cx;

                ProfileSVG.reset(fig)

                fig.rects.forEach(function(rect, i){
                    rect.dblclick(function(){
                        bbox = rect.getBBox();
                        ProfileSVG.move_and_zoom(bbox.cx, bbox.cx, fig.clip_width/bbox.w, fig);
                    })
                    .mouseover(function(){
                        fig.details.nodeValue = rect.node.getAttribute("data-info");
                    })
                    .mouseout(function(){
                        fig.details.nodeValue = "";
                    });

                })
                fig.texts.forEach(function(text, i){
                    text.dblclick(function(){
                        bbox = fig.rects[i].getBBox();
                        ProfileSVG.move_and_zoom(bbox.cx, bbox.cx, fig.clip_width/bbox.w, fig);
                    })
                    .mouseover(function(){
                        fig.details.nodeValue = fig.rects[i].node.getAttribute("data-info");
                    })
                    .mouseout(function(){
                        fig.details.nodeValue = "";
                    });
                })
                Snap.selectAll("#$fig_id-frame .pvbackground").forEach(function(bg){
                   bg.dblclick(function(e){
                       ProfileSVG.reset(fig);
                    });
                });

                function throttle(delay, callback) {
                    var previousCall = new Date().getTime();
                    return function() {
                        var time = new Date().getTime();

                        if ((time - previousCall) >= delay) {
                            previousCall = time;
                            callback.apply(null, arguments);
                        }
                        else {
                            arguments[0].preventDefault();
                        }
                    };
                }

                var MouseWheelHandler = throttle(400, function(e) {
                    e.preventDefault();
                    var e = window.event || e;
                    var delta = Math.max(-1, Math.min(1, (e.wheelDelta || -e.detail)));
                    pt.x = e.clientX;
                    pt.y = e.clientY;

                    pt.matrixTransform(fig.viewport.node.getScreenCTM().inverse());
                    var targetScale = fig.scale + 0.2*delta;
                    ProfileSVG.move_and_zoom(fig.shift, pt.x, targetScale, fig, 400);
                    return false;
                })
                frame = fig.frame.node;
                if (frame.addEventListener) {
                    frame.addEventListener("mousewheel", MouseWheelHandler, false);
                    frame.addEventListener("DOMMouseScroll", MouseWheelHandler, false);
                }
                else {
                    frame.attachEvent("onmousewheel", MouseWheelHandler);
                }

            fig.viewport.drag();
    }); ]]></script>
    </svg>
    """)
end


text:
# ENV["PYCALL_JL_RUNTIME_PYTHON"] = "venv/bin/python"

import Pkg; Pkg.activate(joinpath(@__DIR__, "../../../../../"))
include(joinpath(@__DIR__, "../../../../PlotUtils/imports.jl"))
using DataFrames
using CSV
using PyCall
TSNE = PyCall.pyimport("sklearn.manifold").TSNE
const rfam = CytofRepFAM.Model

function impute(y; yQuantiles=[.0, .25, .5], pBounds=[.05, .8, .05], seed=1)
  I = length(y)

  # Missing mechanism parameters
  beta = hcat([rfam.gen_beta_est(vec(y[i]), yQuantiles, pBounds)
               for i in 1:I]...)

  # Mean of imputed values
  missMean = [beta[2, i] / (-2 * beta[1, i]) for i in 1:I]

  # Preimpute values
  y_imputed = deepcopy(y)
  for i in 1:I
    rfam.preimpute!(y_imputed[i], missMean[i])
  end

  return y_imputed
end

function get_best_lam(samples)
  n = length(samples)
  Ws = [s[:theta__W] for s in samples]
  Zs = [s[:theta__Z] for s in samples]
  lams = [s[:theta__lam] for s in samples]
  I = length(lams[1])

  idx_best = [estimate_ZWi_index(Zs, Ws, i) for i in 1:I]
  return [lams[idx_best[i]][i] for i in 1:I]
end

function compute_tsne(path_to_simdat; use_complete_data=true, seed=0, verbose=2)
  simdat = BSON.load(path_to_simdat)[:simdat]

  if use_complete_data
    y = Matrix{Float64}.(simdat[:y_complete])
  else
    y = impute(simdat[:y])
  end

  # number of samples in y
  I = length(y)

  # Create TSNE objects
  tsne = [TSNE(verbose=verbose, random_state=seed) for _ in 1:I]

  # Fit TSNE for each sample
  foreach(i -> tsne[i].fit(y[i]), 1:I)

  return tsne
end

function compute_combined_tsne(path_to_simdat; use_complete_data=true, seed=0,
                               verbose=2, digits=5)
  # Load simulated data
  simdat = BSON.load(path_to_simdat)[:simdat]

  # Missing data indicator
  m = [isnan.(yi) for yi in simdat[:y]]

  # Determine what to do if missing data is present.
  if use_complete_data
    y = Matrix{Float64}.(simdat[:y_complete])
  else
    y = impute(simdat[:y])
  end

  # Missing indicator combined
  M = vcat(m...)

  # Combine samples into one matrix
  Y = vcat(y...)

  # number of samples in y
  I = length(y)

  # Sample sizes
  N = size.(y, 1)

  # Sample indicators
  sample_ind = vcat([fill(i, N[i]) for i in 1:I]...)

  # Create TSNE object
  tsne = TSNE(verbose=verbose, random_state=seed)

  # Fit TSNE for samples combined
  tsne.fit(Y)

  # Get true labels
  true_labels = vcat(simdat[:lam]...)

  return (round.(tsne.embedding_, digits=digits), sample_ind,
          round.(Y, digits=digits), Int.(M), true_labels)
end

function gen_imgdir(path_to_output, suffix)
  d = splitdir(path_to_output)[1]
  d1, d2 = splitdir(d)
  return joinpath(d1, "tsne", d2, suffix)
end


text:
"""
# `IASP91`

1-D model of velocity.

**N.B.** iasp91 does not contain density.

This version of the model uses the original coefficients and is a
`PREMPolyModel`  as published.  Other implementations available
may use a layered piecewise linear model and/or add density from
other sources.

## Reference

Kennett  B.L.N.  Engdahl  E.  1991. Traveltimes for global earthquake location
and phase identification. Geophys J Int 105  429–465.
https://doi.org/10.1111/j.1365-246X.1991.tb06724.x
"""
const IASP91 = PREMPolyModel(
	r = [   1217.1,     3482,     3631,     5611,     5711,     5961,     6161,     6251,      6336, 6351, 6371],
	vp = [11.24094  10.03904  14.49470  25.14860  25.96984  29.38896  30.78765  25.41389   8.787541   6.5   5.8;
	       0.0       3.75665  -1.47089 -41.15380 -16.93412 -21.40656 -23.25415 -17.69722  -0.749530   0.0   0.0;
	      -4.09689 -13.67046   0.0      51.99320   0.0       0.0       0.0       0.0       0.0        0.0   0.0;
	       0.0       0.0       0.0     -26.60830   0.0       0.0       0.0       0.0       0.0        0.0   0.0],
	vs = [ 3.56454   0.0       8.16616  12.93030  20.76890  17.70732  15.24213   5.75020   6.706231   3.75  3.36;
	       0.0       0.0      -1.58206 -21.25900 -16.53147 -13.50652 -11.08552  -1.27420  -2.248585   0.0   0.0;
	      -3.45241   0.0       0.0      27.89880   0.0       0.0       0.0       0.0       0.0        0.0   0.0;
	       0.0       0.0       0.0     -14.10800   0.0       0.0       0.0       0.0       0.0        0.0   0.0],
	  )


text:
"""
    Pretty

Tablize Vector{<:NamedTuple}
"""
module Pretty # Octo

settings = Dict{Symbol, Any}(
    :pretty => true,
    :nrows => 20,
    :colsize => 30,
)

const pretty_settings_keys = (:nrows, :colsize)

"""
    Pretty.set(pretty::Bool=true; kwargs...)

Set the display options for `Vector{<:NamedTuple}` rows.

```julia-repl
julia> Pretty.set(nrows = 10)     # limit number of the rows

julia> Pretty.set(colsize = 10)   # limit column size

julia> Pretty.set(false)          # do or don't use pretty

```
"""
function set(pretty::Bool=true; kwargs...)
    for (k, v) in kwargs
        if k in pretty_settings_keys && v isa Int && v > 0
            settings[k] = v
        end
    end
    settings[:pretty] = pretty
    settings
end

"""
    Pretty.table(nts::Vector{<:NamedTuple})::String
"""
function table(nts::Vector{<:NamedTuple})::String
    buf = IOBuffer()
    show(buf, MIME"text/plain"(), nts)
    String(take!(buf))
end

function _regularize_text(str::String, padding::Int)::String
    s = escape_string(str)
    if textwidth(s) < padding
        padded_str = rpad(s, padding)
        if textwidth(padded_str) > padding
        else
            return padded_str
        end
    end
    n = 0
    a = []
    for (idx, x) in enumerate(s)
        n += textwidth(x)
        push!(a, x)
        if n > padding - 1
            break
        end
    end
    newstr = join(a)
    newdiff = textwidth(s) - textwidth(newstr)
    if length(s) == length(a)
        news = newstr
        npad = padding - textwidth(newstr)
    else
        if newdiff > 0 && length(a) >= 2
            newstr = join(a[1:end-2])
            newpad = padding - textwidth(newstr)
            news = string(newstr, fill('.', newpad)...)
        else
            news = newstr
        end
        npad = padding - textwidth(news)
    end
    string(news, npad > 0 ? join(fill(' ', npad)) : "")
end

function _print_named_tuple_vector(io::IO, nts::Vector{<:NamedTuple})
    limit_nrows = settings[:nrows]
    limit_colsize = settings[:colsize]
    header_spike(tree, bold) = printstyled(io, tree, bold=bold)
    row_spike(tree, bold)    = printstyled(io, tree, bold=bold)
    colnames = collect(first(typeof(nts).parameters).names)
    ncols = length(colnames)
    real_nrows = length(nts)
    function print_header(pad_, paddings_)
        (boldone, boldtwo) = (false, false)
        header_spike("| ", boldone)
        for colidx in 1:ncols
            printstyled(io, pad_(colidx, colnames[colidx]), color=:cyan)
            ncols != colidx && header_spike(" | ", boldone)
        end
        header_spike(" |\n", boldone)
        header_spike("| ", boldtwo)
        for colidx in 1:ncols
            header_spike(join(fill('-', paddings_[colidx])), boldtwo)
            ncols != colidx && header_spike(" | ", boldtwo)
        end
        header_spike(" |\n", boldtwo)
    end
    function fetched_info()
        printstyled(io, "\n")
        printstyled(io, real_nrows, color=:cyan)
        printstyled(io, " row", real_nrows == 1 ? "" : 's', '.',  real_nrows > limit_nrows ? ".." :  "")
    end
    if isempty(nts)
        paddings = (length ∘ string).(colnames) .+ 2
        print_header((colidx, el) -> rpad(el, paddings[colidx]), paddings)
        fetched_info()
        return
    end
    uno = first(nts)
    nrows = min(limit_nrows, real_nrows)
    A = vcat(map(v -> vcat(v...), nts[1:nrows])...)
    rt = reshape(A, ncols, nrows)
    paddings = maximum((length ∘ string).(rt), dims=2)
    paddings = [maximum(x) for x in zip(paddings, (length ∘ string).(colnames))] .+ 2
    for (colidx, el) in enumerate(values(uno))
        if el isa String
            paddings[colidx] = min(max(limit_colsize, (length ∘ string)(colnames[colidx])), paddings[colidx])
        end
    end
    padfuncs = (x -> x isa Number ? lpad : rpad).(values(uno))
    function pad(colidx, el)
        f = padfuncs[colidx]
        padding = paddings[colidx]
        if el isa String
            _regularize_text(el, padding)
        else
            f(el, padding)
        end
    end
    print_header(pad, paddings)
    for rowidx in 1:nrows
        row = nts[rowidx]
        bold = isodd(rowidx)
        row_spike("| ", bold)
        for (colidx, el) in enumerate(values(row))
            print(io, pad(colidx, el))
            ncols != colidx && row_spike(" | ", false)
        end
        row_spike(" |", false)
        nrows != rowidx && println(io)
    end
    fetched_info()
end

function Base.show(io::IO, mime::MIME"text/plain", nts::Vector{<:NamedTuple})
   if settings[:pretty]
       _print_named_tuple_vector(io, nts)
   else
       Base.show(io, nts)
   end
end

end # module Octo.Pretty


text:
module SatelliteToolbox

export JD_J2000, R0, Rm, m0, J2, Rs, ne, au2m, sunRad
export a_wgs84, b_wgs84, f_wgs84, e_wgs84, el_wgs84

import Base: asin, atan, copy, cos, deepcopy, getindex, mod, setindex!, sin,
       show

using Crayons
using Dates
using DelimitedFiles
using Interpolations
using LinearAlgebra
using OptionalData
using Parameters
using PolynomialRoots
using PrettyTables
using Printf
using ReferenceFrameRotations
using RemoteFiles
using StaticArrays
using SparseArrays
using Statistics
using Reexport

# Re-exporting symbols from ReferenceFrameRotations.jl.
export DCM
export Quaternion

################################################################################
#                                  Submodules
################################################################################

include("submodules/SatelliteToolboxTLE/SatelliteToolboxTLE.jl")
@reexport using .SatelliteToolboxTLE

include("submodules/SGP4/SGP4.jl")
@reexport using .SGP4

################################################################################
#                             Types and Structures
################################################################################

include("types.jl")

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

include("constants.jl")

# Pre-defined crayons.
const _reset_crayon = Crayon(reset = true)
const _crayon_bold  = crayon"bold"
const _crayon_g     = crayon"bold green"
const _crayon_u     = crayon"bold blue"
const _crayon_y     = crayon"bold yellow"

# Escape sequences related to the crayons.
const _b = _crayon_bold
const _d = _reset_crayon
const _g = _crayon_g
const _y = _crayon_y
const _u = _crayon_u

################################################################################
#                                  Exceptions
################################################################################

################################################################################
#                                    Files
################################################################################

include("analysis/beta_angle.jl")
include("analysis/eclipse_time.jl")
include("analysis/lighting_conditions.jl")
include("analysis/ground_stations.jl")
include("analysis/ground_trace.jl")
include("analysis/raan.jl")
include("analysis/payload_optical_analysis.jl")
include("analysis/satellite_position_countries.jl")
include("analysis/sun_angle.jl")
include("analysis/sun_radiation.jl")

include("earth/atmospheric_models/expatmosphere/expatmosphere.jl")
include("earth/atmospheric_models/jb2008/jb2008.jl")
include("earth/atmospheric_models/jr1971/jr1971.jl")
include("earth/atmospheric_models/nrlmsise00/nrlmsise00.jl")
include("earth/gravity_models/embedded_gravity_models.jl")
include("earth/gravity_models/gravity_model.jl")
include("earth/geomagnetic_field_models/dipole/geomag_dipole.jl")
include("earth/geomagnetic_field_models/dipole/geomag_dipole_constants.jl")
include("earth/geomagnetic_field_models/igrf/igrf.jl")
include("earth/geomagnetic_field_models/igrf/igrf12_coefs.jl")
include("earth/geomagnetic_field_models/igrf/igrf12syn_coefs.jl")
include("earth/space_indices/space_indices.jl")

include("./misc/crossing.jl")
include("./misc/legendre.jl")
include("./misc/dlegendre.jl")
include("./misc/icgem.jl")

include("./moon/moon_position.jl")

include("sun/equation_of_time.jl")
include("sun/sun_position.jl")
include("sun/sun_velocity.jl")

include("orbit/general.jl")
include("orbit/anomalies.jl")
include("orbit/orbit_sun_sync.jl")
include("orbit/orbit_sun_sync_ground_reap.jl")
include("orbit/state_vector.jl")
include("orbit/propagators/j2.jl")
include("orbit/propagators/j4.jl")
include("orbit/propagators/twobody.jl")
include("orbit/propagators/api/epoch.jl")
include("orbit/propagators/api/init_orbit_propagator.jl")
include("orbit/propagators/api/propagate.jl")
include("orbit/propagators/api/propagate_to_epoch.jl")
include("orbit/propagators/api/step.jl")

include("transformations/eop.jl")
include("transformations/ecef_to_ecef.jl")
include("transformations/ecef_to_eci.jl")
include("transformations/eci_to_ecef.jl")
include("transformations/eci_to_eci.jl")
include("transformations/geodetic_geocentric.jl")
include("transformations/gmst.jl")
include("transformations/misc.jl")
include("transformations/orbit_elements.jl")
include("transformations/position.jl")
include("transformations/sv_ecef_to_ecef.jl")
include("transformations/sv_ecef_to_eci.jl")
include("transformations/sv_eci_to_ecef.jl")
include("transformations/sv_eci_to_eci.jl")

include("transformations/fk5/fk5.jl")
include("transformations/fk5/nutation.jl")
include("transformations/fk5/precession.jl")

include("transformations/iau2006/iau2006_const.jl")
include("transformations/iau2006/iau2006.jl")
include("transformations/iau2006/precession_nutation_iau2006.jl")

include("transformations/teme/teme.jl")

include("time/julian_day.jl")
include("time/time.jl")

include("deprecations.jl")

end # module


text:
function label(colors::AbstractArray)

    # Set the color-number converter list
    color_key = ["black", "brown", "red", "orange", "yellow", "green", "blue", "violet", "grey", "white"]

    # Initialize the array that hold the converted numbers
    converted_colors = collect(1:3)

    # For every color, put the converted number in the converted_colors array
    for color_index in 1:3
        converted_colors[color_index] = (findfirst(x -> x == colors[color_index], color_key) - 1)
    end

    # If the second digit ends with zero, end the first two digits at the first digit, and increment the power of ten
    # else, keep the two digits and convert the third number into a power of 10
    if converted_colors[2] == 0
        first_two_digits = converted_colors[1]
        converted_colors[3] += 1
        power_of_ten = 10^converted_colors[3]
    else
        first_two_digits = 10*converted_colors[1] + converted_colors[2]
        power_of_ten = 10^converted_colors[3]
    end

    # If the resistance is an integral multiple of 1000 ohms, write the answer in terms of kiloohms
    # Else, write the answer in terms of ohms
    if power_of_ten >= 1000
        power_of_ten = 10^(converted_colors[3]-3)
        suffix = "kiloohms"
        first_two_digits *= power_of_ten
    else
        first_two_digits *= power_of_ten
        suffix = "ohms"
    end

    # print the final anwer with the appropriate suffix
    return ("$first_two_digits" * " " * "$suffix")

end