text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Tests related to TEME transformations.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
# ==============================================================================
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# File: ./src/transformations/teme/teme.jl
# ========================================

# Functions r_teme_to_tod and r_tod_to_teme
# -----------------------------------------

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-15: Performing IAU-76/FK5 reduction.
#
# According to this example and Table 3-6, using:
#
#   JD_TT = 2453101.828154745
#   δΔϵ   = -0.003875"
#   δΔψ   = -0.052195"
#   r_teme = 5094.18016210   i + 6127.64465950   j + 6380.34453270   k [km]
#   v_teme =   -4.7461314870 i +    0.7858180410 j +    5.5319312880 k [km/s]
#
# one gets the following data:
#
#   r_tod = 5094.51620300   i + 6127.36527840   j + 6380.34453270   k [km]
#   v_tod =   -4.7460883850 i +    0.7860783240 j +    5.5319312880 k [km/s]
#
# Furthermore, using:
#
#   JD_TT = 2453101.828154745
#   δΔϵ   = 0"
#   δΔψ   = 0"
#   r_teme = 5094.18016210   i + 6127.64465950   j + 6380.34453270   k [km]
#   v_teme =   -4.7461314870 i +    0.7858180410 j +    5.5319312880 k [km/s]
#
# one gets the following data:
#
#   r_tod = 5094.51478040   i + 6127.36646120   j + 6380.34453270   k [km]
#   v_tod =   -4.7460885670 i +    0.7860772220 j +    5.5319312880 k [km/s]
#
################################################################################

@testset "Functions r_teme_to_tod and r_tod_to_teme" begin
    JD_TT   = 2453101.828154745

    ## r_teme_to_tod
    ## =============

    ## -------------------------------------------------------------------------
    ##                             First Test
    ## -------------------------------------------------------------------------

    r_teme = [5094.18016210; 6127.64465950; 6380.34453270]
    v_teme = [-4.7461314870; 0.7858180410; 5.5319312880]

    ## DCM
    ## ---

    D_TOD_TEME = r_teme_to_tod(
        JD_TT,
        -0.003875 * π / (180 * 3600),
        -0.052195 * π / (180 * 3600)
    )

    r_tod = D_TOD_TEME*r_teme
    v_tod = D_TOD_TEME*v_teme

    @test r_tod[1] ≈ +5094.51620300 atol=1e-7
    @test r_tod[2] ≈ +6127.36527840 atol=1e-7
    @test r_tod[3] ≈ +6380.34453270 atol=1e-7

    @test v_tod[1] ≈ -4.7460883850  atol=1e-9
    @test v_tod[2] ≈ +0.7860783240  atol=1e-9
    @test v_tod[3] ≈ +5.5319312880  atol=1e-9

    ## Quaternion
    ## ----------

    q_TOD_TEME = r_teme_to_tod(
        Quaternion,
        JD_TT,
        -0.003875 * π / (180 * 3600),
        -0.052195 * π / (180 * 3600)
    )

    r_tod = vect(conj(q_TOD_TEME) * r_teme * q_TOD_TEME)
    v_tod = vect(conj(q_TOD_TEME) * v_teme * q_TOD_TEME)

    @test r_tod[1] ≈ +5094.51620300 atol=1e-7
    @test r_tod[2] ≈ +6127.36527840 atol=1e-7
    @test r_tod[3] ≈ +6380.34453270 atol=1e-7

    @test v_tod[1] ≈ -4.7460883850  atol=1e-9
    @test v_tod[2] ≈ +0.7860783240  atol=1e-9
    @test v_tod[3] ≈ +5.5319312880  atol=1e-9

    ## -------------------------------------------------------------------------
    ##                            Second Test
    ## -------------------------------------------------------------------------

    r_teme = [5094.18016210; 6127.64465950; 6380.34453270]
    v_teme = [-4.7461314870; 0.7858180410; 5.5319312880]

    ## DCM
    ## ---

    D_TOD_TEME = r_teme_to_tod(JD_TT)

    r_tod = D_TOD_TEME * r_teme
    v_tod = D_TOD_TEME * v_teme

    @test r_tod[1] ≈ +5094.51478040 atol=1e-7
    @test r_tod[2] ≈ +6127.36646120 atol=1e-7
    @test r_tod[3] ≈ +6380.34453270 atol=1e-7

    @test v_tod[1] ≈ -4.7460885670  atol=1e-9
    @test v_tod[2] ≈ +0.7860772220  atol=1e-9
    @test v_tod[3] ≈ +5.5319312880  atol=1e-9

    ## Quaternion
    ## ----------

    q_TOD_TEME = r_teme_to_tod(Quaternion, JD_TT)

    r_tod = vect(conj(q_TOD_TEME) * r_teme * q_TOD_TEME)
    v_tod = vect(conj(q_TOD_TEME) * v_teme * q_TOD_TEME)

    @test r_tod[1] ≈ +5094.51478040 atol=1e-7
    @test r_tod[2] ≈ +6127.36646120 atol=1e-7
    @test r_tod[3] ≈ +6380.34453270 atol=1e-7

    @test v_tod[1] ≈ -4.7460885670  atol=1e-9
    @test v_tod[2] ≈ +0.7860772220  atol=1e-9
    @test v_tod[3] ≈ +5.5319312880  atol=1e-9

    ## r_tod_to_teme
    ## ==========

    ## -------------------------------------------------------------------------
    ##                             First Test
    ## -------------------------------------------------------------------------

    r_tod = [5094.51620300; 6127.36527840; 6380.34453270]
    v_tod = [-4.7460883850; 0.7860783240; 5.5319312880]

    ## DCM
    ## ---

    D_TEME_TOD = r_tod_to_teme(
        JD_TT,
        -0.003875 * π / (180 * 3600),
        -0.052195 * π / (180 * 3600)
    )

    r_teme = D_TEME_TOD * r_tod
    v_teme = D_TEME_TOD * v_tod

    @test r_teme[1] ≈ +5094.18016210 atol=1e-7
    @test r_teme[2] ≈ +6127.64465950 atol=1e-7
    @test r_teme[3] ≈ +6380.34453270 atol=1e-7

    @test v_teme[1] ≈ -4.7461314870  atol=1e-9
    @test v_teme[2] ≈ +0.7858180410  atol=1e-9
    @test v_teme[3] ≈ +5.5319312880  atol=1e-9

    ## Quaternion
    ## ----------

    q_TEME_TOD = r_tod_to_teme(
        Quaternion,
        JD_TT,
        -0.003875 * π / (180 * 3600),
        -0.052195 * π / (180 * 3600)
    )

    r_teme = vect(conj(q_TEME_TOD) * r_tod * q_TEME_TOD)
    v_teme = vect(conj(q_TEME_TOD) * v_tod * q_TEME_TOD)

    @test r_teme[1] ≈ +5094.18016210 atol=1e-7
    @test r_teme[2] ≈ +6127.64465950 atol=1e-7
    @test r_teme[3] ≈ +6380.34453270 atol=1e-7

    @test v_teme[1] ≈ -4.7461314870  atol=1e-9
    @test v_teme[2] ≈ +0.7858180410  atol=1e-9
    @test v_teme[3] ≈ +5.5319312880  atol=1e-9

    ## -------------------------------------------------------------------------
    ##                            Second Test
    ## -------------------------------------------------------------------------

    r_tod = [5094.51478040; 6127.36646120; 6380.34453270]
    v_tod = [-4.7460885670; 0.7860772220; 5.5319312880]

    ## DCM
    ## ---

    D_TEME_TOD = r_tod_to_teme(JD_TT)

    r_teme = D_TEME_TOD*r_tod
    v_teme = D_TEME_TOD*v_tod

    @test r_teme[1] ≈ +5094.18016210 atol=1e-7
    @test r_teme[2] ≈ +6127.64465950 atol=1e-7
    @test r_teme[3] ≈ +6380.34453270 atol=1e-7

    @test v_teme[1] ≈ -4.7461314870  atol=1e-9
    @test v_teme[2] ≈ +0.7858180410  atol=1e-9
    @test v_teme[3] ≈ +5.5319312880  atol=1e-9

    ## Quaternion
    ## ----------

    q_TEME_TOD = r_tod_to_teme(Quaternion, JD_TT)

    r_teme = vect(conj(q_TEME_TOD) * r_tod * q_TEME_TOD)
    v_teme = vect(conj(q_TEME_TOD) * v_tod * q_TEME_TOD)

    @test r_teme[1] ≈ +5094.18016210 atol=1e-7
    @test r_teme[2] ≈ +6127.64465950 atol=1e-7
    @test r_teme[3] ≈ +6380.34453270 atol=1e-7

    @test v_teme[1] ≈ -4.7461314870  atol=1e-9
    @test v_teme[2] ≈ +0.7858180410  atol=1e-9
    @test v_teme[3] ≈ +5.5319312880  atol=1e-9
end

# Functions r_teme_to_mod and r_mod_to_teme
# -----------------------------------------

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-15: Performing IAU-76/FK5 reduction.
#
# According to this example and Table 3-6, using:
#
#   JD_TT = 2453101.828154745
#   δΔϵ   = -0.003875"
#   δΔψ   = -0.052195"
#   r_teme = 5094.18016210   i + 6127.64465950   j + 6380.34453270   k [km]
#   v_teme =   -4.7461314870 i +    0.7858180410 j +    5.5319312880 k [km/s]
#
# one gets the following data:
#
#   r_mod = 5094.02837450   i + 6127.87081640   j + 6380.24851640   k [km]
#   v_mod =   -4.7462630520 i +    0.7860140450 j +    5.5317905620 k [km/s]
#
# Furthermore, using:
#
#   JD_TT = 2453101.828154745
#   δΔϵ   = 0"
#   δΔψ   = 0"
#   r_teme = 5094.18016210   i + 6127.64465950   j + 6380.34453270   k [km]
#   v_teme =   -4.7461314870 i +    0.7858180410 j +    5.5319312880 k [km/s]
#
# one gets the following data:
#
#   r_mod = 5094.02901670   i + 6127.87093630   j + 6380.24788850   k [km]
#   v_mod =   -4.7462624950 i +    0.7860141490 j +    5.5317910250 k [km/s]
#
################################################################################

@testset "Functions r_teme_to_mod and r_mod_to_teme" begin
    JD_TT   = 2453101.828154745

    ## r_teme_to_mod
    ## ==========

    ## -------------------------------------------------------------------------
    ##                             First Test
    ## -------------------------------------------------------------------------

    r_teme = [5094.18016210; 6127.64465950; 6380.34453270]
    v_teme = [-4.7461314870; 0.7858180410; 5.5319312880]

    ## DCM
    ## ---

    D_MOD_TEME = r_teme_to_mod(
        JD_TT,
        -0.003875 * π / (180 * 3600),
        -0.052195 * π / (180 * 3600)
    )

    r_mod = D_MOD_TEME * r_teme
    v_mod = D_MOD_TEME * v_teme

    @test r_mod[1] ≈ +5094.02837450 atol=1e-7
    @test r_mod[2] ≈ +6127.87081640 atol=1e-7
    @test r_mod[3] ≈ +6380.24851640 atol=1e-7

    @test v_mod[1] ≈ -4.7462630520  atol=1e-9
    @test v_mod[2] ≈ +0.7860140450  atol=1e-9
    @test v_mod[3] ≈ +5.5317905620  atol=1e-9

    ## Quaternion
    ## ----------

    q_MOD_TEME = r_teme_to_mod(
        Quaternion,
        JD_TT,
        -0.003875 * π / (180*3600),
        -0.052195 * π / (180*3600)
    )

    r_mod = vect(conj(q_MOD_TEME) * r_teme * q_MOD_TEME)
    v_mod = vect(conj(q_MOD_TEME) * v_teme * q_MOD_TEME)

    @test r_mod[1] ≈ +5094.02837450 atol=1e-7
    @test r_mod[2] ≈ +6127.87081640 atol=1e-7
    @test r_mod[3] ≈ +6380.24851640 atol=1e-7

    @test v_mod[1] ≈ -4.7462630520  atol=1e-9
    @test v_mod[2] ≈ +0.7860140450  atol=1e-9
    @test v_mod[3] ≈ +5.5317905620  atol=1e-9

    ## -------------------------------------------------------------------------
    ##                            Second Test
    ## -------------------------------------------------------------------------

    r_teme = [5094.18016210; 6127.64465950; 6380.34453270]
    v_teme = [-4.7461314870; 0.7858180410; 5.5319312880]

    ## DCM
    ## ---

    D_MOD_TEME = r_teme_to_mod(JD_TT)

    r_mod = D_MOD_TEME * r_teme
    v_mod = D_MOD_TEME * v_teme

    @test r_mod[1] ≈ +5094.02901670 atol=1e-7
    @test r_mod[2] ≈ +6127.87093630 atol=1e-7
    @test r_mod[3] ≈ +6380.24788850 atol=1e-7

    @test v_mod[1] ≈ -4.7462624950  atol=1e-9
    @test v_mod[2] ≈ +0.7860141490  atol=1e-9
    @test v_mod[3] ≈ +5.5317910250  atol=1e-9

    ## Quaternion
    ## ----------

    q_MOD_TEME = r_teme_to_mod(Quaternion, JD_TT)

    r_mod = vect(conj(q_MOD_TEME) * r_teme * q_MOD_TEME)
    v_mod = vect(conj(q_MOD_TEME) * v_teme * q_MOD_TEME)

    @test r_mod[1] ≈ +5094.02901670 atol=1e-7
    @test r_mod[2] ≈ +6127.87093630 atol=1e-7
    @test r_mod[3] ≈ +6380.24788850 atol=1e-7

    @test v_mod[1] ≈ -4.7462624950  atol=1e-9
    @test v_mod[2] ≈ +0.7860141490  atol=1e-9
    @test v_mod[3] ≈ +5.5317910250  atol=1e-9

    ## r_mod_to_teme
    ## ==========

    ## -------------------------------------------------------------------------
    ##                             First Test
    ## -------------------------------------------------------------------------

    r_mod = [5094.02837450; 6127.87081640; 6380.24851640]
    v_mod = [-4.7462630520; 0.7860140450; 5.5317905620]

    ## DCM
    ## ---

    D_TEME_MOD = r_mod_to_teme(
        JD_TT,
        -0.003875 * π / (180 * 3600),
        -0.052195 * π / (180 * 3600)
    )

    r_teme = D_TEME_MOD * r_mod
    v_teme = D_TEME_MOD * v_mod

    @test r_teme[1] ≈ +5094.18016210 atol=1e-7
    @test r_teme[2] ≈ +6127.64465950 atol=1e-7
    @test r_teme[3] ≈ +6380.34453270 atol=1e-7

    @test v_teme[1] ≈ -4.7461314870  atol=1e-9
    @test v_teme[2] ≈ +0.7858180410  atol=1e-9
    @test v_teme[3] ≈ +5.5319312880  atol=1e-9

    ## Quaternion
    ## ----------

    q_TEME_MOD = r_mod_to_teme(
        Quaternion,
        JD_TT,
        -0.003875 * π / (180 * 3600),
        -0.052195 * π / (180 * 3600)
    )

    r_teme = vect(conj(q_TEME_MOD) * r_mod * q_TEME_MOD)
    v_teme = vect(conj(q_TEME_MOD) * v_mod * q_TEME_MOD)

    @test r_teme[1] ≈ +5094.18016210 atol=1e-7
    @test r_teme[2] ≈ +6127.64465950 atol=1e-7
    @test r_teme[3] ≈ +6380.34453270 atol=1e-7

    @test v_teme[1] ≈ -4.7461314870  atol=1e-9
    @test v_teme[2] ≈ +0.7858180410  atol=1e-9
    @test v_teme[3] ≈ +5.5319312880  atol=1e-9

    ## -------------------------------------------------------------------------
    ##                            Second Test
    ## -------------------------------------------------------------------------

    r_mod = [5094.02901670; 6127.87093630; 6380.24788850]
    v_mod = [-4.7462624950; 0.7860141490; 5.5317910250]

    ## DCM
    ## ---

    D_TEME_MOD = r_mod_to_teme(JD_TT)

    r_teme = D_TEME_MOD * r_mod
    v_teme = D_TEME_MOD * v_mod

    @test r_teme[1] ≈ +5094.18016210 atol=1e-7
    @test r_teme[2] ≈ +6127.64465950 atol=1e-7
    @test r_teme[3] ≈ +6380.34453270 atol=1e-7

    @test v_teme[1] ≈ -4.7461314870  atol=1e-9
    @test v_teme[2] ≈ +0.7858180410  atol=1e-9
    @test v_teme[3] ≈ +5.5319312880  atol=1e-9

    ## Quaternion
    ## ----------

    q_TEME_MOD = r_mod_to_teme(Quaternion, JD_TT)

    r_teme = vect(conj(q_TEME_MOD) * r_mod * q_TEME_MOD)
    v_teme = vect(conj(q_TEME_MOD) * v_mod * q_TEME_MOD)

    @test r_teme[1] ≈ +5094.18016210 atol=1e-7
    @test r_teme[2] ≈ +6127.64465950 atol=1e-7
    @test r_teme[3] ≈ +6380.34453270 atol=1e-7

    @test v_teme[1] ≈ -4.7461314870  atol=1e-9
    @test v_teme[2] ≈ +0.7858180410  atol=1e-9
    @test v_teme[3] ≈ +5.5319312880  atol=1e-9
end

# Functions r_teme_to_gcrf and r_gcrf_to_teme
# -------------------------------------------

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-15: Performing IAU-76/FK5 reduction.
#
# According to this example and Table 3-6, using:
#
#   JD_TT = 2453101.828154745
#   δΔϵ   = -0.003875"
#   δΔψ   = -0.052195"
#   r_teme = 5094.18016210   i + 6127.64465950   j + 6380.34453270   k [km]
#   v_teme =   -4.7461314870 i +    0.7858180410 j +    5.5319312880 k [km/s]
#
# one gets the following data:
#
#   r_gcrf = 5102.50895790  i + 6123.01140070   j + 6378.13692820   k [km]
#   v_gcrf =  -4.7432201570 i + 0.7905364970    j + 5.5337557270    k [km/s]
#
# Furthermore, using:
#
#   JD_TT = 2453101.828154745
#   δΔϵ   = 0"
#   δΔψ   = 0"
#   r_teme = 5094.18016210   i + 6127.64465950   j + 6380.34453270   k [km]
#   v_teme =   -4.7461314870 i +    0.7858180410 j +    5.5319312880 k [km/s]
#
# one gets the following data:
#
#   r_j2000 = 5102.50960000  i + 6123.01152000   j + 6378.13630000   k [km]
#   v_j2000 =  -4.7432196000 i + 0.7905366000    j + 5.5337561900    k [km/s]
#
################################################################################

@testset "Functions r_teme_to_gcrf and r_gcrf_to_teme" begin
    JD_TT   = 2453101.828154745

    ## r_teme_to_gcrf
    ## ===========

    ## -------------------------------------------------------------------------
    ##                             First Test
    ## -------------------------------------------------------------------------

    r_teme = [5094.18016210; 6127.64465950; 6380.34453270]
    v_teme = [-4.7461314870; 0.7858180410; 5.5319312880]

    ## DCM
    ## ---

    D_GCRF_TEME = r_teme_to_gcrf(
        JD_TT,
        -0.003875 * π / (180 * 3600),
        -0.052195 * π / (180 * 3600)
    )

    r_gcrf = D_GCRF_TEME * r_teme
    v_gcrf = D_GCRF_TEME * v_teme

    @test r_gcrf[1] ≈ +5102.50895790 atol=1e-7
    @test r_gcrf[2] ≈ +6123.01140070 atol=1e-7
    @test r_gcrf[3] ≈ +6378.13692820 atol=1e-7

    @test v_gcrf[1] ≈ -4.7432201570  atol=1e-9
    @test v_gcrf[2] ≈ +0.7905364970  atol=1e-9
    @test v_gcrf[3] ≈ +5.5337557270  atol=1e-9

    ## Quaternion
    ## ----------

    q_GCRF_TEME = r_teme_to_gcrf(
        Quaternion,
        JD_TT,
        -0.003875 * π / (180 * 3600),
        -0.052195 * π / (180 * 3600)
    )

    r_gcrf = vect(conj(q_GCRF_TEME) * r_teme * q_GCRF_TEME)
    v_gcrf = vect(conj(q_GCRF_TEME) * v_teme * q_GCRF_TEME)

    @test r_gcrf[1] ≈ +5102.50895790 atol=1e-7
    @test r_gcrf[2] ≈ +6123.01140070 atol=1e-7
    @test r_gcrf[3] ≈ +6378.13692820 atol=1e-7

    @test v_gcrf[1] ≈ -4.7432201570  atol=1e-9
    @test v_gcrf[2] ≈ +0.7905364970  atol=1e-9
    @test v_gcrf[3] ≈ +5.5337557270  atol=1e-9

    ## -------------------------------------------------------------------------
    ##                            Second Test
    ## -------------------------------------------------------------------------

    r_teme = [5094.18016210; 6127.64465950; 6380.34453270]
    v_teme = [-4.7461314870; 0.7858180410; 5.5319312880]

    ## DCM
    ## ---

    D_J2000_TEME = r_teme_to_gcrf(JD_TT)

    r_j2000 = D_J2000_TEME * r_teme
    v_j2000 = D_J2000_TEME * v_teme

    @test r_j2000[1] ≈ +5102.50960000 atol=1e-7
    @test r_j2000[2] ≈ +6123.01152000 atol=1e-7
    @test r_j2000[3] ≈ +6378.13630000 atol=1e-7

    @test v_j2000[1] ≈ -4.7432196000  atol=1e-9
    @test v_j2000[2] ≈ +0.7905366000  atol=1e-9
    @test v_j2000[3] ≈ +5.5337561900  atol=1e-9

    ## Quaternion
    ## ----------

    q_J2000_TEME = r_teme_to_gcrf(Quaternion, JD_TT)

    r_j2000 = vect(conj(q_J2000_TEME) * r_teme * q_J2000_TEME)
    v_j2000 = vect(conj(q_J2000_TEME) * v_teme * q_J2000_TEME)

    @test r_j2000[1] ≈ +5102.50960000 atol=1e-7
    @test r_j2000[2] ≈ +6123.01152000 atol=1e-7
    @test r_j2000[3] ≈ +6378.13630000 atol=1e-7

    @test v_j2000[1] ≈ -4.7432196000  atol=1e-9
    @test v_j2000[2] ≈ +0.7905366000  atol=1e-9
    @test v_j2000[3] ≈ +5.5337561900  atol=1e-9

    ## r_gcrf_to_teme
    ## ===========

    ## -------------------------------------------------------------------------
    ##                             First Test
    ## -------------------------------------------------------------------------

    r_gcrf = [5102.50895790; 6123.01140070; 6378.13692820]
    v_gcrf = [-4.7432201570; 0.7905364970; 5.5337557270]

    ## DCM
    ## ---

    D_TEME_GCRF = r_gcrf_to_teme(
        JD_TT,
        -0.003875 * π / (180 * 3600),
        -0.052195 * π / (180 * 3600)
    )

    r_teme = D_TEME_GCRF * r_gcrf
    v_teme = D_TEME_GCRF * v_gcrf

    @test r_teme[1] ≈ +5094.18016210 atol=1e-7
    @test r_teme[2] ≈ +6127.64465950 atol=1e-7
    @test r_teme[3] ≈ +6380.34453270 atol=1e-7

    @test v_teme[1] ≈ -4.7461314870  atol=1e-9
    @test v_teme[2] ≈ +0.7858180410  atol=1e-9
    @test v_teme[3] ≈ +5.5319312880  atol=1e-9

    ## Quaternion
    ## ----------

    q_TEME_GCRF = r_gcrf_to_teme(
        Quaternion,
        JD_TT,
        -0.003875 * π / (180 * 3600),
        -0.052195 * π / (180 * 3600)
    )

    r_teme = vect(conj(q_TEME_GCRF) * r_gcrf * q_TEME_GCRF)
    v_teme = vect(conj(q_TEME_GCRF) * v_gcrf * q_TEME_GCRF)

    @test r_teme[1] ≈ +5094.18016210 atol=1e-7
    @test r_teme[2] ≈ +6127.64465950 atol=1e-7
    @test r_teme[3] ≈ +6380.34453270 atol=1e-7

    @test v_teme[1] ≈ -4.7461314870  atol=1e-9
    @test v_teme[2] ≈ +0.7858180410  atol=1e-9
    @test v_teme[3] ≈ +5.5319312880  atol=1e-9

    ## -------------------------------------------------------------------------
    ##                            Second Test
    ## -------------------------------------------------------------------------

    r_j2000 = [5102.50960000; 6123.01152000; 6378.13630000]
    v_j2000 = [-4.7432196000; 0.7905366000; 5.5337561900]

    ## DCM
    ## ---

    D_TEME_J2000 = r_gcrf_to_teme(JD_TT)

    r_teme = D_TEME_J2000*r_j2000
    v_teme = D_TEME_J2000*v_j2000

    @test r_teme[1] ≈ +5094.18016210 atol=1e-7
    @test r_teme[2] ≈ +6127.64465950 atol=1e-7
    @test r_teme[3] ≈ +6380.34453270 atol=1e-7

    @test v_teme[1] ≈ -4.7461314870  atol=1e-9
    @test v_teme[2] ≈ +0.7858180410  atol=1e-9
    @test v_teme[3] ≈ +5.5319312880  atol=1e-9

    ## Quaternion
    ## ----------

    q_TEME_J2000 = r_gcrf_to_teme(Quaternion, JD_TT)

    r_teme = vect(conj(q_TEME_J2000) * r_j2000 * q_TEME_J2000)
    v_teme = vect(conj(q_TEME_J2000) * v_j2000 * q_TEME_J2000)

    @test r_teme[1] ≈ +5094.18016210 atol=1e-7
    @test r_teme[2] ≈ +6127.64465950 atol=1e-7
    @test r_teme[3] ≈ +6380.34453270 atol=1e-7

    @test v_teme[1] ≈ -4.7461314870  atol=1e-9
    @test v_teme[2] ≈ +0.7858180410  atol=1e-9
    @test v_teme[3] ≈ +5.5319312880  atol=1e-9
end

# Functions r_teme_to_pef and r_pef_to_teme
# -----------------------------------------

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-15: Performing IAU-76/FK5 reduction.
#
# According to this example and Table 3-6, using:
#
#   JD_UT1 = date_to_jd(2004,4,6,7,51,28.386009) - 0.4399619/86400
#   LOD    = 0.0015563 s
#   r_teme = 5094.18016210   i + 6127.64465950   j + 6380.34453270   k [km]
#   v_teme =   -4.7461314870 i +    0.7858180410 j +    5.5319312880 k [km/s]
#
# one gets the following data:
#
#   r_pef  = -1033.47503130   i + 7901.30558560   j + 6380.34453270   k [km]
#   v_pef  =    -3.2256327470 i -    2.8724425110 j +    5.5319312880 k [km/s]
#
################################################################################

@testset "Functions r_teme_to_pef and r_pef_to_teme" begin
    JD_UT1 = date_to_jd(2004,4,6,7,51,28.386009) - 0.4399619 / 86400
    LOD    = 0.0015563
    w      = 7.292115146706979e-5 * (1 - LOD / 86400)

    ## r_teme_to_pef
    ## ==========

    r_teme = [5094.18016210; 6127.64465950; 6380.34453270]
    v_teme = [-4.7461314870; 0.7858180410; 5.5319312880]

    ## DCM
    ## ---

    D_TEME_PEF = r_teme_to_pef(JD_UT1)

    r_pef = D_TEME_PEF * r_teme
    v_pef = D_TEME_PEF * v_teme - [0, 0, w] × r_pef

    @test r_pef[1] ≈ -1033.47503130 atol=1e-7
    @test r_pef[2] ≈ +7901.30558560 atol=1e-7
    @test r_pef[3] ≈ +6380.34453270 atol=1e-7

    @test v_pef[1] ≈ -3.2256327470  atol=1e-9
    @test v_pef[2] ≈ -2.8724425110  atol=1e-9
    @test v_pef[3] ≈ +5.5319312880  atol=1e-9

    ## Quaternion
    ## ----------

    q_TEME_PEF = r_teme_to_pef(Quaternion, JD_UT1)

    r_pef = vect(conj(q_TEME_PEF) * r_teme * q_TEME_PEF)
    v_pef = vect(conj(q_TEME_PEF) * v_teme * q_TEME_PEF) - [0, 0, w] × r_pef

    @test r_pef[1] ≈ -1033.47503130 atol=1e-7
    @test r_pef[2] ≈ +7901.30558560 atol=1e-7
    @test r_pef[3] ≈ +6380.34453270 atol=1e-7

    @test v_pef[1] ≈ -3.2256327470  atol=1e-9
    @test v_pef[2] ≈ -2.8724425110  atol=1e-9
    @test v_pef[3] ≈ +5.5319312880  atol=1e-9

    ## r_pef_to_teme
    ## ==========

    r_pef  = [-1033.47503130; 7901.30558560; 6380.34453270]
    v_pef  = [-3.2256327470; -2.8724425110; +5.5319312880]

    ## DCM
    ## ---

    D_PEF_TEME = r_pef_to_teme(JD_UT1)

    r_teme = D_PEF_TEME * r_pef
    v_teme = D_PEF_TEME * (v_pef + [0, 0, w] × r_pef)

    @test r_teme[1] ≈ +5094.18016210 atol=1e-7
    @test r_teme[2] ≈ +6127.64465950 atol=1e-7
    @test r_teme[3] ≈ +6380.34453270 atol=1e-7

    @test v_teme[1] ≈ -4.7461314870  atol=1e-9
    @test v_teme[2] ≈ +0.7858180410  atol=1e-9
    @test v_teme[3] ≈ +5.5319312880  atol=1e-9

    ## Quaternion
    ## ----------

    q_PEF_TEME = r_pef_to_teme(Quaternion, JD_UT1)

    r_teme = vect(conj(q_PEF_TEME) * r_pef * q_PEF_TEME)
    v_teme = vect(conj(q_PEF_TEME) * (v_pef + [0, 0, w] × r_pef) * q_PEF_TEME)

    @test r_teme[1] ≈ +5094.18016210 atol=1e-7
    @test r_teme[2] ≈ +6127.64465950 atol=1e-7
    @test r_teme[3] ≈ +6380.34453270 atol=1e-7

    @test v_teme[1] ≈ -4.7461314870  atol=1e-9
    @test v_teme[2] ≈ +0.7858180410  atol=1e-9
    @test v_teme[3] ≈ +5.5319312880  atol=1e-9
end


text:
using Devectorize

macro test_te(expr, expect_ty, expect_mode)
    te = texpr(expr)
    if !isa(te, eval(expect_ty))
        error("texpr failed on: texpr_type of ($expr) == $(eval(expect_ty)) | got $(typeof(te))")
    end
    if !(tmode(te) == eval(expect_mode))
        error("texpr failed on: tmode of ($expr) == $(eval(expect_mode)) | got $(tmode(te))")
    end
end

macro expect_error(expr, errty)
    last_err = nothing
    try
        texpr(expr)
    catch err
        last_err = err
    end
    if last_err == nothing
        error("expect exception from $expr, but it raises nothing (unexpectedly)")
    elseif !isa(last_err, eval(errty))
        error("The expression $expr raises unexpected exception of $(typeof(last_err))")
    end
end

# terminal

@test_te 1      TNum{Int}       ScalarMode()
@test_te 2.5    TNum{Float64}   ScalarMode()
@test_te a      TVar            EWiseMode{0}()
@test_te x.b    TQVar           EWiseMode{0}()
@test_te x.y.b  TQVar           EWiseMode{0}()

# reference

@test_te a[1]       TScalarRef1     EWiseMode{0}()
@test_te a[i]       TGeneralRef1    EWiseMode{0}()
@test_te a.b[i]     TGeneralRef1    EWiseMode{0}()
@test_te a[:]       TRef1D          EWiseMode{1}()
@test_te a.b[:]     TRef1D          EWiseMode{1}()
@test_te a[:,1]     TRefCol         EWiseMode{1}()
@test_te a[:,j]     TRefCol         EWiseMode{1}()
@test_te a[1,:]     TRefRow         EWiseMode{1}()
@test_te a[i,:]     TRefRow         EWiseMode{1}()
@test_te a[:,:]     TRef2D          EWiseMode{2}()

# unary expressions

@test_te -a             TMap    EWiseMode{0}()
@test_te sin(1.2)       TMap    ScalarMode()
@test_te sin(a)         TMap    EWiseMode{0}()
@test_te sin(a[0])      TMap    EWiseMode{0}()
@test_te sin(a[:])      TMap    EWiseMode{1}()
@test_te sin(a[:,:])    TMap    EWiseMode{2}()
@test_te sin(a.b[:])    TMap    EWiseMode{1}()

# binary expressions

@test_te a + 1          TMap    EWiseMode{0}()
@test_te a + b          TMap    EWiseMode{0}()
@test_te a + b[1]       TMap    EWiseMode{0}()
@test_te a + b[:]       TMap    EWiseMode{1}()
@test_te a + b[:,:]     TMap    EWiseMode{2}()
@test_te a + b.c        TMap    EWiseMode{0}()

@test_te 1 + 2          TMap    ScalarMode()
@test_te 1 + b          TMap    EWiseMode{0}()
@test_te 1 + b[1]       TMap    EWiseMode{0}()
@test_te 1 + b[:]       TMap    EWiseMode{1}()
@test_te 1 + b[:,:]     TMap    EWiseMode{2}()

@test_te a[1] + 2           TMap    EWiseMode{0}()
@test_te a[1] + b           TMap    EWiseMode{0}()
@test_te a[1] + b[1]        TMap    EWiseMode{0}()
@test_te a[1] + b[:]        TMap    EWiseMode{1}()
@test_te a[1] + b[:,:]      TMap    EWiseMode{2}()

@test_te a[:] + 1           TMap    EWiseMode{1}()
@test_te a[:] + b           TMap    EWiseMode{1}()
@test_te a[:] + b[1]        TMap    EWiseMode{1}()
@test_te a[:] + b[:]        TMap    EWiseMode{1}()
@expect_error a[:] + b[:,:]   DeError

@test_te a[:,:] + 1         TMap    EWiseMode{2}()
@test_te a[:,:] + b         TMap    EWiseMode{2}()
@test_te a[:,:] + b[1]      TMap    EWiseMode{2}()
@test_te a[:,:] + b[:,:]    TMap    EWiseMode{2}()
@expect_error a[:,:] + b[:]   DeError

# ternary expressions

@test_te clamp(a, b, c)         TMap    EWiseMode{0}()
@test_te clamp(1, b, c)         TMap    EWiseMode{0}()
@test_te clamp(a, 2, c)         TMap    EWiseMode{0}()
@test_te clamp(a, b, 3)         TMap    EWiseMode{0}()
@test_te clamp(1, 2, c)         TMap    EWiseMode{0}()
@test_te clamp(1, b, 3)         TMap    EWiseMode{0}()
@test_te clamp(a, 2, 3)         TMap    EWiseMode{0}()
@test_te clamp(1, 2, 3)         TMap    ScalarMode()

@test_te clamp(a[:], b, c)      TMap    EWiseMode{1}()
@test_te clamp(a, b[:], c)      TMap    EWiseMode{1}()
@test_te clamp(a, b, c[:])      TMap    EWiseMode{1}()
@test_te clamp(a[:], b[:], c)   TMap    EWiseMode{1}()
@test_te clamp(a[:], b, c[:])   TMap    EWiseMode{1}()
@test_te clamp(a, b[:], c[:])   TMap    EWiseMode{1}()
@test_te clamp(a[:], b[:], c[:])    TMap    EWiseMode{1}()

@test_te clamp(a[:], 1, 2)      TMap    EWiseMode{1}()
@test_te clamp(1, b[:], 2)      TMap    EWiseMode{1}()
@test_te clamp(1, 2, c[:])      TMap    EWiseMode{1}()
@test_te clamp(a[:], b[:], 1)   TMap    EWiseMode{1}()
@test_te clamp(a[:], 1, c[:])   TMap    EWiseMode{1}()
@test_te clamp(1, b[:], c[:])   TMap    EWiseMode{1}()

# comparison & logical expressions

@test_te a .== b        TMap    EWiseMode{0}()
@test_te 1 .== a        TMap    EWiseMode{0}()
@test_te a .== 1        TMap    EWiseMode{0}()
@test_te 1 .== 2        TMap    ScalarMode()

@test_te a .!= b        TMap    EWiseMode{0}()
@test_te a .< b         TMap    EWiseMode{0}()
@test_te a .> b         TMap    EWiseMode{0}()
@test_te a .<= b        TMap    EWiseMode{0}()
@test_te a .>= b        TMap    EWiseMode{0}()

@test_te a & b          TMap    EWiseMode{0}()
@test_te a & true       TMap    EWiseMode{0}()
@test_te true & a       TMap    EWiseMode{0}()
@test_te true & false   TMap    ScalarMode()

@test_te a | b          TMap    EWiseMode{0}()
@test_te a | true       TMap    EWiseMode{0}()
@test_te true | a       TMap    EWiseMode{0}()
@test_te true | false   TMap    ScalarMode()

# compound ewise expressions

@test_te a + b .* c     TMap    EWiseMode{0}()
@test_te 1 + 2 .* 3     TMap    ScalarMode()
@test_te a + 2 .* 3     TMap    EWiseMode{0}()
@test_te 1 + 2 .* b     TMap    EWiseMode{0}()

@test_te a[:] + b .* 2      TMap    EWiseMode{1}()
@test_te a + b[:] .* c      TMap    EWiseMode{1}()
@test_te a + 3 .* c[:,:]    TMap    EWiseMode{2}()
@test_te a[0] + 2 .* c[0]   TMap    EWiseMode{0}()
@test_te a.x[:] + 2 .* c.y  TMap    EWiseMode{1}()

@test_te blend(a .> 1, c, d)    TMap    EWiseMode{0}()

# reduction expressions

@test_te sum(a)     TReduc  ReducMode()
@test_te mean(a)    TReduc  ReducMode()
@test_te minimum(a) TReduc  ReducMode()
@test_te maximum(a) TReduc  ReducMode()
@test_te dot(a, b)  TReduc  ReducMode()

@test_te sum(1)     TReduc  ScalarMode()
@test_te mean(1)    TReduc  ScalarMode()
@test_te minimum(1) TReduc  ScalarMode()
@test_te maximum(1) TReduc  ScalarMode()
@test_te dot(1, 2)  TReduc  ScalarMode()

@test_te sum(a[:])      TReduc  ReducMode()
@test_te mean(a[:])     TReduc  ReducMode()
@test_te minimum(a[:])  TReduc  ReducMode()
@test_te maximum(a[:])  TReduc  ReducMode()
@test_te dot(a, b[:])   TReduc  ReducMode()
@test_te sum(a.x)       TReduc  ReducMode()

@test_te sum(a[:,:])        TReduc  ReducMode()
@test_te mean(a[:,:])       TReduc  ReducMode()
@test_te minimum(a[:,:])    TReduc  ReducMode()
@test_te maximum(a[:,:])    TReduc  ReducMode()
@test_te dot(a, b[:,:])     TReduc  ReducMode()

@test_te sum(a + b)     TReduc  ReducMode()
@test_te mean(a + b)    TReduc  ReducMode()
@test_te minimum(a + b)     TReduc  ReducMode()
@test_te maximum(a + b)     TReduc  ReducMode()
@test_te dot(a + b, c .* d)     TReduc  ReducMode()

# scalar expressions

@test_te scalar(x)          TScalarVar      ScalarMode()
@test_te scalar(x[1])       TGeneralScalar  ScalarMode()
@test_te scalar(x[i,j])     TGeneralScalar  ScalarMode()
@test_te scalar(1)          TNum            ScalarMode()
@test_te scalar(sin(1))     TGeneralScalar  ScalarMode()
@test_te scalar(a.b)        TGeneralScalar  ScalarMode()
@test_te scalar(a.b.c)      TGeneralScalar  ScalarMode()

# assignment expressions

@test_te a = 1          TAssign{TVar, TNum{Int}}    ScalarMode()
@test_te a = x          TAssign{TVar, TVar}         EWiseMode{0}()
@test_te a = sin(x)     TAssign{TVar, TMap}         EWiseMode{0}()
@test_te a = x + y      TAssign{TVar, TMap}         EWiseMode{0}()
@test_te a = x + y[:]   TAssign{TVar, TMap}         EWiseMode{1}()
@test_te a = x[:,:]     TAssign{TVar, TRef2D}       EWiseMode{2}()

@test_te a[:] = 1           TAssign{TRef1D, TNum{Int}}  EWiseMode{1}()
@test_te a[:] = x           TAssign{TRef1D, TVar}       EWiseMode{1}()
@test_te a[:] = sin(x)      TAssign{TRef1D, TMap}       EWiseMode{1}()
@test_te a[:] = x + y       TAssign{TRef1D, TMap}       EWiseMode{1}()
@test_te a[:] = x + y[:]    TAssign{TRef1D, TMap}       EWiseMode{1}()
@expect_error a[:] = x[:,:]     DeError

@test_te a[:,i] = 1             TAssign{TRefCol, TNum{Int}} EWiseMode{1}()
@test_te a[:,i] = x             TAssign{TRefCol, TVar}      EWiseMode{1}()
@test_te a[:,i] = sin(x)        TAssign{TRefCol, TMap}      EWiseMode{1}()
@test_te a[:,i] = x + y         TAssign{TRefCol, TMap}      EWiseMode{1}()
@test_te a[:,i] = x + y[:]      TAssign{TRefCol, TMap}      EWiseMode{1}()
@expect_error a[:,i] = x[:,:]   DeError

@test_te a[i,:] = 1             TAssign{TRefRow, TNum{Int}} EWiseMode{1}()
@test_te a[i,:] = x             TAssign{TRefRow, TVar}      EWiseMode{1}()
@test_te a[i,:] = sin(x)        TAssign{TRefRow, TMap}      EWiseMode{1}()
@test_te a[i,:] = x + y         TAssign{TRefRow, TMap}      EWiseMode{1}()
@test_te a[i,:] = x + y[:]      TAssign{TRefRow, TMap}      EWiseMode{1}()
@expect_error a[i,:] = x[:,:]   DeError

@test_te a[:,:] = 1             TAssign{TRef2D, TNum{Int}}  EWiseMode{2}()
@test_te a[:,:] = x             TAssign{TRef2D, TVar}       EWiseMode{2}()
@test_te a[:,:] = sin(x)        TAssign{TRef2D, TMap}       EWiseMode{2}()
@test_te a[:,:] = x + y         TAssign{TRef2D, TMap}       EWiseMode{2}()
@test_te a[:,:] = x[:,:]        TAssign{TRef2D, TRef2D}     EWiseMode{2}()
@expect_error a[:,:] = x + y[:] DeError

ext: jl


text:
"""
    TopHatFilterType

Represents different fitting filter models.  A fitting filter is a linear operator that is 1) symmetric about the
center, and 2) the sum of the elements must be zero.  The standard filter is a top-hat shape, although Gaussian,
triangular, Savitsky-Golay and other shapes are also possible.
"""
abstract type TopHatFilterType end

"""
    VariableWidthFilter

A top-hat filter that varies in width with the FWHM of the detector.
"""
struct VariableWidthFilter <: TopHatFilterType end # Variable width filter

"""
    ConstantWidthFilter

A top-hat filter that has constant width determined by FWHM at Mn Kα for all channels.
"""
struct ConstantWidthFilter <: TopHatFilterType end # Constant width filter

"""
    GaussianFilter

A Gaussian-shaped filter that varies in width with the FWHM of the detector.  The Gaussian is offset to ensure the
sum of the filter elements is zero.
"""
struct GaussianFilter <: TopHatFilterType end # A variable width Gaussian filter

"""
The TopHatFilter struct represents a zero-sum symmetric second-derivative-like filter that when applied to
spectral data has the property of suppressing constant and slowly varying signals (like the continuum) while
retaining a linear signal for faster changing signals like the characteristic peaks.

See
  * F. H. Schamber Proc Symposium of "X-ray Fluorscence Analysis on Environmental Samples" Chapel Hill 1976 T Dzubay Ed.
  * P. Statham Anal Chem 49 no 14 Dec 1977

The TopHatFilter struct optimizes the memory and CPU use when applying top-hat filters to spectrum data.

The easiest way to implement a top-hat filter is as matrix F.  The rows represent the filters.  The product of the
filter and the data vector is the filtered spectrum.  The product of the filter times a diagnonal matrix constructed
from the data times the transpose of the filter is the covariance of the filtered data.  The diagonal matrix constructed
from the spectrum data is the covariance matrix associated with the spectrum data because the channels in the spectrum
data are independent (thus the matrix is diagnonal) and the magnitude equals the counts in each channels because the
spectrum data is nominally Poissonian and in the large number limit, the variance of a Poissonian random variable is
the number itself (σ=sqrt(N) => Var = N)

Notes on memory and code optimization:
The filter matrix is banded diagonal.  Approximately, 2.5% of the elements are non-zero.  This suggest use of the
BandedMatrix type.  The most expensive operation is calculating F⋅D⋅Fᵀ, the covariance matrix of the filtered data.
D is a diagonal matrix and so computing each element in F⋅D⋅Fᵀ reduces to a sum over a single variable.
Furthermore, the weighted least squares fit doesn't require the full F⋅D⋅Fᵀ, just diag(F⋅D⋅Fᵀ).  However, it turns out
that we can do better implementing our own banded matrix type largely because D is fully diagonal and the matrix
product F⋅D⋅Fᵀ reduces down to a sum over a single variable.  The product F⋅d and F⋅D⋅Fᵀ are readily
implemented as element-by-element multiplies and sums.  Thus storing the filter as offsets and row filters is
efficient in both memory and CPU use.
"""
struct TopHatFilter
    filttype::Type{<:TopHatFilterType}
    detector::Detector
    offsets::Vector{Int}  # offset to start of filter in row
    filters::Vector{Vector{Float64}} # Filter data
    weights::Vector{Float64} # Correlation compensation weights

    function TopHatFilter(
        ty::Type{<:TopHatFilterType},
        det::Detector,
        filt::AbstractMatrix{Float64},
        wgts::AbstractVector{Float64},
    )
        # Extract the contiguous non-zero elements
        dim = size(filt)[1]
        offsets = zeros(Int, dim)
        filts = fill(Vector{Float64}(), dim)
        for (r, row) in enumerate(eachrow(filt))
            start = findfirst(i -> i ≠ 0.0, row)
            if !isnothing(start)
                stop = findlast(i -> i ≠ 0.0, row)
                offsets[r] = start
                filts[r] = [row[start:stop]...]
            end
        end
        return new(ty, det, offsets, filts, wgts)
    end
end

function filterdata(filt::TopHatFilter, row::Int)::Vector{Float64}
    res = zeros(Float64, length(filt.filters))
    res[filt.offsets[row]:filt.offsets[row]+length(filt.filters[row])-1] = filt.filters[row]
    return res
end

function filterdata(filt::TopHatFilter, region::AbstractUnitRange{Int})::Matrix{Float64}
    res = zeros(Float64, length(region), length(region))
    foreach(r -> res[r-first(region)+1, :] = filterdata(filt, r)[region], region)
    return res
end

filterdata(filt::TopHatFilter) = #
    filterdata(filt, eachindex(filt.filters))

"""
    filteredcovar(filt::TopHatFilter, specdata::Vector{Float64}, row::Int, col::Int)::Float64

Compute the covariance matrix entry for the specified row and column.  The resulting covariance matrix is equal
to F⋅Ω⋅Fᵀ where F is the filter and Ω is the covariance matrix of the spectrum data.  Since each channel in the spectrum
is 1) independent; 2) Poisson distributed, Ωᵢⱼ = specdata[i] if i==j, 0 otherwise. Because Ω is diagonal, the matrix
multiplication reduces to sum over a single index.  Often this sum is zero because the non-zero portions of the
`row` and `col` filters don't intersect.

Note:  `specdata` should be preprocessed so that no element is less than or equal to zero.
"""
function filteredcovar(
    filt::TopHatFilter,
    specdata::Vector{Float64},
    i::Int,
    l::Int,
)::Float64
    function dot3(a, b, c) # vmapreduce((ai,bi,ci)->ai*bi*ci, +, a, b, c)
        sum = zero(eltype(a))
        @avx for i in eachindex(a)  # @avx takes overall time from 670 μs down to 430 μs
            sum += a[i] * b[i] * c[i]
        end
        return sum
    end
    fi, fl, oi, ol = filt.filters[i], filt.filters[l], filt.offsets[i], filt.offsets[l]
    roi = max(oi, ol):min(oi + length(fi) - 1, ol + length(fl) - 1) # The ROI over which both filters are non-zero.
    return length(roi) > 0 ? #
           dot3(
        view(fi, roi.start-oi+1:roi.stop-oi+1),
        view(specdata, roi),
        view(fl, roi.start-ol+1:roi.stop-ol+1),
    ) : 0.0
end

"""
    filtereddatum(filt::TopHatFilter, specdata::Vector{Float64}, ch::Int)::Float64 =

Compute a single channel in the filtered spectrum.
"""
filtereddatum(filt::TopHatFilter, specdata::Vector{Float64}, i::Int)::Float64 = dot(
    filt.filters[i],
    view(specdata, filt.offsets[i]:filt.offsets[i]+length(filt.filters[i])-1),
)

Base.size(filt::TopHatFilter) = (length(filt.filters), length(filt.filters))
Base.length(filt::TopHatFilter) = length(filt.filters)

Base.show(io::IO, thf::TopHatFilter) = print(io, "$(thf.filttype)[$(thf.detector)]")

"""
    buildfilter(det::Detector, a::Float64=1.0, b::Float64=2.0)::TopHatFilter

Build the default top-hat filter for the specified detector with the specified top and base parameters.
"""
buildfilter(det::Detector, a::Float64 = 1.0, b::Float64 = 2.0)::TopHatFilter =
    buildfilter(VariableWidthFilter, det, a, b)

"""
    buildfilter(::Type{<:TopHatFilterType}, det::Detector, a::Float64=1.0, b::Float64=1.0)::TopHatFilter

Build a top-hat-style filter for the specified detector with the specified top and base parameters. The
VariableWidthFilter and ConstantWidthFilter types are currently supported.
"""
function buildfilter(
    ty::Type{<:TopHatFilterType},
    det::Detector,
    a::Float64 = 1.0, # Top
    b::Float64 = 2.0,  # Base
)::TopHatFilter
    resf(::Type{VariableWidthFilter}, center, det) = resolution(center, det) # FWHM at center
    resf(::Type{ConstantWidthFilter}, center, det) = resolution(energy(n"Mn K-L3"), det) # FWHM at Mn Ka
    intersect(aa, bb, cc, dd) = max(0.0, min(bb, dd) - max(aa, cc)) # Length of intersection [a,b) with [c,d)
    filtint(e0, e1, minb, mina, maxa, maxb) =
        intersect(minb, mina, e0, e1) * (-0.5 / (mina - minb)) +
        intersect(mina, maxa, e0, e1) / (maxa - mina) +
        intersect(maxa, maxb, e0, e1) * (-0.5 / (maxb - maxa))
    M(astop, astart) = (channel(astop, det) - channel(astart, det) - 1.0) / 2.0   # base length in channels
    N(bstart, astart) = channel(astart, det) - channel(bstart, det) # Top length in channels
    cc = channelcount(det)
    filt, wgts = zeros(Float64, (cc, cc)), zeros(Float64, cc)
    for ch1 in eachindex(wgts)
        center = 0.5 * (energy(ch1, det) + energy(ch1 + 1, det)) # midpoint of channel
        res = resf(ty, center, det)
        ea = (center - 0.5 * a * res, center + 0.5 * a * res)
        eb = (ea[1] - 0.5 * b * res, ea[2] + 0.5 * b * res)
        chmin, chmax = channel(eb[1], det), channel(eb[2], det)
        if (chmin >= 1) && (chmax <= cc)
            for i = 0:(chmax-chmin)÷2
                filt[ch1, chmax-i] = (
                    filt[ch1, i+chmin] = filtint(
                        energy(i + chmin, det),
                        energy(i + chmin + 1, det),
                        eb[1],
                        ea[1],
                        ea[2],
                        eb[2],
                    )
                )
            end
            filt[ch1, chmax] =
                (filt[ch1, chmin] -= sum(filt[ch1, chmin:min(cc, chmax)]) / 2.0)
            @assert abs(sum(filt[ch1, :])) < 1.0e-12 "Filter $ch1 does not sum to zero."
            @assert all(i -> filt[ch1, i] == filt[ch1, chmax-(i-chmin)], chmin:chmax) "The $ch1-th filter is not symmetric - O"
        end
        m, n = M(ea[2], ea[1]), N(eb[1], ea[1])
        @assert m > 0.0 "m=$(m)!!!!"
        @assert n > 0.0 "n=$(n)!!!!"
        wgts[ch1] = 2.87 + 1.758e-4 * energy(ch1, det)
        # wgts[ch1] = 2.0*n*m/(n+2.0*m)  # Schamber's formula (see Schamber 1997, note the formula in Statham 1977, Anal Chem 49, 14 doesn't seem to work.)
    end
    @assert all(r -> abs(sum(r)) < 1.0e-12, eachrow(filt)) "A filter does not sum to zero."
    return TopHatFilter(ty, det, filt, wgts)
end

"""
    buildfilter(::Type{GaussianFilter}, det::Detector, a::Float64=1.0, b::Float64=5.0)::TopHatFilter

Build a top-hat filter with Gaussian shape whose width varies with the detector's resolution as a function of X-ray
energy for the specified detector with the specified top and base parameters. The `a` parameter corresponds
to the filter width relative to the detector resolution expressed as Gaussian width.  So `a=1` is a filter
whose width equals the detector resolution at each energy.  The `b` parameter is the extent of the
filter in Gaussian widths.  The default `a=1, b=5` corresponds to a  filter that has the same resolution
as the detector and an extent of 2.5 Gaussian widths above and below the center channel.
"""
function buildfilter(
    ty::Type{GaussianFilter},
    det::Detector,
    a::Float64 = 1.0,
    b::Float64 = 6.0, # Width
)::TopHatFilter
    filtint(center, ee, gw) = exp(-0.5 * ((ee - center) / gw)^2)
    cc = channelcount(det)
    filt, wgts = zeros(Float64, (cc, cc)), zeros(Float64, cc)
    for ch1 in eachindex(wgts)
        center = energy(ch1, det) # midpoint of channel
        res = a * gaussianwidth(resolution(center, det))
        ex = (center - 0.5 * b * res, center + 0.5 * b * res)
        chmin, chmax = channel(ex[1], det), channel(ex[2], det)
        if (chmin >= 1) && (chmax <= cc)
            for i = 0:(chmax-chmin)÷2 # Ensure that it is symmetric
                filt[ch1, chmax-i] =
                    (filt[ch1, chmin+i] = filtint(center, energy(chmin + i, det), res))
            end
            # Offset the Gaussian to ensure the sum is zero.
            filt[ch1, chmin:chmax] .-= sum(filt[ch1, chmin:chmax]) / length(chmin:chmax)
            @assert abs(sum(filt[ch1, :])) < 1.0e-12 "Filter $ch1 does not sum to zero."
            @assert all(i -> filt[ch1, i] == filt[ch1, chmax-(i-chmin)], chmin:chmax) "The $ch1-th filter is not symmetric - G"
        end
        wgts[ch1] = 2.87 + 1.758e-4 * energy(ch1, det)
    end
    @info "The uncertainty estimates will be about a factor of three low for the Gaussian filter."
    return TopHatFilter(ty, det, filt, wgts)
end

abstract type FilteredDatum end

"""
    FilteredReference

Represents the filtered reference spectrum over an ROI.
Carries the minimal data necessary to support filter-fitting a single
region-of-interest (continguous range of channles) and computing
useful output statistics.
"""
struct FilteredReference <: FilteredDatum
    label::ReferenceLabel # A way of identifying this filtered datum
    scale::Float64 # A dose or other scale correction factor
    roi::UnitRange{Int} # ROI for the raw data
    ffroi::UnitRange{Int} # ROI for the filtered data
    data::Vector{Float64} # Spectrum data over ffroi
    filtered::Vector{Float64} # Filtered spectrum data over ffroi
    charonly::Vector{Float64} # Background corrected intensity data over roi
    covscale::Float64  # Schamber's variance scale factor
end

Base.show(io::IO, fd::FilteredReference) = print(io, "Reference[$(fd.label)]")

"""
    _filter(
      spec::Spectrum,
      roi::UnitRange{Int},
      thf::TopHatFilter,
      tol::Float64
    )::FilteredReference

For filtering an ROI on a reference spectrum. Process a portion of a Spectrum with the specified filter.
"""
function _filter(spec::Spectrum, roi::UnitRange{Int}, thf::TopHatFilter, tol::Float64)
    # Extract the spectrum data as Float64 to match the filter
    data = counts(spec, Float64, true)
    # Determine tangents to the two background end points
    tangents = map(st -> estimatebackground(data, st, 5, 2), (roi.start, roi.stop))
    # Replace the non-ROI channels with extensions of the tangent functions
    data[1:roi.start-1] = map(tangents[1], 1-roi.start:-1)
    data[roi.stop+1:end] = map(tangents[2], 1:length(data)-roi.stop)
    # Compute the filtered data
    filtered = [filtereddatum(thf, data, i) for i in eachindex(data)]
    maxval = maximum(filtered)
    ff = findfirst(f -> abs(f) > tol * maxval, filtered)
    if isnothing(ff)
        error("There doesn't appear to be any data in $(spec[:Name]))")
    end
    roiff = ff:findlast(f -> abs(f) > tol * maxval, filtered)
    return (roiff, data[roiff], filtered[roiff])
end

function tophatfilter(
    escLabel::EscapeLabel,
    thf::TopHatFilter,
    scale::Float64,
    tol::Float64 = 1.0e-6,
)::FilteredReference
    spec, roi = spectrum(escLabel), escLabel.roi
    charonly = spec[roi] - modelBackground(spec, roi)
    f = _filter(spec, roi, thf, tol)
    return FilteredReference(
        escLabel,
        scale,
        roi,
        f...,
        charonly,
        thf.weights[(roi.start+roi.stop)÷2],
    )
end

"""
    function labeledextents(
        elm::Element,  # All CharXRay for this element
        det::Detector,
        ampl::Float64,
        maxE::Float64=energy(det.channelcount+1, det) # full detector range
    )::Vector{Tuple{Vector{CharXRay},UnitRange{Int}}}

Creates a vector containing pairs containing a vector of CharXRay and an interval. The interval represents a
contiguous interval over which all the X-rays in the interval are sufficiently close in energy that they will
interfere with each other on the specified detector.
"""
labeledextents(elm::Element, det::Detector, ampl::Float64, maxE::Float64 = 1.0e6) =
    labeledextents(
        isvisible(characteristic(elm, alltransitions, ampl, maxE), det),
        det,
        ampl,
    )

"""
    charXRayLabels(#
      spec::Spectrum, #
      elm::Element, #
      allElms::AbstractSet{Element}, #
      det::Detector, #
      ampl::Float64, #
      maxE::Float64=1.0e6)::Vector{SpectrumFeature}

Creates a vector of CharXRayLabel objects associated with 'elm' for a spectrum containing the elements
'allElms' assuming that it was collected on 'det'.  ROIs in which other elements from 'allElms'
interfere with 'elm' will not be included.
"""
function charXRayLabels(#
    spec::Spectrum, #
    elm::Element, #
    allElms::AbstractSet{Element}, #
    det::Detector, #
    maxE::Float64 = 1.0e6;
    ampl::Float64 = 1.0e-5, #
)::Vector{ReferenceLabel}
    suitable = suitablefor(elm, allElms, det, maxE, ampl)
    length(suitable) > 0 || @error "The spectrum $(name(spec)) provides no suitable ROIs for the element $(symbol(elm))."
    # Find all the ROIs associated with other elements
    return [ CharXRayLabel(spec, roi, xrays) for (xrays, roi) in suitable ]
end
charXRayLabels(#
    spec::Spectrum, #
    elm::Element, #
    allElms::AbstractVector{Element}, #
    det::Detector, #
    maxE::Float64 = 1.0e6;
    ampl::Float64 = 1.0e-5, #
)::Vector{ReferenceLabel} = charXRayLabels(spec, elm, Set{Element}(allElms), det, maxE, ampl=ampl)

"""
    suitablefor(elm::Element, allElms::AbstractSet{Element}, det::Detector, maxE::Float64 = 1.0e6, ampl::Float64 = 1.0e-5)::Vector{Tuple{Vector{CharXRay}, UnitRange{Int64}}}

Given a material containing `allElms` which ROIs for element `elm` is the material a suitable reference on the detector `det`?
This function provides informational, warning and error messages depending upon the suitability of the material.
"""
function suitablefor( #
    elm::Element, #
    allElms::AbstractSet{Element}, #
    det::Detector, #
    maxE::Float64 = 1.0e6,
    ampl::Float64 = 1.0e-5, #
    warnme::Bool = true #
)::Vector{Tuple{Vector{CharXRay}, UnitRange{Int64}}}
    # Find all the ROIs associated with other elements
    otherElms = Set(allElms)
    if elm in allElms
        delete!(otherElms, elm)
        res = if length(otherElms) > 0
            urs = mapreduce(ae -> extents(ae, det, ampl), append!, collect(otherElms))
            # Find elm's ROIs that don't intersect another element's ROI
            filter(labeledextents(elm, det, ampl, maxE)) do lx
                furs = filter(ur -> length(intersect(ur, lx[2])) > 0, urs)
                warnme && (!isempty(furs)) && @info "A material containing $(join(symbol.(allElms),", "," and ")) is not a suitable reference for \"$(lx[1])\" due to $(length(furs)==1 ? "a peak interference." : "$(length(furs)) peak interferences.")"
                isempty(furs)
            end
        else
            labeledextents(elm, det, ampl, maxE)
        end
        length(res) == 0 && warnme && @warn "A material containing $(join(symbol.(allElms),", "," and ")) provides no references for $(symbol(elm))."
    else
        @error "The element $(symbol(elm)) is not contained within the elements $(join(symbol.(allElms),", "," and "))."
        res = Tuple{Vector{CharXRay}, UnitRange{Int64}}[]
    end
    return res
end
function suitablefor( #
    elm::Element, #
    mat::Material, #
    det::Detector, #
    maxE::Float64 = 1.0e6,
    ampl::Float64 = 1.0e-5, #
    warnme::Bool = true #
)
    suitablefor(elm, keys(mat), det, maxE, ampl, warnme)
end

"""
    suitability(elm::Element, mats::AbstractSet{<:Material}, det::Detector; maxE=30.0e3)
    suitability(elm::Element, det::Detector; maxE=30.0e3, minC=0.1)

Tabulates the characteristic X-ray peaks for the `Element` for which there are suitable materials
in `mats` for the specified detector.  The second form uses a default set of Materials in the
file NeXLCore "standards.txt".

This function is helpful for determining which `Material`s are suitable to act as
fitting standards for the specified Element.  It shows how NeXLSpectrum will break up the
characteristic peaks associated with `elm` into contiguous regions each of which will be
fit independently. NeXLSpectrum attempts to break each element into as many independent
regions as possible dependent on the resolution of the specified `EDSDetector`.  If there
is an interference between one of the other elements in the `Material` and `elm` then
this peak will not be suitable as a fitting standard.  However, it can be used as a
similar standard.
"""
function suitability(elm::Element, mats::AbstractSet{<:Material}, det::EDSDetector; maxE=30.0e3, latex=false)
    wrap(matname, ltx) = ltx ? "\\ce{$matname}" : matname
    ex, chk = latex ? ( "\\xmark", "\\checkmark" ) : ("✗", "✓")
    mats = filter(m->value(m[elm])>0.0, mats)
    # Find all elemental ROIs
    rois = Dict{Vector{CharXRay}, Vector{Material}}()
    for (cxrs, _) in NeXLSpectrum.suitablefor(elm, Set( (elm, ) ), det, maxE, 1.0e-5, false)
      rois[cxrs] = Material[]
    end
    for mat in mats
      for (cxrs, _) in NeXLSpectrum.suitablefor(elm, mat, det, maxE, 1.0e-5, false)
        push!(rois[cxrs], mat)
      end
    end
    cxrss = collect(keys(rois))
    sort!(cxrss, lt = (a,b) -> energy(brightest(a))<energy(brightest(b)))
    res = DataFrame(:Material=>String[], :MassFrac =>Float64[], :Count=>Int[], ( Symbol(cxrs)=>String[] for cxrs in cxrss )...)
    for m in mats
      push!(res, [wrap(name(m), latex), round(value(m[elm]), digits=3), count(cxrs->m in rois[cxrs], cxrss), (m in rois[cxrs] ? chk : ex for cxrs in cxrss)...])
    end
    sort!(res, [ :Count, :MassFrac ], rev=true)
end
function suitability(elm::Element, det::EDSDetector; maxE=30.0e3, minC=0.1, latex=false)
    suitability(elm, getstandards(elm, minC), det, maxE=maxE, latex=latex)
end

function escapeextents(elm::Element, det::Detector, ampl::Float64, maxE::Float64)
    vis = isvisible(characteristic(elm, alltransitions, ampl, maxE), det)
    escapeextents(vis, det, ampl, maxE)
end

"""
    escapeLabels(#
      spec::Spectrum, #
      elm::Element, #
      allElms::AbstractSet{Element}, #
      det::Detector, #
      ampl::Float64, #
      maxE::Float64=1.0e6
    )::Vector{EscapeLabel}

Creates a vector EscapeLabel objects associated with 'elm' for a spectrum containing the elements
'allElms' assuming that it was collected on 'det'.  ROIs in which other elements from 'allElms'
interfere with 'elm' will not be included.
"""
function escapeLabels(#
    spec::Spectrum, #
    elm::Element, #
    allElms::AbstractSet{Element}, #
    det::Detector, #
    maxE::Float64 = 1.0e6;
    ampl::Float64 = 1.0e-5, #
)::Vector{ReferenceLabel}
    @assert elm in allElms "$elm must be in $allElms."
    intersects(urs, roi) = !isnothing(findfirst(ur -> length(intersect(ur, roi)) > 0, urs))
    # Find all the ROIs associated with the elements
    urs = collect(Iterators.flatten(extents(ae, det, ampl) for ae in allElms))
    # Find elm's ROIs that don't intersect an element's primary ROI
    lxs = filter(lx -> !intersects(urs, lx[2]), escapeextents(elm, det, ampl, maxE))
    isempty(lxs) || @warn "There are no escape peaks for $elm that don't interfere with one or more of $allElms."
    return ReferenceLabel[EscapeLabel(spec, roi, xrays) for (xrays, roi) in lxs]
end

"""
    tophatfilter(
        charLabel::CharXRayLabel,
        thf::TopHatFilter,
        scale::Float64 = 1.0,
        tol::Float64 = 1.0e-6,
    )::FilteredReference

For filtering an ROI on a reference spectrum. Process a portion of a Spectrum with the specified filter.  Use a simple
edge-based background model.
"""
function tophatfilter(
    charLabel::CharXRayLabel,
    filt::TopHatFilter,
    scale::Float64 = 1.0,
    tol::Float64 = 1.0e-6,
)::FilteredReference
    ashellof(xrays) = inner(xrays[argmax(jumpratio.(inner.(xrays)))])
    spec, roi, ashell = spectrum(charLabel), charLabel.roi, ashellof(charLabel.xrays)
    return FilteredReference(
        charLabel,
        scale,
        roi,
        _filter(spec, roi, filt, tol)...,
        spec[roi] - modelBackground(spec, roi, ashell),
        filt.weights[(roi.start+roi.stop)÷2],
    )
end

function tophatfilter(
    labels::AbstractVector{ReferenceLabel},
    filt::TopHatFilter,
    scale::Float64 = 1.0,
    tol::Float64 = 1.0e-6,
)::Vector{FilteredReference}
    lbls = filter(labels) do lbl
        res = (lbl.roi.start >= 1) && (lbl.roi.stop <= length(filt))
        (!res) && @warn "The ROI $(lbl.roi) not fully contained on [1, $(length(filt))] for $lbl."
        w = (weight(brightest(lbl.xrays)) > 1.0e-3)
        (!w) && @warn "No sufficiently bright X-rays for $lbl."
        res && w
    end
    return [ tophatfilter(lbl, filt, scale, tol) for lbl in lbls ]
end

"""
    tophatfilter(
      reflabel::ReferenceLabel,
      roi::UnitRange{Int},
      thf::TopHatFilter,
      scale = 1.0,
      tol = 1.0e-6
    )::FilteredReference

For filtering an ROI on a reference spectrum. Process a portion of a Spectrum with the specified filter. Use a naive
linear background model.
"""
function tophatfilter(
    reflabel::ReferenceLabel,
    thf::TopHatFilter,
    scale::Float64 = 1.0,
    tol::Float64 = 1.0e-6,
)::FilteredReference
    spec, roi = spectrum(reflabel), reflabel.roi
    return FilteredReference(
        reflabel,
        scale,
        roi,
        _filter(spec, roi, thf, tol)...,
        spec[roi] - modelBackground(spec, roi),
        thf.weights[(roi.start+roi.stop)÷2],
    )
end

"""
    FilteredUnknown

A FilteredDatum representing the unknown (currently only FilteredUnknownW).
"""
abstract type FilteredUnknown <: FilteredDatum end

Base.show(io::IO, fd::FilteredUnknown) = print(io, fd.label)

"""
    extract(fd::FilteredReference, roi::UnitRange{Int})

Extract the filtered data representing the specified range.  `roi` must fully encompass the filtered
data in `fd`.
"""
function NeXLUncertainties.extract(fd::FilteredReference, roi::UnitRange{Int})
    @assert fd.ffroi.start >= roi.start "$(fd.ffroi.start) < $(roi.start) in $(fd)"
    @assert fd.ffroi.stop <= roi.stop "$(fd.ffroi.stop) > $(roi.stop) in $(fd)"
    data = zeros(Float64, length(roi))
    nz = fd.ffroi.start-roi.start+1:fd.ffroi.stop-roi.start+1
    data[nz] = fd.filtered
    return data
end

"""
    extract(fd::FilteredUnknown, roi::UnitRange{Int})::AbstractVector{Float64}

Extract the filtered data representing the specified range.  `roi` must be fully contained within the
filtered data in `fd`.
"""
NeXLUncertainties.extract(
    fd::FilteredUnknown,
    roi::UnitRange{Int},
)::AbstractVector{Float64} = view(fd.filtered, roi)

_buildlabels(ffs::AbstractVector{FilteredReference}) = collect(ff.label for ff in ffs)
_buildscale(unk::FilteredUnknown, ffs::AbstractVector{FilteredReference}) =
    Diagonal([unk.scale / ffs[i].scale for i in eachindex(ffs)])

# Internal: Computes the residual spectrum based on the fit k-ratios
function _computeResidual(
    unk::FilteredUnknown,
    ffs::AbstractVector{FilteredReference},
    kr::UncertainValues,
)
    res = copy(unk.data)
    for ff in ffs
        res[ff.roi] -= (value(kr, ff.label) * ff.scale / unk.scale) * ff.charonly
    end
    return res
end

# Internal: Computes the peak and background count based on the fit k-ratios
function _computecounts( #
    unk::FilteredUnknown, #
    ffs::AbstractVector{FilteredReference}, #
    kr::UncertainValues, #
)::Dict{<:ReferenceLabel,NTuple{2,Float64}}
    res = Dict{ReferenceLabel,NTuple{2,Float64}}()
    for ff in ffs
        su = sum(unk.data[ff.roi])
        res[ff.label] = (
            su,
            su - (value(kr, ff.label) * ff.scale / unk.scale) * sum(ff.charonly)
        )
    end
    return res
end

"""
Ordinary least squares for either FilteredUnknown[G|W]
"""
function fitcontiguouso(
    unk::FilteredUnknown,
    ffs::AbstractVector{FilteredReference},
    chs::UnitRange{Int},
)::UncertainValues
    x, lbls, scale = _buildmodel(ffs, chs), _buildlabels(ffs), _buildscale(unk, ffs)
    return scale * olspinv(extract(unk, chs), x, 1.0, lbls)
end

"""
    selectBestReferences(refs::AbstractVector{FilteredReference})::Vector{FilteredReference}

For each roi in each element, pick the best FilteredReference with the highest
intensity in the characteristic-only data.
"""
function selectBestReferences(
    refs::AbstractVector{FilteredReference},
)::Vector{FilteredReference}
    rois = Dict{UnitRange,Tuple{FilteredReference,Float64}}()
    elms = unique(element(ref.label) for ref in refs)
    for elm in elms
        for ref in filter(ref -> element(ref.label) == elm, refs)
            comp = composition(ref.label)
            # Pick the reference with the largest charonly value but prefer pure elements over compounds
            mrf = maximum(ref.charonly) * (ismissing(comp) ? 0.01 : normalized(comp, elm)^2)
            if (!haskey(rois, ref.roi)) || (mrf > rois[ref.roi][2])
                rois[ref.roi] = (ref, mrf)
            end
        end
    end
    return [v[1] for v in values(rois)]
end

function filterreference(
    filt::TopHatFilter,
    spec::Spectrum,
    elm::Element,
    allElms::AbstractSet{Element};
    props::Dict{Symbol,<:Any} = Dict{Symbol,Any}(),
    withEsc::Bool = false,
)::Vector{FilteredReference}
    @assert elm in allElms "$elm not in $allElms"
    cprops = merge(spec.properties, props)
    if !haskey(cprops, :Detector)
        cprops[:Detector] = filt.detector
    end
    # Creates a list of unobstructed ROIs for elm as CharXRayLabel objects
    lbls = charXRayLabels(spec, elm, allElms, cprops[:Detector], cprops[:BeamEnergy])
    if withEsc
        append!(lbls, escapeLabels(spec, elm, allElms, cprops[:Detector], cprops[:BeamEnergy]))
    end
    # Filters the spectrum and returns a list of FilteredReference objects, one per ROI
    return tophatfilter(lbls, filt, 1.0 / (cprops[:ProbeCurrent] * cprops[:LiveTime]))
end

function filterreference(
    filt::TopHatFilter,
    spec::Spectrum,
    elm::Element,
    comp::Material;
    props::Dict{Symbol,<:Any} = Dict{Symbol,Any}(),
    withEsc::Bool = false,
)::Vector{FilteredReference}
    spec[:Composition] = get(spec, :Composition, comp)
    filterreference(filt, spec, elm, keys(comp), props = props, withEsc = withEsc)
end

filterreferences(
    filt::TopHatFilter,
    refs::Tuple{Spectrum,Element,Material}...;
    props::Dict{Symbol,<:Any} = Dict{Symbol,Any}(),
    withEsc::Bool = false,
) = mapreduce(
    ref -> filterreference(filt, ref..., props = props, withEsc = withEsc),
    append!,
    refs,
)


text:
### A Pluto.jl notebook ###
# v0.16.4

using Markdown
using InteractiveUtils

# ╔═╡ 8952762c-3810-11ec-33d4-9f44a17217c4
begin
	using Pkg
	Pkg.activate("..")
end

# ╔═╡ 19e810a5-da40-4fec-98cc-6c838bb0ad53
begin
	using JLD2
	using DataFrames
	using Plots
	using Dates
	using PrettyTables
	using PlutoUI
	using Statistics
	using StatsBase
end

# ╔═╡ 22370722-2285-43c2-97f5-998540a610bd
md"""
# Looking at Julia Development
Notebook to analize julia development over the years. Used GitHub.jl to gather PR information stored as a dataframe in data/ directory. Did not save actual objects and that seemed a bit risky from a PII standpoint.
"""

# ╔═╡ 375fa75c-5261-49f7-ba5d-d7dc1321e303
md"""
## Load Packages
Using outside enviornment instead of Pluto Pkg manager because GitHub.jl is yet to merge my PR
"""

# ╔═╡ e2fa7a94-2dc9-4f6f-8029-065d58fd696e
md"""
## Data Engineering
"""

# ╔═╡ 3505b1c1-1c67-4606-b647-e6b83dc1799e
df = load("../data/df.jld2")["df"];

# ╔═╡ 54c8fd9f-e6ca-4a61-8e07-ae91af05786c
md"""
Round off dates to days
"""

# ╔═╡ 78ddf957-ca20-41a1-9837-812b3b5c3bf3
safedate(x) = isnothing(x) ? x : Date(x)

# ╔═╡ 8d6ea0a9-2f81-44bf-8c84-7a06f6f9da5a
begin
	datecols = [:closed_at,:created_at,:updated_at, :merged_at]
	transform!(df, datecols .=> (x->safedate.(x)) .=> datecols)
end

# ╔═╡ da8dad6f-74bb-467b-bec3-4a3de25181c7
df_closed = df[df.state.=="closed",:];

# ╔═╡ 4ecc8970-0075-4a6f-97aa-74eed9f42f02
df_closed_per_month = combine(groupby(transform(df_closed, :closed_at .=> (x->floor.(x, Month)) .=> :closed_at), :closed_at), nrow=>:count);

# ╔═╡ 919edb83-990e-49d4-8d66-65e6e63fac5d
begin
	bar(df_closed_per_month.closed_at, df_closed_per_month.count; legend=false)
	vline!([Date(2018,8,9)])
	annotate!(Date(2018,8,9),350, text("v1.0.0",10,:left))
end

# ╔═╡ 7c0d339e-07cb-4f3c-ad1b-327dfab1f52c
with_terminal() do
	pretty_table(
		first(sort(combine(groupby(df, :owner), nrow=>:count), :count, rev=true), 25),
		tf=tf_markdown
	)
end

# ╔═╡ dd1e9e73-60f4-48ca-9bbf-7144c81bf38c
df2 = combine(x->only(x).changed_files,
	groupby(df[Not(isempty.(df.changed_files)),:],Not(:changed_files)))

# ╔═╡ baefdcf1-8006-484b-adb9-cba1ef4ea351
cfdf= combine(groupby(df2, :file),:number=>length=>:prs,:changes=>sum=>:changes)

# ╔═╡ 69c993a3-fb0a-4d2c-9afe-7610de9568d2
with_terminal() do
	pretty_table(
		first(sort(cfdf, :prs, rev=true),50),
		tf=tf_markdown);
end

# ╔═╡ 2c82ed83-3aba-4226-af2c-6cf6db9059e1
bar(filter(p->(p.second>40) && (p.first !="") ,countmap(last.(splitext.(cfdf.file))));legend=false)

# ╔═╡ 6c44e520-1b84-41f9-94ad-306a0be35078
bar(filter(p->p.second>50,countmap(first.(splitpath.(cfdf.file))));legend=false,xrotation=90)

# ╔═╡ 8113d6f1-aa26-4f4f-995f-f9ff21c2e181
days_open = ceil.(df_closed.closed_at .- df_closed.created_at, Dates.Day)

# ╔═╡ c8b105ce-1f47-4255-8c65-d669cd68aa8e
μ_days_open = mean(Dates.value.(days_open))

# ╔═╡ c23948c9-43fc-4245-964b-54bc61d6228f
median_days_open = median(Dates.value.(days_open))

# ╔═╡ cf96819f-8641-4821-af6c-e3a925d7762f
std_days_open = std(Dates.value.(days_open))

# ╔═╡ 3f682cd1-07de-4cf6-927f-9e2cb2e90ce2
begin
	hdo = histogram(days_open; xlims=[0,50],legend=false)
	sdo = scatter(df_closed.created_at,days_open;xrotation=45)
	plot(hdo, sdo;legend=false)
end

# ╔═╡ cb369955-fa62-4948-9348-d8a6380fdd77
scatter(df_closed.created_at,days_open;legend=false)

# ╔═╡ a9fe8e66-b98c-428f-9fbe-5b2c9ff96f27
begin
	hc = histogram(size.(df.commits,1);xlims=[-1,50],title="# Commits")
	hcf = histogram(size.(df.changed_files,1), title="# Changed files")
	plot(hc, hcf;legend=false)
end

# ╔═╡ Cell order:
# ╟─22370722-2285-43c2-97f5-998540a610bd
# ╠═375fa75c-5261-49f7-ba5d-d7dc1321e303
# ╠═8952762c-3810-11ec-33d4-9f44a17217c4
# ╠═19e810a5-da40-4fec-98cc-6c838bb0ad53
# ╠═e2fa7a94-2dc9-4f6f-8029-065d58fd696e
# ╠═3505b1c1-1c67-4606-b647-e6b83dc1799e
# ╟─54c8fd9f-e6ca-4a61-8e07-ae91af05786c
# ╠═78ddf957-ca20-41a1-9837-812b3b5c3bf3
# ╠═8d6ea0a9-2f81-44bf-8c84-7a06f6f9da5a
# ╠═da8dad6f-74bb-467b-bec3-4a3de25181c7
# ╠═4ecc8970-0075-4a6f-97aa-74eed9f42f02
# ╠═919edb83-990e-49d4-8d66-65e6e63fac5d
# ╠═7c0d339e-07cb-4f3c-ad1b-327dfab1f52c
# ╠═dd1e9e73-60f4-48ca-9bbf-7144c81bf38c
# ╠═baefdcf1-8006-484b-adb9-cba1ef4ea351
# ╠═69c993a3-fb0a-4d2c-9afe-7610de9568d2
# ╠═2c82ed83-3aba-4226-af2c-6cf6db9059e1
# ╠═6c44e520-1b84-41f9-94ad-306a0be35078
# ╠═8113d6f1-aa26-4f4f-995f-f9ff21c2e181
# ╠═c8b105ce-1f47-4255-8c65-d669cd68aa8e
# ╠═c23948c9-43fc-4245-964b-54bc61d6228f
# ╠═cf96819f-8641-4821-af6c-e3a925d7762f
# ╠═3f682cd1-07de-4cf6-927f-9e2cb2e90ce2
# ╠═cb369955-fa62-4948-9348-d8a6380fdd77
# ╠═a9fe8e66-b98c-428f-9fbe-5b2c9ff96f27


text:
function Well512a()
    return Well512a(())
end

function Well512a(arg0::Vector{jint})
    return Well512a((Vector{jint},), arg0)
end

function Well512a(arg0::jint)
    return Well512a((jint,), arg0)
end

function Well512a(arg0::jlong)
    return Well512a((jlong,), arg0)
end

function next_int(obj::Well512a)
    return jcall(obj, "nextInt", jint, ())
end


text:
# To ensure that tests can run on travis we have to do a little
# hackadoodle here. The tests require a license file. We include
# a license file that is only valid for one day (the day when
# change is submitted).
# If there is no valid license file, we default to that file.

using MathOptInterfaceMosek

using Test

using MathOptInterface
const MOI = MathOptInterface
const MOIT = MOI.Test
const MOIU = MOI.Utilities
const MOIB = MOI.Bridges

const optimizer = Mosek.Optimizer(QUIET = true, service = "mosek://solve.mosek.com:30080")
# 1e-3 needed for rotatedsoc3 test
const config = MOIT.TestConfig(atol=1e-3, rtol=1e-3, query=false)

# Mosek does not support names
MOIU.@model(Model,
            (MOI.Integer,),
            (MOI.EqualTo, MOI.LessThan, MOI.GreaterThan),
            (MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives, MOI.RotatedSecondOrderCone),
            (),
            (MOI.SingleVariable,),
            (MOI.ScalarAffineFunction,),
            (MOI.VectorOfVariables,),
            (MOI.VectorAffineFunction,))

@testset "SolverName" begin
    @test MOI.get(optimizer, MOI.SolverName()) == "Mosek"
end

@testset "Copy" begin
    # Currently does not work because get is missing for ConstraintSet
    # and ConstraintFunction, see https://github.com/JuliaOpt/MathOptInterfaceMosek.jl/issues/50
    #MOIT.copytest(optimizer, Model{Float64}())
end

@testset "Unit" begin
    # Mosek does not support names
    cached = MOIU.CachingOptimizer(Model{Float64}(), optimizer)
    MOIT.unittest(MOIB.QuadtoSOC{Float64}(MOIB.SplitInterval{Float64}(cached)),
                  config,
                  [# Does not support quadratic objective yet, needs
                   # https://github.com/JuliaOpt/MathOptInterface.jl/issues/529
                   "solve_qp_edge_cases",
                   # Find objective bound of 0.0 which is lower than 4.0
                   "solve_objbound_edge_cases",
                   # Cannot put multiple bound sets of the same type on a variable
                   "solve_integer_edge_cases"
                  ])
end

@testset "Continuous linear problems" begin
    # linear1 is failing for two reasons
    # * it does not remove constraints using a variable if this variable is deleted, see https://github.com/JuliaOpt/MathOptInterface.jl/issues/511
    # * it does not support duplicated terms, see https://github.com/JuliaOpt/MathOptInterfaceMosek.jl/issues/41
    MOIT.contlineartest(optimizer, config, ["linear1"])
end

# include("contquadratic.jl")
# @testset "Continuous quadratic problems" begin
#     # contquadratictest(GurobiSolver())
# end

@testset "Continuous conic problems" begin
    MOIT.contconictest(MOIB.SquarePSD{Float64}(MOIB.RootDet{Float64}(MOIB.GeoMean{Float64}(optimizer))),
                       config, ["exp", "rootdets", "logdet"])
end

@testset "Mixed-integer linear problems" begin
    MOIT.intlineartest(optimizer, config, ["int2"])
end

#include("test_jump.jl")


text:
################ Introduce the model structures ###############

abstract type DSGEModel end
abstract type ModelPrimatives end

struct REModelPrimatives{Q <: AbstractString} <: ModelPrimatives

    # This structure contains the critial model information extracted from a
    # model file.

    states::Array{Q,1}
    jumps::Array{Q,1}
    shocks::Array{Q,1}
    variables::Array{Q,1}
    parameters::Array{Q,1}
    parametervalues::Array{Q,1}
    equations::Array{Q,1}
    unassigned_parameters::Array{Q,1}

end

struct REModel{S <: Integer, Q <: AbstractString} <: DSGEModel

    # This structure contains information about a model in a form that the
    # model-solvers can work with.

    number_states::S
    number_jumps::S
    number_shocks::S
    number_variables::S
    number_equations::S
    jumps_approximated::Array{S,1}
    variables::Array{Q,1}
    nlsolve_static_function::Function
    static_function::Function
    dynamic_function::Function
    each_eqn_function::Array{Function,1}
    closure_function::Function
    closure_function_piecewise::Function

end

struct REModelPartial{S <: Integer, Q <: AbstractString} <: DSGEModel

    # This structure is similar to REModel, but relates to the case
    # where some parameter values have yet to be assigned.

    number_states::S
    number_jumps::S
    number_shocks::S
    number_variables::S
    number_equations::S
    jumps_approximated::Array{S,1}
    variables::Array{Q,1}
    nlsolve_static_function::Function
    static_function::Function
    dynamic_function::Function
    each_eqn_function::Array{Function,1}
    closure_function::Function
    closure_function_piecewise::Function
    unassigned_parameters::Array{Q,1}

end

############### Introduce the solutionscheme structures #######################

abstract type SolutionScheme end
abstract type ProjectionScheme <: SolutionScheme end
abstract type ChebyshevScheme <: ProjectionScheme end
abstract type SmolyakScheme <: ProjectionScheme end
abstract type PiecewiseLinearScheme <: ProjectionScheme end

struct PerturbationScheme{T <: Real, Q <: AbstractString} <: SolutionScheme

    steady_state::Array{T,1}
    cutoff::T
    order::Q

end

struct ChebyshevSchemeStoch{T <: Real, S <: Integer} <: ChebyshevScheme

    initial_guess::Union{T,Array{T,1}}
    node_generator::Function
    node_number::Union{S,Array{S,1}}
    num_quad_nodes::S
    order::Union{S,Array{S,1}}
    domain::Union{Array{T,1},Array{T,2}}
    tol_fix_point_solver::T
    tol_variables::T
    maxiters::S

end

struct ChebyshevSchemeDet{T <: Real, S <: Integer} <: ChebyshevScheme

    initial_guess::Union{T,Array{T,1}}
    node_generator::Function
    node_number::Union{S,Array{S,1}}
    order::Union{S,Array{S,1}}
    domain::Union{Array{T,1},Array{T,2}}
    tol_fix_point_solver::T
    tol_variables::T
    maxiters::S

end

struct SmolyakSchemeStoch{T <: Real, S <: Integer} <: SmolyakScheme

    initial_guess::Union{T,Array{T,1}}
    node_generator::Function
    num_quad_nodes::S
    layer::Union{S,Array{S,1}}
    domain::Union{Array{T,1},Array{T,2}}
    tol_fix_point_solver::T
    tol_variables::T
    maxiters::S

end

struct SmolyakSchemeDet{T <: Real, S <: Integer} <: SmolyakScheme

    initial_guess::Union{T,Array{T,1}}
    node_generator::Function
    layer::Union{S,Array{S,1}}
    domain::Union{Array{T,1},Array{T,2}}
    tol_fix_point_solver::T
    tol_variables::T
    maxiters::S

end

struct PiecewiseLinearSchemeStoch{T <: Real, S <: Integer} <: PiecewiseLinearScheme

    initial_guess::Union{T,Array{T,1}}
    node_number::Union{S,Array{S,1}}
    num_quad_nodes::S
    domain::Union{Array{T,1},Array{T,2}}
    tol_fix_point_solver::T
    tol_variables::T
    maxiters::S

end

struct PiecewiseLinearSchemeDet{T <: Real, S <: Integer} <: PiecewiseLinearScheme

    initial_guess::Union{T,Array{T,1}}
    node_number::Union{S,Array{S,1}}
    domain::Union{Array{T,1},Array{T,2}}
    tol_fix_point_solver::T
    tol_variables::T
    maxiters::S

end

##################### Introduce the solution structures ########################

abstract type ModelSolution end

abstract type PerturbationSolution <: ModelSolution end
abstract type PerturbationSolutionDet <: PerturbationSolution end
abstract type PerturbationSolutionStoch <: PerturbationSolution end

abstract type ProjectionSolution <: ModelSolution end
abstract type ProjectionSolutionDet <: ProjectionSolution end
abstract type ProjectionSolutionStoch <: ProjectionSolution end

struct FirstOrderSolutionStoch{T <: Real,S <: Integer} <: PerturbationSolutionStoch

    # x(t+1) = hx*x(t) + k*v(t+1)
    #   y(t) = gx*x(t)

    hbar::Union{T,Array{T,1}}           # steady state values for predetermined variables
    hx::Union{Array{T,2},Array{T,1}}    # Transition matrix for predetermined variables
    k::Union{Array{T,2},Array{T,1}}     # Innovation loading matrix
    gbar::Union{T,Array{T,1}}           # steady state values for nonpredetermined variables
    gx::Union{Array{T,2},Array{T,1}}    # Decision rule matrix linking nonpredetermined variables to predetermined variables
    sigma::Union{Array{T,2},Array{T,1}} # Innovation variance-covariance matrix
    grc::S                              # Number of eigenvalues greater than cutoff
    soln_type::String                   # "determinate", "indeterminate", or "explosive"

end

struct FirstOrderSolutionDet{T <: Real,S <: Integer} <: PerturbationSolutionDet

    # x(t+1) = hx*x(t)
    #   y(t) = gx*x(t)

    hbar::Union{T,Array{T,1}}           # steady state values for predetermined variables
    hx::Union{Array{T,2},Array{T,1}}    # Transition matrix for predetermined variables
    gbar::Union{T,Array{T,1}}           # steady state values for nonpredetermined variables
    gx::Union{Array{T,2},Array{T,1}}    # Decision rule matrix linking nonpredetermined variables to predetermined variables
    grc::S                              # Number of eigenvalues greater than cutoff
    soln_type::String                   # "determinate", "indeterminate", or "explosive"

end

struct SecondOrderSolutionStoch{T <: Real,S <: Integer} <: PerturbationSolutionStoch

    # x(t+1) = hx*x(t) + (1/2)*hss + (1/2)*[kron(I,x(t))]'hxx*[kron(I,x(t))] + k*v(t+1)
    #   y(t) = gx*x(t) + (1/2)*gss + (1/2)*[kron(I,x(t))]'gxx*[kron(I,x(t))]

    hbar::Union{T,Array{T,1}}           # steady state values for predetermined variables
    hx::Union{Array{T,2},Array{T,1}}    # Linear component in predetermined block
    hss::Array{T,1}                     # Intercepts in predetermined block
    hxx::Array{T,2}                     # Quadratic component in predetermined block
    k::Union{Array{T,2},Array{T,1}}     # Innovation loading matrix
    gbar::Union{T,Array{T,1}}           # steady state values for nonpredetermined variables
    gx::Union{Array{T,2},Array{T,1}}    # Linear component in non-predetermined block
    gss::Array{T,1}                     # Intercepts in predetermined block
    gxx::Array{T,2}                     # Quadratic component in non-predetermined block
    sigma::Union{Array{T,2},Array{T,1}} # Innovation variance-covariance matrix
    grc::S                              # Number of eigenvalues greater than cutoff
    soln_type::String                   # "determinate", "indeterminate", or "explosive"

end

struct SecondOrderSolutionDet{T <: Real,S <: Integer} <: PerturbationSolutionDet

    # x(t+1) = hx*x(t) + (1/2)*[kron(I,x(t))]'hxx*[kron(I,x(t))]
    #   y(t) = gx*x(t) + (1/2)*[kron(I,x(t))]'gxx*[kron(I,x(t))]

    hbar::Union{T,Array{T,1}}           # steady state values for predetermined variables
    hx::Union{Array{T,2},Array{T,1}}    # Linear component in predetermined block
    hxx::Array{T,2}                     # Quadratic component in predetermined block
    gbar::Union{T,Array{T,1}}           # steady state values for nonpredetermined variables
    gx::Union{Array{T,2},Array{T,1}}    # Linear component in non-predetermined block
    gxx::Array{T,2}                     # Quadratic component in non-predetermined block
    grc::S                              # Number of eigenvalues greater than cutoff
    soln_type::String                   # "determinate", "indeterminate", or "explosive"

end

struct ThirdOrderSolutionStoch{T <: Real, S <: Integer} <: PerturbationSolutionStoch

    # x(t+1) = hx*x(t) + (1/2)*hss + (1/2)*hxx*[kron(x(t),x(t)]
    #        + (1/6)*hsss + (1/6)*hssx*[x(t)] + (1/6)*hxxx*[kron(x(t),x(t),x(t))]
    #        + k*v(t+1)

    #   y(t) = gx*x(t) + (1/2)*gss + (1/2)*gxx*[kron(x(t),x(t)]
    #        +(1/6)*gsss + (1/6)*gssx*[x(t)] + (1/6)*gxxx*[kron(x(t),x(t),x(t))]

    hbar::Union{T,Array{T,1}}              # steady state values for predetermined variables
    hx::Union{Array{T,2},Array{T,1}}       # Linear component in predetermined block
    hss::Array{T,1}                        # Intercepts in predetermined block
    hxx::Array{T,2}                        # Quadratic component in predetermined block
    hsss::Array{T,1}                       # Third-order component in predetermined block
    hssx::Array{T,2}                       # Third-order component in predetermined block
    hxxx::Array{T,2}                       # Third-order component in predetermined block
    k::Union{Array{T,2},Array{T,1}}        # Innovation loading matrix
    gbar::Union{T,Array{T,1}}              # steady state values for nonpredetermined variables
    gx::Union{Array{T,2},Array{T,1}}       # Linear component in non-predetermined block
    gss::Array{T,1}                        # Intercepts in predetermined block
    gxx::Array{T,2}                        # Quadratic component in non-predetermined block
    gsss::Array{T,1}                       # Third-order component in non-predetermined block
    gssx::Array{T,2}                       # Third-order component in non-predetermined block
    gxxx::Array{T,2}                       # Third-order component in non-predetermined block
    sigma::Union{Array{T,2},Array{T,1}}    # Innovation variance-covariance matrix
    skewness::Union{Array{T,1},Array{T,2}} # Innovation skewness matrix
    grc::S                                 # Number of eigenvalues greater than cutoff
    soln_type::String                      # "determinate", "indeterminate", or "explosive"

end

struct ThirdOrderSolutionDet{T <: Real, S <: Integer} <: PerturbationSolutionDet

    # x(t+1) = hx*x(t) + (1/2)*hxx*[kron(x(t),x(t)] + (1/6)*hxxx*[kron(x(t),x(t),x(t))]

    #   y(t) = gx*x(t) + (1/2)*gxx*[kron(x(t),x(t)] + (1/6)*gxxx*[kron(x(t),x(t),x(t))]

    hbar::Union{T,Array{T,1}}              # steady state values for predetermined variables
    hx::Union{Array{T,2},Array{T,1}}       # Linear component in predetermined block
    hxx::Array{T,2}                        # Quadratic component in predetermined block
    hxxx::Array{T,2}                       # Third-order component in predetermined block
    gbar::Union{T,Array{T,1}}              # steady state values for nonpredetermined variables
    gx::Union{Array{T,2},Array{T,1}}       # Linear component in non-predetermined block
    gxx::Array{T,2}                        # Quadratic component in non-predetermined block
    gxxx::Array{T,2}                       # Third-order component in non-predetermined block
    grc::S                                 # Number of eigenvalues greater than cutoff
    soln_type::String                      # "determinate", "indeterminate", or "explosive"

end

struct ChebyshevSolutionStoch{T <: Real, S <: Integer, N} <: ProjectionSolutionStoch

    variables::Array{Array{T,N},1}       # Variables
    weights::Array{Array{T,N},1}         # Chebyshev polynomials
    nodes::Array{Array{T,1},1}           # Chebyshev nodes
    order::Union{S,Array{S,1}}           # Complete polynomial / tensor-product
    domain::Union{Array{T,2},Array{T,1}} # Domain for state variables / state variable
    sigma::Union{Array{T,2},Array{T,1}}  # Innovation variance-covariance matrix
    iteration_count::S                   # Number of iterations needed for convergence
    node_generator::Function             # Function to generate the nodes

end

struct ChebyshevSolutionDet{T <: Real, S <: Integer, N} <: ProjectionSolutionDet

    variables::Array{Array{T,N},1}       # Variables
    weights::Array{Array{T,N},1}         # Chebyshev polynomials
    nodes::Array{Array{T,1},1}           # Chebyshev nodes
    order::Union{S,Array{S,1}}           # Complete polynomial / tensor-product
    domain::Union{Array{T,2},Array{T,1}} # Domain for state variables / state variable
    iteration_count::S                   # Number of iterations needed for convergence
    node_generator::Function             # Function to generate the nodes

end

struct SmolyakSolutionStoch{T <: Real, S <: Integer} <: ProjectionSolutionStoch

    variables::Array{Array{T,1},1}       # Variables
    weights::Array{Array{T,1},1}         # Smolyak polynominals
    grid::Array{T,2}                     # Smolyak grid
    multi_index::Array{S,2}              # Smolyak multi index
    layer::Union{S,Array{S,1}}           # Isotropic / anisotropic
    domain::Union{Array{T,2},Array{T,1}} # Domain for state variables / state variable
    sigma::Union{Array{T,2},Array{T,1}}  # Innovation variance-covariance matrix
    iteration_count::S                   # Number of iterations needed for convergence
    node_generator::Function             # Function to generate the nodes

end

struct SmolyakSolutionDet{T <: Real, S <: Integer} <: ProjectionSolutionDet

    variables::Array{Array{T,1},1}       # Variables
    weights::Array{Array{T,1},1}         # Smolyak polynominals
    grid::Array{T,2}                     # Smolyak grid
    multi_index::Array{S,2}              # Smolyak multi index
    layer::Union{S,Array{S,1}}           # Isotropic / anisotropic
    domain::Union{Array{T,2},Array{T,1}} # Domain for state variables / state variable
    iteration_count::S                   # Number of iterations needed for convergence
    node_generator::Function             # Function to generate the nodes

end

struct PiecewiseLinearSolutionStoch{T <: Real, S <: Integer, N} <: ProjectionSolutionStoch

    variables::Array{Array{T,N},1}       # Variables
    nodes::Array{Array{T,1},1}           # Nodes
    domain::Union{Array{T,2},Array{T,1}} # Domain for state variables / state variable
    sigma::Union{Array{T,2},Array{T,1}}  # Innovation variance-covariance matrix
    iteration_count::S                   # Number of iterations needed for convergence

end

struct PiecewiseLinearSolutionDet{T <: Real, S <: Integer, N} <: ProjectionSolutionDet

    variables::Array{Array{T,N},1}       # Variables
    nodes::Array{Array{T,1},1}           # Nodes
    domain::Union{Array{T,2},Array{T,1}} # Domain for state variables / state variable
    iteration_count::S                   # Number of iterations needed for convergence

end


text:
vline(x::AbstractVector{<:Real}; args...) = plot([x'; x'], ylim(); args...)
hline(y::AbstractVector{<:Real}; args...) = plot(xlim(), [y'; y']; args...)

function grid(ax; args...)
    x = UnitRange{Int}(ax.get_xlim()...)
    y = UnitRange{Int}(ax.get_ylim()...)
    vline(x; args...)
    hline(y; args...)
end

@inline colorname(rgb::NTuple{4, Float64}) = colorname(rgb[1:3])
function colorname(rgb::NTuple{3, Float64})
    if haskey(COLORS, rgb)
        return COLORS[rgb]
    else
        rgbs = collect(keys(COLORS))
        (mn, k) = findmin(map(x->norm(x .- rgb), rgbs))
        return mn < 1.0 ? COLORS[rgbs[k]] : ""
    end
end

function decomment(x::Vector{<:AbstractString})
    filter!(a->a[1] != '#', x)
    return x
end

function clear!(x::PatchDict)
    for k in keys(x)
        x[k].remove()
        delete!(x, k)
    end
    return x
end

function message_box(b::Board, str::String, color::String="orange")
    ht = b.ax.text(0.5, 0.5, str, fontsize=36, horizontalalignment="center",
        verticalalignment="center", transform=b.ax.transAxes,
        bbox=Dict{String,Any}("facecolor"=>color, "alpha"=>0.8))
    return ht
end

function blank_figure(hf::PyPlot.Figure)
    if plt.fignum_exists(hf.number)
        h = hf
    else
        h = figure(num="Robo-puzzle")
    end

    for ax in h.axes
        h.delaxes(ax)
    end

    ax = h.add_axes([0, 0, 1, 1])

    ax.set_frame_on(false)
    ax.set_xticks([])
    ax.set_yticks([])
    return h, ax
end

function init_board(ncol::Integer, nrow::Integer)
    if isassigned(BOARD, 1)
        clear!(BOARD[1])
        h, ax = blank_figure(BOARD[1].hf)
    else
        h, ax = blank_figure(figure(num="Robo-puzzle"))
    end

    h.set_size_inches(ncol, nrow)

    ax.set_xlim(1, ncol+1)
    ax.set_ylim(1, nrow+1)

    grid(ax; color="black")

    return Board(h, ax)
end

function load_level(n::Integer)
    name = @sprintf("%02d.pzl", n)
    pth = joinpath(@__DIR__, "..", "levels", name)
    return load_level(pth)
end

function load_level(pth::AbstractString)
    if isfile(pth)
        make_board(pth)
    end
    return nothing
end


text:
### A Pluto.jl notebook ###
# v0.12.10

using Markdown
using InteractiveUtils

# ╔═╡ 7eeb28ae-f2e1-11ea-2933-955cee17a9bd
using Pkg, DrWatson

# ╔═╡ 7eeb666e-f2e1-11ea-1237-db0c2818822c
begin
	using Distributions
	using StatsPlots, Plots
	using Turing
	using Logging
	using LaTeXStrings
end

# ╔═╡ 6688dfb8-f2df-11ea-2511-2d9736b27265
md"## Clip-02-08t.jl"

# ╔═╡ 7eebee5e-f2e1-11ea-237e-89059c0b9245
md"### snippet 2.8"

# ╔═╡ 7ef86274-f2e1-11ea-25be-5f2b1538b966
md"##### Simple Metropolis algorithm."

# ╔═╡ 7ef91c84-f2e1-11ea-1ed4-c3ef76963769
begin
	n_samples = 10000
	a3d = ones(n_samples,1,1)
	w = 6; l = 3; n = w +l
	p = [0.5]
	for i in 2:n_samples
  		p_new = rand(Normal(p[i-1], 0.1), 1)[1]
  		if  p_new < 0
    		p_new = abs(p_new)
  		end
  		if p_new > 1
    		p_new = 2 - p_new
  		end
  		q0 = pdf(Binomial(n, p[i-1]), w)
  		q1 = pdf(Binomial(n, p_new), w)
  		append!(p, [rand(Uniform(0, 1), 1)[1] < q1/q0 ? p_new : p[i-1]])
	end
end

# ╔═╡ 7f05a540-f2e1-11ea-0d4a-f7c33efa0832
md"##### Create an MCMCChains.Chains object. This Chains object has 10000 samples, one variable and a single chain."

# ╔═╡ 7f0651fe-f2e1-11ea-00f7-79322fafaaf9
begin
	a3d[:, 1, 1] = p
	chns = MCMCChains.Chains(a3d, ["p",])
	Text(sprint(show, "text/plain", summarize(chns)))
end

# ╔═╡ 7f0bac38-f2e1-11ea-0817-093aba68676d
md"##### Show density and computed conjugate solution."

# ╔═╡ 7f16b116-f2e1-11ea-14fb-1ffcca62978a
begin
	x = 0:0.01:1
	density(chns, lab="Samples")
	plot!( x, pdf.(Beta( w+1 , n-w+1 ) , x ), lab="Conjugate solution")
end

# ╔═╡ 7f175f44-f2e1-11ea-31bf-0fbdc4545246
md"## End of clip-02-08t.jl"

# ╔═╡ Cell order:
# ╠═6688dfb8-f2df-11ea-2511-2d9736b27265
# ╠═7eeb28ae-f2e1-11ea-2933-955cee17a9bd
# ╠═7eeb666e-f2e1-11ea-1237-db0c2818822c
# ╟─7eebee5e-f2e1-11ea-237e-89059c0b9245
# ╟─7ef86274-f2e1-11ea-25be-5f2b1538b966
# ╠═7ef91c84-f2e1-11ea-1ed4-c3ef76963769
# ╟─7f05a540-f2e1-11ea-0d4a-f7c33efa0832
# ╠═7f0651fe-f2e1-11ea-00f7-79322fafaaf9
# ╟─7f0bac38-f2e1-11ea-0817-093aba68676d
# ╠═7f16b116-f2e1-11ea-14fb-1ffcca62978a
# ╟─7f175f44-f2e1-11ea-31bf-0fbdc4545246