text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Functions related with the model IAU-76/FK5.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
#   [2] Gontier, A. M., Capitaine, N (1991). High-Accuracy Equation of Equinoxes
#       and VLBI Astrometric Modelling. Radio Interferometry: Theory, Techniques
#       and Applications, IAU Coll. 131, ASP Conference Series, Vol. 19.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

export rITRFtoPEF_fk5, rPEFtoITRF_fk5
export rPEFtoTOD_fk5,  rTODtoPEF_fk5
export rTODtoMOD_fk5,  rMODtoTOD_fk5
export rMODtoGCRF_fk5, rGCRFtoMOD_fk5

export rITRFtoGCRF_fk5, rGCRFtoITRF_fk5
export rPEFtoMOD_fk5,   rMODtoPEF_fk5

################################################################################
#                            IAU-76/FK5 Reductions
################################################################################
#
# The conversion between the Geocentric Celestial Reference Frame (GCRF) to the
# International Terrestrial Reference Frame (ITRF) is done by means of:
#
#                    GCRF <=> MOD <=> TOD <=> PEF <=> ITRF
#
# in which:
#   - MOD: Mean of Date frame.
#   - TOD: True of Date frame.
#   - PEF: Pseudo-Earth fixed frame.
#
# Every rotation will be coded as a function using the IAU-76/FK5 theory.
# Additionally, composed rotations will also available. In general, the API is:
#
#   function r<Origin Frame>to<Destination Frame>_fk5
#
# The arguments vary depending on the origin and destination frame and should be
# verified using the function documentation.
#
################################################################################

################################################################################
#                               Single Rotations
################################################################################

#                                 ITRF <=> PEF
# ==============================================================================

"""
    function rITRFtoPEF_fk5([T,] x_p::Number, y_p::Number)

Compute the rotation that aligns the International Terrestrial Reference Frame
(ITRF) with the Pseudo-Earth Fixed (PEF) frame considering the polar motion
represented by the angles `x_p` [rad] and `y_p` [rad] that are obtained from
IERS EOP Data (see `get_iers_eop`).

`x_p` is the polar motion displacement about X-axis, which is the IERS Reference
Meridian direction (positive south along the 0˚ longitude meridian). `y_p` is
the polar motion displacement about Y-axis (90˚W or 270˚E meridian).

The rotation type is described by the optional variable `T`. If it is `DCM`,
then a DCM will be returned. Otherwise, if it is `Quaternion`, then a Quaternion
will be returned. In case this parameter is omitted, then it falls back to
`DCM`.

# Returns

The rotation that aligns the ITRF frame with the PEF frame. The rotation
representation is selected by the optional parameter `T`.

# Remarks

The ITRF is defined based on the International Reference Pole (IRP), which is
the location of the terrestrial pole agreed by international committees [1]. The
Pseudo-Earth Fixed, on the other hand, is defined based on the Earth axis of
rotation, or the Celestial Intermediate Pole (CIP). Hence, PEF XY-plane contains
the True Equator. Furthermore, since the recovered latitude and longitude are
sensitive to the CIP, then it should be computed considering the PEF frame.

"""
rITRFtoPEF_fk5(x_p::Number, y_p::Number) = rITRFtoPEF_fk5(DCM, x_p, y_p)

rITRFtoPEF_fk5(T::Type, x_p::Number, y_p::Number) =
    # Notice that `x_p` and `y_p` are displacements in X and Y directions and
    # **not** rotation angles. Hence, a displacement in X is a rotation in Y and
    # a displacement in Y is a rotation in X.
    smallangle_to_rot(T, +y_p, +x_p, 0)

"""
    function rPEFtoITRF_fk5([T,] x_p::Number, y_p::Number)

Compute the rotation that aligns the Pseudo-Earth Fixed (PEF) with the
International Terrestrial Reference Frame (ITRF) considering the polar motion
represented by the angles `x_p` [rad] and `y_p` [rad] that are obtained from
IERS EOP Data (see `get_iers_eop`).

`x_p` is the polar motion displacement about X-axis, which is the IERS Reference
Meridian direction (positive south along the 0˚ longitude meridian). `y_p` is
the polar motion displacement about Y-axis (90˚W or 270˚E meridian).

The rotation type is described by the optional variable `T`. If it is `DCM`,
then a DCM will be returned. Otherwise, if it is `Quaternion`, then a Quaternion
will be returned. In case this parameter is omitted, then it falls back to
`DCM`.

# Returns

The rotation that aligns the PEF frame with the ITRF. The rotation
representation is selected by the optional parameter `T`.

# Remarks

The ITRF is defined based on the International Reference Pole (IRP), which is
the location of the terrestrial pole agreed by international committees [1]. The
Pseudo-Earth Fixed, on the other hand, is defined based on the Earth axis of
rotation, or the Celestial Intermediate Pole (CIP). Hence, PEF XY-plane contains
the True Equator. Furthermore, since the recovered latitude and longitude are
sensitive to the CIP, then it should be computed considering the PEF frame.

"""
rPEFtoITRF_fk5(x_p::Number, y_p::Number) = rPEFtoITRF_fk5(DCM, x_p, y_p)

rPEFtoITRF_fk5(T::Type, x_p::Number, y_p::Number) =
    # Notice that `x_p` and `y_p` are displacements in X and Y directions and
    # **not** rotation angles. Hence, a displacement in X is a rotation in Y and
    # a displacement in Y is a rotation in X.
    smallangle_to_rot(T, -y_p, -x_p, 0)

#                                 PEF <=> TOD
# ==============================================================================

"""
    function rPEFtoTOD_fk5([T,] JD_UT1::Number, JD_TT::Number [, δΔψ_1980::Number])

Compute the rotation that aligns the Pseudo-Earth Fixed (PEF) frame with the
True of Date (TOD) frame at the Julian Day `JD_UT1` [UT1] and `JD_TT`
[Terrestrial Time]. This algorithm uses the IAU-76/FK5 theory. Notice that one
can provide correction for the nutation in longitude (`δΔψ_1980`) \\[rad] that
is usually obtained from IERS EOP Data (see `get_iers_eop`).

The Julian Day in UT1 is used to compute the Greenwich Mean Sidereal Time (GMST)
(see `JDtoGMST`), whereas the Julian Day in Terrestrial Time is used to compute
the nutation in the longitude. Notice that the Julian Day in UT1 and in
Terrestrial Time must be equivalent, i.e. must be related to the same instant.
This function **does not** check this.

The rotation type is described by the optional variable `T`. If it is `DCM`,
then a DCM will be returned. Otherwise, if it is `Quaternion`, then a Quaternion
will be returned. In case this parameter is omitted, then it falls back to
`DCM`.

# Returns

The rotation that aligns the PEF frame with the TOD frame. The rotation
representation is selected by the optional parameter `T`.

# Remarks

The Pseudo-Earth Fixed (PEF) frame is rotated into the True of Date (TOD) frame
considering the 1980 IAU Theory of Nutation. The IERS EOP corrections must be
added if one wants to make the rotation consistent with the Geocentric Celestial
Reference Systems (GCRS).

"""
rPEFtoTOD_fk5(JD_UT1::Number, JD_TT::Number, δΔψ_1980::Number = 0) =
    rPEFtoTOD_fk5(DCM, JD_UT1, JD_TT, δΔψ_1980)

function rPEFtoTOD_fk5(T::Type,
                       JD_UT1::Number,
                       JD_TT::Number,
                       δΔψ_1980::Number = 0)
    # Compute the nutation in the Julian Day (Terrestrial Time) `JD_TT`.
    (mϵ_1980, Δϵ_1980, Δψ_1980) = nutation_fk5(JD_TT)

    # Add the corrections to the nutation in obliquity and longitude.
    Δψ_1980 += δΔψ_1980

    # Evaluate the Delaunay parameters associated with the Moon in the interval
    # [0,2π]°.
    #
    # The parameters here were updated as stated in the errata [2].
    T_TT = (JD_TT - JD_J2000)/36525
    r    = 360
    Ω_m  = @evalpoly(T_TT, + 125.04452222,
                           - (5r + 134.1362608),
                           + 0.0020708,
                           + 2.2e-6)
    Ω_m  = mod(Ω_m, 360)*pi/180

    # Compute the equation of Equinoxes.
    #
    # According to [2], the constant unit before `sin(2Ω_m)` is also in [rad].
    Eq_equinox1982 = Δψ_1980*cos(mϵ_1980) +
        ( 0.002640*sin(1Ω_m) + 0.000063*sin(2Ω_m) )*pi/648000

    # Compute the Mean Greenwich Sidereal Time.
    θ_gmst = JDtoGMST(JD_UT1)

    # Compute the Greenwich Apparent Sidereal Time (GAST).
    #
    # TODO: Should GAST be moved to a new function as the GMST?
    θ_gast = θ_gmst + Eq_equinox1982

    # Compute the rotation matrix.
    angle_to_rot(T, -θ_gast, 0, 0, :ZYX)
end

"""
    function rTODtoPEF_fk5([T,] JD_UT1::Number, JD_TT::Number [, δΔψ_1980::Number])

Compute the rotation that aligns the True of Date (TOD) frame with the
Pseudo-Earth Fixed (PEF) frame at the Julian Day `JD_UT1` [UT1] and `JD_TT`
[Terrestrial Time]. This algorithm uses the IAU-76/FK5 theory. Notice that one
can provide correction for the nutation in longitude (`δΔψ_1980`) \\[rad] that
is usually obtained from IERS EOP Data (see `get_iers_eop`).

The Julian Day in UT1 is used to compute the Greenwich Mean Sidereal Time (GMST)
(see `JDtoGMST`), whereas the Julian Day in Terrestrial Time is used to compute
the nutation in the longitude. Notice that the Julian Day in UT1 and in
Terrestrial Time must be equivalent, i.e. must be related to the same instant.
This function **does not** check this.

The rotation type is described by the optional variable `T`. If it is `DCM`,
then a DCM will be returned. Otherwise, if it is `Quaternion`, then a Quaternion
will be returned. In case this parameter is omitted, then it falls back to
`DCM`.

# Returns

The rotation that aligns the TOD frame with the PEF frame. The rotation
representation is selected by the optional parameter `T`.

# Remarks

The True of Date (TOD) frame is rotated into the Pseudo-Earth Fixed (PEF) frame
considering the 1980 IAU Theory of Nutation. The IERS EOP corrections must be
added if one wants to make the rotation consistent with the Geocentric Celestial
Reference Systems (GCRS).

"""
rTODtoPEF_fk5(JD_UT1::Number, JD_TT::Number, δΔψ_1980::Number = 0) =
    rTODtoPEF_fk5(DCM, JD_UT1, JD_TT, δΔψ_1980)

rTODtoPEF_fk5(T::T_ROT, JD_UT1::Number, JD_TT::Number, δΔψ_1980::Number = 0) =
    inv_rotation(rPEFtoTOD_fk5(T, JD_UT1, JD_TT, δΔψ_1980))

#                                 TOD <=> MOD
# ==============================================================================

"""
    function rTODtoMOD_fk5([T,] JD_TT::Number [, δΔϵ_1980::Number, δΔψ_1980::Number])

Compute the rotation that aligns the True of Date (TOD) frame with the Mean of
Date (MOD) frame at the Julian Day `JD_TT` [Terrestrial Time]. This algorithm
uses the IAU-76/FK5 theory. Notice that one can provide corrections for the
nutation in obliquity (`δΔϵ_1980`) \\[rad] and in longitude (`δΔψ_1980`)
\\[rad] that are usually obtained from IERS EOP Data (see `get_iers_eop`).

The rotation type is described by the optional variable `T`. If it is `DCM`,
then a DCM will be returned. Otherwise, if it is `Quaternion`, then a Quaternion
will be returned. In case this parameter is omitted, then it falls back to
`DCM`.

# Returns

The rotation that aligns the TOD frame with the MOD frame. The rotation
representation is selected by the optional parameter `T`.

# Remarks

The True of Date (TOD) frame is rotated into the Mean of Date (MOD) frame
considering the 1980 IAU Theory of Nutation. The IERS EOP corrections must be
added if one wants to make the rotation consistent with the Geocentric Celestial
Reference Systems (GCRS).

"""
rTODtoMOD_fk5(JD_TT::Number, δΔϵ_1980::Number = 0, δΔψ_1980::Number = 0) =
    rTODtoMOD_fk5(DCM, JD_TT, δΔϵ_1980, δΔψ_1980)

function rTODtoMOD_fk5(T::Type,
                       JD_TT::Number,
                       δΔϵ_1980::Number = 0,
                       δΔψ_1980::Number = 0)

    # Compute the nutation in the Julian Day (Terrestrial Time) `JD_TT`.
    (mϵ_1980, Δϵ_1980, Δψ_1980) = nutation_fk5(JD_TT)

    # Add the corrections to the nutation in obliquity and longitude.
    Δϵ_1980 += δΔϵ_1980
    Δψ_1980 += δΔψ_1980

    # Compute the obliquity.
    ϵ_1980 = mϵ_1980 + Δϵ_1980

    # Compute and return the Direction Cosine DCM.
    angle_to_rot(T, ϵ_1980, Δψ_1980, -mϵ_1980, :XZX)
end

"""
    function rMODtoTOD_fk5([T,] JD_TT::Number [, δΔϵ_1980::Number, δΔψ_1980::Number])

Compute the rotation that aligns the Mean of Date (MOD) frame with the True of
Date (TOD) frame at the Julian Day `JD_TT` [Terrestrial Time]. This algorithm
uses the IAU-76/FK5 theory. Notice that one can provide corrections for the
nutation in obliquity (`δΔϵ_1980`) \\[rad] and in longitude (`δΔψ_1980`) \\[rad]
that are usually obtained from IERS EOP Data (see `get_iers_eop`).

The rotation type is described by the optional variable `T`. If it is `DCM`,
then a DCM will be returned. Otherwise, if it is `Quaternion`, then a Quaternion
will be returned. In case this parameter is omitted, then it falls back to
`DCM`.

# Returns

The rotation that aligns the MOD frame with the TOD frame. The rotation
representation is selected by the optional parameter `T`.

# Remarks

The Mean of Date (MOD) frame is rotated into the True of Date (TOD) frame
considering the 1980 IAU Theory of Nutation. The IERS EOP corrections must be
added if one wants to make the rotation consistent with the Geocentric Celestial
Reference Systems (GCRS).

"""
rMODtoTOD_fk5(JD_TT::Number, δΔϵ_1980::Number = 0, δΔψ_1980::Number = 0) =
    rMODtoTOD_fk5(DCM, JD_TT, δΔϵ_1980, δΔψ_1980)

rMODtoTOD_fk5(T::T_ROT, JD_TT::Number, δΔϵ_1980::Number = 0, δΔψ_1980::Number = 0) =
    inv_rotation(rTODtoMOD_fk5(T, JD_TT, δΔϵ_1980, δΔψ_1980))

#                                 MOD <=> GCRF
# ==============================================================================

"""
    function rMODtoGCRF_fk5([T,] JD_TT::Number)

Compute the rotation that aligns the Mean of Date (MOD) frame with the
Geocentric Celestial Reference Frame (GCRF) at the Julian Day [Terrestrial Time]
`JD_TT`. This algorithm uses the IAU-76/FK5 theory.

The rotation type is described by the optional variable `T`. If it is `DCM`,
then a DCM will be returned. Otherwise, if it is `Quaternion`, then a Quaternion
will be returned. In case this parameter is omitted, then it falls back to
`DCM`.

# Returns

The rotation that aligns the MOD frame with the GCRF frame. The rotation
representation is selected by the optional parameter `T`.

# Remarks

The Mean of Date (MOD) frame is rotated into the Geocentric Celestial Reference
Frame (GCRF) considering the IAU 1976 Precession model.

Notice that if the conversion `TOD => MOD` is performed **without** considering
the EOP corrections, then the GCRF obtained by this rotation is what is usually
called the J2000 reference frame.

"""
rMODtoGCRF_fk5(JD_TT::Number) = rMODtoGCRF_fk5(DCM,JD_TT)

function rMODtoGCRF_fk5(T::Type, JD_TT::Number)
    (ζ,Θ,z) = precession_fk5(JD_TT)
    angle_to_rot(T, z, -Θ, ζ, :ZYZ)
end

"""
    function rGCRFtoMOD_fk5([T,] JD_TT::Number)

Compute the rotation that aligns the Geocentric Celestial Reference Frame (GCRF)
with the Mean of Date (MOD) frame at the Julian Day [Terrestrial Time] `JD_TT`.
This algorithm uses the IAU-76/FK5 theory.

The rotation type is described by the optional variable `T`. If it is `DCM`,
then a DCM will be returned. Otherwise, if it is `Quaternion`, then a Quaternion
will be returned. In case this parameter is omitted, then it falls back to
`DCM`.

# Returns

The rotation that aligns the GCRF frame with the MOD frame. The rotation
representation is selected by the optional parameter `T`.

# Remarks

The Geocentric Celestial Reference Frame (GCRF) is rotated into the Mean of Date
(MOD) frame considering the IAU 1976 Precession model.

Notice that if the conversion `MOD => TOD` is performed **without** considering
the EOP corrections, then the GCRF in this rotation is what is usually called
the J2000 reference frame.

"""
rGCRFtoMOD_fk5(JD_TT::Number) = rGCRFtoMOD_fk5(DCM,JD_TT)

rGCRFtoMOD_fk5(T::T_ROT,JD_TT::Number) = inv_rotation(rMODtoGCRF_fk5(T, JD_TT))

################################################################################
#                              Multiple Rotations
################################################################################

# The functions with multiple rotations must be added only in two cases:
#
#   * ITRF <=> GCRF (Full rotation between ECI and ECEF).
#   * When the it will decrease the computational burden compared to
#     calling the functions with the single rotations.
#

#                                ITRF <=> GCRF
# ==============================================================================

"""
    function rITRFtoGCRF_fk5([T,] JD_UT1::Number, JD_TT::Number, x_p::Number, y_p::Number [, δΔϵ_1980::Number, δΔψ_1980::Number])

Compute the rotation that aligns the International Terrestrial Reference Frame
(ITRF) with the Geocentric Celestial Reference Frame (GCRF) at the Julian Day
`JD_UT1` [UT1] and `JD_TT` [Terrestrial Time], and considering the IERS EOP Data
`x_p` [rad], `y_p` [rad], `δΔϵ_1980` [rad], and `δΔψ_1980` [rad] \\(see
`get_iers_eop`). This algorithm uses the IAU-76/FK5 theory.

`x_p` is the polar motion displacement about X-axis, which is the IERS Reference
Meridian direction (positive south along the 0˚ longitude meridian). `y_p` is
the polar motion displacement about Y-axis (90˚W or 270˚E meridian). `δΔϵ_1980`
is the nutation in obliquity. `δΔψ_1980` is the nutation in longitude.

The Julian Day in UT1 is used to compute the Greenwich Mean Sidereal Time (GMST)
(see `JDtoGMST`), whereas the Julian Day in Terrestrial Time is used to compute
the nutation in the longitude. Notice that the Julian Day in UT1 and in
Terrestrial Time must be equivalent, i.e. must be related to the same instant.
This function **does not** check this.

The rotation type is described by the optional variable `T`. If it is `DCM`,
then a DCM will be returned. Otherwise, if it is `Quaternion`, then a Quaternion
will be returned. In case this parameter is omitted, then it falls back to
`DCM`.

# Returns

The rotation that aligns the ITRF frame with the GCRF frame. The rotation
representation is selected by the optional parameter `T`.

# Remarks

The EOP data related to the polar motion (`x_p` and `y_p`) is required, since
this is the only way available to compute the conversion ITRF <=> PEF (the
models are highly imprecise since the motion is still not very well understood
[1]). However, the EOP data related to the nutation of the obliquity
(`δΔϵ_1980`) and the nutation of the longitude (`δΔψ_1980`) can be omitted. In
this case, the GCRF frame is what is usually called J2000 reference frame.

"""
rITRFtoGCRF_fk5(JD_UT1::Number,
                JD_TT::Number,
                x_p::Number,
                y_p::Number,
                δΔϵ_1980::Number = 0,
                δΔψ_1980::Number = 0) =
    rITRFtoGCRF_fk5(DCM, JD_UT1, JD_TT, x_p, y_p, δΔϵ_1980, δΔψ_1980)

function rITRFtoGCRF_fk5(T::Type,
                         JD_UT1::Number,
                         JD_TT::Number,
                         x_p::Number,
                         y_p::Number,
                         δΔϵ_1980::Number = 0,
                         δΔψ_1980::Number = 0)
    # Compute the rotation ITRF => PEF.
    r_PEF_ITRF = rITRFtoPEF_fk5(T, x_p, y_p)

    # Compute the rotation PEF => MOD.
    r_MOD_PEF = rPEFtoMOD_fk5(T, JD_UT1, JD_TT, δΔϵ_1980, δΔψ_1980)

    # Compute the rotation MOD => GCRF.
    r_GCRF_MOD = rMODtoGCRF_fk5(T, JD_TT)

    # Return the full rotation.
    compose_rotation(r_PEF_ITRF, r_MOD_PEF, r_GCRF_MOD)
end

"""
    function rGCRFtoITRF_fk5([T,] JD_UT1::Number, JD_TT::Number, x_p::Number, y_p::Number [, δΔϵ_1980::Number, δΔψ_1980::Number])

Compute the rotation that aligns the Geocentric Celestial Reference Frame (GCRF)
with the International Terrestrial Reference Frame (ITRF) at the Julian Day
`JD_UT1` [UT1] and `JD_TT` [Terrestrial Time], and considering the IERS EOP Data
`x_p` [rad], `y_p` [rad], `δΔϵ_1980` [rad], and `δΔψ_1980` [rad] \\(see
`get_iers_eop`). This algorithm uses the IAU-76/FK5 theory.

`x_p` is the polar motion displacement about X-axis, which is the IERS Reference
Meridian direction (positive south along the 0˚ longitude meridian). `y_p` is
the polar motion displacement about Y-axis (90˚W or 270˚E meridian). `δΔϵ_1980`
is the nutation in obliquity. `δΔψ_1980` is the nutation in longitude.

The Julian Day in UT1 is used to compute the Greenwich Mean Sidereal Time (GMST)
(see `JDtoGMST`), whereas the Julian Day in Terrestrial Time is used to compute
the nutation in the longitude. Notice that the Julian Day in UT1 and in
Terrestrial Time must be equivalent, i.e. must be related to the same instant.
This function **does not** check this.

The rotation type is described by the optional variable `T`. If it is `DCM`,
then a DCM will be returned. Otherwise, if it is `Quaternion`, then a Quaternion
will be returned. In case this parameter is omitted, then it falls back to
`DCM`.

# Returns

The rotation that aligns the GCRF frame with the ITRF frame. The rotation
representation is selected by the optional parameter `T`.

# Remarks

The EOP data related to the polar motion (`x_p` and `y_p`) is required, since
this is the only way available to compute the conversion ITRF <=> PEF (the
models are highly imprecise since the motion is still not very well understood
[1]). However, the EOP data related to the nutation of the obliquity
(`δΔϵ_1980`) and the nutation of the longitude (`δΔψ_1980`) can be omitted. In
this case, the GCRF frame is what is usually called J2000 reference frame.

"""
rGCRFtoITRF_fk5(JD_UT1::Number,
                JD_TT::Number,
                x_p::Number,
                y_p::Number,
                δΔϵ_1980::Number = 0,
                δΔψ_1980::Number = 0) =
    rGCRFtoITRF_fk5(DCM, JD_UT1, JD_TT, x_p, y_p, δΔϵ_1980, δΔψ_1980)

rGCRFtoITRF_fk5(T::T_ROT,
                JD_UT1::Number,
                JD_TT::Number,
                x_p::Number,
                y_p::Number,
                δΔϵ_1980::Number = 0,
                δΔψ_1980::Number = 0) =
    inv_rotation(rITRFtoGCRF_fk5(T, JD_UT1, JD_TT, x_p, y_p, δΔϵ_1980, δΔψ_1980))

#                                 PEF <=> MOD
# ==============================================================================

"""
    function rPEFtoMOD_fk5([T,] JD_UT1::Number, JD_TT::Number [, δΔϵ_1980::Number, δΔψ_1980::Number])

Compute the rotation that aligns the Pseudo-Earth Fixed (PEF) frame with the
Mean of Date (MOD) at the Julian Day `JD_UT1` [UT1] and `JD_TT` [Terrestrial
Time]. This algorithm uses the IAU-76/FK5 theory. Notice that one can provide
corrections for the nutation in obliquity (`δΔϵ_1980`) \\[rad] and in longitude
(`δΔψ_1980`) \\[rad] that are usually obtained from IERS EOP Data (see
`get_iers_eop`).

The Julian Day in UT1 is used to compute the Greenwich Mean Sidereal Time (GMST)
(see `JDtoGMST`), whereas the Julian Day in Terrestrial Time is used to compute
the nutation in the longitude. Notice that the Julian Day in UT1 and in
Terrestrial Time must be equivalent, i.e. must be related to the same instant.
This function **does not** check this.

The rotation type is described by the optional variable `T`. If it is `DCM`,
then a DCM will be returned. Otherwise, if it is `Quaternion`, then a Quaternion
will be returned. In case this parameter is omitted, then it falls back to
`DCM`.

# Returns

The rotation that aligns the PEF frame with the TOD frame. The rotation
representation is selected by the optional parameter `T`.

"""
rPEFtoMOD_fk5(JD_UT1::Number,
              JD_TT::Number,
              δΔϵ_1980::Number = 0,
              δΔψ_1980::Number = 0) =
    rPEFtoMOD_fk5(DCM, JD_UT1, JD_TT, δΔϵ_1980, δΔψ_1980)

function rPEFtoMOD_fk5(T::Type,
                       JD_UT1::Number,
                       JD_TT::Number,
                       δΔϵ_1980::Number = 0,
                       δΔψ_1980::Number = 0)

    # Notice that, in this case, we will not use `rPEFtoTOD` and `rTODtoMOD`
    # because this would call the function `nutation` twice, leading to a huge
    # performance drop. Hence, the code of those two functions is almost
    # entirely rewritten here.

    # Compute the nutation in the Julian Day (Terrestrial Time) `JD_TT`.
    (mϵ_1980, Δϵ_1980, Δψ_1980) = nutation_fk5(JD_TT)

    # Add the corrections to the nutation in obliquity and longitude.
    Δϵ_1980 += δΔϵ_1980
    Δψ_1980 += δΔψ_1980

    # Compute the obliquity.
    ϵ_1980 = mϵ_1980 + Δϵ_1980

    # Evaluate the Delaunay parameters associated with the Moon in the interval
    # [0,2π]°.
    #
    # The parameters here were updated as stated in the errata [2].
    T_TT = (JD_TT - JD_J2000)/36525
    r    = 360
    Ω_m  = @evalpoly(T_TT, + 125.04452222,
                           - (5r + 134.1362608),
                           + 0.0020708,
                           + 2.2e-6)
    Ω_m  = mod(Ω_m, 360)*pi/180

    # Compute the equation of Equinoxes.
    #
    # According to [2], the constant unit before `sin(2Ω_m)` is also in [rad].
    Eq_equinox1982 = Δψ_1980*cos(mϵ_1980) +
        ( 0.002640*sin(1Ω_m) + 0.000063*sin(2Ω_m) )*pi/648000

    # Compute the Mean Greenwich Sidereal Time.
    θ_gmst = JDtoGMST(JD_UT1)

    # Compute the Greenwich Apparent Sidereal Time (GAST).
    #
    # TODO: Should GAST be moved to a new function as the GMST?
    θ_gast = θ_gmst + Eq_equinox1982

    # Compute the rotation PEF => TOD.
    r_TOD_PEF = angle_to_rot(T, -θ_gast, 0, 0, :ZYX)

    # Compute the rotation TOD => MOD.
    r_MOD_TOD = angle_to_rot(T, ϵ_1980, Δψ_1980, -mϵ_1980, :XZX)

    compose_rotation(r_TOD_PEF, r_MOD_TOD)
end

"""
    function rMODtoPEF_fk5([T,] JD_UT1::Number, JD_TT::Number [, δΔϵ_1980::Number, δΔψ_1980::Number])

Compute the rotation that aligns the Mean of Date (MOD) reference frame with the
Pseudo-Earth Fixed (PEF) frame at the Julian Day `JD_UT1` [UT1] and `JD_TT`
[Terrestrial Time]. This algorithm uses the IAU-76/FK5 theory. Notice that one
can provide corrections for the nutation in obliquity (`δΔϵ_1980`) \\[rad] and
in longitude (`δΔψ_1980`) \\[rad] that are usually obtained from IERS EOP Data
(see `get_iers_eop`).

The Julian Day in UT1 is used to compute the Greenwich Mean Sidereal Time (GMST)
(see `JDtoGMST`), whereas the Julian Day in Terrestrial Time is used to compute
the nutation in the longitude. Notice that the Julian Day in UT1 and in
Terrestrial Time must be equivalent, i.e. must be related to the same instant.
This function **does not** check this.

The rotation type is described by the optional variable `T`. If it is `DCM`,
then a DCM will be returned. Otherwise, if it is `Quaternion`, then a Quaternion
will be returned. In case this parameter is omitted, then it falls back to
`DCM`.

# Returns

The rotation that aligns the MOD frame with the PEF frame. The rotation
representation is selected by the optional parameter `T`.

"""
rMODtoPEF_fk5(JD_UT1::Number,
              JD_TT::Number,
              δΔϵ_1980::Number = 0,
              δΔψ_1980::Number = 0) =
    rMODtoPEF_fk5(DCM, JD_UT1, JD_TT, δΔϵ_1980, δΔψ_1980)

rMODtoPEF_fk5(T::T_ROT,
              JD_UT1::Number,
              JD_TT::Number,
              δΔϵ_1980::Number = 0,
              δΔψ_1980::Number = 0) =
    inv_rotation(rPEFtoMOD_fk5(T, JD_UT1, JD_TT, δΔϵ_1980, δΔψ_1980))


text:
using MDP
using Test

include("./valueIteration.jl")
include("./badmodels.jl")

ext: jl


text:
module LittleGibbs

using AbstractMCMC
import AbstractMCMC: step!, transition_type
using Random

export Gibbs, BlockModel, Transition

struct Gibbs{B, T<:NamedTuple} <: AbstractSampler
    init_θ::T

    Gibbs(init_θ::NamedTuple{B, D}) where {B, D} = new{B, NamedTuple{B, D}}(init_θ)
end

struct BlockModel{B, C<:NamedTuple} <: AbstractModel
    conditionals::C

    BlockModel(conditionals::NamedTuple{B, D}) where {B, D} = new{B, NamedTuple{B, D}}(conditionals)
end

struct Transition{B, T<:NamedTuple} <: AbstractTransition
    params::T

    Transition(params::NamedTuple{B, D}) where {B, D} = new{B, NamedTuple{B, D}}(params)
end

transition_type(model::BlockModel{B, T}, spl::Gibbs{B}) where {B, T} = Transition{B, T}

function step!(rng::AbstractRNG, model::BlockModel{B}, spl::Gibbs{B},
               ::Integer; kwargs...) where {B}
    return Transition(spl.init_θ)
end

# perform one determinstic scan update (i.e., all conditionals, once, in order of storage)
function step!(rng::AbstractRNG, model::BlockModel{B}, spl::Gibbs{B, T},
               ::Integer, θ_prev::Transition{B, T}; kwargs...) where {B, T}
    params = θ_prev.params
    conditionals = model.conditionals
    return Transition(updated(rng, conditionals, params))
end

"""
    updated(rng, conditionals::NamedTuple, params::NamedTuple)

Update all `params` using the `conditionals` after each other, i.e.
```
(θ_{1}, θ_{2}, ..., θ_{N}) <- (conditionals[1](θ_{2},...,θ_{N}), θ_{2}, ..., θ_{N})
(θ_{1}, θ_{2}, θ_{3:end}) <- (θ_{1}, conditionals[2](θ_{1}, θ_{3}, ..., θ_{N}), θ_{3}, ..., θ_{N})
...
(θ_{1}, ..., θ_{N-1}, θ_{N}) <- (θ_{1}, ..., θ_{N-1}, conditionals[N](θ_{1}, ..., θ_{N-1}))
```

(Thanks David Widmann for this most elegant implementation!)
"""
function updated end

updated(rng, conditionals::NamedTuple{()}, params::NamedTuple) = params

function updated(rng, conditionals::NamedTuple{B}, params::NamedTuple) where B
    conditional_dist = first(conditionals)
    conditioned_values = Base.tail(params)
    return updated(rng, Base.tail(conditionals),
                   merge(conditioned_values,
                        (; first(B) => rand(rng, conditional_dist(conditioned_values)))))
end

end # module


text:
using GHOST
setup()
@unpack conn, schema, pat = GHOST.PARALLELENABLER

users = execute(conn,
                """
                SELECT author
                FROM $schema.usr_to_get_login_full A
                LEFT JOIN $schema.test_usr B
                ON A.author = B.id
                WHERE B.id is null
                ;
                """) |>
    (obj -> getproperty.(obj, :author))
users = [ data[i:min(i + 4_999, lastindex(users))] for i in 1:5_000:length(users) ]

function query_users(users::AbstractVector{<:AbstractString})
    nodes = [ users[i:min(i + 99, lastindex(users))] for i in 1:100:length(users) ]
    vars = Dict(zip(string.("x", eachindex(nodes)), nodes))
    query = string("fragment a on Node {id __typename",
                   "...on User{login} ",
                   "...on Organization {login} ",
                   "...on Bot {login}}",
                   "query A(",
                   join(("\$x$i:[ID!]!" for i in eachindex(nodes)), ','),
                   "){",
                   join(("_$i:nodes(ids:\$x$i){...a}" for i in eachindex(nodes)), ' '),
                   "}")
    result = graphql(query, vars = vars, max_retries = 2)
    json = JSON3.read(result.Data)
    output = reduce(vcat, DataFrame(r for r in values(elem) if ~isnothing(r)) for elem in values(json.data))
    execute(conn, "BEGIN;")
    GHOST.load!(output, conn, "INSERT INTO $schema.test_usr VALUES(\$1,\$2,\$3) ON CONFLICT DO NOTHING;")
    execute(conn, "COMMIT;")
    nothing
end

for i in 1:lastindex(users)
    println(i)
    query_users(users[i])
    sleep(0.25)
end


text:
function record_camera_calibration()
    picsfolder = joinpath(tempdir(), constants["calibration"]["camera"]["picsfolder"])
    if !isdir(picsfolder)
        mkdir(picsfolder)
    else
        map(rm, readdir(picsfolder, join = true))
    end

    takeapicture(joinpath(picsfolder, "img%04d.png"), options = `--keypress --fullpreview`)
end

function assemble_camera_calibration()
    datafolder = constants["calibration"]["camera"]["datafolder"]
    file = joinpath(homedir(), datafolder, "calibration_matrices.jld2")

    @load file tform itform

    M = SVector{2, Float64}.(zip(tform[:x], tform[:y]))
    _tform = interpolate(M, BSpline(Linear()))
    tform = Base.splat(extrapolate(scale(_tform, tform[:cols], tform[:rows]), Interpolations.Flat()))
    M = SVector{2, Float64}.(zip(itform[:x], itform[:y]))
    _itform = interpolate(M, BSpline(Linear()))
    itform = Base.splat(extrapolate(scale(_itform, itform[:xs], itform[:ys]), Interpolations.Flat()))

    file = joinpath(homedir(), datafolder, "calibration.jld2")
    @save file tform itform
end

function save_mask()
    mktemp() do file, _
        takeapicture(file)
        img = FileIO.load(file)
        w = constants["raspioptions"]["width"]
        h = constants["raspioptions"]["height"]
        h2 = h ÷ 2
        w2 = w ÷ 2
        m = 10
        seeds = vcat(((CartesianIndex(h2 + i, w2 + j), 1) for i in -m:m for j in -m:m)..., ((CartesianIndex(i, j), 2) for i in vcat(1:m,h-m + 1:h) for j in vcat(1:m,w-m + 1:w))...)
        segments = seeded_region_growing(img, seeds)
        mask = labels_map(segments) .== 1
        file = joinpath(homedir(), "mask.png")
        FileIO.save(file, mask)
    end
    nothing
end

function textlayer!(ax::LAxis)
    pxa = lift(AbstractPlotting.zero_origin, ax.scene.px_area)
    Scene(ax.scene, pxa, raw = true, camera = campixel!)
end
function AbstractPlotting.annotations!(textlayer::Scene, ax::LAxis, texts, positions; kwargs...)
    positions = positions isa Observable ? positions : Observable(positions)
    screenpositions = lift(positions, ax.scene.camera.projectionview, ax.scene.camera.pixel_space) do positions, pv, pspace
        p4s = to_ndim.(Vec4f0, to_ndim.(Vec3f0, positions, 0.0), 1.0)
        p1m1s = [pv *  p for p in p4s]
        projected = [inv(pspace) * p1m1 for p1m1 in p1m1s]
        pdisplay = [Point2(p[1:2]...) for p in projected]
    end
    annotations!(textlayer, texts, screenpositions; kwargs...)
end

function horizon(θ, ps)
    rm = RotMatrix(θ)
    std([last(rm*p) for p in ps])
end

function get_ends(labels, tform)
    map(1:maximum(labels)) do i
        j = findall(==(i), labels)
        xy = SVector{2, Float64}.(reverse.(Tuple.(j)))
        μ = mean(xy)
        xy .-= Ref(μ)
        o = optimize(θ -> horizon(θ, xy), -π/2, π)
        θ = -o.minimizer
        rm = RotMatrix(-θ)
        xy2 = [rm*p for p in xy]
        l = -(reverse(extrema(first, xy2))...)
        p1 = μ .+ reverse(sincos(θ)).*l./2
        p2 = μ .+ reverse(sincos(θ + π)).*l./2
        p1 = SVector{2, Float64}(tform(p1))
        p2 = SVector{2, Float64}(tform(p2))
        Δ = p2 - p1
        (length = norm(Δ), ends = p1 => p2, θ = atan(reverse(Δ)...))
    end
end

function test_camera_calibration()
    mktemp() do imgfile, _
        takeapicture(imgfile)
        img = FileIO.load(imgfile)

        maskfile = joinpath(homedir(), "mask.png")
        maskimg = FileIO.load(maskfile)
        mask = maskimg .== zero(eltype(maskimg))

        bw = Gray.(img)
        tf = bw .< otsu_threshold(bw)
        markers = label_components(tf)
        for l in unique(markers)
            i = findall(==(l), markers)
            if any(mask[i])
                markers[i] .= 0
            end
        end
        labels = replace(markers, (u => i - 1 for (i, u) in enumerate(unique(markers)))...)

        calibrationfolder = constants["calibration"]["camera"]["datafolder"]
        file = joinpath(homedir(), calibrationfolder, "calibration.jld2")
        @load file tform itform

        lxy = get_ends(labels, tform)

        indices = ImageTransformations.autorange(img, tform)
        imgw = parent(warp(img, itform, indices))

        unit = constants["calibration"]["camera"]["board"]["unit"]
        scene, layout = layoutscene(0, resolution = size(imgw))
        ax = layout[1, 1] = LAxis(scene, aspect = DataAspect(), xlabel = "X ($unit)", ylabel = "Y ($unit)")
        image!(ax, indices..., imgw)
        textlayer = textlayer!(ax)
        annotations!(textlayer, ax, [string(round(i.length, digits = 2), " ", unit) for i in lxy], [mean(i.ends) for i in lxy], color = :yellow, align = (:center, :bottom),  rotation = [i.θ for i in lxy])
        linesegments!(ax, [i.ends for i in lxy], color = :cyan)
        tightlimits!(ax)
        handler(a,b) = scene
        JSServe.Application(handler, "0.0.0.0", constants["web"]["port"])
    end
end


text:
include("scaling.jl")
include("dimspecs.jl")
include("nointerp.jl")
include("withextrap.jl")


text:
module path

end

max_stars_repo_name: tmaone/xsh

max_issues_repo_name: turing-complete/xsh

max_forks_repo_name: turing-complete/xsh


text:
using MPI
using CLIMA
using CLIMA.Mesh.Topologies
using CLIMA.Mesh.Grids
using CLIMA.Mesh.Filters
using CLIMA.DGBalanceLawDiscretizations
using CLIMA.MPIStateArrays
using Printf
using LinearAlgebra
using Logging

@static if haspkg("CuArrays")
  using CUDAdrv
  using CUDAnative
  using CuArrays
  CuArrays.allowscalar(false)
  const ArrayTypes = (CuArray, )
else
  const ArrayTypes = (Array, )
end

using Test
function run(mpicomm, dim, ArrayType, Ne, FT)
  N = 3

  brickrange = ntuple(j->range(FT(-1); length=Ne[j]+1, stop=1), dim)
  topl = BrickTopology(mpicomm, brickrange, periodicity=ntuple(j->true, dim))

  grid = DiscontinuousSpectralElementGrid(topl,
                                          FloatType = FT,
                                          DeviceArray = ArrayType,
                                          polynomialorder = N,
                                         )

  spacedisc = DGBalanceLaw(grid = grid,
                           length_state_vector = 4,
                           flux! = (x...) -> (),
                           numerical_flux! = (x...) -> ())

  filter = CutoffFilter(grid, 2)

  # Legendre Polynomials
  l0(r) = 1
  l1(r) = r
  l2(r) = (3*r^2-1)/2
  l3(r) = (5*r^3-3r)/2

  low(x, y, z)  = l0(x) * l0(y) + 4 * l1(x) * l1(y) + 5 * l1(z) +
                  6 * l1(z) * l1(x)

  high(x, y, z) = l2(x) * l3(y) + l3(x) + l2(y) + l3(z) * l1(y)

  filtered_both(x, y, z) = high(x, y, z)

  filtered_vertical(x, y, z) = (dim == 2) ? l2(x) * l3(y) + l2(y) :
                                            l3(z) * l1(y)
  filtered_horizontal(x, y, z) = (dim == 2) ? l2(x) * l3(y) + l3(x) :
                                              l2(x) * l3(y) + l3(x) + l2(y)

  for horizontal = false:true
    for vertical = false:true

      if horizontal && vertical
        filtered = filtered_both
      elseif horizontal
        filtered = filtered_horizontal
      elseif vertical
        filtered = filtered_vertical
      else
        filtered = (x...) -> zero(x[1])
      end

      Q = MPIStateArray(spacedisc) do Q, x, y, z, _...
        @inbounds begin
          Q[1] = low(x, y, z) + high(x, y, z)
          Q[2] = low(x, y, z) + high(x, y, z)
          Q[3] = low(x, y, z) + high(x, y, z)
          Q[4] = low(x, y, z) + high(x, y, z)
        end
      end
      P = MPIStateArray(spacedisc) do P, x, y, z, _...
        @inbounds begin
          P[1] = low(x, y, z) + high(x, y, z) - filtered(x, y, z)
          P[2] = low(x, y, z) + high(x, y, z)
          P[3] = low(x, y, z) + high(x, y, z) - filtered(x, y, z)
          P[4] = low(x, y, z) + high(x, y, z)
        end
      end

      Filters.apply!(Q, (1,3), spacedisc.grid, filter;
                     horizontal=horizontal, vertical=vertical)

      @test Array(Q.data) ≈ Array(P.data)
    end
  end
end

let
  MPI.Initialized() || MPI.Init()

  mpicomm = MPI.COMM_WORLD
  ll = uppercase(get(ENV, "JULIA_LOG_LEVEL", "INFO"))
  loglevel = ll == "DEBUG" ? Logging.Debug :
             ll == "WARN"  ? Logging.Warn  :
             ll == "ERROR" ? Logging.Error : Logging.Info
  logger_stream = MPI.Comm_rank(mpicomm) == 0 ? stderr : devnull
  global_logger(ConsoleLogger(logger_stream, loglevel))
  @static if haspkg("CUDAnative")
    device!(MPI.Comm_rank(mpicomm) % length(devices()))
  end

  numelem = (1, 1, 1)
  lvls = 1

  @testset "$(@__FILE__)" for ArrayType in ArrayTypes
    for FT in (Float64, Float32)
      for dim = 2:3
        err = zeros(FT, lvls)
        for l = 1:lvls
          @info (ArrayType, FT, dim)
          run(mpicomm, dim, ArrayType, ntuple(j->2^(l-1) * numelem[j], dim),
              FT)
        end
      end
    end
  end
end

nothing

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
"""
plot_Trajectory(t::Trajectory; filtering=0, control=true)
"""
plot_Trajectory
@recipe function plot_Trajectory(t::Trajectory; filtering=0, control=true)
    layout --> (control ? (2,1) : 1)
    show --> false
    xlabel --> "Time"
    @series begin
        title --> "States"
        control && (subplot --> 1)
        filtering > 0 ? filt(ones(filtering),[filtering], t.x') : t.x'
    end
    control && @series begin
        title --> "Control signal"
        subplot --> 2
        filtering > 0 ? filt(ones(filtering),[filtering], t.u') : t.u'
    end
    delete!(plotattributes, :filtering)
    delete!(plotattributes, :control)
    nothing
end

@userplot EigvalPlot

"""
eigvalplot(models, traj; ds=10, cont=false, onlyat=0)
"""
eigvalplot
@recipe function eigvalplot(h::EigvalPlot; ds=10, cont=false, onlyat=0)
    ms = h.args[1]
    t = h.args[2]
    T = length(t)
    title --> "Eigenvalue spectrum"
    show --> false
    seriestype := :scatter
    @unpack xu,nu,nx = t
    lim = 1e-3
    truejac = length(h.args) >= 3
    ħ = 1
    inds = onlyat == 0 ? (1:ds:T) : onlyat
    cmap = colormap("Blues", T)
    cmapt = colormap("Reds", T÷ds+1)
    delete!(plotattributes, :onlyat)
    for evalpoint = inds
        J = jacobian(ms, xu[:,evalpoint])[1]
        e = eigvals(J[1:nx,1:nx])
        if cont
            e = log.(Complex.(e))./ħ
        end
        lim = max(maximum(imag.(e)), lim)
        @series begin
            c --> (onlyat == 0 ? cmap[evalpoint] : (:blue))
            real.(e), imag.(e)
        end
    end
    if truejac
        truejacfun = h.args[3]
        Jtrue = map(x->truejacfun(x...), (t.x[:,i],t.u[:,i]) for i=inds)
        for (i,J) in enumerate(Jtrue)
            e = eigvals(J[1:nx,1:nx])
            if cont
                e = log.(Complex.(e))./ħ
            end
            @series begin
                c --> (onlyat == 0 ? cmapt[i] : (:red))
                markersize --> 4
                real.(e), imag.(e)
            end
        end

    end
    delete!(plotattributes, :ds)
    phi = range(-π/4, stop=π/4, length=150)
    seriestype := :path
    if !cont
        ratio --> :equal
        @series begin
            legend --> false; color --> :black; linestyle --> :dash
            real.(exp.(phi.*im)), imag.(exp.(phi.*im))
        end
    else
        @series (legend := false; linestyle := :dash; color := :black; ([0,0],[-lim,lim]))
    end
    delete!(plotattributes, :cont)
    nothing
end

@userplot PredPlot
"""
predplot(models, traj; filtering=1)
"""
PredPlot
@recipe function predplot(h::PredPlot; filtering=1)
    ms = h.args[1]
    t = h.args[2]
    lab --> "Prediction"
    ms isa Vector && Flux.reset!.(getfield.(ms, :m))
    if ms isa Vector || ms isa ADDiffSystem
        r = filt(ones(filtering),[filtering], predict(ms,t.xu)[1]')
    else
        r = filt(ones(filtering),[filtering], LTVModelsBase.predict(ms,t)')
    end
    ms isa Vector && Flux.reset!.(getfield.(ms, :m))
    r
end

@userplot SimPlot
"""
simplot(models, traj; filtering=1)
"""
SimPlot
@recipe function simplot(h::SimPlot; filtering=1)
    ms = h.args[1]
    t = h.args[2]
    lab --> "Simulation"
    simfun = ms isa LTVModelsBase.AbstractModel ? LTVModelsBase.simulate : simulate
    ms isa Vector && Flux.reset!.(getfield.(ms, :m))
    r = filt(ones(filtering),[filtering], simfun(ms,t)')
    ms isa Vector && Flux.reset!.(getfield.(ms, :m))
    r
end

@userplot PredSimPlot
"""
predsimplot(models, traj; filtering=1)
"""
predsimplot
@recipe function predsimplot(h::PredSimPlot; filtering=1)
    ms = h.args[1]
    t = h.args[2]
    layout --> size(t.x,1)
    @series begin
        label --> "True"
        filt(ones(filtering),[filtering], t.x')
    end
    ms isa Vector && Flux.reset!.(getfield.(ms, :m))
    @series begin
        lab --> "Prediction"
        filt(ones(filtering),[filtering], predict(ms,t)[1]')
    end
    ms isa Vector && Flux.reset!.(getfield.(ms, :m))
    @series begin
        lab --> "Simulation"
        filt(ones(filtering),[filtering], simulate(ms,t)')
    end
    ms isa Vector && Flux.reset!.(getfield.(ms, :m))
    nothing
end

@userplot JacPlot
"""
jacplot(models, traj [, true_jac::Function]; ds=10, cont=false)
"""
jacplot
@recipe function jacplot(h::JacPlot; ds=10, cont=false)
    ms = h.args[1]
    t = h.args[2]
    truejac = length(h.args) >= 3
    if truejac
        truejacfun = h.args[3]
        Jtrue = map(x->truejacfun(x...), [(t.x[:,i],t.u[:,i]) for i=1:ds:size(t.x,2)])
    end
    show --> false
    if ms isa LTVModelsBase.LTVStateSpaceModel
        Jm = reshape(cat(2,ms.At, ms.Bt), :, length(ms))[:,1:ds:end]
    else
        Jm,Js = jacobians(ms, t, ds)
    end
    N = size(Jm,1)
    layout --> N
    legend := false
    for i = 1:N
        seriestype --> :line
        subplot := i
        truejac && @series begin
            linestyle --> :dash
            getindex.(Jtrue,i)
        end
        @series begin
            # Js == nothing || (ribbon --> 2Js[i,:])
            label --> "Estimated"
            fillalpha --> 0.4
            Jm[i,:]
        end
        @series [0 0]
    end
    delete!(plotattributes, :ds)
    delete!(plotattributes, :cont)
    nothing
end

@userplot MutRegPlot
"""
mutregplot(modeltrainer, trajs, truejacfun; useprior = true, showltv = true)
mutregplot(modeltrainer, truejacfun; useprior = true, showltv = true)
"""
mutregplot
@recipe function mutregplot(h::MutRegPlot; useprior = true, showltv = true)
    @assert length(h.args) >= 2 "Call with (mt::ModelTrainer, t::Trajectory, true_jacobian::(x,u)->J)\n or (mt::ModelTrainer, true_jacobian::(x,u)->J"
    mt         = h.args[1]
    manytrajs  = length(h.args) == 2
    if manytrajs
        t          = mt.trajs
        truejacfun = h.args[2]
    else
        t          = [h.args[2]]
        truejacfun = h.args[3]
    end

    errorhistory = map(mt.modelhistory) do ms
        map(t) do t
            showltv && (ltvmodel = KalmanModel(mt, t, ms, useprior=useprior))
            error_nn = 0.
            error_ltv = 0.
            for (i,xu) in enumerate(t)
                xi,ui = xu
                Jtrue = truejacfun(xi,ui)
                Jm, Js = jacobian(ms, [xi;ui])
                error_nn += eval_jac(Jm, Jtrue)
                showltv && (error_ltv += eval_jac([ltvmodel.At[:,:,i] ltvmodel.Bt[:,:,i]], Jtrue))
            end
            error_nn/length(t), error_ltv/length(t)
        end
    end
    title --> "Jacobian error"
    xlabel --> "Number of training sessions"
    ylims --> (0, Inf)
    N = length(errorhistory)
    errorhistory = [errorhistory[mi][ti] for mi in 1:N, ti in eachindex(errorhistory[1])]
    seriestype --> (showltv ? (:scatter) : (:bar))
    c1,c2 = :red, :blue
    @series (color --> c1;label := "NN"; ((1:N) .- 0.1, getindex.(errorhistory,1)))
    showltv && @series (color --> c2;label := "LTV"; ((1:N) .+ 0.1, getindex.(errorhistory,2)))
end

@userplot TracePlot
"""
traceplot(h::TracePlot)
"""
traceplot
@recipe function traceplot(h::TracePlot)
    trace         = h.args[1]
    if !(trace isa ValueHistory)
        trace = trace.trace
    end
    yscale --> :log10
    trace
end

# function display_modeltrainer(mt::ModelTrainer; kwargs...)
#     ms = mt.models
#     modeldict = OrderedDict("All models" => 0,
#                 ["Model "*string(i) => i for i = 1:length(mt.models)]...)
#
#     hist = length(mt.modelhistory)
#     ui = @manipulate for  t  = togglebuttons(1:length(mt.trajs), selected=1),
#                           mn = togglebuttons(modeldict, selected=1),
#                           modelversion = slider(1:hist,  value=hist, label="Model state"),
#                           f       = slider(1:100, label="Filtering", value=0),
#                           ds      = slider(1:100, value=4, label="Downsampling"),
#                           eigvals = [true,false],
#                           cont    = [true,false],
#                           trajplot = togglebuttons(OrderedDict( "Prediction"=>1,
#                                                                 "Simulation"=>2,
#                                                                 "LTV pred"=>3,
#                                                                 "LTV sim"=>4),
#                                                                 multiselect=true)
#
#         ms        = mt.modelhistory[modelversion]
#         modelinds = (mn <= 0 || mn > length(ms)) ? (1:length(ms)) : (mn:mn)
#         ms = ms[modelinds]
#         if eigvals
#             eigvalplot(ms, mt.trajs[t]; ds=ds, cont=cont, kwargs...)
#         else
#             fig = plot(mt.trajs[t]; filtering=f, lab="True", kwargs...)
#             1 ∈ trajplot && predplot!(ms,mt.trajs[t]; filtering=f, l=:dash, subplot=1)
#             2 ∈ trajplot && simplot!(ms,mt.trajs[t]; filtering=f, l=:dash, subplot=1)
#             3 ∈ trajplot && predplot!(KalmanModel(mt,mt.trajs[t]),mt.trajs[t]; filtering=f, l=:dash, subplot=1)
#             4 ∈ trajplot && simplot!(KalmanModel(mt,mt.trajs[t]),mt.trajs[t]; filtering=f, l=:dash, subplot=1)
#             fig
#         end
#     end
# end


text:
using VoronoiDelaunay, Gadfly

points = [
 Point2D(-1.0563841812533212, -1.4606363138997696)
 Point2D(0.06312982975128989, -0.48031801152366027)
 Point2D(0.1624918689993189, -0.19919450833195906)
 Point2D(-1.5293344878962758, -0.7657808444340142)
 Point2D(0.5319064220493406, 0.6107808132476504)
 Point2D(-0.3670342825169435, 0.8207427582546951)
 Point2D(-1.9797019290444364, -0.5066353099040788)
 Point2D(-1.5811584920674324, 1.0346409888830976)
 Point2D(1.2185165319349451, 1.4177209167374605)
 Point2D(-1.5991536318191626, -1.3063986775765466) ];

tess, ranges = DelaunayTessellation( points )
xy = getplotxy( delaunayedges( tess, ranges ) )
l1 = layer( x=getx.(points), y=gety.(points), Theme(default_color=colorant"orange"), Geom.point )
l2 = layer( x=xy[1], y=xy[2], Geom.path )
plot( l1, l2, Coord.cartesian(fixed=true) )


text:
using SmallFixedSizeSets
using Base.Test

using Combinatorics

function test_issubset0()
    Empty = SFSSFromOrderedVec(Int32[])
    # Size-0 subset, is subset
    A = SFSSFromOrderedVec(Int32[])
    @test SFSS_issubset(A, Empty)
    @test SFSS_issubset(A, SFSSFromOrderedVec([11]))
    @test SFSS_issubset(A, SFSSFromOrderedVec([11, 111]))
    @test SFSS_issubset(A, SFSSFromOrderedVec([11, 111, 1111]))
    @test SFSS_issubset(A, SFSSFromOrderedVec([11, 111, 1111, 11111]))
    @test SFSS_issubset(A, SFSSFromOrderedVec([11, 111, 1111, 11111, 111111]))
end

function test_issubset1()
    Empty = SFSSFromOrderedVec(Int32[])
    # Size-1 subset
    A = SFSSFromOrderedVec([11])
    # Size-1, size-0, smaller
    @test !SFSS_issubset(A, Empty)
    # Size-1, size-1, is subset
    @test SFSS_issubset(A, SFSSFromOrderedVec([11]))
    # Size-1, size-2, is subset
    @test SFSS_issubset(A, SFSSFromOrderedVec([11, 12]))
    @test SFSS_issubset(A, SFSSFromOrderedVec([10, 11]))
    # Size-1, size-3, is subset
    @test SFSS_issubset(A, SFSSFromOrderedVec([8, 9, 11]))
    @test SFSS_issubset(A, SFSSFromOrderedVec([9, 11, 13]))
    @test SFSS_issubset(A, SFSSFromOrderedVec([11, 12, 13]))
    # Size-1, size-4, is subset
    @test SFSS_issubset(A, SFSSFromOrderedVec([11, 12, 13, 14]))
    @test SFSS_issubset(A, SFSSFromOrderedVec([9, 11, 12, 13]))
    @test SFSS_issubset(A, SFSSFromOrderedVec([8, 9, 11, 13]))
    @test SFSS_issubset(A, SFSSFromOrderedVec([7, 8, 9, 11]))
    # Size-1, size-5, is subset
    @test SFSS_issubset(A, SFSSFromOrderedVec([11, 12, 13, 14, 15]))
    @test SFSS_issubset(A, SFSSFromOrderedVec([9, 11, 12, 13, 14]))
    @test SFSS_issubset(A, SFSSFromOrderedVec([8, 9, 11, 12, 13]))
    @test SFSS_issubset(A, SFSSFromOrderedVec([7, 8, 9, 11, 12]))
    @test SFSS_issubset(A, SFSSFromOrderedVec([6, 7, 8, 9, 11]))

    # Size 1, not subset
    @test !SFSS_issubset(A, SFSSFromOrderedVec([21]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([21, 22]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([21, 22, 23]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([21, 22, 23, 24]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([21, 22, 23, 24, 25]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([7]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([7, 22]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([7, 22, 23]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([7, 22, 23, 24]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([7, 22, 23, 24, 25]))
end

function test_issubset2()
    Empty = SFSSFromOrderedVec(Int32[])
    # Size-2 subset
    A = SFSSFromOrderedVec([11, 111])

    # Size-2, Size-2, is subset
    @test SFSS_issubset(A, SFSSFromOrderedVec([11, 111]))
    # Size-2, Size-3, is subset
    @test SFSS_issubset(A, SFSSFromOrderedVec([11, 111, 112]))
    @test SFSS_issubset(A, SFSSFromOrderedVec([9, 11, 111]))
    @test SFSS_issubset(A, SFSSFromOrderedVec([11, 109, 111]))
    # Size-2, Size-4, is subset
    @test SFSS_issubset(A, SFSSFromOrderedVec([11, 111, 1111, 11111]))
    @test SFSS_issubset(A, SFSSFromOrderedVec([11, 109, 111, 1111]))
    @test SFSS_issubset(A, SFSSFromOrderedVec([11, 108, 109, 111]))
    @test SFSS_issubset(A, SFSSFromOrderedVec([9, 11, 111, 1111]))
    @test SFSS_issubset(A, SFSSFromOrderedVec([9, 11, 109, 111]))
    @test SFSS_issubset(A, SFSSFromOrderedVec([8, 9, 11, 111]))
    # Size-2, Size-5, is subset
    @test SFSS_issubset(A, SFSSFromOrderedVec([11, 111, 1111, 11111, 111111]))
    @test SFSS_issubset(A, SFSSFromOrderedVec([11, 109, 111, 1111, 11111]))
    @test SFSS_issubset(A, SFSSFromOrderedVec([11, 108, 109, 111, 1111]))
    @test SFSS_issubset(A, SFSSFromOrderedVec([11, 107, 108, 109, 111]))
    @test SFSS_issubset(A, SFSSFromOrderedVec([9, 11, 111, 1111, 11111]))
    @test SFSS_issubset(A, SFSSFromOrderedVec([9, 11, 109, 111, 1111]))
    @test SFSS_issubset(A, SFSSFromOrderedVec([9, 11, 108, 109, 111]))
    @test SFSS_issubset(A, SFSSFromOrderedVec([8, 9, 11, 111, 1111]))
    @test SFSS_issubset(A, SFSSFromOrderedVec([8, 9, 11, 109, 111]))
    @test SFSS_issubset(A, SFSSFromOrderedVec([7, 8, 9, 11, 111]))

    # Size-2, smaller, not subset
    @test !SFSS_issubset(A, Empty)
    @test !SFSS_issubset(A, SFSSFromOrderedVec([50]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([111]))
    # Size-2, size-2, not subset
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 112]))
    # Size-2, size-3, not subset
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 111, 113]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 112, 113]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 112, 113]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 12, 111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 11, 112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 12, 112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 109, 111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 109, 112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 109, 112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 111, 1111, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 112, 1111, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 112, 1111, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 109, 111, 1111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 109, 112, 1111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 109, 112, 1111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 108, 109, 111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 108, 109, 112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 108, 109, 112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 12, 111, 1111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 11, 112, 1111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 12, 112, 1111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 12, 109, 111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 11, 109, 112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 12, 109, 112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([8, 9, 12, 111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([8, 9, 11, 112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([8, 9, 12, 112]))
    # Size-2, size-4, not subset
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 111, 1111, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 112, 1111, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 112, 1111, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 109, 111, 1111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 109, 112, 1111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 109, 112, 1111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 108, 109, 111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 108, 109, 112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 108, 109, 112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 12, 111, 1111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 11, 112, 1111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 12, 112, 1111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 12, 109, 111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 11, 109, 112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 12, 109, 112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([8, 9, 12, 111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([8, 9, 11, 112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([8, 9, 12, 112]))
    # Size-2, Size-5, not subset
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 111, 1111, 11111, 111111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 112, 1111, 11111, 111111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 112, 1111, 11111, 111111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 109, 111, 1111, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 109, 112, 1111, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 109, 112, 1111, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 108, 109, 111, 1111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 108, 109, 112, 1111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 108, 109, 112, 1111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 107, 108, 109, 111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 107, 108, 109, 112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 107, 108, 109, 112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 12, 111, 1111, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 11, 112, 1111, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 12, 112, 1111, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 12, 109, 111, 1111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 11, 109, 112, 1111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 12, 109, 112, 1111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 12, 108, 109, 111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 11, 108, 109, 112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 12, 108, 109, 112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([8, 9, 12, 111, 1111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([8, 9, 11, 112, 1111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([8, 9, 12, 112, 1111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([8, 9, 12, 109, 111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([8, 9, 11, 109, 112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([8, 9, 12, 109, 112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([7, 8, 9, 12, 111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([7, 8, 9, 11, 112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([7, 8, 9, 12, 112]))
end

function test_issubset3()
    Empty = SFSSFromOrderedVec(Int32[])
    A = SFSSFromOrderedVec([11, 111, 1111])
    # Size-3, Size-3, is subset
    @test SFSS_issubset(A, SFSSFromOrderedVec([11, 111, 1111]))
    # Size-3, Size-4, is subset
    @test SFSS_issubset(A, SFSSFromOrderedVec([11, 111, 1111, 11111]))
    @test SFSS_issubset(A, SFSSFromOrderedVec([9, 11, 111, 1111]))
    @test SFSS_issubset(A, SFSSFromOrderedVec([11, 109, 111, 1111]))
    @test SFSS_issubset(A, SFSSFromOrderedVec([11, 111, 1109, 1111]))
    @test SFSS_issubset(A, SFSSFromOrderedVec([11, 111, 1111, 1112]))
    # Size-3, Size-5, is subset
    @test SFSS_issubset(A, SFSSFromOrderedVec([11, 111, 1111, 11111, 111111]))
    @test SFSS_issubset(A, SFSSFromOrderedVec([11, 109, 111, 1111, 11111]))
    @test SFSS_issubset(A, SFSSFromOrderedVec([11, 109, 111, 1109, 1111]))
    @test SFSS_issubset(A, SFSSFromOrderedVec([11, 108, 109, 111, 1111]))
    @test SFSS_issubset(A, SFSSFromOrderedVec([9, 11, 111, 1111, 11111]))
    @test SFSS_issubset(A, SFSSFromOrderedVec([9, 11, 109, 111, 1111]))
    @test SFSS_issubset(A, SFSSFromOrderedVec([8, 9, 11, 111, 1111]))
    # Size-3 subset, smaller, not subset
    @test !SFSS_issubset(A, Empty)
    @test !SFSS_issubset(A, SFSSFromOrderedVec([50]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([50, 51]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 50]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([111, 150]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([50, 1111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([1111, 1150]))
    # Size-3, Size-3, not subset
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 111, 1111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 112, 1111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 111, 1112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 112, 1111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 112, 1112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 112, 1112]))
    # Size-3, Size-4, not subset
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 111, 1111, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 112, 1111, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 111, 1112, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 112, 1111, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 111, 1112, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 112, 1112, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 112, 1112, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 12, 111, 1111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 11, 112, 1111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 11, 111, 1112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 12, 112, 1111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 12, 111, 1112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 11, 112, 1112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 12, 112, 1112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 109, 111, 1111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 109, 112, 1111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 109, 111, 1112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 109, 112, 1111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 109, 111, 1112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 109, 112, 1112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 109, 112, 1112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 111, 1109, 1111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 112, 1109, 1111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 111, 1109, 1112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 112, 1109, 1111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 111, 1109, 1112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 112, 1109, 1112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 112, 1109, 1112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 111, 1111, 1113]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 112, 1111, 1113]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 111, 1112, 1113]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 112, 1111, 1113]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 111, 1112, 1113]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 112, 1112, 1113]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 112, 1112, 1113]))
    # Size-3, Size-5, not subset
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 111, 1111, 11111, 111111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 112, 1111, 11111, 111111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 111, 1112, 11111, 111111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 112, 1111, 11111, 111111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 111, 1112, 11111, 111111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 112, 1112, 11111, 111111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 112, 1112, 11111, 111111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 109, 111, 1111, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 109, 112, 1111, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 109, 111, 1112, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 109, 112, 1111, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 109, 111, 1112, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 109, 112, 1112, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 109, 112, 1112, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 109, 111, 1109, 1111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 109, 112, 1109, 1111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 109, 111, 1109, 1112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 109, 112, 1109, 1111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 109, 111, 1109, 1112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 109, 112, 1109, 1112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 109, 112, 1109, 1112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 108, 109, 111, 1111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 108, 109, 112, 1111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 108, 109, 111, 1112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 108, 109, 112, 1111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 108, 109, 111, 1112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 108, 109, 112, 1112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 108, 109, 112, 1112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 12, 111, 1111, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 11, 112, 1111, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 11, 111, 1112, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 12, 112, 1111, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 12, 111, 1112, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 11, 112, 1112, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 12, 112, 1112, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 12, 109, 111, 1111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 11, 109, 112, 1111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 11, 109, 111, 1112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 12, 109, 112, 1111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 12, 109, 111, 1112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 11, 109, 112, 1112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 12, 109, 112, 1112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([8, 9, 12, 111, 1111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([8, 9, 11, 112, 1111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([8, 9, 11, 111, 1112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([8, 9, 12, 112, 1111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([8, 9, 12, 111, 1112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([8, 9, 11, 112, 1112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([8, 9, 12, 112, 1112]))
end

function test_issubset4()
    Empty = SFSSFromOrderedVec(Int32[])
    A = SFSSFromOrderedVec([11, 111, 1111, 11111])

    # Size-4, Size-4, is subset
    @test SFSS_issubset(A, SFSSFromOrderedVec([11, 111, 1111, 11111]))
    # Size-4, Size-5, is subset
    @test SFSS_issubset(A, SFSSFromOrderedVec([11, 111, 1111, 11111, 111111]))
    @test SFSS_issubset(A, SFSSFromOrderedVec([9, 11, 111, 1111, 11111]))
    @test SFSS_issubset(A, SFSSFromOrderedVec([11, 109, 111, 1111, 11111]))
    @test SFSS_issubset(A, SFSSFromOrderedVec([11, 111, 1109, 1111, 11111]))
    @test SFSS_issubset(A, SFSSFromOrderedVec([11, 111, 1111, 11109, 11111]))
    @test SFSS_issubset(A, SFSSFromOrderedVec([11, 111, 1111, 11111, 11112]))

    # Size-4 subset, smaller, not subset
    @test !SFSS_issubset(A, Empty)
    @test !SFSS_issubset(A, SFSSFromOrderedVec([50]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([1111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([50, 51]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 50]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([111, 150]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([1111, 150]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([50, 1111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([1111, 1150]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([50, 51, 52]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 51, 52]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([1111, 51, 52]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([51, 1111, 11112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 111, 112]))
    # Size-4, Size-4, not subset
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 111, 1111, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 112, 1111, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 111, 1112, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 111, 1111, 11112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 112, 1111, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 111, 1112, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 111, 1111, 11112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 112, 1112, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 112, 1111, 11112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 111, 1112, 11112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 112, 1112, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 112, 1111, 11112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 111, 1112, 11112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 112, 1112, 11112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 112, 1112, 11112]))
    # Size-4, Size-5, not subset
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 111, 1111, 11111, 111111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 112, 1111, 11111, 111111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 111, 1112, 11111, 111111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 111, 1111, 11112, 111111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 112, 1111, 11111, 111111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 111, 1112, 11111, 111111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 111, 1111, 11112, 111111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 112, 1112, 11111, 111111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 112, 1111, 11112, 111111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 111, 1112, 11112, 111111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 112, 1112, 11111, 111111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 112, 1111, 11112, 111111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 111, 1112, 11112, 111111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 112, 1112, 11112, 111111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 112, 1112, 11112, 111111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 12, 111, 1111, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 11, 112, 1111, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 11, 111, 1112, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 11, 111, 1111, 11112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 12, 112, 1111, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 12, 111, 1112, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 12, 111, 1111, 11112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 11, 112, 1112, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 11, 112, 1111, 11112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 11, 111, 1112, 11112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 12, 112, 1112, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 12, 112, 1111, 11112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 12, 111, 1112, 11112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 11, 112, 1112, 11112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([9, 12, 112, 1112, 11112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 109, 111, 1111, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 109, 112, 1111, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 109, 111, 1112, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 109, 111, 1111, 11112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 109, 112, 1111, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 109, 111, 1112, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 109, 111, 1111, 11112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 109, 112, 1112, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 109, 112, 1111, 11112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 109, 111, 1112, 11112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 109, 112, 1112, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 109, 112, 1111, 11112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 109, 111, 1112, 11112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 109, 112, 1112, 11112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 109, 112, 1112, 11112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 111, 1109, 1111, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 112, 1109, 1111, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 111, 1109, 1112, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 111, 1109, 1111, 11112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 112, 1109, 1111, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 111, 1109, 1112, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 111, 1109, 1111, 11112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 112, 1109, 1112, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 112, 1109, 1111, 11112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 111, 1109, 1112, 11112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 112, 1109, 1112, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 112, 1109, 1111, 11112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 111, 1109, 1112, 11112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 112, 1109, 1112, 11112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 112, 1109, 1112, 11112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 111, 1111, 11109, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 112, 1111, 11109, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 111, 1112, 11109, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 111, 1111, 11109, 11112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 112, 1111, 11109, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 111, 1112, 11109, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 111, 1111, 11109, 11112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 112, 1112, 11109, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 112, 1111, 11109, 11112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 111, 1112, 11109, 11112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 112, 1112, 11109, 11111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 112, 1111, 11109, 11112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 111, 1112, 11109, 11112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 112, 1112, 11109, 11112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 112, 1112, 11109, 11112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 111, 1111, 11111, 11113]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 112, 1111, 11111, 11113]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 111, 1112, 11111, 11113]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 111, 1111, 11112, 11113]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 112, 1111, 11111, 11113]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 111, 1112, 11111, 11113]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 111, 1111, 11112, 11113]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 112, 1112, 11111, 11113]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 112, 1111, 11112, 11113]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 111, 1112, 11112, 11113]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 112, 1112, 11111, 11113]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 112, 1111, 11112, 11113]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 111, 1112, 11112, 11113]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 112, 1112, 11112, 11113]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 112, 1112, 11112, 11113]))
end

function test_issubset5()
    Empty = SFSSFromOrderedVec(Int32[])
    A = SFSSFromOrderedVec([11, 111, 1111, 11111, 111111])
    # Size-5, Size-5 is subset
    @test SFSS_issubset(A, SFSSFromOrderedVec([11, 111, 1111, 11111, 111111]))

    # Size-5 subset, smaller, not subset
    @test !SFSS_issubset(A, Empty)
    @test !SFSS_issubset(A, SFSSFromOrderedVec([50]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([1111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([50, 51]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 50]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([111, 150]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([1111, 150]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([50, 1111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([1111, 1150]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([50, 51, 52]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 51, 52]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([1111, 51, 52]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([51, 1111, 11112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 111, 112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 111, 150, 1111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 111, 1111, 1112]))
    # Size-5, Size-5, not subset
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 111, 1111, 11111, 111111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 112, 1111, 11111, 111111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 111, 1112, 11111, 111111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 111, 1111, 11112, 111111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 111, 1111, 11111, 111112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 112, 1111, 11111, 111111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 111, 1112, 11111, 111111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 111, 1111, 11112, 111111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 111, 1111, 11111, 111112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 112, 1112, 11111, 111111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 112, 1111, 11112, 111111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 112, 1111, 11111, 111112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 111, 1112, 11112, 111111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 111, 1112, 11111, 111112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 111, 1111, 11112, 111112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 112, 1112, 11111, 111111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 112, 1111, 11112, 111111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 112, 1111, 11111, 111112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 111, 1112, 11112, 111111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 111, 1112, 11111, 111112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 111, 1111, 11112, 111112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 112, 1112, 11112, 111111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 112, 1112, 11111, 111112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 112, 1111, 11112, 111112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 111, 1112, 11112, 111112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 112, 1112, 11112, 111111]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 112, 1112, 11111, 111112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 112, 1111, 11112, 111112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 111, 1112, 11112, 111112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([11, 112, 1112, 11112, 111112]))
    @test !SFSS_issubset(A, SFSSFromOrderedVec([12, 112, 1112, 11112, 111112]))
end

function test_intersect()
    function run_test(Av::Vector{Int64}, Bv::Vector{Int64}, Cv::Vector{Int64})
        A = SFSSFromOrderedVec(Av)
        B = SFSSFromOrderedVec(Bv)
        C = SFSSFromOrderedVec(Cv)
        @test SFSS_intersect(A, B) == C
        @test SFSS_intersect(B, A) == C
    end

    E = Int64[]

    run_test(E, [11], E)
    run_test(E, [11, 12], E)
    run_test(E, [11, 12, 13], E)
    run_test(E, [11, 12, 13, 14], E)
    run_test(E, [11, 12, 13, 14, 15], E)

    run_test([11], [12], E)
    run_test([11], [12, 13], E)
    run_test([11], [12, 13], E)
    run_test([11], [12, 13, 14], E)
    run_test([11], [12, 13, 14, 15], E)
    run_test([11], [12, 13, 14, 15, 16], E)

    run_test([11], [11], [11])
    run_test([11], [11, 13], [11])
    run_test([11], [11, 13, 14], [11])
    run_test([11], [11, 13, 14, 15], [11])
    run_test([11], [11, 13, 14, 15, 16], [11])
    run_test([13], [11, 13], [13])
    run_test([14], [11, 13, 14], [14])
    run_test([15], [11, 13, 14, 15], [15])
    run_test([16], [11, 13, 14, 15, 16], [16])

    run_test([11, 12], [21, 22], E)
    run_test([11, 12], [11, 22], [11])
    run_test([11, 12], [11, 12], [11, 12])
    run_test([11, 12], [9, 11], [11])
    run_test([11, 12], [12, 13], [12])

    run_test([11, 111], [11, 22, 23], [11])
    run_test([11, 111], [9, 11, 12], [11])
    run_test([11, 111], [9, 10, 11], [11])
    run_test([11, 111], [111, 112, 113], [111])
    run_test([11, 111], [110, 111, 113], [111])
    run_test([11, 111], [109, 110, 111], [111])

    run_test([11, 111], [21, 22, 23, 24], E)
    run_test([11, 111], [11, 22, 23, 24], [11])
    run_test([11, 111], [9, 11, 23, 24], [11])
    run_test([11, 111], [8, 9, 11, 24], [11])
    run_test([11, 111], [8, 9, 10, 11], [11])
    run_test([11, 111], [111, 122, 123, 124], [111])
    run_test([11, 111], [19, 111, 123, 124], [111])
    run_test([11, 111], [8, 9, 111, 124], [111])
    run_test([11, 111], [8, 9, 10, 111], [111])

    run_test([11, 111], [11, 111, 122, 133], [11, 111])
    run_test([11, 111], [11, 12, 111, 133], [11, 111])
    run_test([11, 111], [11, 12, 13, 111], [11, 111])
    run_test([11, 111], [9, 11, 111, 122], [11, 111])
    run_test([11, 111], [9, 11, 110, 111], [11, 111])
    run_test([11, 111], [9, 10, 11, 111], [11, 111])

    run_test([11, 111, 2222], [12, 112, 123], E)
    run_test([11, 111, 2222], [9, 11, 123], [11])
    run_test([11, 111, 2222], [9, 10, 11], [11])
    run_test([11, 111, 2222], [111, 112, 122], [111])
    run_test([11, 111, 2222], [10, 111, 122], [111])
    run_test([11, 111, 2222], [10, 110, 111], [111])
    run_test([11, 111, 2222], [2222, 2223, 2224], [2222])
    run_test([11, 111, 2222], [2221, 2222, 2224], [2222])
    run_test([11, 111, 2222], [2220, 2221, 2222], [2222])
    run_test([11, 111, 2222], [11, 111, 123], [11, 111])
    run_test([11, 111, 2222], [9, 11, 111], [11, 111])
    run_test([11, 111, 2222], [11, 12, 111], [11, 111])
    run_test([11, 111, 2222], [11, 112, 2222], [11, 2222])
    run_test([11, 111, 2222], [11, 2222, 2223], [11, 2222])
    run_test([11, 111, 2222], [9, 11, 2222], [11, 2222])
    run_test([11, 111, 2222], [11, 111, 2222], [11, 111, 2222])

    run_test([10, 20, 30], [100, 200, 300, 400], E)
    run_test([10, 20, 30], [10, 200, 300, 400], [10])
    run_test([10, 20, 30], [5, 20, 300, 400], [20])
    run_test([10, 20, 30], [11, 21, 30, 400], [30])
    run_test([10, 20, 30], [10, 20, 300, 400], [10, 20])
    run_test([10, 20, 30], [5, 6, 10, 20], [10, 20])
    run_test([10, 20, 30], [10, 30, 300, 400], [10, 30])
    run_test([10, 20, 30], [5, 6, 10, 30], [10, 30])
    run_test([10, 20, 30], [10, 20, 30, 40], [10, 20, 30])
    run_test([10, 20, 30], [10, 20, 25, 30], [10, 20, 30])
    run_test([10, 20, 30], [10, 15, 20, 30], [10, 20, 30])
    run_test([10, 20, 30], [5, 10, 20, 30], [10, 20, 30])
    run_test([10, 20, 30], [100, 200, 300, 400, 500], E)
    run_test([10, 20, 30], [10, 20, 30, 400, 500], [10, 20, 30])
    run_test([10, 20, 30], [5, 6, 10, 20, 30], [10, 20, 30])
    run_test([10, 20, 30], [5, 10, 16, 20, 30], [10, 20, 30])
    run_test([10, 20, 30], [5, 10, 20, 26, 30], [10, 20, 30])
    run_test([10, 20, 30], [5, 10, 20, 30, 36], [10, 20, 30])
    run_test([10, 20, 30], [5, 6, 11, 20, 30], [20, 30])
    run_test([10, 20, 30], [5, 11, 16, 20, 30], [20, 30])
    run_test([10, 20, 30], [5, 9, 20, 26, 30], [20, 30])
    run_test([10, 20, 30], [5, 9, 20, 30, 36], [20, 30])
    run_test([10, 20, 30], [5, 6, 10, 21, 30], [10, 30])
    run_test([10, 20, 30], [5, 10, 16, 21, 30], [10, 30])
    run_test([10, 20, 30], [5, 10, 21, 26, 30], [10, 30])
    run_test([10, 20, 30], [5, 10, 21, 30, 36], [10, 30])

    run_test([10, 20, 30, 40], [100, 200, 300, 400, 500], E)
    run_test([10, 20, 30, 40], [10, 20, 30, 40, 50], [10, 20, 30, 40])
    run_test([10, 20, 30, 40], [10, 20, 25, 30, 40], [10, 20, 30, 40])
    run_test([10, 20, 30, 40], [10, 12, 20, 30, 40], [10, 20, 30, 40])
    run_test([10, 20, 30, 40], [8, 10, 20, 30, 40], [10, 20, 30, 40])
    run_test([10, 20, 30, 40], [10, 21, 25, 30, 40], [10, 30, 40])
    run_test([10, 20, 30, 40], [10, 20, 25, 31, 40], [10, 20, 40])
    run_test([10, 20, 30, 40], [9, 20, 25, 30, 40], [20, 30, 40])
    run_test([10, 20, 30, 40], [10, 21, 25, 30, 41], [10, 30])
    run_test([10, 20, 30, 40], [10, 21, 25, 31, 40], [10, 40])
    run_test([10, 20, 30, 40], [9, 10, 21, 25, 40], [10, 40])
    run_test([10, 20, 30, 40], [9, 11, 20, 30, 41], [20, 30])
    run_test([10, 20, 30, 40], [9, 11, 20, 31, 41], [20])
    run_test([10, 20, 30, 40], [40, 41, 42, 43, 44], [40])

    run_test([10, 20, 30, 40, 50], [100, 200, 300, 400, 500], E)
    run_test([10, 20, 30, 40, 50], [10, 20, 30, 40, 50], [10, 20, 30, 40, 50])
    run_test([10, 20, 30, 40, 50], [10, 20, 30, 40, 51], [10, 20, 30, 40])
    run_test([10, 20, 30, 40, 50], [9, 10, 20, 30, 40], [10, 20, 30, 40])
    run_test([10, 20, 30, 40, 50], [10, 21, 30, 40, 50], [10, 30, 40, 50])
    run_test([10, 20, 30, 40, 50], [10, 30, 32, 40, 50], [10, 30, 40, 50])
    run_test([10, 20, 30, 40, 50], [10, 30, 40, 42, 50], [10, 30, 40, 50])
    run_test([10, 20, 30, 40, 50], [9, 10, 30, 40, 50], [10, 30, 40, 50])
    run_test([10, 20, 30, 40, 50], [10, 30, 40, 50, 60], [10, 30, 40, 50])
    run_test([10, 20, 30, 40, 50], [11, 20, 30, 40, 51], [20, 30, 40])
    run_test([10, 20, 30, 40, 50], [20, 21, 30, 40, 51], [20, 30, 40])
    run_test([10, 20, 30, 40, 50], [20, 30, 31, 40, 51], [20, 30, 40])
    run_test([10, 20, 30, 40, 50], [20, 30, 40, 41, 51], [20, 30, 40])
    run_test([10, 20, 30, 40, 50], [9, 20, 30, 40, 51], [20, 30, 40])
    run_test([10, 20, 30, 40, 50], [10, 21, 31, 40, 50], [10, 40, 50])
    run_test([10, 20, 30, 40, 50], [9, 10, 31, 40, 50], [10, 40, 50])
    run_test([10, 20, 30, 40, 50], [8, 9, 10, 40, 50], [10, 40, 50])
    run_test([10, 20, 30, 40, 50], [8, 9, 20, 33, 50], [20, 50])
    run_test([10, 20, 30, 40, 50], [8, 20, 21, 33, 50], [20, 50])
    run_test([10, 20, 30, 40, 50], [8, 20, 21, 50, 51], [20, 50])
    run_test([10, 20, 30, 40, 50], [8, 9, 20, 50, 51], [20, 50])
    run_test([10, 20, 30, 40, 50], [8, 9, 20, 53, 55], [20])
    run_test([10, 20, 30, 40, 50], [8, 20, 51, 53, 55], [20])
    run_test([10, 20, 30, 40, 50], [10, 21, 51, 53, 55], [10])
    run_test([10, 20, 30, 40, 50], [8, 9, 10, 53, 55], [10])
    run_test([10, 20, 30, 40, 50], [2, 3, 4, 5, 10], [10])
    run_test([10, 20, 30, 40, 50], [50, 51, 52, 53, 54], [50])
    run_test([10, 20, 30, 40, 50], [9, 50, 52, 53, 54], [50])
    run_test([10, 20, 30, 40, 50], [9, 11, 50, 52, 53], [50])
    run_test([10, 20, 30, 40, 50], [9, 11, 49, 50, 53], [50])
    run_test([10, 20, 30, 40, 50], [9, 11, 48, 49, 50], [50])
end

function test_setdiff()
    function run_test(Av::Vector{Int64}, Bv::Vector{Int64}, Cv::Vector{Int64})
        A = SFSSFromOrderedVec(Av)
        B = SFSSFromOrderedVec(Bv)
        C = SFSSFromOrderedVec(Cv)
        @test SFSS_setdiff(A, B) == C
    end

    function run_from_base(base::Vector{Int64})
        for r = 1:length(base)
            for inds in combinations(1:length(base), r)
                keep = ones(Bool, length(base))
                keep[inds] = false
                start = base[inds]
                result = base[keep]
                for additional = 0:(5 - length(start))
                    rest = [60 + j for j = 1:additional]
                    other_set = [start; rest]
                    run_test(base, other_set, result)
                end
            end
        end
    end
    run_from_base([10, 20, 30, 40, 50])
    run_from_base([10, 20, 30, 40])
    run_from_base([10, 20, 30])
    run_from_base([10, 20])
    run_from_base([10])
    run_from_base(Int64[])
end

function test_union()
    function run_test(Av::Vector{Int64}, Bv::Vector{Int64})
        A = SFSSFromOrderedVec(Av)
        B = SFSSFromOrderedVec(Bv)
        C = SFSSFromVec(union(Av, Bv))
        @test SFSS_union(A, B) == C
        @test SFSS_union(B, A) == C
    end

    S = [10, 20, 30, 40, 50]
    nS = length(S)
    for r1 = 0:nS, r2 = 0:nS
        for c1 in combinations(S, r1)
            for c2 in combinations(S, r2)
                run_test(c1, c2)
            end
        end
    end
end

function test_random(num_trials::Int64)
    function run_test(A::Set{Int64}, B::Set{Int64},
                      sA::SmallFixedSizeSet, sB::SmallFixedSizeSet)
        @test (A ⊆ B) == SFSS_issubset(sA, sB)
        C = A ∩ B
        sC = SFSS_intersect(sA, sB)
        @test SFSSFromSet(C) == sC
        D = setdiff(A, B)
        sD = SFSS_setdiff(sA, sB)
        @test SFSSFromSet(D) == sD
        U = C ∪ D
        sU = SFSSFromSet(U)
        @test SFSS_union(sC, sD) == sU
        @test SFSS_union(sC, sA) == sA
        @test SFSS_union(sA, sC) == sA
        @test SFSS_union(sC, sB) == sB
        @test SFSS_union(sB, sC) == sB
    end

    base = collect(1:15)
    for _ in 1:num_trials
        A = Set{Int64}(rand(base, rand(0:5)))
        B = Set{Int64}(rand(base, rand(0:5)))
        sA = SFSSFromSet(A)
        sB = SFSSFromSet(B)
        run_test(A, B, sA, sB)
    end
end

function test_all()
    test_issubset0()
    test_issubset1()
    test_issubset2()
    test_issubset3()
    test_issubset4()
    test_issubset5()
    test_intersect()
    test_setdiff()
    test_union()
    test_random(5000)
end

test_all()


text:
using Documenter
using NCTSSOS

makedocs(
sitename = "NCTSSOS",
pages = ["Home" => "index.md",
         "Noncommutative Polynomial Optimization" => "ncpop.md"],
modules = [NCTSSOS],
format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true")
)

deploydocs(
    repo = "github.com/wangjie212/NCTSSOS.git"
)