text:
function preprocess_index(natom::Int, index::AbstractVector)::AbstractVector
    if isempty(index)
        index2 = collect(1:natom)
    elseif eltype(index) <: Bool
        index2 = findall(index)
    elseif typeof(index) <: UnitRange{Int}
        index2 = collect(index)
    else
        index2 = index
    end
    if natom < maximum(index2)
        @printf "Warning: some indices are larger than # of atoms, truncated.\n"
        id = index2 .<= natom
        index2 = index2[id]
    end
    if isempty(index2)
        error("no atom can be selected by the given index")
    end
    return index2
end

function preprocess_weight(natom::Int, isweight::Bool, mass::Vector{T}) where {T}
    if isweight
        if isempty(mass)
            #@printf "Warning: not weighted because masses are not provided.\n"
            weight = ones(T, natom)
        else
            weight = mass
        end
    else
        weight = ones(T, natom)
    end
    return weight
end

"""
    centerofmass(ta::TrjArray; isweight=true, index::Union{UnitRange,Vector,BitArray}) -> com::TrjArray

Calculates the center of mass coordinates (COM) of the given trajectory, a `TrjArray` object `ta`.
If `isweight` is `true` (default), coordinates are weighted by `ta.mass` (as long as `ta.mass` is not empty).
Uses can also specify column indices for the COM calculation by giving `index` which should be a Vector of integers, a UnitRange, or a BiArray.

Returns the center of mass coordinates as a TrjArray object, which have virtual single "atom" whose coordinates are COMs.

# Example
```julia-repl
julia> ta = mdload("ak.pdb")
julia> com = centersofmass(ta)
```
"""
function centerofmass(ta::TrjArray{T, U};
    isweight::Bool=true, index::AbstractVector=Vector{U}(undef, 0))::TrjArray{T, U} where {T, U}
    #start = time()
    nframe = ta.nframe
    natom = ta.natom
    index2 = preprocess_index(natom, index)
    natom_sub = length(index2)
    index3 = to3(index2)
    weight = preprocess_weight(natom, isweight, ta.mass)
    weight = weight[index2]
    wsum_inv = one(T) / sum(weight)
    weight = reshape(weight, 1, natom_sub)
    #elapsed = time() - start; println("elapsed1 = $(elapsed)")

    #start = time()
    xyz = zeros(T, nframe, 3)
    if natom_sub == 1
        return ta[:, index2]
    else
        xyz[:, 1:1] .= sum(weight .* view(ta.xyz, :, index3[1:3:end]), dims=2) .* wsum_inv
        xyz[:, 2:2] .= sum(weight .* view(ta.xyz, :, index3[2:3:end]), dims=2) .* wsum_inv
        xyz[:, 3:3] .= sum(weight .* view(ta.xyz, :, index3[3:3:end]), dims=2) .* wsum_inv
        #x = view(ta.xyz, :, index3[1:3:end])
        #y = view(ta.xyz, :, index3[2:3:end])
        #z = view(ta.xyz, :, index3[3:3:end])
        #@time xyz[:, 1:1] .= sum(weight .* x, dims=2) .* wsum_inv
        #@time xyz[:, 2:2] .= sum(weight .* y, dims=2) .* wsum_inv
        #@time xyz[:, 3:3] .= sum(weight .* z, dims=2) .* wsum_inv
    end
    #elapsed = time() - start; println("elapsed2 = $(elapsed)")

    TrjArray{T, U}(xyz=xyz)
end

"""
    decenter(ta::TrjArray; isweight=true, index::Union{UnitRange,Vector,BitArray}) -> ta_decentered::TrjArray, com::TrjArray

Calculates the centers of mass coordinates (COM) of the given trajectory, a TrjArray object ta.
And translates the coordinates so that the COM of the output is identical to the origin (x=y=z=0).
If `isweight` is `true` (default), coordinates are weighted by `ta.mass` (as long as `ta.mass` is not empty).
Uses can also specify column indices for the COM calculation by giving `index` which should be a Vector of integers, a UnitRange, or a BiArray.

Returns a TrjArray object whose COM is the origin (x=y=z=0), and the COM of the given TrjArray `ta`.

# Example
```julia-repl
julia> ta = mdload("ak.pdb")
julia> ta_decentered = decenter(ta)
```
"""
function decenter(ta::TrjArray{T, U}; isweight::Bool=true,
    index::AbstractVector=Vector{Int64}(undef, 0))::Tuple{TrjArray{T, U}, TrjArray{T, U}} where {T, U}
    com = centerofmass(ta, isweight=isweight, index=index)
    xyz = similar(ta.xyz)
    xyz[:, 1:3:end] .= ta.xyz[:, 1:3:end] .- com.xyz[:, 1:1]
    xyz[:, 2:3:end] .= ta.xyz[:, 2:3:end] .- com.xyz[:, 2:2]
    xyz[:, 3:3:end] .= ta.xyz[:, 3:3:end] .- com.xyz[:, 3:3]
    TrjArray(ta, xyz=xyz), com
end

"""
    decenter!(ta::TrjArray; isweight=true, index::Union{UnitRange,Vector,BitArray}) -> com::TrjArray

Translates the coordinates of the given trajectory, a TrjArray `ta`,
so that the COM of it become identical to the origin (x=y=z=0).
If `isweight` is `true` (default), coordinates are weighted by `ta.mass` (as long as `ta.mass` is not empty).
Uses can also specify column indices for the COM calculation by giving `index` which should be an Array of integers.

Returns the center of mass coordinates of the given trajectory.

# Example
```julia-repl
julia> ta = mdload("ak.pdb")
julia> decenter!(ta)
```
"""
function decenter!(ta::TrjArray{T, U}; isweight::Bool=true,
    index::AbstractVector=Vector{Int64}(undef, 0)) where {T, U}
    com = centerofmass(ta, isweight=isweight, index=index)
    natom3 = ta.natom*3
    x = view(ta.xyz, :, 1:3:natom3)
    y = view(ta.xyz, :, 2:3:natom3)
    z = view(ta.xyz, :, 3:3:natom3)
    x .= x .- com.xyz[:, 1:1]
    y .= y .- com.xyz[:, 2:2]
    z .= z .- com.xyz[:, 3:3]
    return com
end

"""
    orient!(ta::TrjArray, index::Union{UnitRange,Vector,BitArray})

Orient the molecule using its principal axes of inertia.

# Example
```julia-repl
julia> ta = mdload("ak.pdb")
julia> orient!(ta)
```
"""
function orient!(ta::TrjArray{T, U}; index::AbstractVector=Vector{U}(undef, 0)) where {T, U}
    natom = ta.natom
    decenter!(ta, index=index)

    index2 = preprocess_index(natom, index)
    ta_sub = ta[:, index2]
    natom_sub = ta_sub.natom
    natom3_sub = 3*natom_sub
    if isempty(ta_sub.mass)
        mass = ones(T, ta_sub.natom)
    else
        mass = ta_sub.mass
    end

    I = zeros(T, 3, 3)
    for iframe = 1:ta.nframe
        x = view(ta_sub.xyz, iframe, 1:3:natom3_sub)
        y = view(ta_sub.xyz, iframe, 2:3:natom3_sub)
        z = view(ta_sub.xyz, iframe, 3:3:natom3_sub)

        I[1, 1] = sum(mass.*(y.^2 + z.^2));
        I[2, 2] = sum(mass.*(x.^2 + z.^2));
        I[3, 3] = sum(mass.*(x.^2 + y.^2));

        I[1, 2] = - sum(mass.*(x.*y));
        I[2, 1] = I[1, 2]

        I[1, 3] = - sum(mass.*(x.*z));
        I[3, 1] = I[1, 3]

        I[2, 3] = - sum(mass.*(y.*z));
        I[3, 2] = I[2, 3]

        F = svd(I)
        p_axis = F.Vt[end:-1:1, :]; #z-axis has the largest inertia

        # check reflection
        if det(p_axis) < 0.0
            p_axis[1, :] .= - p_axis[1, :]
        end

        # project onto the principal axis of inertia
        proj = p_axis * reshape(ta.xyz[iframe, :], 3, natom)
        ta.xyz[iframe, 1:3:end] .= proj[1, :]
        ta.xyz[iframe, 2:3:end] .= proj[2, :]
        ta.xyz[iframe, 3:3:end] .= proj[3, :]
    end
end

###### superimpose #################
function innerproduct(iframe::U, ref_xyz::Matrix{T}, ta_xyz::Matrix{T},
                       index::Vector{U}, isweight::Bool, mass::Vector{T}) where {T, U}
    A = zeros(T, 9)
    # A = Vector{Float64}(undef, 9)
    # A[:] .= 0.0
    G1 = G2 = zero(T)
    if isweight
        for i in index
            ref_x = ref_xyz[1, 3*(i-1)+1]
            ref_y = ref_xyz[1, 3*(i-1)+2]
            ref_z = ref_xyz[1, 3*(i-1)+3]
            x1 = mass[i] * ref_x
            y1 = mass[i] * ref_y
            z1 = mass[i] * ref_z
            G1 += x1 * ref_x + y1 * ref_x + z1 * ref_z
            x2 = ta_xyz[iframe, 3*(i-1)+1]
            y2 = ta_xyz[iframe, 3*(i-1)+2]
            z2 = ta_xyz[iframe, 3*(i-1)+3]
            G2 += mass[i] * (x2 * x2 + y2 * y2 + z2 * z2)
            A[1] +=  (x1 * x2)
            A[2] +=  (x1 * y2)
            A[3] +=  (x1 * z2)
            A[4] +=  (y1 * x2)
            A[5] +=  (y1 * y2)
            A[6] +=  (y1 * z2)
            A[7] +=  (z1 * x2)
            A[8] +=  (z1 * y2)
            A[9] +=  (z1 * z2)
        end
    else
        for i in index
            ref_x = ref_xyz[1, 3*(i-1)+1]
            ref_y = ref_xyz[1, 3*(i-1)+2]
            ref_z = ref_xyz[1, 3*(i-1)+3]
            x1 = ref_x
            y1 = ref_y
            z1 = ref_z
            G1 += x1 * ref_x + y1 * ref_y + z1 * ref_z
            x2 = ta_xyz[iframe, 3*(i-1)+1]
            y2 = ta_xyz[iframe, 3*(i-1)+2]
            z2 = ta_xyz[iframe, 3*(i-1)+3]
            G2 += (x2 * x2 + y2 * y2 + z2 * z2)
            A[1] +=  (x1 * x2)
            A[2] +=  (x1 * y2)
            A[3] +=  (x1 * z2)
            A[4] +=  (y1 * x2)
            A[5] +=  (y1 * y2)
            A[6] +=  (y1 * z2)
            A[7] +=  (z1 * x2)
            A[8] +=  (z1 * y2)
            A[9] +=  (z1 * z2)
        end
    end
    return A, (G1 + G2) * T(0.5)
end

"""
this code is licensed under the BSD license (Copyright (c) 2009-2016 Pu Liu and Douglas L. Theobald), see LICENSE.md
"""
function fastCalcRMSDAndRotation(A::Vector{T}, E0::T, wsum_inv::T) where {T}
    oldg = zero(T)
    C = zeros(T, 3)
    rot = zeros(T, 9)
    evecprec = T(1e-6)
    evalprec = T(1e-11)

    Sxx = A[1]; Sxy = A[2]; Sxz = A[3]
    Syx = A[4]; Syy = A[5]; Syz = A[6]
    Szx = A[7]; Szy = A[8]; Szz = A[9]

    Sxx2 = Sxx * Sxx
    Syy2 = Syy * Syy
    Szz2 = Szz * Szz

    Sxy2 = Sxy * Sxy
    Syz2 = Syz * Syz
    Sxz2 = Sxz * Sxz

    Syx2 = Syx * Syx
    Szy2 = Szy * Szy
    Szx2 = Szx * Szx

    SyzSzymSyySzz2 = 2*(Syz*Szy - Syy*Szz)
    Sxx2Syy2Szz2Syz2Szy2 = Syy2 + Szz2 - Sxx2 + Syz2 + Szy2

    C[3] = -2 * (Sxx2 + Syy2 + Szz2 + Sxy2 + Syx2 + Sxz2 + Szx2 + Syz2 + Szy2)
    C[2] = 8 * (Sxx*Syz*Szy + Syy*Szx*Sxz + Szz*Sxy*Syx - Sxx*Syy*Szz - Syz*Szx*Sxy - Szy*Syx*Sxz)

    SxzpSzx = Sxz + Szx
    SyzpSzy = Syz + Szy
    SxypSyx = Sxy + Syx
    SyzmSzy = Syz - Szy
    SxzmSzx = Sxz - Szx
    SxymSyx = Sxy - Syx
    SxxpSyy = Sxx + Syy
    SxxmSyy = Sxx - Syy
    Sxy2Sxz2Syx2Szx2 = Sxy2 + Sxz2 - Syx2 - Szx2

    C[1] = Sxy2Sxz2Syx2Szx2 * Sxy2Sxz2Syx2Szx2 +
         (Sxx2Syy2Szz2Syz2Szy2 + SyzSzymSyySzz2) * (Sxx2Syy2Szz2Syz2Szy2 - SyzSzymSyySzz2) +
         (-(SxzpSzx)*(SyzmSzy)+(SxymSyx)*(SxxmSyy-Szz)) * (-(SxzmSzx)*(SyzpSzy)+(SxymSyx)*(SxxmSyy+Szz)) +
         (-(SxzpSzx)*(SyzpSzy)-(SxypSyx)*(SxxpSyy-Szz)) * (-(SxzmSzx)*(SyzmSzy)-(SxypSyx)*(SxxpSyy+Szz)) +
         (+(SxypSyx)*(SyzpSzy)+(SxzpSzx)*(SxxmSyy+Szz)) * (-(SxymSyx)*(SyzmSzy)+(SxzpSzx)*(SxxpSyy+Szz)) +
         (+(SxypSyx)*(SyzmSzy)+(SxzmSzx)*(SxxmSyy-Szz)) * (-(SxymSyx)*(SyzpSzy)+(SxzmSzx)*(SxxpSyy-Szz))

    # Newton-Raphson
    mxEigenV = E0
    icount = 0
    for i in 1:50
        icount += 1
        oldg = mxEigenV
        x2 = mxEigenV*mxEigenV
        b = (x2 + C[3])*mxEigenV
        a = b + C[2]
        delta = ((a*mxEigenV + C[1])/(T(2)*x2*mxEigenV + b + a));
        mxEigenV -= delta;
        # printf("\n diff[%3d]: %16g %16g %16g", i, mxEigenV - oldg, evalprec*mxEigenV, mxEigenV)
        if abs(mxEigenV - oldg) < abs(evalprec*mxEigenV)
            break
        end
    end

    if icount == 50
        print("More than 50 iterations needed!")
    end

    # the abs() is to guard against extremely small, but *negative* numbers due to floating point error */
    r1 = sqrt(abs(2 * (E0 - mxEigenV) * wsum_inv))
    # printf("\n\n %16g %16g %16g \n", r1, E0, 2.0 * (E0 - mxEigenV)/len)

    a11 = SxxpSyy + Szz-mxEigenV; a12 = SyzmSzy; a13 = - SxzmSzx; a14 = SxymSyx
    a21 = SyzmSzy; a22 = SxxmSyy - Szz-mxEigenV; a23 = SxypSyx; a24= SxzpSzx
    a31 = a13; a32 = a23; a33 = Syy-Sxx-Szz - mxEigenV; a34 = SyzpSzy
    a41 = a14; a42 = a24; a43 = a34; a44 = Szz - SxxpSyy - mxEigenV
    a3344_4334 = a33 * a44 - a43 * a34; a3244_4234 = a32 * a44-a42*a34
    a3243_4233 = a32 * a43 - a42 * a33; a3143_4133 = a31 * a43-a41*a33
    a3144_4134 = a31 * a44 - a41 * a34; a3142_4132 = a31 * a42-a41*a32
    q1 =  a22*a3344_4334-a23*a3244_4234+a24*a3243_4233
    q2 = -a21*a3344_4334+a23*a3144_4134-a24*a3143_4133
    q3 =  a21*a3244_4234-a22*a3144_4134+a24*a3142_4132
    q4 = -a21*a3243_4233+a22*a3143_4133-a23*a3142_4132

    qsqr = q1 * q1 + q2 * q2 + q3 * q3 + q4 * q4

    # The following code tries to calculate another column in the adjoint matrix when the norm of the
    # current column is too small.
    # Usually this block will never be activated.  To be absolutely safe this should be
    # uncommented, but it is most likely unnecessary.
    while qsqr < evecprec
        q1 =  a12*a3344_4334 - a13*a3244_4234 + a14*a3243_4233
        q2 = -a11*a3344_4334 + a13*a3144_4134 - a14*a3143_4133
        q3 =  a11*a3244_4234 - a12*a3144_4134 + a14*a3142_4132
        q4 = -a11*a3243_4233 + a12*a3143_4133 - a13*a3142_4132
        qsqr = q1*q1 + q2 *q2 + q3*q3+q4*q4

        if qsqr < evecprec
            a1324_1423 = a13 * a24 - a14 * a23; a1224_1422 = a12 * a24 - a14 * a22
            a1223_1322 = a12 * a23 - a13 * a22; a1124_1421 = a11 * a24 - a14 * a21
            a1123_1321 = a11 * a23 - a13 * a21; a1122_1221 = a11 * a22 - a12 * a21

            q1 =  a42 * a1324_1423 - a43 * a1224_1422 + a44 * a1223_1322
            q2 = -a41 * a1324_1423 + a43 * a1124_1421 - a44 * a1123_1321
            q3 =  a41 * a1224_1422 - a42 * a1124_1421 + a44 * a1122_1221
            q4 = -a41 * a1223_1322 + a42 * a1123_1321 - a43 * a1122_1221
            qsqr = q1*q1 + q2 *q2 + q3*q3+q4*q4

            if qsqr < evecprec
                q1 =  a32 * a1324_1423 - a33 * a1224_1422 + a34 * a1223_1322
                q2 = -a31 * a1324_1423 + a33 * a1124_1421 - a34 * a1123_1321
                q3 =  a31 * a1224_1422 - a32 * a1124_1421 + a34 * a1122_1221
                q4 = -a31 * a1223_1322 + a32 * a1123_1321 - a33 * a1122_1221
                qsqr = q1*q1 + q2 *q2 + q3*q3 + q4*q4;

                if qsqr < evecprec
                    # if qsqr is still too small, return the identity matrix.
                    rot[1] = rot[5] = rot[9] = one(T)
                    rot[2] = rot[3] = rot[4] = rot[6] = rot[7] = rot[8] = zero(T)
                    return r1, rot
                end
            end
        end
    end

    normq = sqrt(qsqr)
    q1 /= normq
    q2 /= normq
    q3 /= normq
    q4 /= normq

    a2 = q1 * q1
    x2 = q2 * q2
    y2 = q3 * q3
    z2 = q4 * q4

    xy = q2 * q3
    az = q1 * q4
    zx = q4 * q2
    ay = q1 * q3
    yz = q3 * q4
    ax = q1 * q2

    rot[1] = a2 + x2 - y2 - z2
    rot[2] = 2 * (xy + az)
    rot[3] = 2 * (zx - ay)
    rot[4] = 2 * (xy - az)
    rot[5] = a2 - x2 + y2 - z2
    rot[6] = 2 * (yz + ax)
    rot[7] = 2 * (zx + ay)
    rot[8] = 2 * (yz - ax)
    rot[9] = a2 - x2 - y2 + z2
    return r1, rot
end

#function applyrotation!(iframe::U, x::AbstractArray{T, 2}, y::AbstractArray{T, 2}, z::AbstractArray{T, 2}, xyz::Matrix{T}, rot::Vector{T}) where {T, U}
#    natom = size(x, 2)
#    for iatom = 1:natom
#        x[iframe, iatom] = rot[1] * xyz[iframe, 3*(iatom-1)+1] + rot[2] * xyz[iframe, 3*(iatom-1)+2] + rot[3] * xyz[iframe, 3*(iatom-1)+3]
#        y[iframe, iatom] = rot[4] * xyz[iframe, 3*(iatom-1)+1] + rot[5] * xyz[iframe, 3*(iatom-1)+2] + rot[6] * xyz[iframe, 3*(iatom-1)+3]
#        z[iframe, iatom] = rot[7] * xyz[iframe, 3*(iatom-1)+1] + rot[8] * xyz[iframe, 3*(iatom-1)+2] + rot[9] * xyz[iframe, 3*(iatom-1)+3]
#    end
#end

function applyrotation!(iframe::U, x::AbstractArray{T, 2}, y::AbstractArray{T, 2}, z::AbstractArray{T, 2}, xyz::Matrix{T}, rot::Vector{T}) where {T, U}
    natom = size(x, 2)
    @simd for iatom = 1:natom
        x[iframe, iatom] = rot[1] * xyz[iframe, 3*(iatom-1)+1] + rot[2] * xyz[iframe, 3*(iatom-1)+2] + rot[3] * xyz[iframe, 3*(iatom-1)+3]
        y[iframe, iatom] = rot[4] * xyz[iframe, 3*(iatom-1)+1] + rot[5] * xyz[iframe, 3*(iatom-1)+2] + rot[6] * xyz[iframe, 3*(iatom-1)+3]
        z[iframe, iatom] = rot[7] * xyz[iframe, 3*(iatom-1)+1] + rot[8] * xyz[iframe, 3*(iatom-1)+2] + rot[9] * xyz[iframe, 3*(iatom-1)+3]
    end
end

function applyrotation!(iframe::U, xyz::Matrix{T}, rot::Vector{T}) where {T, U}
    natom = length(xyz[iframe, 1:3:end])
    @simd for iatom = 1:natom
        x = rot[1] * xyz[iframe, 3*(iatom-1)+1] + rot[2] * xyz[iframe, 3*(iatom-1)+2] + rot[3] * xyz[iframe, 3*(iatom-1)+3]
        y = rot[4] * xyz[iframe, 3*(iatom-1)+1] + rot[5] * xyz[iframe, 3*(iatom-1)+2] + rot[6] * xyz[iframe, 3*(iatom-1)+3]
        z = rot[7] * xyz[iframe, 3*(iatom-1)+1] + rot[8] * xyz[iframe, 3*(iatom-1)+2] + rot[9] * xyz[iframe, 3*(iatom-1)+3]
        xyz[iframe, 3*(iatom-1)+1] = x
        xyz[iframe, 3*(iatom-1)+2] = y
        xyz[iframe, 3*(iatom-1)+3] = z
    end
end

#function applyrotation!(iframe::U, x::AbstractArray{T, 2}, y::AbstractArray{T, 2}, z::AbstractArray{T, 2}, xyz::Matrix{T}, rot::Vector{T}) where {T, U}
#    @simd @inbounds x[iframe, :] .= rot[1] .* xyz[iframe, 1:3:end] .+ rot[2] .* xyz[iframe, 2:3:end] .+ rot[3] .* xyz[iframe, 3:3:end]
#    @simd @inbounds y[iframe, :] .= rot[4] .* xyz[iframe, 1:3:end] .+ rot[5] .* xyz[iframe, 2:3:end] .+ rot[6] .* xyz[iframe, 3:3:end]
#    @simd @inbounds z[iframe, :] .= rot[7] .* xyz[iframe, 1:3:end] .+ rot[8] .* xyz[iframe, 2:3:end] .+ rot[9] .* xyz[iframe, 3:3:end]
#end

#function applyrotation!(iframe::U, x::AbstractArray{T, 2}, y::AbstractArray{T, 2}, z::AbstractArray{T, 2}, x_ta::AbstractArray{T, 2}, y_ta::AbstractArray{T, 2}, z_ta::AbstractArray{T, 2}, rot::Vector{T}) where {T, U}
#    x[iframe, :] .= rot[1] .* x_ta[iframe, :] .+ rot[2] .* y_ta[iframe, :] .+ rot[3] .* z_ta[iframe, :]
#    y[iframe, :] .= rot[4] .* x_ta[iframe, :] .+ rot[5] .* y_ta[iframe, :] .+ rot[6] .* z_ta[iframe, :]
#    z[iframe, :] .= rot[7] .* x_ta[iframe, :] .+ rot[8] .* y_ta[iframe, :] .+ rot[9] .* z_ta[iframe, :]
#end

"""
    rotate(ta::TrjArray, quaternion::Vector) -> ta_rotated::TrjArray

Performs rotation of coordinates of the given TrjArray object `ta`.
Rotation matrix is constructed fromt the given `quaternion` vector.

Returns the trajecotry with rotated coordinates.

# Example
```julia-repl
julia> ta = mdload("ak.dcd")
julia> quaternion = [0.80951112,  0.10657412,  0.35146912,  0.45804312]
julia> ta_rotated = rotate(ta, quaternion)
```
"""
function rotate(ta::TrjArray{T, U}, quater::AbstractVector{T})::TrjArray{T, U} where {T, U}
    xyz = similar(ta.xyz)
    rot = similar(quater, 9)
    rot[1] = 1 - 2 * quater[2] * quater[2] - 2 * quater[3] * quater[3]
    rot[2] = 2 * (quater[1] * quater[2] + quater[3] * quater[4])
    rot[3] = 2 * (quater[1] * quater[3] - quater[2] * quater[4])
    rot[4] = 2 * (quater[1] * quater[2] - quater[3] * quater[4])
    rot[5] = 1 - 2 * quater[1] * quater[1] - 2 * quater[3] * quater[3]
    rot[6] = 2 * (quater[2] * quater[3] + quater[1] * quater[4])
    rot[7] = 2 * (quater[1] * quater[3] + quater[2] * quater[4])
    rot[8] = 2 * (quater[2] * quater[3] - quater[1] * quater[4])
    rot[9] = 1 - 2 * quater[1] * quater[1] - 2 * quater[2] * quater[2]
    x = rot[1] .* ta.xyz[:, 1:3:end] .+ rot[2] .* ta.xyz[:, 2:3:end] .+ rot[3] .* ta.xyz[:, 3:3:end]
    y = rot[4] .* ta.xyz[:, 1:3:end] .+ rot[5] .* ta.xyz[:, 2:3:end] .+ rot[6] .* ta.xyz[:, 3:3:end]
    z = rot[7] .* ta.xyz[:, 1:3:end] .+ rot[8] .* ta.xyz[:, 2:3:end] .+ rot[9] .* ta.xyz[:, 3:3:end]
    xyz[:, 1:3:end] .= x
    xyz[:, 2:3:end] .= y
    xyz[:, 3:3:end] .= z
    return TrjArray(ta, xyz=xyz)
end

"""
    rotate!(ta::TrjArray, quaternion::Vector)

Performs rotation of coordinates of the given TrjArray object `ta`.
Rotation matrix is constructed fromt the given `quaternion` vector.

The coordinates of `ta` is replaced with the rotated ones.

# Example
```julia-repl
julia> ta = mdload("ak.dcd")
julia> quaternion = [0.80951112,  0.10657412,  0.35146912,  0.45804312]
julia> rotate!(ta, quaternion)
```
"""
function rotate!(ta::TrjArray{T, U}, quater::AbstractVector{T}) where {T, U}
    rot = similar(quater, 9)
    rot[1] = 1 - 2 * quater[2] * quater[2] - 2 * quater[3] * quater[3]
    rot[2] = 2 * (quater[1] * quater[2] + quater[3] * quater[4])
    rot[3] = 2 * (quater[1] * quater[3] - quater[2] * quater[4])
    rot[4] = 2 * (quater[1] * quater[2] - quater[3] * quater[4])
    rot[5] = 1 - 2 * quater[1] * quater[1] - 2 * quater[3] * quater[3]
    rot[6] = 2 * (quater[2] * quater[3] + quater[1] * quater[4])
    rot[7] = 2 * (quater[1] * quater[3] + quater[2] * quater[4])
    rot[8] = 2 * (quater[2] * quater[3] - quater[1] * quater[4])
    rot[9] = 1 - 2 * quater[1] * quater[1] - 2 * quater[2] * quater[2]
    x = rot[1] .* ta.xyz[:, 1:3:end] .+ rot[2] .* ta.xyz[:, 2:3:end] .+ rot[3] .* ta.xyz[:, 3:3:end]
    y = rot[4] .* ta.xyz[:, 1:3:end] .+ rot[5] .* ta.xyz[:, 2:3:end] .+ rot[6] .* ta.xyz[:, 3:3:end]
    z = rot[7] .* ta.xyz[:, 1:3:end] .+ rot[8] .* ta.xyz[:, 2:3:end] .+ rot[9] .* ta.xyz[:, 3:3:end]
    ta.xyz[:, 1:3:end] .= x
    ta.xyz[:, 2:3:end] .= y
    ta.xyz[:, 3:3:end] .= z
    return nothing
end

"""
    rotate(ta::TrjArray, quaternions::Matrix) -> ta_rotated::TrjArray

Performs multiple rotations of the coordinates of the 1st frame of the given TrjArray object `ta`.
Rotation matrices are construcuted from the 2-dimensional Array `quaternions`.
One set of quaternions for rotaton should be given in a row vector in `quaternions`.

Returns the trajecotry with rotated coordinates from the 1st frame of `ta`.

# Example
```julia-repl
julia> ta = mdload("ak.pdb")
julia> quaternions = [0.80951112  0.10657412  0.35146912  0.45804312;
>                     0.10657412  0.80951112  -0.35146912  0.45804312]
julia> ta_rotated = rotate(ta, quaternions)
```
"""
function rotate(ta_single::TrjArray{T, U}, quater::AbstractMatrix{T})::TrjArray{T, U} where {T, U}
    rot = similar(quater, size(quater, 1), 9)
    rot[:, 1:1] .= 1 .- 2 .* quater[:, 2:2] .* quater[:, 2:2] .- 2.0 .* quater[:, 3:3] .* quater[:, 3:3]
    rot[:, 2:2] .= 2 .* (quater[:, 1:1] .* quater[:, 2:2] .+ quater[:, 3:3] .* quater[:, 4:4])
    rot[:, 3:3] .= 2 .* (quater[:, 1:1] .* quater[:, 3:3] .- quater[:, 2:2] .* quater[:, 4:4])
    rot[:, 4:4] .= 2 .* (quater[:, 1:1] .* quater[:, 2:2] .- quater[:, 3:3] .* quater[:, 4:4])
    rot[:, 5:5] .= 1 .- 2 .* quater[:, 1:1] .* quater[:, 1:1] .- 2 .* quater[:, 3:3] .* quater[:, 3:3]
    rot[:, 6:6] .= 2 .* (quater[:, 2:2] .* quater[:, 3:3] .+ quater[:, 1:1] .* quater[:, 4:4])
    rot[:, 7:7] .= 2 .* (quater[:, 1:1] .* quater[:, 3:3] .+ quater[:, 2:2] .* quater[:, 4:4])
    rot[:, 8:8] .= 2 .* (quater[:, 2:2] .* quater[:, 3:3] .- quater[:, 1:1] .* quater[:, 4:4])
    rot[:, 9:9] .= 1 .- 2 .* quater[:, 1:1] .* quater[:, 1:1] .- 2 .* quater[:, 2:2] .* quater[:, 2:2]
    x = rot[:, 1:1] .* ta_single.xyz[1:1, 1:3:end] .+ rot[:, 2:2] .* ta_single.xyz[1:1, 2:3:end] .+ rot[:, 3:3] .* ta_single.xyz[1:1, 3:3:end]
    y = rot[:, 4:4] .* ta_single.xyz[1:1, 1:3:end] .+ rot[:, 5:5] .* ta_single.xyz[1:1, 2:3:end] .+ rot[:, 6:6] .* ta_single.xyz[1:1, 3:3:end]
    z = rot[:, 7:7] .* ta_single.xyz[1:1, 1:3:end] .+ rot[:, 8:8] .* ta_single.xyz[1:1, 2:3:end] .+ rot[:, 9:9] .* ta_single.xyz[1:1, 3:3:end]
    xyz = similar(ta_single.xyz[1:1, :])
    xyz[:, 1:3:end] .= x # fixme, todo reduce memory
    xyz[:, 2:3:end] .= y
    xyz[:, 3:3:end] .= z
    return TrjArray(ta_single, xyz=xyz)
end

"""
    rotate_with_matrix(ta::TrjArray, R::Matrix) -> ta_rotated::TrjArray

Performs rotation of coordinates of the given TrjArray object `ta` with the given rotation matrix.

The coordinates of `ta` is replaced with the rotated ones.

# Example
```julia-repl
julia> ta = mdload("ak.dcd")
julia> R = [0.36 0.48 -0.8;
            -0.8 0.60 0.0;
            -.48 0.64 0.60]
julia> ra_rotated = rotate(ta, R)
```
"""
function rotate_with_matrix(ta::TrjArray{T, U}, R::AbstractMatrix{T})::TrjArray{T, U} where {T, U}
    x = R[1, 1].*ta.xyz[:, 1:3:end] .+ R[1, 2].*ta.xyz[:, 2:3:end] .+ R[1, 3].*ta.xyz[:, 3:3:end]
    y = R[2, 1].*ta.xyz[:, 1:3:end] .+ R[2, 2].*ta.xyz[:, 2:3:end] .+ R[2, 3].*ta.xyz[:, 3:3:end]
    z = R[3, 1].*ta.xyz[:, 1:3:end] .+ R[3, 2].*ta.xyz[:, 2:3:end] .+ R[3, 3].*ta.xyz[:, 3:3:end]
    xyz = similar(ta.xyz) # fixme, todo reduce memory
    xyz[:, 1:3:end] .= x
    xyz[:, 2:3:end] .= y
    xyz[:, 3:3:end] .= z
    return TrjArray(ta, xyz=xyz)
end

"""
    superimpose(ref::TrjArray, ta::TrjArray; isweight=true, index::Union{UnitRange,Vector,BitArray}, isdecenter=false) -> ta_superimposed::TrjArray

Performs the least-squares fitting of the given trajectory (TrjArray object `ta`)
to the reference structure (1st frame of TrjArray object `ref`).
If `isweight` is `true` (default), coordinates are weighted by `ta.mass` (as long as `ta.mass` is not empty).
Uses can specify column indices for the COM calculation by giving `index` which should be a Vector of integers, a UnitRange, or a BiArray.
When `isdecenter` is `true`, the function assumes the COMs of both structures are located at the origin (default is `isdecenter=false`).

Returns the superimposed trajecotry.

The code of this function is licensed under the BSD license (Copyright (c) 2009-2016 Pu Liu and Douglas L. Theobald), see LICENSE.md

# Example
```julia-repl
julia> ref = mdload("ak.pdb")
julia> ta = mdload("ak.dcd", top=pdb)
julia> ta_superimposed = superimpose(ref, ta)
```
# References
```
The algorithm of this function is based on
P. Liu, D. K. Agrafiotis, and D. L. Theobald,
J. Comput. Chem. 31, 1561-1563 (2010).
```
"""
function superimpose(ref::TrjArray{T, U}, ta::TrjArray{T, U};
    isweight::Bool=true, index::AbstractVector=Vector{U}(undef, 0), isdecenter::Bool=false)::TrjArray{T, U} where {T, U}
    nframe = ta.nframe
    natom = max(ref.natom, ta.natom)
    natom3_ta = ta.natom*3

    index2 = preprocess_index(natom, index)
    if min(ref.natom, ta.natom) < maximum(index2)
        @printf "Warning: some atom indices are larger than # of atoms, truncated.\n"
        id = index2 .<= min(ref.natom, ta.natom)
        index2 = index2[id]
    end

    weight = preprocess_weight(min(ref.natom, ta.natom), isweight, ta.mass)
    wsum_inv = one(T) / sum(weight)

    ref2 = deepcopy(ref[1, :])
    ta2 = deepcopy(ta)

    #println("step4")
    if !isdecenter
        com = decenter!(ref2, isweight=isweight, index=index2)
        decenter!(ta2, isweight=isweight, index=index2)
    end

    xyz = Matrix{T}(undef, nframe, natom3_ta)
    x = view(xyz, :, 1:3:natom3_ta)
    y = view(xyz, :, 2:3:natom3_ta)
    z = view(xyz, :, 3:3:natom3_ta)

    #println("step5")
    Threads.@threads for iframe in 1:nframe
        A, E0 = innerproduct(iframe, ref2.xyz, ta2.xyz, index2, isweight, weight)
        rmsd, rot = fastCalcRMSDAndRotation(A, E0, wsum_inv)
        applyrotation!(iframe, x, y, z, ta2.xyz, rot)
        #applyrotation!(iframe, x, y, z, x_ta, y_ta, z_ta, rot)
    end

    #println("step6")

    if !isdecenter
        x .= x .+ com.xyz[:, 1]
        y .= y .+ com.xyz[:, 2]
        z .= z .+ com.xyz[:, 3]
        #@time xyz .+= com.xyz
    end

    #@time xyz = Matrix{T}(undef, nframe, natom*3)
    #@time xyz[:, 1:3:end] .= x
    #@time xyz[:, 2:3:end] .= y
    #@time xyz[:, 3:3:end] .= z

    return TrjArray(ta, xyz=xyz)
end

function superimpose!(ref::TrjArray{T, U}, ta::TrjArray{T, U};
    isweight::Bool=true, index::AbstractVector=Vector{U}(undef, 0), isdecenter::Bool=false)::Vector{T} where {T, U}
    nframe = ta.nframe
    natom = max(ref.natom, ta.natom)

    index2 = preprocess_index(natom, index)
    if min(ref.natom, ta.natom) < maximum(index2)
        @printf "Warning: some atom indices are larger than # of atoms, truncated.\n"
        id = index2 .<= min(ref.natom, ta.natom)
        index2 = index2[id]
    end

    weight = preprocess_weight(min(ref.natom, ta.natom), isweight, ta.mass)
    wsum_inv = one(T) / sum(weight)

    if !isdecenter
        com = decenter!(ref, isweight=isweight, index=index2)
        decenter!(ta, isweight=isweight, index=index2)
    end

    rmsd = zeros(T, nframe)
    Threads.@threads for iframe in 1:nframe
        A, E0 = innerproduct(iframe, ref.xyz, ta.xyz, index2, isweight, weight)
        rmsd[iframe], rot = fastCalcRMSDAndRotation(A, E0, wsum_inv)
        applyrotation!(iframe, ta.xyz, rot)
    end

    if !isdecenter
        ta.xyz[:, 1:3:end] .= ta.xyz[:, 1:3:end] .+ com.xyz[:, 1]
        ta.xyz[:, 2:3:end] .= ta.xyz[:, 2:3:end] .+ com.xyz[:, 2]
        ta.xyz[:, 3:3:end] .= ta.xyz[:, 3:3:end] .+ com.xyz[:, 3]
    end

    return rmsd
end

function superimpose_serial(ref::TrjArray{T, U}, ta::TrjArray{T, U};
    isweight::Bool=true, index::AbstractVector=Vector{U}(undef, 0), isdecenter::Bool=false)::TrjArray{T, U} where {T, U}
    nframe = ta.nframe
    natom = ta.natom

    if isempty(index)
        index2 = collect(1:natom)
    else
        index2 = index
    end

    if isweight && length(ref.mass) == natom && length(ta.mass) == natom
        isweight2 = true
        weight = ta.mass
    else
        isweight2 = false
    end

    if isweight2
        weight2 = reshape(weight[index2], length(index2))
        wsum_inv = one(T) / sum(weight2)
    else
        wsum_inv = one(T) / T(length(index2))
    end

    if isdecenter
        ta2 = copy(ta)
        ref2 = ref[1, :]
    else
        ta2, = decenter(ta, isweight=isweight2, index=index)
        ref2, com = decenter(ref[1, :], isweight=isweight2, index=index)
    end

    x = Matrix{T}(undef, nframe, natom)
    y = Matrix{T}(undef, nframe, natom)
    z = Matrix{T}(undef, nframe, natom)
    for iframe in 1:nframe
        A, E0 = innerproduct(iframe, ref2, ta2, index2, isweight2)
        rmsd, rot = fastCalcRMSDAndRotation(A, E0, wsum_inv)
        applyrotation!(iframe, x, y, z, ta2.xyz, rot)
    end

    if !isdecenter
        x = x .+ com.x
        y = y .+ com.y
        z = z .+ com.z
    end

    xyz = Matrix{T}(undef, nframe, natom*3)
    xyz[:, 1:3:end] .= x
    xyz[:, 2:3:end] .= y
    xyz[:, 3:3:end] .= z

    return TrjArray(ta, xyz=xyz)
end

"""
    compute_rmsd(ref::TrjArray, ta::TrjArray; isweight=true, index=index::Union{UnitRange,Vector,BitArray}) -> rmsd

Calculates the root mean square deviations (RMSD) of the given trajectory (TrjArray object `ta`)
from the reference structure (1st frame of TrjArray object `ref`).
If `isweight` is `true` (default), coordinates are weighted by `ta.mass` (as long as `ta.mass` is not empty).
Uses can specify column indices for the RMSD calculation by giving `index` which should be a Vector of integers, a UnitRange, or a BiArray.

Returns RMSD values.

# Example
```julia-repl
julia> ref = mdload("ak.pdb")
julia> ta = mdload("ak.nc", top=pdb)
julia> ta_superimposed = superimpose(ref, ta)
julia> rmsd = compute_rmsd(ref, ta_superimposed)
```
"""
function compute_rmsd(ref::TrjArray{T, U}, ta::TrjArray{T, U};
    isweight::Bool=true, index::AbstractVector=Vector{U}(undef, 0))::Vector{T} where {T, U}
    nframe = ta.nframe
    natom = ta.natom

    index2 = preprocess_index(natom, index)
    index3 = to3(index2)
    weight = preprocess_weight(natom, isweight, ref.mass)
    wsum_inv = one(T) / sum(weight[index2])
    weight = reshape(weight[index2], 1, length(index2))

    ref_x = view(ref.xyz, 1:1, index3[1:3:end])
    ref_y = view(ref.xyz, 1:1, index3[2:3:end])
    ref_z = view(ref.xyz, 1:1, index3[3:3:end])
    ta_x = view(ta.xyz, :, index3[1:3:end])
    ta_y = view(ta.xyz, :, index3[2:3:end])
    ta_z = view(ta.xyz, :, index3[3:3:end])
    d =  sum(weight .* ((ta_x .- ref_x).^2 .+ (ta_y .- ref_y).^2 .+ (ta_z .- ref_z).^2), dims=2)
    d = d .* wsum_inv
    d = sqrt.(d)
    reshape(d, nframe)
end

"""
    meanstructure(ta::TrjArray; isweight=true, index=index::Union{UnitRange,Vector,BitArray}) -> ta_mean::TrjArray, ta_superimposed::TrjArray

Calculates the mean structure of the given trajectory (TrjArray object `ta`)
by iteratively fitting the trajecotry to tentative mean structures until the mean structure converges.
If `isweight` is `true` (default), coordinates are weighted by `ta.mass` (as long as `ta.mass` is not empty).
Uses can specify column indices for the fitting (superimpose function) by giving `index` which should be a Vector of integers, a UnitRange, or a BiArray.

Returns the mean structure as TrjArray object, and the supserimposed trajectory to the mean structure.

# Example
```julia-repl
julia> ta = mdload("ak.nc", top=pdb)
julia> ta_mean, ta_superimposed = meanstructure(ta)
```
"""
function meanstructure(ta::TrjArray{T, U};
    isweight::Bool=true, index::Vector{U}=Vector{U}(undef, 0))::Tuple{TrjArray{T, U},TrjArray{T, U}} where {T, U}
    nframe = ta.nframe
    natom = ta.natom
    ta2 = ta
    ref = ta2[1, :]

    r = [one(T)]
    tolerance = tolerance = T(10^(-6))
    while r[1] > tolerance
        ref_old = ref;
        ta2 = superimpose(ref, ta2, isweight=isweight, index=index)
        ref = TrjArray{T, U}(xyz=mean(ta2.xyz, dims=1)) # TODO: mean(ta) should be available in the futre
        r .= compute_rmsd(ref_old, ref, isweight=isweight, index=index)
        println("rmsd from the previous mean structure: ", r[1])
    end

    ta2 = superimpose(ref, ta2, isweight=isweight, index=index)
    ref, ta2
end

"""
    compute_rmsf(ta::TrjArray{T, U}; isweight::Bool=true) -> rmsf

rmsf (root mean square fluctuation)

Calculates the root mean square fluctuations (RMSFs) of the atoms in the given trajectory (TrjArray object `ta`).
from the average coordinates of the trajectory.
If `isweight` is `true` (default), coordinates are weighted by `ta.mass` (as long as `ta.mass` is not empty).

Returns RMSF values.

# Example
```julia-repl
julia> ta = mdload("ak.nc", top=pdb)
julia> ta_mean, ta_superimposed = meanstructure(ta)
julia> rmsf = compute_rmsf(ta)
```
"""
function compute_rmsf(ta::TrjArray{T, U}; isweight::Bool=true)::Vector{T} where {T, U}
    nframe = ta.nframe
    natom = ta.natom
    if isweight && length(ta.mass) == natom
        weight = ta.mass
    else
        weight = ones(T, natom)
    end
    wsum_inv = 1.0 / sum(weight)
    weight_row = reshape(weight, 1, length(weight))

    ref_x = sum(ta.xyz[1:1, 1:3:end] .* weight_row * wsum_inv, dims=1)
    ref_y = sum(ta.xyz[1:1, 2:3:end] .* weight_row * wsum_inv, dims=1)
    ref_z = sum(ta.xyz[1:1, 3:3:end] .* weight_row * wsum_inv, dims=1)
    ta_x = ta.xyz[:, 1:3:end]
    ta_y = ta.xyz[:, 2:3:end]
    ta_z = ta.xyz[:, 3:3:end]
    d =  sum(weight_row .* ((ta_x .- ref_x).^2 .+ (ta_y .- ref_y).^2 .+ (ta_z .- ref_z).^2), dims=2)
    d = d .* wsum_inv
    d = sqrt.(d)
    rmsf = mean(d, dims=1)
    reshape(rmsf, natom)
end

"""
    compute_distance(ta::TrjArray, index::Matrix)

Calculates distances between two atom pairs specified by the Matrix object `index`.
Each row vector in `index` contains two column indices for calculating distance.
Default of `index` is `[1 2]`.

Returns distances specified pairs in `index`.

# Example
```julia-repl
julia> ta = mdload("ak.dcd")
julia> d = compute_distance(ta, [1 2; 1; 3])
```
"""
function compute_distance(ta::TrjArray{T, U}, index=[1 2]::Matrix{U})::Matrix{T} where {T, U}
    xyz1 = view(ta.xyz, :, to3(index[:, 1]))
    xyz2 = view(ta.xyz, :, to3(index[:, 2]))
    dx = xyz1[:, 1:3:end] .- xyz2[:, 1:3:end]
    dy = xyz1[:, 2:3:end] .- xyz2[:, 2:3:end]
    dz = xyz1[:, 3:3:end] .- xyz2[:, 3:3:end]
    if !isempty(ta.boxsize)
        wrap_dx!(dx, ta.boxsize[:, 1])
        wrap_dx!(dy, ta.boxsize[:, 2])
        wrap_dx!(dz, ta.boxsize[:, 3])
    end
    d = sqrt.(dx.^2 .+ dy.^2 .+ dz.^2)
    return d
end

"""
    compute_distance(ta1::TrjArray, ta2::TrjArray, index::Matrix)

Calculates distances between two atom pairs in `ta1` and `ta2` specified by the Matrix object `index`.
Each row vector in `index` contains a column index of `ta1` and column index of `ta2` for calculating distance.
Default of `index` is `[1 1]`.

Returns distances specified pairs in `index`.

# Example
```julia-repl
julia> ta = mdload("ak.pdb")
julia> ta = mdload("ak.dcd", top=ta)
julia> d = compute_distance(ta[:, "atomid 1"], ta[:, "atomid 9"])
```
"""
function compute_distance(ta1::TrjArray{T, U}, ta2::TrjArray{T, U}, index=[1 1]::Matrix{U})::Matrix{T} where {T, U}
    xyz1 = view(ta1.xyz, :, to3(index[:, 1]))
    xyz2 = view(ta2.xyz, :, to3(index[:, 2]))
    dx = xyz1[:, 1:3:end] .- xyz2[:, 1:3:end]
    dy = xyz1[:, 2:3:end] .- xyz2[:, 2:3:end]
    dz = xyz1[:, 3:3:end] .- xyz2[:, 3:3:end]
    if !isempty(ta1.boxsize)
        wrap_dx!(dx, ta1.boxsize[:, 1])
        wrap_dx!(dy, ta1.boxsize[:, 2])
        wrap_dx!(dz, ta1.boxsize[:, 3])
    end
    d = sqrt.(dx.^2 .+ dy.^2 .+ dz.^2)
    return d
end

"""
    compute_distancemap(ta::TrjArray; kneighbor=3)

Calculates distance-map vectors from the given TrjArray `ta`.
If `atomid=j` and `atomid=i` are closer to each other than the specified `kneighbor`, i.e., `abs(i-j) < kneighbor`,
the calculation of the distance is ignored, considering that they are not nonbonded pairs.
The default value of `kneighbor` is `3`.

Returns distance-map vectors of all frames. The distance-map of each frame is contained in a row vector of the output.

# Example
```julia-repl
julia> ta = mdload("ak.pdb")
julia> ta = mdload("ak.dcd", top=ta)
julia> d = compute_distancemap(ta["atomname CA"])
```
"""
function compute_distancemap(ta::TrjArray{T, U}; kneighbor=3)::Matrix{T} where {T, U}
    natom = ta.natom
    nframe = ta.nframe
    count = 0
    for i = 1:natom
        for j = (i+kneighbor):natom
            count += 1
        end
    end
    index_all = zeros(U, count, 2)
    count = 0
    for i = 1:natom
        for j = (i+kneighbor):natom
            count += 1
            index_all[count, 1] = i
            index_all[count, 2] = j
        end
    end
    d = compute_distance(ta, index_all)
    return d
end

"""
    compute_contactmap(ta::TrjArray; rcut=8.5, kneighbor=3)

Calculates contact-map vectors from the given TrjArray `ta`.
Cut-off value for contact can be given in `rcut` (default is `rcut=8.5`).
If `atomid=j` and `atomid=i` are closer to each other than the specified `kneighbor`, i.e., `abs(i-j) < kneighbor`,
the calculation of the distance is ignored, considering that they are not nonbonded pairs.
The default value of `kneighbor` is `3`.

Returns contact-map vectors of all frames. The contact-map of each frame is contained in a row vector of the output.

# Example
```julia-repl
julia> ta = mdload("ak.pdb")
julia> ta = mdload("ak.dcd", top=ta)
julia> d = compute_distancemap(ta["atomname CA"])
```
"""
function compute_contactmap(ta::TrjArray{T, U}; rcut=8.5, kneighbor=3)::Matrix{T} where {T, U}
    d = compute_distancemap(ta, kneighbor=kneighbor)
    c = T.(d .< rcut)
    return c
end

"""
    compute_angle(ta1::TrjArray, ta2::TrjArray, ta3::TrjArray)

Calculates angles for the atom coordinates or the COMs of groups from the triplet of `ta1`, `ta2` and `ta3`.

Returns angles.

# Example
```julia-repl
julia> ta = mdload("ak.pdb")
julia> ta = mdload("ak.dcd", top=ta)
julia> a = compute_distance(ta[:, "atomid 1"], ta[:, "atomid 9"], ta[:, "atomid 11"])
```
"""
function compute_angle(ta1::TrjArray{T, U}, ta2::TrjArray{T, U}, ta3::TrjArray{T, U})::Vector{T} where {T, U}
    nframe = ta1.nframe
    com1 = centerofmass(ta1, isweight=true)
    com2 = centerofmass(ta2, isweight=true)
    com3 = centerofmass(ta3, isweight=true)
    a = zeros(T, nframe)
    for iframe in 1:nframe
        d1 = [com1.xyz[iframe, 1] - com2.xyz[iframe, 1]; com1.xyz[iframe, 2] - com2.xyz[iframe, 2]; com1.xyz[iframe, 3] - com2.xyz[iframe, 3]]
        d2 = [com3.xyz[iframe, 1] - com2.xyz[iframe, 1]; com3.xyz[iframe, 2] - com2.xyz[iframe, 2]; com3.xyz[iframe, 3] - com2.xyz[iframe, 3]]
        a[iframe] = acos(dot(d1, d2)/(norm(d1)*norm(d2)))
    end
    a = (a ./ pi) .* T(180)
end

"""
    compute_dihedral(ta1::TrjArray, ta2::TrjArray, ta3::TrjArray, ta4::TrjArray)

Calculates dihedral angles for the atom coordinates or the COMs of groups from the quadruplet of `ta1`, `ta2`, `ta3` and `ta4`.

Returns dihedral angles.

# Example
```julia-repl
julia> ta = mdload("ak.pdb")
julia> ta = mdload("ak.dcd", top=ta)
julia> d = compute_dihedral(ta[:, "atomid 1"], ta[:, "atomid 9"], ta[:, "atomid 11", ta[:, "atom 13"]])
```
"""
function compute_dihedral(ta1::TrjArray{T, U}, ta2::TrjArray{T, U}, ta3::TrjArray{T, U}, ta4::TrjArray{T, U})::Vector{T} where {T, U}
    nframe = ta1.nframe
    com1 = centerofmass(ta1, isweight=true)
    com2 = centerofmass(ta2, isweight=true)
    com3 = centerofmass(ta3, isweight=true)
    com4 = centerofmass(ta4, isweight=true)
    a = zeros(T, nframe)
    # Threads.@threads for iframe in 1:nframe
    for iframe in 1:nframe
        d1 = [com1.xyz[iframe, 1] - com2.xyz[iframe, 1]; com1.xyz[iframe, 2] - com2.xyz[iframe, 2]; com1.xyz[iframe, 3] - com2.xyz[iframe, 3]]
        d2 = [com3.xyz[iframe, 1] - com2.xyz[iframe, 1]; com3.xyz[iframe, 2] - com2.xyz[iframe, 2]; com3.xyz[iframe, 3] - com2.xyz[iframe, 3]]
        d3 = [com3.xyz[iframe, 1] - com4.xyz[iframe, 1]; com3.xyz[iframe, 2] - com4.xyz[iframe, 2]; com3.xyz[iframe, 3] - com4.xyz[iframe, 3]]
        m1 = cross(d1, d2)
        m2 = cross(d2, d3)
        a[iframe] = acos(dot(m1, m2)/(norm(m1)*norm(m2)))
        rotdirection = dot(d2,cross(m1,m2))
        if rotdirection < zero(T)
            a[iframe] = -a[iframe]
        end
    end
    a .= (a ./ pi) .* T(180)
end

"""
    compute_dihedral(ta1::TrjArray, ta2::TrjArray, ta3::TrjArray, ta4::TrjArray)

Calculates dihedral angles for the atom coordinates or the COMs of groups from the quadruplet of `ta1`, `ta2`, `ta3` and `ta4`.

Returns dihedral angles.

# Example
```julia-repl
julia> ta = mdload("ak.pdb")
julia> ta = mdload("ak.dcd", top=ta)
julia> d = compute_dihedral(ta[:, "atomid 1"], ta[:, "atomid 9"], ta[:, "atomid 11", ta[:, "atom 13"]])
```
"""
function compute_dihedral(ta::TrjArray{T, U}, array_index) where {T, U}
    nframe = ta.nframe
    ntorsion = length(array_index)
    natom3 = ta.natom*3
    a = zeros(T, nframe, ntorsion)
    x = view(ta.xyz, :, 1:3:natom3)
    y = view(ta.xyz, :, 2:3:natom3)
    z = view(ta.xyz, :, 3:3:natom3)
    for itorsion = 1:ntorsion
        id = array_index[itorsion]
        for iframe in 1:nframe
            d1 = [x[iframe, id[1]] - x[iframe, id[2]]; y[iframe, id[1]] - y[iframe, id[2]]; z[iframe, id[1]] - z[iframe, id[2]]]
            d2 = [x[iframe, id[3]] - x[iframe, id[2]]; y[iframe, id[3]] - y[iframe, id[2]]; z[iframe, id[3]] - z[iframe, id[2]]]
            d3 = [x[iframe, id[3]] - x[iframe, id[4]]; y[iframe, id[3]] - y[iframe, id[4]]; z[iframe, id[3]] - z[iframe, id[4]]]
            m1 = cross(d1, d2)
            m2 = cross(d2, d3)
            a[iframe, itorsion] = acos(dot(m1, m2)/(norm(m1)*norm(m2)))
            rotdirection = dot(d2,cross(m1,m2))
            if rotdirection < zero(T)
                a[iframe, itorsion] = -a[iframe, itorsion]
            end
        end
    end
    a .= (a ./ pi) .* T(180)
end

function compute_phi(ta::TrjArray{T, U}) where {T, U}
    return 0
end

"""
    compute_qscore(native::TrjArray, ta::TrjArray) -> qscore::Array

Calculates all-atom based Q-score from given heavy atom coordinates.

Returns Q-scores.

# Example
```julia-repl
julia> native = mdload("ak.pdb")
julia> ta = mdload("ak.dcd", top=native)
julia> qscore = compute_qscore(native["not hydrogen"], ta["not hydrogen"])
```

# References
```
Definition of Q-score from heavy atoms is given in
R. B. Best, G. Hummer, and W. A. Eaton, PNAS 110, 17874 (2013).
```
"""
function compute_qscore(native::TrjArray{T, U}, ta::TrjArray{T, U})::Vector{T} where {T, U}
    natom = native.natom
    index_all = zeros(U, U(natom*(natom-1)/2), 2)
    count = 1
    for i = 1:natom
        for j = (i+1):natom
            index_all[count, 1] = i
            index_all[count, 2] = j
            count += 1
        end
    end
    id = abs.(native.resid[index_all[:, 1]] .- native.resid[index_all[:, 2]]) .> 3
    index_all = index_all[id, :]
    d = compute_distance(native[1, :], index_all)
    id = (d .< 4.5)[:]
    index_contact = index_all[id, :]
    d_contact = d[:, id]

    d = compute_distance(ta, index_contact)
    BETA = 5.0
    LAMBDA = 1.8
    # qscore = mean(1.0./(1.0 + exp(BETA * bsxfun(@minus, dist, LAMBDA*dist0))), 2);
    qscore = mean(1.0 ./ (1.0 .+ exp.(BETA .* (d .- (LAMBDA .* d_contact)))), dims=2)[:]
    return qscore
end

"""
    compute_drms(native1::TrjArray, native2::TrjArray, ta1::TrjArray, ta2::TrjArray) -> drms

Calculates distance-based root mean square displacements (DRMS) from the given trajectories.
First, native contacts between molecule1 and molecule2 are identified from the two native structures `native1` and `native2` where cutoff distance of 6.0 Angstrom is used.
Then, the DRMS are calculated from the trajectories of molecule1 `ta1` and molecule2 `ta2`.

Returns distance-based root mean square displacements.

# Example
```julia-repl
julia> native1 = mdload("protein.pdb")["not hydrogen"]
julia> native2 = mdload("ligand.pdb")["not hydrogen"]
julia> ta1 = mdload("protein.dcd")["not hydrogen"]
julia> ta2 = mdload("ligand.dcd")["not hydrogen"]
julia> drms = compute_drms(native1, native2, ta1, ta2)
```

# References
```
Definition of DRMS is given in
J. Domański, G. Hedger, R. B. Best, P. J. Stansfeld, and M. S. P. Sansom,
Convergence and Sampling in Determining Free Energy Landscapes for Membrane Protein Association,
J. Phys. Chem. B 121, 3364 (2017).
```
"""
function compute_drms(native1::TrjArray{T, U}, native2::TrjArray{T, U}, ta1::TrjArray{T, U}, ta2::TrjArray{T, U})::Vector{T} where {T, U}
    natom1 = native1.natom
    natom2 = native2.natom
    index_all = zeros(U, natom1*natom2, 2)
    count = 1
    for i = 1:natom1
        for j = 1:natom2
            index_all[count, 1] = i
            index_all[count, 2] = j
            count += 1
        end
    end
    d = compute_distance(native1[1, :], native2[1, :], index_all)
    id = (1.0 .< d)[:] .& (d .< 6.0)[:]
    index_contact = index_all[id, :]
    d_contact = d[:, id]

    d = compute_distance(ta1, ta2, index_contact)
    ncontact = T(size(index_contact, 1))
    drms = sum((d_contact .- d).^2, dims=2)[:]
    drms .= drms ./ ncontact
    return drms
end

function wrap_x!(x, boxsize)
    x .= x .- floor.(x./boxsize).*boxsize
end

function wrap_dx!(dx, boxsize)
    dx .= dx .- round.(dx./boxsize).*boxsize
end

minimum_image(n::Int, N::Int) = n - N * sign(n) * floor(typeof(n), (Float64(abs(n)) + 0.5*Float64(N)) / Float64(N))

function pairwise_distance(x::AbstractVector, y::AbstractVector, z::AbstractVector, index::AbstractVector, rcut2::Real)
    n = length(index)
    npair = Int(n*(n-1)/2)
    pair = zeros(eltype(index), npair, 2)
    dist = zeros(eltype(x), npair)
    count = 0
    @inbounds for i = 1:n
        ii = index[i]
        for j = (i + 1):n
            jj = index[j]
            r2 = (x[ii] - x[jj])^2 + (y[ii] - y[jj])^2 + (z[ii] - z[jj])^2
            count += 1
            pair[count, 1] = ii
            pair[count, 2] = jj
            dist[count] = r2
        end
    end
    logical_index = dist .< rcut2
    pair[logical_index, :], sqrt.(dist[logical_index])
end

function pairwise_distance(x::AbstractVector, y::AbstractVector, z::AbstractVector, index::AbstractVector, rcut2::Real, boxsize::AbstractVector)
    n = length(index)
    npair = Int(n*(n-1)/2)
    pair = zeros(eltype(index), npair, 2)
    dist = zeros(eltype(x), npair)
    count = 0
    for i = 1:n
        ii = index[i]
        for j = (i + 1):n
            jj = index[j]
            dx = x[ii] - x[jj]
            dy = y[ii] - y[jj]
            dz = z[ii] - z[jj]
            dx = dx - round(dx/boxsize[1])*boxsize[1]
            dy = dy - round(dy/boxsize[2])*boxsize[2]
            dz = dz - round(dz/boxsize[3])*boxsize[3]
            r2 = dx^2 + dy^2 + dz^2
            count += 1
            pair[count, 1] = ii
            pair[count, 2] = jj
            dist[count] = r2
        end
    end
    logical_index = dist .< rcut2
    pair[logical_index, :], sqrt.(dist[logical_index])
end

function pairwise_distance(x::AbstractVector, y::AbstractVector, z::AbstractVector, index1::AbstractVector, index2::AbstractVector, rcut2::Real)
    n1 = length(index1)
    n2 = length(index2)
    npair = n1*n2
    pair = zeros(eltype(index1), npair, 2)
    dist = zeros(eltype(x), npair)
    count = 0
    @inbounds for i = 1:n1
        ii = index1[i]
        for j = 1:n2
            jj = index2[j]
            r2 = (x[ii] - x[jj])^2 + (y[ii] - y[jj])^2 + (z[ii] - z[jj])^2
            count += 1
            pair[count, 1] = ii
            pair[count, 2] = jj
            dist[count] = r2
        end
    end
    logical_index = dist .< rcut2
    pair[logical_index, :], sqrt.(dist[logical_index])
end

function pairwise_distance(x::AbstractVector, y::AbstractVector, z::AbstractVector, index1::AbstractVector, index2::AbstractVector, rcut2::Real, boxsize::AbstractVector)
    n1 = length(index1)
    n2 = length(index2)
    npair = n1*n2
    pair = zeros(eltype(index1), npair, 2)
    dist = zeros(eltype(x), npair)
    count = 0
    @inbounds for i = 1:n1
        ii = index1[i]
        for j = 1:n2
            jj = index2[j]
            dx = x[ii] - x[jj]
            dy = y[ii] - y[jj]
            dz = z[ii] - z[jj]
            dx = dx - round(dx/boxsize[1])*boxsize[1]
            dy = dy - round(dy/boxsize[2])*boxsize[2]
            dz = dz - round(dz/boxsize[3])*boxsize[3]
            r2 = dx^2 + dy^2 + dz^2
            count += 1
            pair[count, 1] = ii
            pair[count, 2] = jj
            dist[count] = r2
        end
    end
    logical_index = dist .< rcut2
    pair[logical_index, :], sqrt.(dist[logical_index])
end

"""
    compute_pairlist(ta::TrjArray, rcut; iframe=1::Int) -> F

Make a pairlist for the given strcuture `ta1` by searching pairs within a cutoff distance `rcut`.
By default, the 1st frame of `ta1` is used. Users can specify the frame by `iframe`.

Returns a NamedTuple object `F` which contains the pair lists in `F.pair`,
and the distances of corresponding pairs in `F.dist`.

# Example
```julia-repl
julia> ta = mdload("ak.pdb")
julia> F = compute_pairlist(ta, 8.0)
```

# References
```
The algorithm of this function is based on
T.N. Heinz, and P.H. Hünenberger,
A fast pairlist-construction algorithm for molecular simulations under periodic boundary conditions.
J Comput Chem 25, 1474–1486 (2004).
```
"""
function compute_pairlist(ta::TrjArray{T, U}, rcut::T; iframe=1::Int) where {T, U}
    # TODO: PBC support
    natom = ta.natom
    rcut2 = rcut^2
    is_pbc = !isempty(ta.boxsize)
    x = ta.xyz[iframe, 1:3:end]
    y = ta.xyz[iframe, 2:3:end]
    z = ta.xyz[iframe, 3:3:end]
    boxsize = zeros(T, 3)

    if is_pbc
        boxsize .= ta.boxsize[iframe, :]
        wrap_x!(x, boxsize[1])
        wrap_x!(y, boxsize[2])
        wrap_x!(z, boxsize[3])
    else
        x .= x .- minimum(x)
        y .= y .- minimum(y)
        z .= z .- minimum(z)
        boxsize[1] = maximum(x) + one(T)
        boxsize[2] = maximum(y) + one(T)
        boxsize[3] = maximum(z) + one(T)
    end

    if any(boxsize .< (2.0*rcut))
        println("working with brute-force search...")
        return compute_pairlist_bruteforce(ta, rcut, iframe=iframe)
    end

    rcell = map(x -> max(x, rcut), T.([8.0, 8.0, 8.0]))
    ncell = floor.(U, boxsize./T.(rcell))
    rcell .= boxsize./T.(ncell)

    nx = floor.(U, x./rcell[1])
    ny = floor.(U, y./rcell[2])
    nz = floor.(U, z./rcell[3])

    m = (ncell[1]*ncell[2]).*nz .+ ncell[1].*ny .+ nx .+ one(U)
    M = prod(ncell)

    # calculate cell mask
    dm = collect(U, 1:(M-1))
    dmx = mod.(dm, ncell[1])
    dmy = floor.(U, mod.(dm, ncell[1]*ncell[2]) ./ ncell[1])
    dmz = floor.(U, dm ./ (ncell[1]*ncell[2]))
    dnx = abs.(minimum_image.(dmx, ncell[1]))
    dny = similar(dmy)
    logical_index = iszero.(dmx) .| (isequal.(dmy, ncell[2] .- one(U)) .& isequal.(dmz, ncell[3] .- one(U)))
    dny[logical_index] .= abs.(minimum_image.(dmy[logical_index], ncell[2]))
    dny[.~logical_index] .= min.(abs.(minimum_image.(dmy[.~logical_index], ncell[2])), abs.(minimum_image.(dmy[.~logical_index] .+ 1, ncell[2])))
    dnz = similar(dmz)
    logical_index = isequal.(dmz, ncell[3] .- one(U)) .| (iszero.(dmx) .& iszero.(dmy))
    dnz[logical_index] .= abs.(minimum_image.(dmz[logical_index], ncell[3]))
    dnz[.~logical_index] .= min.(abs.(minimum_image.(dmz[.~logical_index], ncell[3])), abs.(minimum_image.(dmz[.~logical_index] .+ 1, ncell[3])))
    mask = T.((max.(dnx, one(U)) .- one(U)).^2) .* rcell[1].^2 .+
           T.((max.(dny, one(U)) .- one(U)).^2) .* rcell[2].^2 .+
           T.((max.(dnz, one(U)) .- one(U)).^2) .* rcell[3].^2 .< rcut2
    mask_index = findall(mask)

    # calculate cell mask pointer
    K = round(U, (natom/M) * 20) ######### TODO, 20 is really OK?
    cell_pointer = zeros(U, M)
    cell_atom = zeros(U, K, M)
    for iatom = 1:natom
        cell_pointer[m[iatom]] += 1
        cell_atom[cell_pointer[m[iatom]], m[iatom]] = iatom
    end
    if any(cell_pointer .> K)
        error("# of atoms in a cell greater than expected")
    end

    # calculate distances of atoms in masked cells
    pair = Matrix{U}(undef, natom*1000, 2)
    dist = Vector{T}(undef, natom*1000)
    count = 0
    @inbounds for m1 = 1:M
        if iszero(cell_pointer[m1])
            continue
        end
        m1_index = cell_atom[1:cell_pointer[m1], m1]
        if isempty(m1_index)
            continue
        elseif length(m1_index) > 1
            if is_pbc
                pair_local, dist_local = pairwise_distance(x, y, z, m1_index, rcut2, boxsize)
            else
                pair_local, dist_local = pairwise_distance(x, y, z, m1_index, rcut2)
            end
            num = length(dist_local)
            if (num > 0)
                pair[(count+1):(count+num), :] .= pair_local
                dist[(count+1):(count+num)] .= dist_local
                count += num
            end
        end

        for m2 = (m1 .+ mask_index)
            if m2 > M
                break
            end
            if iszero(cell_pointer[m2])
                continue
            end
            m2_index = cell_atom[1:cell_pointer[m2], m2]
            if is_pbc
                pair_local, dist_local = pairwise_distance(x, y, z, m1_index, m2_index, rcut2, boxsize)
            else
                pair_local, dist_local = pairwise_distance(x, y, z, m1_index, m2_index, rcut2)
            end
            num = length(dist_local)
            if (num > 0)
                pair[(count+1):(count+num), :] .= pair_local
                dist[(count+1):(count+num)] .= dist_local
                count += num
            end
        end
    end
    return (pair=pair[1:count, :], dist=dist[1:count])
end

"""
    compute_pairlist_bruteforce(ta::TrjArray, rcut; iframe=1::Int) -> F

Bruteforce search version of `compute_pairlist` function. Mainly used for debugging.

Returns a NamedTuple object `F` which contains the pair lists in `F.pair`,
and the distances of corresponding pairs in `F.dist`.

# Example
```julia-repl
julia> ta = mdload("ak.pdb")
julia> F = compute_pairlist(ta, 8.0)
```
"""
function compute_pairlist_bruteforce(ta::TrjArray{T, U}, rcut::T; iframe=1::Int) where {T, U}
    rcut2 = rcut^2
    natom3 = ta.natom*3
    x = view(ta.xyz, iframe, 1:3:natom3)
    y = view(ta.xyz, iframe, 2:3:natom3)
    z = view(ta.xyz, iframe, 3:3:natom3)
    index = collect(U, 1:ta.natom)
    is_pbc = !isempty(ta.boxsize)
    boxsize = zeros(T, 3)
    if is_pbc
        boxsize .= ta.boxsize[iframe, :]
        pair, dist = pairwise_distance(x, y, z, index, rcut2, boxsize)
    else
        pair, dist = pairwise_distance(x, y, z, index, rcut2)
    end
    return (pair=pair, dist=dist)
end

"""
    compute_skrew(t0::TrjArray, t1::TrjArray) -> S

Compute the skrew motion prameters by interpreting two structural placements
as a rigid body skrew motion

Returns a NamedTuple object `S` which contains
the angle of rotation about the skrew axis in `S.phi`,
the translation along the skrew axis in `S.d`,
a unit vector parallel to the screw axis in `S.omega`
and the shortest vector to the point on the screw axis (perpendicular to the screw axis) in `S.prou`.

# Example
```julia-repl
julia> t0 = mdload("monomer1.pdb")
julia> t1 = mdload("monomer2.pdb")
julia> S = compute_skrew(t1, t2)
```
"""
function compute_skrew(t0::TrjArray{T, U}, t1::TrjArray{T, U}) where {T, U}
    @assert t0.natom == t1.natom
    i = 1
    j = round(Int, t0.natom / 2)
    k = t1.natom
    P0 = t0[1, i].xyz[:]
    Q0 = t0[1, j].xyz[:]
    R0 = t0[1, k].xyz[:]
    P1 = t1[1, i].xyz[:]
    Q1 = t1[1, j].xyz[:]
    R1 = t1[1, k].xyz[:]
    nume = cross((Q1 .- Q0) .- (R1 .- R0), (P1 .- P0) .- (R1 .- R0))
    denom = dot((Q1 .- Q0) .- (R1 .- R0), (P1 .+ P0) .- (R1 .+ R0))
    omega = nume ./ denom
    omega_norm = norm(omega, 2.0)
    omega = omega ./ omega_norm
    phi = 2.0 * atan(omega_norm)
    prou = 0.5 .* (cross(omega, P1 .- P0)./tan(0.5*phi) .- dot(omega, P1 .+ P0) .* omega .+ P0 .+ P1)
    d = dot(omega, P1 .- P0)
    return (omega=omega, phi=phi, prou=prou, d=d)
end

function disp_skrew(S, d=500.0, d0=0.0)
    x1 = S.prou[1] + d0*S.omega[1]
    y1 = S.prou[2] + d0*S.omega[2]
    z1 = S.prou[3] + d0*S.omega[3]
    x2 = S.prou[1] + d0*S.omega[1] + d*S.omega[1]
    y2 = S.prou[2] + d0*S.omega[2] + d*S.omega[2]
    z2 = S.prou[3] + d0*S.omega[3] + d*S.omega[3]
    @printf "pseudoatom pt1, pos=[%f, %f, %f]\n" x1 y1 z1
    @printf "pseudoatom pt2, pos=[%f, %f, %f]\n" x2 y2 z2
    @printf "distance /pt1, /pt2\n"
    @printf "set dash_gap, 0\n"
    @printf "set dash_width, 10\n"
end


text:
using BinaryBuilder, Pkg
using Base.BinaryPlatforms: arch, os

include("../../fancy_toys.jl")

name = "TensorRT"
version = v"8.0.3"

script = raw"""
cd ${WORKSPACE}/srcdir

mkdir -p ${bindir} ${libdir} ${includedir}
cd TensorRT*
mv bin/* ${bindir}
mv include/* ${includedir}
mv lib/*.${dlext}* ${libdir}
if [[ ${target} == x86_64-w64-mingw32 ]]; then
    chmod +x ${bindir}/*.{dll,exe}
    install_license doc/TensorRT-SLA.pdf
else
    install_license doc/pdf/TensorRT-SLA.pdf
fi
"""

lib_names = [
    "nvinfer",
    "nvinfer_plugin",
    "nvonnxparser",
    "nvparsers"
]

products = vcat(
    [LibraryProduct(["lib" * lib_name, lib_name], Symbol("lib" * lib_name); dont_dlopen=true) for lib_name in lib_names],
    [ExecutableProduct("trtexec", :trtexec)]
)

dependencies = [Dependency("CUDNN_jll", v"8.2.1"; compat="8.2")]

cuda_versions = [v"10.2", v"11.0", v"11.1", v"11.2", v"11.3"]
for cuda_version in cuda_versions
    cuda_tag = "$(cuda_version.major).$(cuda_version.minor)"
    include("build_$(cuda_tag).jl")

    for (platform, sources) in platforms_and_sources
        augmented_platform = Platform(arch(platform), os(platform); cuda=cuda_tag)
        should_build_platform(triplet(augmented_platform)) || continue
        build_tarballs(ARGS, name, version, sources, script, [augmented_platform],
                       products, dependencies; lazy_artifacts=true)
    end
end

ext: jl

max_stars_repo_name: cmhamel/Yggdrasil

max_issues_repo_name: cmhamel/Yggdrasil

max_forks_repo_name: cmhamel/Yggdrasil


text:
# conversion of SizedArray from StaticArrays.jl
"""
    SizedAbstractArray{Tuple{dims...}}(array)

Wraps an `AbstractArray` with a static size, so to take advantage of the (faster)
methods defined by the static array package. The size is checked once upon
construction to determine if the number of elements (`length`) match, but the
array may be reshaped.
"""
struct SizedAbstractArray{S<:Tuple,T,N,M,TData<:AbstractArray{T,M}} <: StaticArray{S,T,N}
    data::TData

    function SizedAbstractArray{S,T,N,M,TData}(
        a::TData,
    ) where {S,T,N,M,TData<:AbstractArray}
        if length(a) != StaticArrays.tuple_prod(S)
            error("Dimensions $(size(a)) don't match static size $S")
        end
        return new{S,T,N,M,TData}(a)
    end

    function SizedAbstractArray{S,T,N,1}(::UndefInitializer) where {S,T,N}
        return new{S,T,N,1,Array{T,1}}(Array{T,1}(undef, StaticArrays.tuple_prod(S)))
    end
    function SizedAbstractArray{S,T,N,N}(::UndefInitializer) where {S,T,N}
        return new{S,T,N,N,Array{T,N}}(Array{T,N}(undef, size_to_tuple(S)...))
    end
end

@inline function SizedAbstractArray{S,T,N}(
    a::TData,
) where {S,T,N,M,TData<:AbstractArray{T,M}}
    return SizedAbstractArray{S,T,N,M,TData}(a)
end
@inline function SizedAbstractArray{S,T}(a::TData) where {S,T,M,TData<:AbstractArray{T,M}}
    return SizedAbstractArray{S,T,StaticArrays.tuple_length(S),M,TData}(a)
end
@inline function SizedAbstractArray{S}(a::TData) where {S,T,M,TData<:AbstractArray{T,M}}
    return SizedAbstractArray{S,T,StaticArrays.tuple_length(S),M,TData}(a)
end
@inline function SizedAbstractArray{S,T,N}(::UndefInitializer) where {S,T,N}
    return SizedAbstractArray{S,T,N,N}(undef)
end
@inline function SizedAbstractArray{S,T}(::UndefInitializer) where {S,T}
    return SizedAbstractArray{S,T,StaticArrays.tuple_length(S),StaticArrays.tuple_length(S)}(
        undef,
    )
end
@generated function (::Type{SizedAbstractArray{S,T,N,M,TData}})(
    x::NTuple{L,Any},
) where {S,T,N,M,TData,L}
    if L != StaticArrays.tuple_prod(S)
        error("Dimension mismatch")
    end
    exprs = [:(a[$i] = x[$i]) for i in 1:L]
    return quote
        $(Expr(:meta, :inline))
        a = SizedAbstractArray{S,T,N,M}(undef)
        @inbounds $(Expr(:block, exprs...))
        return a
    end
end
@inline function SizedAbstractArray{S,T,N}(x::Tuple) where {S,T,N}
    return SizedAbstractArray{S,T,N,N,Array{T,N}}(collect(x))
end
@inline function SizedAbstractArray{S,T}(x::Tuple) where {S,T}
    return SizedAbstractArray{S,T,StaticArrays.tuple_length(S)}(x)
end
@inline function SizedAbstractArray{S}(x::NTuple{L,T}) where {S,T,L}
    return SizedAbstractArray{S,T}(x)
end

# Overide some problematic default behaviour
@inline function Base.convert(
    ::Type{SA},
    sa::SizedAbstractArray,
) where {SA<:SizedAbstractArray}
    return SA(sa.data)
end
@inline Base.convert(::Type{SA}, sa::SA) where {SA<:SizedAbstractArray} = sa

# Back to Array (unfortunately need both convert and construct to overide other methods)
@inline function Base.Array(sa::SizedAbstractArray{S}) where {S}
    return Array(reshape(sa.data, size_to_tuple(S)))
end
@inline function Base.Array{T}(sa::SizedAbstractArray{S,T}) where {T,S}
    return Array(reshape(sa.data, size_to_tuple(S)))
end
@inline function Base.Array{T,N}(sa::SizedAbstractArray{S,T,N}) where {T,S,N}
    return Array(reshape(sa.data, size_to_tuple(S)))
end

@inline function Base.convert(::Type{Array}, sa::SizedAbstractArray{S}) where {S}
    return Array(reshape(sa.data, size_to_tuple(S)))
end
@inline function Base.convert(::Type{Array{T}}, sa::SizedAbstractArray{S,T}) where {T,S}
    return Array(reshape(sa.data, size_to_tuple(S)))
end
@inline function Base.convert(
    ::Type{Array{T,N}},
    sa::SizedAbstractArray{S,T,N},
) where {T,S,N}
    return Array(reshape(sa.data, size_to_tuple(S)))
end

Base.@propagate_inbounds Base.getindex(a::SizedAbstractArray, i::Int) = getindex(a.data, i)

Base.@propagate_inbounds function Base.setindex!(a::SizedAbstractArray, v, i::Int)
    return setindex!(a.data, v, i)
end

const SizedAbstractVector{S,T,M} = SizedAbstractArray{Tuple{S},T,1,M}

@inline function SizedAbstractVector{S}(a::TData) where {S,T,M,TData<:AbstractArray{T,M}}
    return SizedAbstractArray{Tuple{S},T,1,M,TData}(a)
end
@inline function SizedAbstractVector{S}(x::NTuple{L,T}) where {S,T,L}
    return SizedAbstractArray{Tuple{S},T,1,1,Vector{T}}(x)
end

const SizedAbstractMatrix{S1,S2,T,M} = SizedAbstractArray{Tuple{S1,S2},T,2,M}

@inline function SizedAbstractMatrix{S1,S2}(
    a::TData,
) where {S1,S2,T,M,TData<:AbstractArray{T,M}}
    return SizedAbstractArray{Tuple{S1,S2},T,2,M,TData}(a)
end
@inline function SizedAbstractMatrix{S1,S2}(x::NTuple{L,T}) where {S1,S2,T,L}
    return SizedAbstractArray{Tuple{S1,S2},T,2,2,Matrix{T}}(x)
end

Base.dataids(sa::SizedAbstractArray) = Base.dataids(sa.data)

function Base.promote_rule(
    ::Type{<:SizedAbstractArray{S,T,N,M,TDataA}},
    ::Type{<:SizedAbstractArray{S,U,N,M,TDataB}},
) where {S,T,U,N,M,TDataA,TDataB}
    TU = promote_type(T, U)
    return SizedAbstractArray{
        S,
        TU,
        N,
        M,
        promote_type(TDataA, TDataB)::Type{<:AbstractArray{TU}},
    }
end

@inline Base.copy(a::SizedAbstractArray) = typeof(a)(copy(a.data))

function Base.similar(::Type{<:SizedAbstractArray{S,T,N,M}}, ::Type{T2}) where {S,T,N,M,T2}
    return SizedAbstractArray{S,T2,N,M}(undef)
end
function Base.similar(::Type{SA}, ::Type{T}, s::Size{S}) where {SA<:SizedAbstractArray,T,S}
    return sizedabstractarray_similar_type(T, s, StaticArrays.length_val(s))(undef)
end

function sizedabstractarray_similar_type(
    ::Type{T},
    s::Size{S},
    ::Type{Val{D}},
) where {T,S,D}
    return SizedAbstractArray{Tuple{S...},T,D,length(s)}
end


text:
@testset "normalized array" begin
    grdt!(buf, _) = (buf .= 0; buf[1] = 1; buf)
    result = optimize(x->x[1], grdt!, randn(2,2), ConjugateGradient(manifold=Sphere()))
    @test result.minimizer ≈ [-1 0; 0 0]
    @test result.minimum ≈ -1
end

@testset "input types" begin
    f(X) = (10 - X[1])^2 + (0 - X[2])^2 + (0 - X[3])^2 + (5 - X[4])^2

    function g!(storage, x)
        storage[1] = -20 + 2 * x[1]
        storage[2] = 2 * x[2]
        storage[3] = 2 * x[3]
        storage[4] = -10 + 2 * x[4]
        return
    end
    @testset "vector" begin
        for m in (AcceleratedGradientDescent, ConjugateGradient, BFGS, LBFGS, NelderMead, GradientDescent, MomentumGradientDescent, NelderMead, ParticleSwarm, SimulatedAnnealing, NGMRES, OACCEL)
            debug_printing && printstyled("Solver: "*string(m); color = :green)
            res = optimize(f, g!, [1., 0., 1., 0.], m())
            @test typeof(Optim.minimizer(res)) <: Vector
            if !(m in (NelderMead, SimulatedAnnealing, ParticleSwarm))
                @test norm(Optim.minimizer(res) - [10.0, 0.0, 0.0, 5.0]) < 10e-8
            end
        end
    end

    @testset "matrix" begin
        for m in (AcceleratedGradientDescent, ConjugateGradient, BFGS, LBFGS, ConjugateGradient,  GradientDescent, MomentumGradientDescent, ParticleSwarm, SimulatedAnnealing, NGMRES, OACCEL)
            res = optimize(f, g!, Matrix{Float64}(I, 2, 2), m())
            @test typeof(Optim.minimizer(res)) <: Matrix
            if !(m in (SimulatedAnnealing, ParticleSwarm))
                @test norm(Optim.minimizer(res) - [10.0 0.0; 0.0 5.0]) < 10e-8
            end
        end
    end

    @testset "tensor" begin
        eye3 = zeros(2,2,1)
        eye3[:,:,1] = Matrix{Float64}(I, 2, 2)
        for m in (AcceleratedGradientDescent, ConjugateGradient, BFGS, LBFGS, ConjugateGradient,  GradientDescent, MomentumGradientDescent, ParticleSwarm, SimulatedAnnealing, NGMRES, OACCEL)
            res = optimize(f, g!, eye3, m())
            _minimizer = Optim.minimizer(res)
            @test typeof(_minimizer) <: Array{Float64, 3}
            @test size(_minimizer) == (2,2,1)
            if !(m in (SimulatedAnnealing, ParticleSwarm))
                @test norm(_minimizer - [10.0 0.0; 0.0 5.0]) < 10e-8
            end
        end
    end
end

using RecursiveArrayTools
@testset "arraypartition input" begin

    function polynomial(x)
            return (10.0 - x[1])^2 + (7.0 - x[2])^4 + (108.0 - x[3])^4
        end

    function polynomial_gradient!(storage, x)
            storage[1] = -2.0 * (10.0 - x[1])
            storage[2] = -4.0 * (7.0 - x[2])^3
            storage[3] = -4.0 * (108.0 - x[3])^3
        end

    function polynomial_hessian!(storage, x)
            storage[1, 1] = 2.0
            storage[1, 2] = 0.0
            storage[1, 3] = 0.0
            storage[2, 1] = 0.0
            storage[2, 2] = 12.0 * (7.0 - x[2])^2
            storage[2, 3] = 0.0
            storage[3, 1] = 0.0
            storage[3, 2] = 0.0
            storage[3, 3] = 12.0 * (108.0 - x[3])^2
        end

    ap = ArrayPartition(rand(1), rand(2))

    optimize(polynomial, polynomial_gradient!, polynomial_hessian!, ap, NelderMead())
    optimize(polynomial, polynomial_gradient!, polynomial_hessian!, ap, ParticleSwarm())
    optimize(polynomial, polynomial_gradient!, polynomial_hessian!, ap, SimulatedAnnealing())

    optimize(polynomial, polynomial_gradient!, polynomial_hessian!, ap, GradientDescent())
    optimize(polynomial, polynomial_gradient!, polynomial_hessian!, ap, AcceleratedGradientDescent())
    optimize(polynomial, polynomial_gradient!, polynomial_hessian!, ap, MomentumGradientDescent())

    optimize(polynomial, polynomial_gradient!, polynomial_hessian!, ap, ConjugateGradient())

    optimize(polynomial, polynomial_gradient!, polynomial_hessian!, ap, BFGS())
    optimize(polynomial, polynomial_gradient!, polynomial_hessian!, ap, LBFGS())

    optimize(polynomial, polynomial_gradient!, polynomial_hessian!, ap, Newton())
    optimize(polynomial, polynomial_gradient!, polynomial_hessian!, ap, NewtonTrustRegion())
end


text:
# based on https://dlmf.nist.gov/11.13.iii, works for real and complex arguments
# fails at z=0 because Struve K and Bessel Y diverge there

using SpecialFunctions
using QuadGK

integrate(f, a, b) = QuadGK.quadgk(f, a, b)[1]

# only valid for real(z) > 0, https://dlmf.nist.gov/11.5.E2
_K0_integral(z::Real) = (2 / pi) * integrate(t -> exp(-z * sinh(t)), 0, Inf)
_K_integral(ν, z::Real) = 2(0.5z)^ν / (sqrt(pi) * gamma(ν + 0.5)) *
    integrate(t -> exp(-z * t) * (1 + t^2)^(ν - 0.5), 0, Inf)
_H0_integral(z::Real) = _K0_integral(z) + bessely0(z)
_H_integral(ν, z::Real) = _K_integral(ν, z) + bessely(ν, z)

# only valid for real(ν) > -0.5, https://dlmf.nist.gov/11.5.E1
_H0_integral(z::Complex{T}) where {T <: Real} = (2 / pi) *
    integrate(ϑ -> sin(z * cos(ϑ)), 0, 0.5pi)
_H_integral(ν, z::Complex{T}) where {T <: Real} = 2(0.5z)^ν /
    (sqrt(pi) * gamma(ν + 0.5)) *
    integrate(t -> (1 - t^2)^(ν - 0.5) * sin(z * t), 0, 1)
_K0_integral(z::Complex{T}) where {T <: Real} = _H0_integral(z) - bessely0(z)
_K_integral(ν, z::Complex{T}) where {T <: Real} = _H_integral(ν, z) - bessely(ν, z)

# only valid for real(ν) > -0.5, https://dlmf.nist.gov/11.5.E4
_M0_integral(z) = -(2 / pi) * integrate(ϑ -> exp(-z * cos(ϑ)), 0, 0.5pi)
_M_integral(ν, z) = -2(0.5z)^ν / (sqrt(pi) * gamma(ν + 0.5)) *
    integrate(t -> exp(-z * t) * (1 - t^2)^(ν - 0.5), 0, 1)

_L0_integral(z)   = besseli(0, z) + _M0_integral(z)
_L_integral(ν, z) = besseli(ν, z) + _M_integral(ν, z)


text:
using Distributed
using Hwloc
if length(workers()) == 1
    addprocs(Hwloc.num_physical_cores())  # get ready for parallelism
end
using Distributions:  Normal, quantile
using MLJ: coerce!, autotype
using CategoricalArrays: CategoricalArray, levels
using DataFrames
using Missings: ismissing
using Printf

function setScientificTypes!(X::DataFrame)
    for c in 1:ncol(X)
        if unique(X[:,c])[1] === missing
            println("setScientificTypes:  Column :", names(X)[c], " has missing values.  They are not allowed.")
            return (missing)
        end
        # remove the missing type that MLJ modeling does not like
        if typeof(X[:,c]) == Array{Union{Missing,Float64},1}
            X[!,c] = Float64.(X[!,c])
        end
        if typeof(X[:,c]) == Array{Union{Missing,Int64},1}
            X[!,c] = Int64.(X[!,c])
        end
        if typeof(X[:,c]) == Array{Union{Missing,String},1}
            X[!,c] = String.(X[!,c])
        end
    end

    coerce!(X, autotype(X, :discrete_to_continuous))
    coerce!(X, autotype(X, :string_to_multiclass))
    return (X)
end # setScientificTypes!

@everywhere function getResiduals(
    E::Array{Int64},
    e::Int64,
    residuals::Array{Float64},
    type::String = "inEnvironment",
)
    if type == "inEnvironment"
        # only itself
        l = [residuals[i] for i in 1:length(E) if E[i] == e]
    else
        # all other excluding itself
        l = [residuals[i] for i in 1:length(E) if E[i] != e]
    end
    return (l)
end # getResiduals

function normalizer(v)
    fMax = partialsort(v, 1, rev = true)
    fMin = partialsort(v, 1, rev = false)
    if ismissing(fMax - fMin)
        n = Vector{Union{Missing,Float64}}(missing, length(v))
    elseif (fMax - fMin) == 0
        n = v
    else
        n = (v .- fMin) ./ (fMax - fMin)
    end
    return n
end

function string_to_float(str)
    try
        convert(Float64, str)
    catch
        return (missing)
    end
end

function getIndexFromPredictor(X::DataFrame)
    indexFromPredictor = Dict{String,Int64}()
    i = 1
    for n in names(X)
        indexFromPredictor[String(n)] = i
        i += 1
    end
    return (indexFromPredictor)
end # getIndexFromPredictor

function getPredictorFromIndex(X::DataFrame)
    predictorFromIndex = Dict{Int64,String}()
    i = 1
    for n in names(X)
        predictorFromIndex[i] = String(n)
        i += 1
    end
    return (predictorFromIndex)
end # getPredictorFromIndex

function getQuantileNormalStd(α::Float64, coefsStdError)
    return ( quantile(Normal(0, 1), 1 - α / 4) .* coefsStdError )
end # getQuantileNormalStd

function StoreAcceptedOuput!(
    r,
    acceptedTest::Array{Int64,1},
    PValuesAccepted::Dict{Int64,Array{Float64}},
    Coeff = missing,
    CoeffVar = missing,
)
    i = 1
    for v in acceptedTest
        if !ismissing(CoeffVar)
            if !haskey(Coeff, v)
                CoeffVar[v] = Float64[]
                Coeff[v] = Float64[]
            end
        end
        PValuesAccepted[v] = Float64[]
        if !ismissing(CoeffVar)
            push!(Coeff[v], r.coefs[i])
            push!(CoeffVar[v], r.coefsStdError[i])
        end
        push!(PValuesAccepted[v], r.pvalue)
        i += 1
    end
end # StoreAcceptedOuput!

function StoreNotAcceptedOuput!(
    r,
    notAcceptedTest::Array{Int64,1},
    PValuesNotAccepted::Dict{Int64,Array{Float64}}
)
    i = 1
    for v in notAcceptedTest
        if !haskey(PValuesNotAccepted, v)
            PValuesNotAccepted[v] = Float64[]
        end
        push!(PValuesNotAccepted[v], r.pvalue)
        i += 1
    end
end # StoreNotAcceptedOuput!

function doOutput(
    X,
    featureImportances,
    acceptedSets,
    PValuesAccepted,
    PValuesNotAccepted,
    Pall,
    gof,
    verbose,
    ConfInt = missing,
    Coeff = missing,
    CoeffVar = missing,
)
    if maximum(Pall) < gof
        println("Goodness Of Fit is bad.  Cut off gof = ", gof, " but min and max p-values are:")
        @printf("%9.6f\t%9.6f\n", minimum(Pall), maximum(Pall))
        return (missing)
    end

    if verbose == true
        if !ismissing(Coeff)
            println("\n\n\nCoefficients:")
            for c = 1:ncol(X)
                @printf("%s \t", names(X)[c])
                println(Coeff[c])
            end

            println("\n\nStadard Error of the point-estimates:")
            for c = 1:ncol(X)
                @printf("%s \t", names(X)[c])
                println(CoeffVar[c])
            end
        end

        if !ismissing(featureImportances)
            println("\nFeature Importances: Mean Effect Percent")
            for f in featureImportances
                @printf("%s \t %9.12f\n", f.feature_name, f.meanEffectPercent)
            end
        end

        println("\n\nP-Values NOT Accepted:")
        for c = 1:ncol(X)
            @printf("%s \t", names(X)[c])
            try
                for i in PValuesNotAccepted[c]
                    @printf("%9.12f\t", i)
                end
            catch
                print("None")
            end
            println()
        end

        println("\n\nP-Values Accepted:")
        for c = 1:ncol(X)
            @printf("%-25s \t", names(X)[c])
            try
                for i in PValuesAccepted[c]
                    @printf("%9.6f\t", i)
                end
            catch
                print("None")
            end
            println()
        end
    end # verbose

    println("\nAccepted Sets:")
    for a in acceptedSets[2:end]
        println(a)
    end

    if !ismissing(ConfInt)
        println("\n\nConfidence Intervals:")
        @printf("%-35s \t  %-9s \t  %-9s \t  %-9s \t\n", "Predictor", "Low", "High", "Difference")
        i = 1
        for n in names(X)
            high = ConfInt[1,i]
            low = ConfInt[2,i]
            @printf("%-35s \t %9.6f\t %9.6f \t %9.6f\n", n, low, high, high - low)
            i += 1
        end
    end

    println("\n\nInvariant Casual Predictors:")
    icp = Set()
    for a in acceptedSets[end]
        push!(icp, a)
    end
    for a in acceptedSets[2:end]
        intersect!(icp, Set(a))
    end
    if length(icp) == 0
        println("Invariant Causal Predictions were not found at the given confidence level. ")
        println("First, run linear ICP.  If all linear models were rejected, then run nonlinear forest ICP .")
        println("If both linear and non linear do not find any casual predictors, there might be hidden variables.  Then, try the linear hiddenICP function in R.\n")
        return (missing)
    else
        casualPredictors = Dict{Int64,String}()
        predictorFromIndex = getPredictorFromIndex(X)
        for i in sort(collect(icp))
            print(i, " = ", predictorFromIndex[i], "\t")
            casualPredictors[i] =  predictorFromIndex[i]
        end
        println("\n")
        return (casualPredictors)
    end

    GC.gc()
end # doOutput

function hotEncoder(X::DataFrame)
    setScientificTypes!(X)
    for n in names(X)
        if typeof(CategoricalArray(["a"])) == typeof(X[:, n]) || typeof(CategoricalArray([1])) == typeof(X[:, n])
            if length(levels(X[:, n])) == 1
                # one level category is useless
                println("hotEncoder: Removing column :", n, " since it has only one category value")
                select!(X, Not([n]))
            elseif length(levels(X[:, n])) == 2
                # 2 levels only do not need a new hot encoder column
                encode = 0.0
				X.newColumn = zeros(length(X[:, n]))
                for l in unique(X[!, n])
                    X[X[:, n] .== l, :newColumn] .= encode
                    encode += 1.0
                end
				select!(X, Not([n]))
                rename!(X, Dict(:newColumn => n))
            else
                for l in levels(X[:, n])[1:end]
                    encoder = string(n, "__", strip(l))
                    X.newColumn = zeros(length(X[:, n]))
                    X[X[:, n] .== l, :newColumn] .= 1.0
                    rename!(X, Dict(:newColumn => encoder))
                end
                select!(X, Not([n]))
            end
        end
    end
    setScientificTypes!(X)
    return (X)
end # hotEncoder

# use to transform Y and ExpInd from strings/categories to Integer
function categoryIntegerMapper!(df, startingNumber=1)
    df.ICP___Mapper = Array{Int64,1}(undef, nrow(df))
    n = startingNumber
    for u in unique(sort(df[:,1]))
        println( u, " --> ", n)
        for r in 1:nrow(df)
            if df[r,1] == u
                df[r,2] = n
            end
        end
        n += 1
    end
    select!(df, :ICP___Mapper)
end

function labelEncoder(X::DataFrame)
    setScientificTypes!(X)
    for n in names(X)
        if typeof(CategoricalArray(["a"])) == typeof(X[:, n]) || typeof(CategoricalArray([1])) == typeof(X[:, n])
            if length(levels(X[:, n])) == 1
                # one level category is useless
                println("labelEncoder: Removing column :", n, " since it has only one category value")
                select!(X, Not([n]))
            else
                encode = 0.0
				X.newColumn = zeros(length(X[:, n]))
                for l in unique(X[!, n])
                    X[X[:, n] .== l, :newColumn] .= encode
                    encode += 1.0
                end
				select!(X, Not([n]))
                rename!(X, Dict(:newColumn => n))
            end
        end
    end
    setScientificTypes!(X)
    return (X)
end # labelEncoder


text:
struct NireODESol{uType,tType}
   t::Array{tType,1}
   u::Array{uType,1}
   iter::Array{Float64,1}
   retcode::Symbol
   f
end

function  IRK8(u0, t0, tf, h, n, m, f, p;  initial_interp=true)
    order=16
    utype = typeof(u0)
    ttype = typeof(t0)
    itermax = 100
    sdt = sign(h)
    (hb, hc, mu, nu) = IRK8Coefficients(h)
    s = length(hb)
    U = Array{utype}(undef, s)
    Uz = Array{utype}(undef, s)
    L = Array{utype}(undef, s)
    L2 = Array{utype}(undef, s)
    Dmin=Array{utype}(undef,s)
    F = Array{utype}(undef, s)

    for i in 1:s
      Dmin[i] = zeros(eltype(u0), size(u0))
      F[i] = zeros(eltype(u0), size(u0))
    end

    ej=zeros(eltype(u0), size(u0))

    for i in 1:s
        U[i] = zeros(eltype(u0), size(u0))
        Uz[i] = zeros(eltype(u0), size(u0))
        L[i] = zeros(eltype(u0), size(u0))
        L2[i] = zeros(eltype(u0), size(u0))
    end
    uu = Array{typeof(u0)}(undef, n+1)
    tt = Array{ttype}(undef, n+1)

    iters = zeros(n+1)
    uu[1] = copy(u0)
    tt[1] = t0
    tj = [t0, zero(t0)]
    uj = copy(u0)
    j = 0
    cont = true
    while cont
        j += 1
        for i in 1:m
           it = IRK8step!(tj,uj,ej,p,f,h,hb,hc,mu,nu,U,Uz,L,F,Dmin,itermax,initial_interp)
           iters[j+1] += it
        end
        iters[j+1] /= m
        uu[j+1] = uj + ej
        tt[j+1] = tj[1] + tj[2]
        cont = (sdt*tt[j+1] < sdt*tf) && (j<n)
    end
    sol = NireODESol(tt[1:(j+1)],uu[1:(j+1)],iters[1:(j+1)],:Successs,f)
    return sol
  end

function  IRK8(u0, t0, tf, n, m, f, p;  initial_interp=true)
    order=16
    utype = typeof(u0)
    ttype = typeof(t0)
    itermax = 100
    h = (tf-t0)/(n*m)
    (hb, hc, mu, nu) = IRK8Coefficients(h)
    s = length(hb)
    U = Array{utype}(undef, s)
    Uz = Array{utype}(undef, s)
    L = Array{utype}(undef, s)
    L2 = Array{utype}(undef, s)
    Dmin=Array{utype}(undef,s)
    F = Array{utype}(undef, s)

    for i in 1:s
      Dmin[i] = zeros(eltype(u0), size(u0))
      F[i] = zeros(eltype(u0), size(u0))
    end

    ej=zeros(eltype(u0), size(u0))

    for i in 1:s
        U[i] = zeros(eltype(u0), size(u0))
        Uz[i] = zeros(eltype(u0), size(u0))
        L[i] = zeros(eltype(u0), size(u0))
        L2[i] = zeros(eltype(u0), size(u0))
    end
    uu = Array{typeof(u0)}(undef, n+1)
    tt = Array{ttype}(undef, n+1)

    iters = zeros(n+1)
    uu[1] = copy(u0)
    tt[1] = t0
    tj = [t0, zero(t0)]
    uj = copy(u0)
    cont = true
    for j in 1:n
        for i in 1:m
           it = IRK8step!(tj,uj,ej,p,f,h,hb,hc,mu,nu,U,Uz,L,F,Dmin,itermax,initial_interp)
           iters[j+1] += it
        end
        iters[j+1] /= m
        uu[j+1] = uj + ej
        tt[j+1] = tj[1] + tj[2]
    end
    sol = NireODESol(tt,uu,iters,:Successs,f)
    return sol
  end

function IRK8step!(ttj,uj,ej,p,f,h,hb,hc,mu,nu,U,Uz,L,F,Dmin,itermax,initial_interp)
       s = length(hb)
       elems = s*length(uj)
       tj = ttj[1]
       te = ttj[2]
       if initial_interp
          for is in 1:s
            for k in eachindex(uj)
                aux=0.
                for js in 1:s
                    aux+=nu[is,js]*L[js][k]
                end
                U[is][k]=(uj[k]+ej[k])+aux
            end
          end
       else
       for is in 1:s
            @. U[is] = uj + ej
          end
       end

       for is in 1:s
            f(F[is], U[is], p, tj + hc[is])
            @. L[is] =  hb[is]*F[is]
            Dmin[is] .= Inf
        end

        iter = true # Initialize iter outside the for loop
        plusIt=true

        j=1

        while (j<itermax && iter)

            j+=1
            iter=false
            D0=0

            for is in 1:s

                Uz[is] .= U[is]
                @. U[is] = uj + (ej + mu[is,1]*L[1] + mu[is,2]*L[2] + mu[is,3]*L[3] + mu[is,4]*L[4] +
                         mu[is,5]*L[5] + mu[is,6]*L[6] + mu[is,7]*L[7] + mu[is,8]*L[8])

            end

            for is in 1:s

                eval=false
                for k in eachindex(uj)
                            DY=abs(U[is][k]-Uz[is][k])
                            if DY>0.
                               eval=true
                               if DY< Dmin[is][k]
                                  Dmin[is][k]=DY
                                  iter=true
                               end
                           else
                               D0+=1
                           end
                end

               if eval==true
                  f(F[is], U[is], p,  tj + hc[is])
                  @. L[is] = hb[is]*F[is]
               end
           end

            if (iter==false && D0<elems && plusIt)
                iter=true
                plusIt=false
            else
                plusIt=true
            end

        end # while

        indices = eachindex(uj)

        for k in indices    #Batura konpentsatuaren parekoa
            e0 = ej[k]
            for is in 1:s
	         e0 += muladd(F[is][k], hb[is], -L[is][k])
            end
            res = Base.TwicePrecision(uj[k], e0)
            for is in 1:s
	       res += L[is][k]
            end
            uj[k] = res.hi
            ej[k] = res.lo
         end

         res = Base.TwicePrecision(tj, te) + h
         ttj[1] = res.hi
         ttj[2] = res.lo
         return  (j)

end

function IRK8Coefficients(h)

"""
Float64

    mu = [0.5 -0.081894963105581419782 0.040042703777945054533 -0.024721345803200289737 0.016976173236371128183 -0.012225914113298097519 0.0087485667691973301174 -0.0054828082532158983753
          1.0818949631055814198 0.5 -0.086958924300832629584 0.044941126302625700184 -0.028759775474749282864 0.020017127636408726943 -0.014074355889166945133 0.0087485667691973301174
          0.95995729622205494547 1.0869589243008326296 0.5 -0.088093838732308249462 0.046165464814799994642 -0.029603873064977914709 0.020017127636408726943 -0.012225914113298097519
          1.0247213458032002897 0.95505887369737429982 1.0880938387323082495 0.5 -0.088347161109827876402 0.046165464814799994642 -0.028759775474749282864 0.016976173236371128183
          0.98302382676362887182 1.0287597754747492829 0.95383453518520000536 1.0883471611098278764 0.5 -0.088093838732308249462 0.044941126302625700184 -0.024721345803200289737
          1.0122259141132980975 0.97998287236359127306 1.0296038730649779147 0.95383453518520000536 1.0880938387323082495 0.5 -0.086958924300832629584 0.040042703777945054533
          0.99125143323080266988 1.0140743558891669451 0.97998287236359127306 1.0287597754747492829 0.95505887369737429982 1.0869589243008326296 0.5 -0.081894963105581419782
          1.0054828082532158984 0.99125143323080266988 1.0122259141132980975 0.98302382676362887182 1.0247213458032002897 0.95995729622205494547 1.0818949631055814198 0.5];

          hc = h*[0.019855071751231884 0.10166676129318664 0.2372337950418355 0.4082826787521751 0.591717321247825 0.7627662049581645 0.8983332387068134 0.9801449282487681]
          hb = h*[0.05061426814518813 0.11119051722668724 0.15685332293894363 0.181341891689181 0.181341891689181 0.15685332293894363 0.11119051722668724 0.05061426814518813]

"""

    mu = convert.(typeof(h), [
               # s=1
               parse(BigFloat,"0.5") parse(BigFloat,"-8.1894963105581497136508164735930892e-02") parse(BigFloat,"+4.0042703777945052339969045601553216e-02") parse(BigFloat,"-2.4721345803200374868044581645082913e-02") parse(BigFloat,"+1.6976173236371093026881708761116103e-02 ") parse(BigFloat,"-1.2225914113298206053942531721943548e-02") parse(BigFloat,"+8.7485667691973688117766530139122287e-03") parse(BigFloat,"-5.4828082532158826793409353214950813e-03")
               # s=2
               parse(BigFloat,"+1.0818949631055814971365081647359309e+00") parse(BigFloat,"0.5") parse(BigFloat,"-8.6958924300832723329070964616248016e-02") parse(BigFloat,"+4.4941126302625688139830943466131237e-02") parse(BigFloat,"-2.8759775474749310978230557041068536e-02")  parse(BigFloat,"+2.0017127636408709173710417097426705e-02") parse(BigFloat,"-1.4074355889166929145973516652599415e-02") parse(BigFloat,"+8.7485667691973688117766530139122287e-03")
               # s=3
               parse(BigFloat,"+9.5995729622205494766003095439844678e-01") parse(BigFloat,"+1.0869589243008327233290709646162480e+00") parse(BigFloat,"0.5") parse(BigFloat,"-8.8093838732308313442213871391320316e-02") parse(BigFloat,"+4.6165464814800034116730885592456977e-02")  parse(BigFloat,"-2.9603873064977937463012598212122325e-02") parse(BigFloat,"+2.0017127636408709173710417097426705e-02") parse(BigFloat,"-1.2225914113298206053942531721943548e-02")
               # s=4
               parse(BigFloat,"+1.0247213458032003748680445816450829e+00") parse(BigFloat,"+9.5505887369737431186016905653386876e-01") parse(BigFloat,"+1.0880938387323083134422138713913203e+00") parse(BigFloat,"0.5") parse(BigFloat,"-8.8347161109827784250707380600804499e-02")  parse(BigFloat,"+4.6165464814800034116730885592456977e-02") parse(BigFloat,"-2.8759775474749310978230557041068536e-02") parse(BigFloat,"+1.6976173236371093026881708761116103e-02")
               # s=5
               parse(BigFloat,"+9.8302382676362890697311829123888390e-01") parse(BigFloat,"+1.0287597754747493109782305570410685e+00 ") parse(BigFloat,"+9.5383453518519996588326911440754302e-01") parse(BigFloat,"+1.0883471611098277842507073806008045e+00") parse(BigFloat,"0.5") parse(BigFloat,"-8.8093838732308313442213871391320316e-02") parse(BigFloat,"+4.4941126302625688139830943466131237e-02") parse(BigFloat,"-2.4721345803200374868044581645082913e-02")
               # s=6
               parse(BigFloat,"+1.0122259141132982060539425317219435e+00") parse(BigFloat,"+9.7998287236359129082628958290257329e-01")  parse(BigFloat,"+1.0296038730649779374630125982121223e+00")  parse(BigFloat,"+9.5383453518519996588326911440754302e-01") parse(BigFloat,"+1.0880938387323083134422138713913203e+00") parse(BigFloat,"0.5") parse(BigFloat,"-8.6958924300832723329070964616248016e-02") parse(BigFloat,"+4.0042703777945052339969045601553216e-02")
               # s=7
               parse(BigFloat,"+9.9125143323080263118822334698608777e-01") parse(BigFloat,"+1.0140743558891669291459735166525994e+00")  parse(BigFloat,"+9.7998287236359129082628958290257329e-01")  parse(BigFloat,"+1.0287597754747493109782305570410685e+00") parse(BigFloat,"+9.5505887369737431186016905653386876e-01 ") parse(BigFloat,"+1.0869589243008327233290709646162480e+00") parse(BigFloat,"0.5") parse(BigFloat,"-8.1894963105581497136508164735930892e-02")
               # s=8
               parse(BigFloat,"+1.0054828082532158826793409353214951e+00") parse(BigFloat,"+9.9125143323080263118822334698608777e-01")  parse(BigFloat,"+1.0122259141132982060539425317219435e+00")  parse(BigFloat,"+9.8302382676362890697311829123888390e-01") parse(BigFloat,"+1.0247213458032003748680445816450829e+00")  parse(BigFloat,"+9.5995729622205494766003095439844678e-01") parse(BigFloat,"+1.0818949631055814971365081647359309e+00") parse(BigFloat,"0.5")
        ])

        hb = convert.(typeof(h),
           h*[ parse(BigFloat,"+5.0614268145188129576265677154981094e-02"),
               parse(BigFloat,"+1.1119051722668723527217799721312045e-01"),
               parse(BigFloat,"+1.5685332293894364366898110099330067e-01"),
               parse(BigFloat,"+1.8134189168918099148257522463859781e-01"),
               parse(BigFloat,"+1.8134189168918099148257522463859781e-01"),
               parse(BigFloat,"+1.5685332293894364366898110099330067e-01"),
               parse(BigFloat,"+1.1119051722668723527217799721312045e-01"),
               parse(BigFloat,"+5.0614268145188129576265677154981094e-02")
             ])

        hc= convert.(typeof(h),
           h*[ parse(BigFloat,"+1.9855071751231884158219565715263505e-02"),
               parse(BigFloat,"+1.0166676129318663020422303176208480e-01"),
               parse(BigFloat,"+2.3723379504183550709113047540537686e-01"),
               parse(BigFloat,"+4.0828267875217509753026192881990801e-01"),
               parse(BigFloat,"+5.9171732124782490246973807118009203e-01"),
               parse(BigFloat,"+7.6276620495816449290886952459462321e-01"),
               parse(BigFloat,"+8.9833323870681336979577696823791522e-01"),
               parse(BigFloat,"+9.8014492824876811584178043428473653e-01")
        ])

""" Interpolate coefficients """
        nu = convert.(typeof(h),
           [#s=1
            parse(BigFloat,"-2.3580434359907869272604448668704418e-02") parse(BigFloat,"+3.7564028204350781120220588381324588e-02") parse(BigFloat,"-5.2311682298576367742803308364795570e-02") parse(BigFloat,"+7.2153031436804597582100216023947545e-02") parse(BigFloat,"-1.0367897445499288427629400186099435e-01") parse(BigFloat,"+1.6283634881117387717390478413870362e-01") parse(BigFloat,"-3.0234060551288978571240477869701592e-01") parse(BigFloat,"+7.6796597320813176288773924435846961e-01")
            #s=2
            parse(BigFloat,"-7.3721915347372876383416365681676195e-01") parse(BigFloat,"+1.1683495316711616559345205630057529e+00") parse(BigFloat,"-1.6094743347141669348961276646308222e+00") parse(BigFloat,"+2.1754011196351946282108735964539601e+00") parse(BigFloat,"-3.0074288183636139831253139808157757e+00") parse(BigFloat,"+4.3553757998443780247769815248837769e+00") parse(BigFloat,"-6.6165263753225524986317973726684325e+00") parse(BigFloat,"+9.1860239215521482142367391478777201e+00")
            #s=3
            parse(BigFloat,"-1.2297700288095738198718883246880343e+01") parse(BigFloat,"+1.9333875551079588083957605043596343e+01") parse(BigFloat,"-2.6197717938782395621819055691791508e+01") parse(BigFloat,"+3.4375513292689120907660371301821534e+01") parse(BigFloat,"-4.5083720476972868489541375557064921e+01") parse(BigFloat,"+5.9246620811109433433154163113334023e+01") parse(BigFloat,"-7.4772488278621580963027410652542431e+01") parse(BigFloat,"+7.4720772807309976713169982344387532e+01")
            #s=4
            parse(BigFloat,"-1.2360519192775641529208291530012872e+02") parse(BigFloat,"+1.9271209302089311449325157212657236e+02") parse(BigFloat,"-2.5678715794267672761631251206680453e+02") parse(BigFloat,"+3.2738013789187108730152371652322038e+02") parse(BigFloat,"-4.0945219322597510368323043549108768e+02") parse(BigFloat,"+4.9820118998766368260622671514051444e+02") parse(BigFloat,"-5.5851667346968918682289190242626758e+02") parse(BigFloat,"+4.8118799386789480845212083059902066e+02")
            #s=5
            parse(BigFloat,"-7.6506359873607050490475404918734150e+02") parse(BigFloat,"+1.1843177154904885692635422843492179e+03") parse(BigFloat,"-1.5561007397261703002242559137962282e+03") parse(BigFloat,"+1.9384516434527721535098498400498934e+03") parse(BigFloat,"-2.3391243164654789426290647731129215e+03") parse(BigFloat,"+2.7004970213121135878376040434322658e+03") parse(BigFloat,"-2.8241828008846508731118675847222545e+03") parse(BigFloat,"+2.2683066185998760812106944011486825e+03")
            #s=6
            parse(BigFloat,"-2.9864681431085356613428555972148677e+03") parse(BigFloat,"+4.5981543451008061740598602621277281e+03") parse(BigFloat,"-5.9789738886452262969037360399834883e+03") parse(BigFloat,"+7.3248683186540351736460833184871867e+03") parse(BigFloat,"-8.6243827952223739681429053998929657e+03") parse(BigFloat,"+9.6280058867130293266015441780974974e+03") parse(BigFloat,"-9.6712897035873639047772333270447980e+03") parse(BigFloat,"+7.4939111643640677935209888348628654e+03")
            #s=7
            parse(BigFloat,"-7.5322560258736694963352803649936405e+03") parse(BigFloat,"+1.1555244400412346782160633692641882e+04") parse(BigFloat,"-1.4922008897475610276386022948789031e+04") parse(BigFloat,"+1.8084657457385087410934304600504492e+04") parse(BigFloat,"-2.0967791943855618026873030594152413e+04") parse(BigFloat,"+2.2942168855337640310062689548583393e+04") parse(BigFloat,"-2.2527726562941111873889414492088629e+04") parse(BigFloat,"+1.7129917709645864090891579107793160e+04")
            #s=8
            parse(BigFloat,"-1.2513853497502627647496026877301948e+04") parse(BigFloat,"+1.9160065004504862609943292716902368e+04") parse(BigFloat,"-2.4651395180933566824098540012109580e+04") parse(BigFloat,"+2.9706130282097634507632738110877420e+04") parse(BigFloat,"-3.4168328879509508253154372762217900e+04") parse(BigFloat,"+3.7008535511486677576566037943315160e+04") parse(BigFloat,"-3.5939320855736553497649468561881539e+04") parse(BigFloat,"+2.7086805366868524974942204129999368e+04")
        ])

s = length(hb)
for i in 1:s
    for j in i+1:s
        mu[i,j] = 1 - mu[j,i]
    end
end

hb1 = (h-sum(hb[2:end-1]))/2
hb[1] = hb1
hb[end] = hb1

return (hb, hc, mu, nu)
end


text:
get_aggregation_topology_accessor(::Type{Area}) = get_area
get_aggregation_topology_accessor(::Type{LoadZone}) = get_load_zone

set_load_zone!(bus::Bus, load_zone::LoadZone) = bus.load_zone = load_zone
set_area!(bus::Bus, area::Area) = bus.area = area

"""
Remove the aggregation topology in a Bus
"""
_remove_aggregration_topology!(bus::Bus, ::LoadZone) = bus.load_zone = nothing
_remove_aggregration_topology!(bus::Bus, ::Area) = bus.area = nothing

"""
Calculate the admittance of AC branches
"""
get_series_susceptance(b::ACBranch) = 1 / get_x(b)

"""
Returns the series susceptance of a controllable transformer following the convention
in power systems to define susceptance as the inverse of the imaginary part of the impedance.
In the case of phase shifter transformers the angle is ignored.
"""
function get_series_susceptance(b::Union{PhaseShiftingTransformer, TapTransformer})
    y = 1 / get_x(b)
    y_a = y / (get_tap(b))
    return y_a
end

get_series_admittance(b::ACBranch) = 1 / (get_r(b) + get_x(b) * 1im)

"""
Return the max active power for a device from get_active_power_limits.max
"""
function get_max_active_power(d::T) where {T <: Device}
    return get_active_power_limits(d).max
end

"""
Return the max reactive power for a device from get_reactive_power_limits.max
"""
function get_max_reactive_power(d::T) where {T <: Device}
    return get_reactive_power_limits(d).max
end

"""
Return the max reactive power for the Renewable Generation calculated as the rating * power_factor if
reactive_power_limits is nothing
"""
function get_max_reactive_power(d::RenewableDispatch)
    reactive_power_limits = get_reactive_power_limits(d)
    if isnothing(reactive_power_limits)
        return get_rating(d) * sin(acos(get_power_factor(d)))
    end
    return reactive_power_limits.max
end

get_active_power_limits(::T) where {T <: Device} =
    throw(ArgumentError("get_active_power_limits not implemented for $T"))
get_reactive_power_limits(::T) where {T <: Device} =
    throw(ArgumentError("get_reactive_power_limits not implemented for $T"))
get_rating(::T) where {T <: Device} =
    throw(ArgumentError("get_rating not implemented for $T"))
get_power_factor(::T) where {T <: Device} =
    throw(ArgumentError("get_power_factor not implemented for $T"))