text:
module EarthAttitude

using ..AstroBase: sec2rad
using ..Bodies
using ReferenceFrameRotations: angle_to_dcm, angleaxis_to_dcm, compose_rotation

export
    J2000,
    bias_precession_matrix_00,
    celestial_to_intermediate,
    earth_rotation_angle,
    equation_of_equinoxes_00,
    equation_of_equinoxes_complementary_terms,
    equation_of_origins,
    fukushima_williams_matrix,
    general_precession_in_longitude,
    greenwich_mean_sidereal_time00,
    greenwich_mean_sidereal_time06,
    greenwich_mean_sidereal_time82,
    mean_elongation,
    mean_longitude,
    mean_longitude_ascending_node,
    mean_longitude_minus_lan,
    mean_obliquity_of_ecliptic,
    numat,
    nutation,
    nutation_00a,
    nutation_00b,
    obliquity_of_ecliptic_06,
    polar_motion,
    precession_fukushima_williams06,
    precession_rate_part_of_nutation,
    s00,
    s06,
    tio_locator,
    xy06,
    nutation_a06,
    greenwich_apparent_sidereal_time06,
    equation_of_equinoxes_94,
    nutation_matrix80,
    precession_nutation00,
    precession_nutation_a00,
    precession_nutation_b00,
    precession_nutation_matrix_a00,
    precession_nutation_matrix_b00,
    precession_nutation_matrix_a06,
    s00a,
    s00b,
    nutation_matrix_day_a00,
    nutation_matrix_day_b00,
    num06a,
    nutation_matrix_day,
    xys00a,
    xys00b,
    s06a,
    equation_of_equinoxes_a00,
    equation_of_equinoxes_b00,
    greenwich_mean_sidereal_time_a06,
    equation_of_equinoxes_a06,
    greenwich_apparent_sidereal_time94,
    greenwich_apparent_sidereal_time_a00,
    greenwich_apparent_sidereal_time_b00,
    greenwich_apparent_sidereal_time_a06,
    celestial_to_intermediate_frame_of_date,
    celestial_to_intermediate_matrix,
    celestial_to_intermediate_matrix_a00,
    celestial_to_intermediate_matrix_b00,
    celestial_to_terrestrial_matrix,
    celestial_to_terrestrial_a00,
    celestial_to_terrestrial_b00,
    c2tpe

const J2000 = 2451545.0
const DAYS_PER_CENTURY = 36525.0
const ARCSECONDS_IN_CIRCLE = 1296000.0
const SECONDS_PER_DAY = 24.0 * 60.0 * 60.0

include(joinpath("constants", "mfals.jl"))
include(joinpath("constants", "nut_const.jl"))
include(joinpath("constants", "cio_locator.jl"))
include(joinpath("constants", "EE00.jl"))
include(joinpath("constants", "NUTATION80.jl"))

"""
    celestial_to_intermediate(x, y, s)

Returns celestial to intermediate-frame-of-date transformation matrix given
the Celestial Intermediate Pole location (`x`, `y` and the CIO locator `s`).

# Example

```jldoctest
julia> celestial_to_intermediate(0.2, 0.2, 0.1)
3×3 RotZYZ{Float64}(0.785398, 0.286757, -0.885398):
  0.976728   0.0774803  0.2
 -0.11811    0.972651   0.2
 -0.179034  -0.218968   0.959166
```
"""
function celestial_to_intermediate(x, y, s)
    r2 = x^2 + y^2
    e = r2 > 0.0 ? atan(y, x) : 0.0
    d = atan(sqrt(r2 / (1.0 - r2)))
    angle_to_dcm(e + s, -d, -e, :ZYZ)
end

"""
    polar_motion(rx, ry, sp)

Form the matrix of polar motion for coordinates of the pole (radians).

# Example

```jldoctest
julia> polar_motion(20, 30, 50)
3×3 RotZYX{Float64}(50.0, -20.0, -30.0):
  0.393785  -0.829946  -0.395124
 -0.10707    0.385514  -0.916469
  0.912945   0.403198   0.0629472
```
"""
function polar_motion(rx, ry, sp)
    angle_to_dcm(ry, rx, -sp, :XYZ)
end

"""
    earth_rotation_angle(jd1, jd2)

Return Earth rotation angle (radians) for a given UT1 2-part Julian Date (jd1, jd2).

# Example

```jldoctest
julia> earth_rotation_angle(2.4578265e6, 0.30434616919175345)
4.912208135094597
```
"""
function earth_rotation_angle(jd1, jd2)
    if jd1 < jd2
        d1 = jd1
        d2 = jd2
    else
        d1 = jd2
        d2 = jd1
    end
    t = d1 + d2 - J2000
    f = mod(d1, 1.0) + mod(d2, 1.0)
    mod2pi(2pi * (f + 0.7790572732640 + 0.00273781191135448 * t))
end

"""
    tio_locator(jd1, jd2)

Returns TIO locator s' position for a given TT 2-part Julian date (jd1, jd2).

# Example

```jldoctest
julia> AstroBase.tio_locator(2.4578265e6, 0.30434616919175345)
-3.9189245827947945e-11
```
"""
function tio_locator(jd1, jd2)
    t = (jd1 - J2000 + jd2) / DAYS_PER_CENTURY
    sec2rad(-47e-6 * t)
end

"""
    mean_obliquity_of_ecliptic(jd1, jd2)

Returns  obliquity of the ecliptic (radians) for a given Julian 2 part date (TT).

# Example

```jldoctest
julia> mean_obliquity_of_ecliptic(2.4578265e6, 0.30434616919175345)
0.40905376936136706
```
"""
function mean_obliquity_of_ecliptic(jd1, jd2)
    t = ((jd1 - J2000) + jd2) / DAYS_PER_CENTURY
    sec2rad(@evalpoly t 84381.448 -46.8150 0.00059 0.001813)
end

"""
    obliquity_of_ecliptic_06(jd1, jd2)

Returns obliquity of ecliptic (radians) for a given 2 part Julain date (TT).

# Example

```jldoctest
julia> obliquity_of_ecliptic_06(2.4578265e6, 0.30434616919175345)
0.409053547482157
```
"""
function obliquity_of_ecliptic_06(jd1, jd2)
    t = ((jd1 - J2000) + jd2) / DAYS_PER_CENTURY
    sec2rad(@evalpoly t 84381.406 -46.836769 -0.0001831 0.00200340 -0.000000576 -0.0000000434)
end

"""
    precession_fukushima_williams06(jd1, jd2)

Returns fukushima angles(radians) for a given 2 part Julian date (TT).

# Example

```jldoctest
julia> precession_fukushima_williams06(2.4578265e6, 0.30434616919175345)
(8.616170933989655e-6, 0.4090536093366178, 0.004201176043952816, 0.409053547482157)
```
"""
function precession_fukushima_williams06(jd1, jd2)
    t = ((jd1 - J2000) + jd2) / DAYS_PER_CENTURY

    sec2rad(@evalpoly t -0.052928 10.556378 0.4932044 -0.00031238 -0.000002788 0.0000000260),
    sec2rad(@evalpoly t 84381.412819 -46.811016 0.0511268 0.00053289 -0.000000440 -0.0000000176),
    sec2rad(@evalpoly t -0.041775 5038.481484 1.5584175 -0.00018522 -0.000026452 -0.0000000148),
    obliquity_of_ecliptic_06(jd1, jd2)
end

"""
    mean_anomaly(::Luna, t)

Returns mean anomaly of Moon for Julian centuries since J2000.0 in TDB.

# Example

```jldoctest
julia> mean_anomaly(moon, 23.0)
0.5891752616281019
```
"""
function mean_anomaly(::Luna, t)
    p = @evalpoly t 485868.249036 1717915923.2178 31.8792 0.051635 -0.00024470
    sec2rad(p % ARCSECONDS_IN_CIRCLE)
end

"""
    mean_anomaly(::Sun, t)

# Example

Returns mean anomaly of the Sun for Julian centuries since J2000.0 in TDB.
julia> mean_anomaly(sun, 23.0)
5.857396217361825
```
"""
function mean_anomaly(::Sun, t)
    p = @evalpoly t 1287104.793048 129596581.0481 -0.5532 0.000136 -0.00001149
    sec2rad(p % ARCSECONDS_IN_CIRCLE)
end

"""
    mean_longitude_minus_lan(::Luna, t)

Returns mean longitude of the Moon minus mean longitude of the ascending node for Julian
centuries since J2000.0 in TDB.

# Example

```jldoctest
julia> mean_longitude_minus_lan(23.0)
3.103138156410118
```
"""
function mean_longitude_minus_lan(::Luna, t)
    p = @evalpoly t 335779.526232 1739527262.8478 -12.7512 -0.001037 0.00000417
    sec2rad(p % ARCSECONDS_IN_CIRCLE)
end

"""
    mean_elongation(::Luna, t)

Returns mean elongation of the Moon from the Sun for Julian centuries since
J2000.0 in TDB.

# Example

```jldoctest
julia> mean_elongation(::Luna, 23.0)
2.8012040574296484
```
"""
function mean_elongation(::Luna, t)
    p = @evalpoly t 1072260.703692 1602961601.2090 -6.3706 0.006593 -0.00003169
    sec2rad(p % ARCSECONDS_IN_CIRCLE)
end

"""
    mean_longitude_ascending_node(::Luna, t)

Return fundamental argument for Julian centuries since J2000.0 in TDB.

# Example
```jldoctest
julia> mean_longitude_ascending_node(luna, 23.0)
4.904897783682109
```
"""
function mean_longitude_ascending_node(::Luna, t)
    p = @evalpoly t 450160.398036 -6962890.5431 7.4722 0.007702 -0.00005939
    sec2rad(p % ARCSECONDS_IN_CIRCLE)
end

"""
    mean_longitude(::Mercury, t)

Returns mean longitude of Mercury for Julian centuries since J2000.0 in TDB.

# Example

```jldoctest
julia> mean_longitude(mercury, 23.0)
2.160150897150834
```
"""
mean_longitude(::Mercury, t) = mod2pi(4.402608842 + 2608.7903141574t)

"""
    mean_longitude(::Venus, t)

Returns mean longitude of Venus for Julian centuries since J2000.0 in TDB.

# Example

```jldoctest
julia> mean_longitude(venus, 23.0)
0.9030394378238363
```
"""
mean_longitude(::Venus, t) = mod2pi(3.176146697 + 1021.3285546211t)

"""
    mean_longitude(::Earth, t)

Returns mean longitude of Earth for Julian centuries since J2000.0 in TDB.

# Example

```jldoctest
julia> mean_longitude(earth, 23.0)
1.501718780251826
```
"""
mean_longitude(::Earth, t) = mod2pi(1.753470314 + 628.3075849991t)

"""
    mean_longitude(::Mars, t)

Returns mean longitude of Mars for Julian centuries since J2000.0 in TDB.

# Example

```jldoctest
julia> mean_longitude(mars, 23.0)
5.276431642365657
```
"""
mean_longitude(::Mars, t) = mod2pi(6.203480913 + 334.0612426700t)

"""
    mean_longitude(::Jupiter, t)

Returns mean longitude of Jupiter for Julian centuries since J2000.0 in TDB.

# Example

```jldoctest
julia> mean_longitude(jupiter, 23.0)
6.233996285639864
```
"""
mean_longitude(::Jupiter, t) = mod2pi(0.599546497 + 52.9690962641t)

"""
    mean_longitude(::Saturn, t)

Returns mean longitude of Saturn for Julian centuries since J2000.0 in TDB.

# Example

```jldoctest
julia> mean_longitude(saturn, 23.0)
1.3735042049922535
```
"""
mean_longitude(::Saturn, t) = mod2pi(0.874016757 + 21.3299104960t)

"""
    mean_longitude(::Uranus, t)

Returns mean longitude of Uranus for Julian centuries since J2000.0 in TDB.

# Example

```jldoctest
julia> mean_longitude(uranus, 23.0)
1.5497819750715893
```
"""
mean_longitude(::Uranus, t) = mod2pi(5.481293872 + 7.4781598567t)

"""
    mean_longitude(::Neptune, t)

Returns mean longitude of Neptune for Julian centuries since J2000.0 in TDB.

# Example

```jldoctest
julia> mean_longitude(neptune, 23.0)
5.053273953885775
```
"""
mean_longitude(::Neptune, t) = mod2pi(5.311886287 + 3.8133035638t)

"""
    general_precession_in_longitude(t)

Returns general accumulated precession in longitude for Julian centuries since J2000.0 in TDB.

# Example

```jldoctest
julia> general_precession_in_longitude(23.0)
0.56362992539
```
"""
general_precession_in_longitude(t) = @evalpoly t 0.00000538691 0.024381750

"""
    xy06(jd1, jd2)

Returns X, Y coordinates of celestial intermediate pole for a given 2-part Julian date in TT.

# Example

```jldoctest
julia> xy06(2.4578265e6, 0.30440190993249416)
(0.0016558850230577835, -3.986943362456243e-5)
```
"""
function xy06(jd1, jd2)
    t = ((jd1 - J2000) + jd2) / DAYS_PER_CENTURY
    # Powers of T.
    pt = [t^i for i = 0:6]

    fa = (mean_anomaly(luna, t),
          mean_anomaly(sun, t),
          mean_longitude_minus_lan(luna, t),
          mean_elongation(luna, t),
          mean_longitude_ascending_node(luna, t),
          mean_longitude(mercury, t),
          mean_longitude(venus, t),
          mean_longitude(earth, t),
          mean_longitude(mars, t),
          mean_longitude(jupiter, t),
          mean_longitude(saturn, t),
          mean_longitude(uranus, t),
          mean_longitude(neptune, t),
          general_precession_in_longitude(t))

    xpr = @evalpoly t x_coeff[1] x_coeff[2] x_coeff[3] x_coeff[4] x_coeff[5] x_coeff[6]
    ypr = @evalpoly t y_coeff[1] y_coeff[2] y_coeff[3] y_coeff[4] y_coeff[5] y_coeff[6]

    xypl = zeros(2)
    xyls = zeros(2)

    ialast = length(amp)
    for ifreq in reverse(eachindex(planetary))
        arg = 0.0
        for i in eachindex(fa)
            m = planetary[ifreq][i]
            arg += float(m) * fa[i]
        end

        sc = sincos(arg)

        ia = pointers_to_amp[ifreq + length(luni_solar)]
        for i in ialast:-1:ia
            j = i - ia + 1
            jxy = jaxy[j]
            jsc = jasc[j]
            jpt = japt[j]
            xypl[jxy] += amp[i] * sc[jsc] * pt[jpt]
        end
        ialast = ia - 1
    end

    for ifreq in reverse(eachindex(luni_solar))
        arg = 0.0
        for i in 1:5
           m = luni_solar[ifreq][i]
           arg += float(m) * fa[i]
        end

        sc = sincos(arg)

        ia = pointers_to_amp[ifreq]
        for i in ialast:-1:ia
            j = i - ia + 1
            jxy = jaxy[j]
            jsc = jasc[j]
            jpt = japt[j]
            xyls[jxy] += amp[i] * sc[jsc] * pt[jpt]
        end
        ialast = ia - 1
    end
    sec2rad((xpr + (xyls[1] + xypl[1]) / 1e6)), sec2rad(ypr + (xyls[2] + xypl[2]) / 1e6)
end

"""
    nutation_00b(jd1, jd2)

Returns luni-solar and planetary nutation for a given 2 part Julian date (TT).

# Example

```jldoctest
julia> nutation_00b(2.4578265e6, 0.30440190993249416)
(-3.7589177912131684e-5, -3.6657431214029895e-5)
```
"""
function nutation_00b(jd1, jd2)
    t = ((jd1 - J2000) + jd2) / DAYS_PER_CENTURY

    el  = sec2rad(mod(485868.249036 + (1717915923.2178) * t, ARCSECONDS_IN_CIRCLE))
    elp = sec2rad(mod(1287104.79305 + (129596581.0481) * t, ARCSECONDS_IN_CIRCLE))
    f   = sec2rad(mod(335779.526232 + (1739527262.8478) * t, ARCSECONDS_IN_CIRCLE))
    d   = sec2rad(mod(1072260.70369 + (1602961601.2090) * t, ARCSECONDS_IN_CIRCLE))
    om  = sec2rad(mod(450160.398036 + (-6962890.5431) * t, ARCSECONDS_IN_CIRCLE))

    dp = 0.0
    de = 0.0

    for i in reverse(eachindex(x_nutation))
        arg = mod(x_nutation[i][1] * el
                  + x_nutation[i][2] * elp
                  + x_nutation[i][3] * f
                  + x_nutation[i][4] * d
                  + x_nutation[i][5] * om, 2pi)
        sarg, carg = sincos(arg)

        dp += (x_nutation[i][6] + x_nutation[i][7] * t) * sarg + x_nutation[i][8] * carg
        de += (x_nutation[i][9] + x_nutation[i][10] * t) * carg + x_nutation[i][11] * sarg
    end

    sec2rad(-0.135e-3) + sec2rad(1e-7dp), sec2rad(0.388e-3) + sec2rad(1e-7de)
end

 """
    fukushima_williams_matrix(gamb, phib, psi, eps)

Returns  obliquity of the ecliptic (radians) for a given Julian 2 part date (TT).

# Example

```jldoctest
julia> fukushima_williams_matrix(0.2,0.3,0.5,0.6)
3×3 RotMatrix{Float64}:
  0.951082   0.21718   0.219716
 -0.274534   0.920305  0.278692
 -0.14168   -0.325378  0.93491
```
"""
function fukushima_williams_matrix(gamb, phib, psi, eps)
    compose_rotation(angleaxis_to_dcm(eps, [1.0, 0.0, 0.0]), angle_to_dcm(psi, -phib, -gamb, :ZXZ))
end

"""
    numat(epsa, dpsi, deps)

Returns nutation matrix for a given epsa(mean obliquity), dpsi and deps nutation.

# Example

```jldoctest
julia> numat(0.7, 1.4, 1.3)
3×3 RotXZX{Float64}(0.7, -1.4, -2.0):
  0.169967  -0.410092   0.896067
 -0.753714   0.531687   0.386296
 -0.634844  -0.741035  -0.218722
```
"""
function numat(epsa, dpsi, deps)
    angle_to_dcm(epsa + deps, dpsi, -epsa, :XZX)
end

"""
    greenwich_mean_sidereal_time82(jd1, jd2)

Returns Greenwich mean sidereal time(radians) for given 2 part Julian dates (UT1).
(consistent with IAU 1982 model)

# Example

```jldoctest
julia> greenwich_mean_sidereal_time82(2.4578265e6, 0.30434616919175345)
4.916054244834956
```
"""
function greenwich_mean_sidereal_time82(jd1, jd2)
    A = 24110.54841  -  SECONDS_PER_DAY / 2.0
    B = 8640184.812866
    C = 0.093104
    D = -6.2e-6

    if jd1 < jd2
        d1 = jd1
        d2 = jd2
    else
        d1 = jd2
        d2 = jd1
    end

    t = (d1 + (d2 - J2000)) / DAYS_PER_CENTURY
    f = SECONDS_PER_DAY * (mod(d1, 1.0) + mod(d2, 1.0))
    mod2pi(7.272205216643039903848712e-5 * (@evalpoly t A + f B C D))
end

"""
    greenwich_mean_sidereal_time00(ut1, ut2, tt1, tt2)

Returns Greenwich mean sidereal time(radians) for given two, 2 part Julian dates (TT and UT1).
(consistent with IAU 2000 precession)

# Example

```jldoctest
julia> greenwich_mean_sidereal_time00(2.4579405e6, 0.0, 2.4579405e6, -0.0007966009351851851)
4.9596733720586075
```
"""
function greenwich_mean_sidereal_time00(ut1, ut2, tt1, tt2)
    t = ((tt1 - J2000) + tt2) / DAYS_PER_CENTURY
    mod2pi(earth_rotation_angle(ut1, ut2) +
           sec2rad(@evalpoly t 0.014506 4612.15739966 1.39667721 -0.00009344 0.00001882))
end

"""
    greenwich_mean_sidereal_time06(ut1, ut2, tt1, tt2)

Returns Greenwich mean sidereal time(radians) for given two, 2 part Julian dates (TT and UT1).
(consistent with IAU 2006 precession)

# Example

```jldoctest
julia> greenwich_mean_sidereal_time06(2.4579405e6, 0.0, 2.4579405e6, -0.0007966009351851851)
4.959673370568533
```
"""
function greenwich_mean_sidereal_time06(ut1, ut2, tt1, tt2)
    t = ((tt1 - J2000) + tt2) / DAYS_PER_CENTURY
    mod2pi(earth_rotation_angle(ut1, ut2) +
           sec2rad(@evalpoly t 0.014506 4612.156534 1.3915817 -0.00000044 -0.000029956 -0.0000000368))
end

"""
    precession_rate_part_of_nutation(jd1, jd2)

Returns precession corrections for a given 2 part Julian date (TT).

# Example

```jldoctest
julia> precession_rate_part_of_nutation(2400000.5, 53736)
-0.8716465172668347629e-7, -0.7342018386722813087e-8
```
"""
function precession_rate_part_of_nutation(jd1, jd2)
    t = ((jd1 - J2000) + jd2) / DAYS_PER_CENTURY
    -sec2rad(0.29965) * t, -sec2rad(0.02524) * t
end

"""
    bias_precession_matrix_00(jd1, jd2)

Returns a tuple of array, rb (frame bias), rp (precession matrix), rbp (bias precession matrix)
given 2 part Julian date (TT).

# Example

julia> bias_precession_matrix_00(2.4578265e6, 0.30434616919175345)
(RotZYX(-7.07828e-8, -8.05622e-8, 3.30604e-8), [0.999991 0.00384587 0.00167106; -0.00384587 0.999993 -3.2343e-6; -0.00167106 -3.1924e-6 0.999999], [0.999991 0.00384594 0.00167098; -0.00384594 0.999993 -3.26705e-6; -0.00167098 -3.15946e-6 0.999999])
```
"""
function bias_precession_matrix_00(jd1, jd2)
    t = ((jd1 - J2000) + jd2) / DAYS_PER_CENTURY

    dpbias = sec2rad(-0.041775)
    debias = sec2rad(-0.0068192)
    dra0 = sec2rad(-0.0146)
    eps0 = sec2rad(84381.448)

    psia77 = sec2rad((@evalpoly t 0 5038.7784 -1.07259 -0.001147))
    oma77  = eps0 + sec2rad(@evalpoly t 0 0 0.05127 -0.007726)
    chia   = sec2rad((@evalpoly t 0 10.5526 -2.38064 -0.001125))

    dpsipr, depspr = precession_rate_part_of_nutation(jd1, jd2)
    psia = psia77 + dpsipr
    oma  = oma77 + depspr

    rbw = angle_to_dcm(debias, -(dpbias * sin(eps0)), -dra0, :XYZ)
    rp = compose_rotation(angleaxis_to_dcm(-chia, [0.0, 0.0, 1.0]),
                          angle_to_dcm(oma, psia, -eps0, :XZX))

    rbw, rp, compose_rotation(rp, rbw)
end

"""
    equation_of_origins(rnpb, s)

Returns the equation of origins(radians) for given nutation-bias-precession matrix and the CIO locator.

 # Example

 ```jldoctest
equation_of_origins(rand(3,3), 0.2)
1.7738370040531068
 ```
"""
function equation_of_origins(rnpb, s)
    x = rnpb[1, 3]
    ax = x / (1.0 + rnpb[3, 3])
    xs = 1.0 - ax * x
    ys = -ax * rnpb[2, 3]
    zs = -x
    p = rnpb[1, 1] * xs + rnpb[2, 1] * ys + rnpb[3, 1] * zs
    q = rnpb[1, 2] * xs + rnpb[2, 2] * ys + rnpb[3, 2] * zs
    p != 0 || q != 0 ? s - atan(q, p) : s
end

"""
    nutation(jd1, jd2)

Returns nutation in longitude(radians) and obliquity(radians) for a given 2 part Julian date (TT format).

# Example

julia> nutation(2.4578265e6, 0.30434616919175345)
(-3.7565297299394694e-5, -3.665617105048724e-5
```
"""
function nutation(jd1, jd2)

    t = ((jd1 - J2000) + jd2) / DAYS_PER_CENTURY

    mean_longitude_moon_minus_mean_longitude_moon_perigee = rem2pi(
        sec2rad(@evalpoly t 485866.733 715922.633 31.310 0.064 )
      + mod(1325.0 * t, 1.0) * 2pi, RoundNearest)

    mean_longitude_sun_minus_mean_longitude_sun_perigee = rem2pi(
        sec2rad(@evalpoly t 1287099.804 1292581.224 -0.577 - 0.012)
      + mod(99.0 * t, 1.0) * 2pi, RoundNearest)

    mean_longitude_moon_minus_mean_longitude_moon_node = rem2pi(
        sec2rad(@evalpoly t 335778.877 295263.137 -13.257 0.011)
      + mod(1342.0 * t, 1.0) * 2pi, RoundNearest)

    mean_elongation_moon_from_sun = rem2pi(
        sec2rad(@evalpoly t 1072261.307 1105601.328 -6.891 0.019)
      + mod(1236.0 * t, 1.0) * 2pi, RoundNearest)

    mean_ascending_node_lunar_orbit_ecliptic_measured_mean_equinox_date = rem2pi(
        sec2rad(@evalpoly t 450160.280 -482890.539 7.455 0.008)
      + mod(-5.0 * t, 1.0) * 2pi, RoundNearest)

    dp = 0.0
    de = 0.0

    for i in reverse(eachindex(multiples_of_coefficients))
        arg = (multiples_of_coefficients[i][1]  * mean_longitude_moon_minus_mean_longitude_moon_perigee
        + multiples_of_coefficients[i][2] * mean_longitude_sun_minus_mean_longitude_sun_perigee
        + multiples_of_coefficients[i][3]  * mean_longitude_moon_minus_mean_longitude_moon_node
        + multiples_of_coefficients[i][4]  * mean_elongation_moon_from_sun
        + multiples_of_coefficients[i][5] * mean_ascending_node_lunar_orbit_ecliptic_measured_mean_equinox_date)

        s = multiples_of_arguments[i][1] + multiples_of_arguments[i][2] * t
        c = multiples_of_arguments[i][3] + multiples_of_arguments[i][4] * t
        sinarg, cosarg = sincos(arg)
        dp += s * sinarg
        de += c * cosarg
    end

    sec2rad(1e-4dp), sec2rad(1e-4de)
end

"""
    equation_of_equinoxes_complementary_terms(jd1, jd2)

Returns complementary terms for a given 2 part Julian date (TT).

# Example

julia> equation_of_equinoxes_complementary_terms(2.4578265e6, 0.30434616919175345)
5.706799075288604e-9
```
"""
function equation_of_equinoxes_complementary_terms(jd1, jd2)
    t = ((jd1 - J2000) + jd2) / DAYS_PER_CENTURY

    fa = (mean_anomaly(luna, t),
          mean_anomaly(sun, t),
          mean_longitude_minus_lan(luna, t),
          mean_elongation(luna, t),
          mean_longitude_ascending_node(luna, t),
          mean_longitude(venus, t),
          mean_longitude(earth, t),
          general_precession_in_longitude(t))

    s0 = 0.0
    s1 = 0.0

    for i in reverse(eachindex(e0_coefficent))
        a = 0.0
        for j in 1:8
            a += e0_coefficent[i][j] * fa[j]
        end
        s0 += e0_arg[i][1] * sin(a) + e0_arg[i][2] * cos(a)
    end

    for i in reverse(eachindex(e1_coefficent))
        a = 0.0;
        for j in 1:8
            a += e1_coefficent[i][j] * fa[j]
        end
        s1 += e1_arg[i][1] * sin(a) + e1_arg[i][2] * cos(a)
    end

    sec2rad(s0 + s1 * t)
end

"""
    equation_of_equinoxes_00(jd1, jd2, epsa, dpsi)

Return equation of equinoxes for given 2 part Julian date (TT), mean obliquity and nutation in longitude.

# Example

```jldoctest
julia> equation_of_equinoxes_00(2.4578265e6, 0.30440190993249416, 1.5, 1.7)
0.12025324854189404
```
"""
function equation_of_equinoxes_00(jd1, jd2, epsa, dpsi)
    dpsi * cos(epsa) + equation_of_equinoxes_complementary_terms(jd1, jd2)
end

function cio_locator(coeffs, terms0, terms1, terms2, terms3, terms4, jd1, jd2, x, y)
    t = ((jd1 - J2000) + jd2) / DAYS_PER_CENTURY

    fa = (
        mean_anomaly(luna, t),
        mean_anomaly(sun, t),
        mean_longitude_minus_lan(luna, t),
        mean_elongation(luna, t),
        mean_longitude_ascending_node(luna,t),
        mean_longitude(venus, t),
        mean_longitude(earth, t),
        general_precession_in_longitude(t),
    )

    w = collect(coeffs)

    for ct in reverse(terms0)
        a = 0.0
        for i in reverse(eachindex(fa))
            a += ct.coeffs[i] * fa[i]
        end
        s, c = sincos(a)
        w[1] += ct.sc * s + ct.cc * c
    end

    for ct in reverse(terms1)
        a = 0.0
        for i in reverse(eachindex(fa))
            a += ct.coeffs[i] * fa[i]
        end
        s, c = sincos(a)
        w[2] += ct.sc * s + ct.cc * c
    end

    for ct in reverse(terms2)
        a = 0.0
        for i in reverse(eachindex(fa))
            a += ct.coeffs[i] * fa[i]
        end
        s, c = sincos(a)
        w[3] += ct.sc * s + ct.cc * c
    end

    for ct in reverse(terms3)
        a = 0.0
        for i in reverse(eachindex(fa))
            a += ct.coeffs[i] * fa[i]
        end
        s, c = sincos(a)
        w[4] += ct.sc * s + ct.cc * c
    end

    for ct in reverse(terms4)
        a = 0.0
        for i in reverse(eachindex(fa))
            a += ct.coeffs[i] * fa[i]
        end
        s, c = sincos(a)
        w[5] += ct.sc * s + ct.cc * c
    end

    sec2rad((@evalpoly t w[1] w[2] w[3] w[4] w[5] w[6])) - x * y / 2.0
end

const s00(jd1, jd2, x, y) = cio_locator(CIO_COEFFS_00,
                                        CIO_TERMS_0TH_00,
                                        CIO_TERMS_1ST_00,
                                        CIO_TERMS_2ND_00,
                                        CIO_TERMS_3RD_00,
                                        CIO_TERMS_4TH_00,
                                        jd1, jd2, x, y)

const s06(jd1, jd2, x, y) = cio_locator(CIO_COEFFS_06,
                                        CIO_TERMS_0TH_06,
                                        CIO_TERMS_1ST_06,
                                        CIO_TERMS_2ND_06,
                                        CIO_TERMS_3RD_06,
                                        CIO_TERMS_4TH_06,
                                        jd1, jd2, x, y)

"""
    s00(jd1, jd2, x, y)

Returns Celestial Intermediate Origin(CIO) for a given 2-part Julian date (jd1, jd2)
with CIP coordinates (x, y). Compatible with IAU-2000 precession-nutation.

# Example

```jldoctest
julia> AstroBase.s00(2.4578265e6, 0.30434616919175345, 20, 50)
-500.00000000383193
```
"""
s00

"""
    s06(jd1, jd2, x, y)

Returns Celestial Intermediate Origin(CIO) for a given 2-part Julian date (jd1, jd2)
with CIP coordinates (x, y)

 # Example

julia> AstroBase.s06(2.4578265e6, 0.30434616919175345, 20, 50)
-500.00000000383193
```
"""
s06

"""
    nutation_00a(jd1, jd2)

Returns luni-solar and planetary nutations for a given 2 part Julian date(TT).

# Example

```jldoctest
julia> nutation_00a(2.4578265e6, 0.30440190993249416)
(-3.7589903391357206e-5, -3.6659049617818334e-5)
```
"""
function nutation_00a(jd1, jd2)
    t = ((jd1 - J2000) + jd2) / DAYS_PER_CENTURY

    el = mean_anomaly(luna, t)
    f  = mean_longitude_minus_lan(luna, t)
    om = mean_longitude_ascending_node(luna, t)

    apa  = general_precession_in_longitude(t)
    alme = mean_longitude(mercury, t)
    alve = mean_longitude(venus, t)
    alea = mean_longitude(earth, t)
    alma = mean_longitude(mars, t)
    alju = mean_longitude(jupiter, t)
    alsa = mean_longitude(saturn, t)
    alur = mean_longitude(uranus, t)

    elp = sec2rad(mod((@evalpoly t 1287104.79305 129596581.0481 -0.5532 0.000136 -0.00001149),
                      ARCSECONDS_IN_CIRCLE))
    d   = sec2rad(mod((@evalpoly t 1072260.70369 1602961601.2090 -6.3706 0.006593 -0.00003169),
                      ARCSECONDS_IN_CIRCLE))

    al   = mod2pi(2.35555598 + 8328.6914269554 * t)
    af   = mod2pi(1.627905234 + 8433.466158131 * t)
    ad   = mod2pi(5.198466741 + 7771.3771468121 * t)
    aom  = mod2pi(2.18243920 - 33.757045 * t)
    alne = mod2pi(5.321159000 + 3.8127774000 * t)

    dpls = 0.0
    dels = 0.0

    for i in reverse(eachindex(xls_nutation))
        arg = mod2pi(xls_nutation[i][1]  * el +
                    xls_nutation[i][2] * elp +
                    xls_nutation[i][3]  * f +
                    xls_nutation[i][4]  * d +
                    xls_nutation[i][5] * om)
        sarg, carg = sin(arg), cos(arg)
        dpls += (xls_nutation[i][6] + xls_nutation[i][7] * t) * sarg + xls_nutation[i][8] * carg
        dels += (xls_nutation[i][9] + xls_nutation[i][10] * t) * carg + xls_nutation[i][11] * sarg
    end

    dppl = 0.0
    depl = 0.0

    for i in reverse(eachindex(xpl_nutation))
        arg = mod2pi(xpl_nutation[i][1] * al +
             xpl_nutation[i][2] * af   +
             xpl_nutation[i][3] * ad   +
             xpl_nutation[i][4] * aom  +
             xpl_nutation[i][5] * alme +
             xpl_nutation[i][6] * alve +
             xpl_nutation[i][7] * alea +
             xpl_nutation[i][8] * alma +
             xpl_nutation[i][9] * alju +
             xpl_nutation[i][10] * alsa +
             xpl_nutation[i][11] * alur +
             xpl_nutation[i][12] * alne +
             xpl_nutation[i][13] * apa)
        sarg, carg = sin(arg), cos(arg)

        dppl += xpl_nutation[i][14] * sarg + xpl_nutation[i][15] * carg
        depl += xpl_nutation[i][16] * sarg + xpl_nutation[i][17] * carg
    end

    sec2rad(1e-7dpls) + sec2rad(1e-7dppl), sec2rad(1e-7dels) + sec2rad(1e-7depl)
end

"""
    nutation_a06(jd1, jd2)

Returns nutation, luni-solar and planetary for a given 2 part Julian date (TT).

# Example

```jldoctest
julia> nutation_a06(2.4578265e6, 0.30434616919175345)
(-3.758986477837638e-5, -3.665903094375903e-5)
```
"""
function nutation_a06(jd1, jd2)
    t = ((jd1 - J2000) + jd2) / DAYS_PER_CENTURY
    dp, de = nutation_00a(jd1, jd2)

    dp + dp * (0.4697e-6 - 2.7774e-6 * t), de + de * (-2.7774e-6 * t)
end

"""
    greenwich_mean_sidereal_time06(uta, utb, tta, ttb, rnpb)

Returns Greenwich mean sidereal time(radians) for given two, 2 part Julian dates (TT and UT1).

# Example

```jldoctest
julia> gst06(2.4579405e6, 0.0, 2.4579405e6, -0.0007966009351851851, rand(3,3))
3.738379033673736
```
"""
function greenwich_apparent_sidereal_time06(uta, utb, tta, ttb, rnpb)
    x, y = rnpb[1, 3], rnpb[2, 3]
    s = s06(tta, ttb, x, y)
    era = earth_rotation_angle(uta, utb)
    eors = equation_of_origins(rnpb,s)
    mod2pi(era - eors)
end

"""
    equation_of_equinoxes_94(jd1, jd2)

Returns equation of equinoxes for a given 2 part Julian date (TT). (IAU 1994 model)

# Example

```jldoctest
julia> equation_of_equinoxes_94(2.4578265e6, 0.30434616919175345)
-3.446039061100115e-5
```
"""
function equation_of_equinoxes_94(jd1, jd2)
    t = ((jd1 - J2000) + jd2) / DAYS_PER_CENTURY

    om = rem2pi(sec2rad(@evalpoly t 450160.280 -482890.539 7.455 0.008)
      + mod(-5.0 * t, 1.0) * 2pi, RoundNearest)

    dpsi, deps = nutation(jd1, jd2)
    eps0 = mean_obliquity_of_ecliptic(jd1, jd2)

    dpsi*cos(eps0) + sec2rad((0.00264 * sin(om) + 0.000063 * sin(om + om)))
end

"""
    nutation_matrix80(jd1, jd2)

Returns nutation matrix for a given 2 part Julian date (TDB).

# Example

```jldoctest
julia> nutation_matrix80(2.4578265e6, 0.30434616919175345)
3×3 Rotations.RotXZX{Float64}(0.409054, 3.75653e-5, -0.409017):
 1.0         -3.44666e-5  -1.494e-5
 3.44661e-5   1.0         -3.66564e-5
 1.49413e-5   3.66559e-5   1.0
 ```
"""
function nutation_matrix80(jd1, jd2)
    dpsi, deps = nutation(jd1, jd2)
    epsa = mean_obliquity_of_ecliptic(jd1, jd2)
    numat(epsa, dpsi, deps)
end

"""
    precession_nutation00(jd1, jd2, dpsi, deps)

Returns
    epsa
    rb
    rp
    rbp
    rn
    rbpn
 for a given 2 part Julian date (TT) and nutations.

# Example

```jldoctest
julia> pn00(2.4578265e6, 0.30434616919175345, 0.2, 0.2)
(0.40905374831590824, [1.0 7.07828e-8 -8.05622e-8; -7.07828e-8 1.0 -3.30604e-8; 8.05622e-8 3.30604e-8 1.0], [0.999991 0.00384587 0.00167106; -0.00384587 0.999993 -3.24146e-6; -0.00167106 -3.18524e-6 0.999999], [0.999991 0.00384594 0.00167098; -0.00384594 0.999993 -3.27421e-6; -0.00167098 -3.1523e-6 0.999999], [0.980067 0.162947 0.113657; -0.182279 0.965066 0.188206; -0.079019 -0.205172 0.975531], [0.979225 0.166314 0.11601; -0.186046 0.964433 0.187765; -0.080656 -0.205447 0.975339])
```
"""
function precession_nutation00(jd1, jd2, dpsi, deps)
    dpsipr, depspr = precession_rate_part_of_nutation(jd1, jd2)
    epsa = mean_obliquity_of_ecliptic(jd1, jd2) + depspr
    rb, rp, rbpw = bias_precession_matrix_00(jd1, jd2)
    rnw = numat(epsa, dpsi, deps)
    epsa, rb, rp, rbpw, rnw, rbpw * rnw
end

"""
    precession_nutation_a00(jd1, jd2, dpsi, deps)

Returns
    dpsi, deps
    epsa
    rb
    rp
    rbp
    rn
    rbpn
for a given 2 part Julian date (TT) and nutations.

# Example

```jldoctest
julia> precession_nutation_a00(2.4578265e6, 0.30434616919175345)
(-3.758986507815969e-5, -3.6659048454874034e-5, (0.40905374831590824, [1.0 7.07828e-8 -8.05622e-8; -7.07828e-8 1.0 -3.30604e-8; 8.05622e-8 3.30604e-8 1.0], [0.999991 0.00384587 0.00167106; -0.00384587 0.999993 -3.24146e-6; -0.00167106 -3.18524e-6 0.999999], [0.999991 0.00384594 0.00167098; -0.00384594 0.999993 -3.27421e-6; -0.00167098 -3.1523e-6 0.999999], [1.0 -3.44892e-5 -1.49498e-5; 3.44886e-5 1.0 -3.66593e-5; 1.4951e-5 3.66588e-5 1.0], [0.999991 0.00381152 0.00165589; -0.00381145 0.999993 -3.98758e-5; -0.00165603 3.35641e-5 0.999999]))
```
"""
function precession_nutation_a00(jd1, jd2)
    dpsi, deps = nutation_00a(jd1, jd2)
    dpsi, deps, precession_nutation00(jd1, jd2, dpsi, deps)
end

"""
    precession_nutation_b00(jd1, jd2, dpsi, deps)

Returns
    dpsi, deps
    epsa
    rb
    rp
    rbp
    rn
    rbpn
for a given 2 part Julian date (TT) and nutations.

# Example

```jldoctest
julia> precession_nutation_b00(2.4578265e6, 0.30434616919175345)
(-3.758913944197044e-5, -3.665743003046017e-5, (0.40905374831590824, [1.0 7.07828e-8 -8.05622e-8; -7.07828e-8 1.0 -3.30604e-8; 8.05622e-8 3.30604e-8 1.0], [0.999991 0.00384587 0.00167106; -0.00384587 0.999993 -3.24146e-6; -0.00167106 -3.18524e-6 0.999999], [0.999991 0.00384594 0.00167098; -0.00384594 0.999993 -3.27421e-6; -0.00167098 -3.1523e-6 0.999999], [1.0 -3.44885e-5 -1.49495e-5; 3.44879e-5 1.0 -3.66577e-5; 1.49508e-5 3.66572e-5 1.0], [0.999991 0.00381152 0.00165589; -0.00381146 0.999993 -3.98741e-5; -0.00165603 3.35625e-5 0.999999]))
```
"""
function precession_nutation_b00(jd1, jd2)
    dpsi, deps = nutation_00b(jd1, jd2)
    dpsi, deps, precession_nutation00(jd1, jd2, dpsi, deps)
end

"""
    precession_nutation_matrix_a00(jd1, jd2)

Returns classical NBP matrix for a given 2 part Julian date(TT)

# Example

```jldoctest
julia> precession_nutation_matrix_a00(2.4578265e6, 0.30434616919175345)
3×3 Rotations.RotMatrix{3,Float64,9}:
  0.999991    0.00381152   0.00165589
 -0.00381145  0.999993    -3.98758e-5
 -0.00165603  3.35641e-5   0.999999
```
"""
function precession_nutation_matrix_a00(jd1, jd2)
    precession_nutation_a00(jd1, jd2)[3][6]
end

"""
    precession_nutation_matrix_b00(jd1, jd2)

Returns classical NBP matrix for a given 2 part Julian date(TT)

# Example

```jldoctest
julia> precession_nutation_matrix_b00(2.4578265e6, 0.30434616919175345)
3×3 Rotations.RotMatrix{3,Float64,9}:
  0.999991    0.00381152   0.00165589
 -0.00381146  0.999993    -3.98741e-5
 -0.00165603  3.35625e-5   0.999999
```
"""
function precession_nutation_matrix_b00(jd1, jd2)
    precession_nutation_b00(jd1, jd2)[3][6]
end

"""
    precession_nutation_matrix_a06(jd1, jd2)

Returns classical NBP matrix for a given 2 part Julian date(TT)

# Example

```jldoctest
julia> precession_nutation_matrix_a06(2.4578265e6, 0.30434616919175345)
3×3 Rotations.RotMatrix{3,Float64,9}:
  0.999991    0.00381151   0.00165589
 -0.00381145  0.999993    -3.98694e-5
 -0.00165602  3.35578e-5   0.999999
```
"""
function precession_nutation_matrix_a06(jd1, jd2)
    gamb, phib, psib, epsa = precession_fukushima_williams06(jd1, jd2)
    dp, de = nutation_a06(jd1, jd2)
    fukushima_williams_matrix(gamb, phib, psib + dp, epsa + de)
end
"""
    s00a(jd1, jd2)

Returns the CIO locator (radians) for a given 2 part Julian date (TT).

# Example

```jldoctest
julia> s00a(2.4578265e6, 0.30434616919175345)
2.9183401976144617e-8
```
"""
function s00a(jd1, jd2)
    rbpn = precession_nutation_matrix_a00(jd1, jd2)
    s00(jd1, jd2, rbpn[1, 3], rbpn[2, 3])
end

"""
    s00b(jd1, jd2)

Returns the CIO locator (radians) for a given 2 part Julian date (TT).

# Example

```jldoctest
julia> s00b(2.4578265e6, 0.30434616919175345)
2.9182068810695328e-8
```
"""
function s00b(jd1, jd2)
    rbpn = precession_nutation_matrix_b00(jd1, jd2)
    s00(jd1, jd2, rbpn[1,3], rbpn[2,3])
end

"""
    nutation_matrix_day_a00(jd1, jd2)

Returns nutation matrix for a given 2 part Julian date(TT).

# Example

```jldoctest
julia> precession_nutation_a00(2.4578265e6, 0.30434616919175345)[3][5]
3×3 Rotations.RotXZX{Float64}(0.409054, 3.75891e-5, -0.409017):
 1.0         -3.44885e-5  -1.49495e-5
 3.44879e-5   1.0         -3.66577e-5
 1.49508e-5   3.66572e-5   1.0
```
"""
function nutation_matrix_day_a00(jd1, jd2)
    precession_nutation_a00(jd1, jd2)[3][5]
end

"""
    nutation_matrix_day_b00(jd1, jd2)

Returns nutation matrix for a given 2 part Julian date(TT).

# Example

```jldoctest
julia> precession_nutation_b00(2.4578265e6, 0.30434616919175345)[3][5]
3×3 Rotations.RotXZX{Float64}(0.409054, 3.75891e-5, -0.409017):
 1.0         -3.44885e-5  -1.49495e-5
 3.44879e-5   1.0         -3.66577e-5
 1.49508e-5   3.66572e-5   1.0
```
"""
function nutation_matrix_day_b00(jd1, jd2)
    precession_nutation_b00(jd1, jd2)[3][5]
end

"""
    nutation_matrix_day(jd1,jd2)

Returns nutation matrix for a given 2 part Julian date(TT).

# Example

```jldoctest
julia> precession_nutation_b00(2.4578265e6, 0.30434616919175345)[3][5]
3×3 Rotations.RotXZX{Float64}(0.409054, 3.75891e-5, -0.409017):
 1.0         -3.44885e-5  -1.49495e-5
 3.44879e-5   1.0         -3.66577e-5
 1.49508e-5   3.66572e-5   1.0
```
"""
function nutation_matrix_day(jd1, jd2)
    eps = mean_obliquity_of_ecliptic(jd1, jd2)
    dp, de = nutation_a06(jd1, jd2)

    numat(eps, dp, de)
end

"""
    xys00a(jd1, jd2)

Returns Celestial intermediate pole (x, y) and the CIO locator for a given 2 part Julian date. (TT)

# Example

```jldoctest
julia> xys00a(2.4578265e6, 0.30434616919175345)
(0.001655885227051451, -3.987617878202287e-5, 2.9183401976144617e-8)
```
"""
function xys00a(jd1, jd2)
    rbpn = precession_nutation_matrix_a00(jd1, jd2)
    rbpn[1,3], rbpn[2,3], s00(jd1, jd2, rbpn[1,3], rbpn[2,3])
end

"""
xys00a(jd1, jd2)

Returns Celestial intermediate pole (x, y) and the CIO locator for a given 2 part Julian date. (TT)

# Example

```jldoctest
julia> xys00b(2.4578265e6, 0.30434616919175345)
(0.001655885521808701, -3.987456146931851e-5, 2.9182068810695328e-8)
```
"""
function xys00b(jd1, jd2)
    rbpn = precession_nutation_matrix_b00(jd1, jd2)
    rbpn[1,3], rbpn[2,3], s00(jd1, jd2, rbpn[1,3], rbpn[2,3])
end

"""
    s06a(jd1, jd2)

Returns the CIO locator for a given 2 part Julian date (TT).

# Example

```jldoctest
julia> s06a(2.4578265e6, 0.30434616919175345)
2.917767075719633e-8
```
"""
function s06a(jd1, jd2)
    rnpb = precession_nutation_matrix_a06(jd1, jd2)
    s06(jd1, jd2, rnpb[1,3], rnpb[2,3])
end

"""
    equation_of_equinoxes_a00(jd1, jd2)

Returns equation of equinoxes for a given 2 part Julian date(TT).

# Example

```jldoctest
julia> equation_of_equinoxes_a00(2.4578265e6, 0.30434616919175345)
-3.4482238309105874e-5
```
"""
function equation_of_equinoxes_a00(jd1, jd2)
    dpsipr, depspr = precession_rate_part_of_nutation(jd1, jd2)
    epsa = mean_obliquity_of_ecliptic(jd1, jd2) + depspr
    dpsi, deps = nutation_00a(jd1, jd2)
    dpsi * cos(epsa) + equation_of_equinoxes_complementary_terms(jd1, jd2)
end

"""
    equation_of_equinoxes_b00(jd1, jd2)

Returns equation of equinoxes for a given 2 part Julian date(TT).

# Example

```jldoctest
julia> equation_of_equinoxes_b00(2.4578265e6, 0.30434616919175345)
-3.4482238309105874e-5
```
"""
function equation_of_equinoxes_b00(jd1, jd2)
    dpsipr, depspr = precession_rate_part_of_nutation(jd1, jd2)
    epsa = mean_obliquity_of_ecliptic(jd1, jd2) + depspr
    dpsi, deps = nutation_00b(jd1, jd2)
    dpsi * cos(epsa) + equation_of_equinoxes_complementary_terms(jd1, jd2)
end

"""
    greenwich_mean_sidereal_time_a06(uta, utb, tta, ttb, rnpb)

Returns greenwich apparent sidereal time for given two, 2 part Julian dates (UT1, TT).

# Example

```jldoctest

```
"""
function greenwich_mean_sidereal_time_a06(uta, utb, tta, ttb)
    rnpb = precession_nutation_matrix_a06(tta, ttb)
    greenwich_apparent_sidereal_time06(uta, utb, tta, ttb, rnpb)
end

"""
    equation_of_equinoxes_a06(jd1, jd2)

Returns equation of equinoxes for a given 2 part Julian date (TT).

# Example

```jldoctest
julia> equation_of_equinoxes_a06(2.4578265e6, 0.30434616919175345)
-3.448290610741367e-5
```
"""
function equation_of_equinoxes_a06(jd1, jd2)
    gst06a = greenwich_mean_sidereal_time_a06(0.0, 0.0, jd1, jd2)
    gmst06 = greenwich_mean_sidereal_time06(0.0, 0.0, jd1, jd2)

    rem2pi(gst06a - gmst06, RoundNearest)
end

"""
    greenwich_apparent_sidereal_time94(uta, utb)

Returns greenwich appart sidereal time for given 2 part Julian date (UT1).

# Example

```jldoctest
julia> greenwich_apparent_sidereal_time94(2.4578265e6, 0.30434616919175345)
4.9160197844443445
```
"""
function greenwich_apparent_sidereal_time94(uta, utb)
    gmst82 = greenwich_mean_sidereal_time82(uta, utb)
    eqeq94 = equation_of_equinoxes_94(uta, utb)
    mod2pi(gmst82 + eqeq94)
end

"""
    greenwich_apparent_sidereal_time_a00(uta, utb, tta, ttb)

Returns Greenwich apparent sidereal time for given two, 2 part Julian date. (UT1, TT)

# Example

```jldoctest
julia> greenwich_apparent_sidereal_time_a00(2.4579405e6, 0.0, 2.4579405e6, -0.0007966009351851851)
4.959632360790386
```
"""
function greenwich_apparent_sidereal_time_a00(uta, utb, tta, ttb)
    gmst00 = greenwich_mean_sidereal_time00(uta, utb, tta, ttb)
    ee00a = equation_of_equinoxes_a00(tta, ttb)
    mod2pi(gmst00 + ee00a)
end

"""
    greenwich_apparent_sidereal_time_b00(uta, utb, tta, ttb)

Returns Greenwich apparent sidereal time for given two, 2 part Julian date. (UT1, TT)

# Example

```jldoctest
julia> greenwich_apparent_sidereal_time_b00(2.4579405e6, 0.0, 2.4579405e6, -0.0007966009351851851)
4.95963236232993
```
"""
function greenwich_apparent_sidereal_time_b00(uta, utb)
    gmst00 = greenwich_mean_sidereal_time00(uta, utb, uta, utb)
    ee00b = equation_of_equinoxes_b00(uta, utb)
    mod2pi(gmst00 + ee00b)
end

"""
    greenwich_apparent_sidereal_time_a06(uta, utb, tta, ttb)

Returns Greenwich apparent sidereal time for given two, 2 part Julian dates. (UT1, TT)

# Example

```jldoctest
julia> greenwich_apparent_sidereal_time_a06(2.4579405e6, 0.0, 2.4579405e6, -0.0007966009351851851)
4.959632359298287
```
"""
function greenwich_apparent_sidereal_time_a06(uta, utb, tta, ttb)
    rnpb = precession_nutation_matrix_a06(tta, ttb)
    greenwich_apparent_sidereal_time06(uta, utb, tta, ttb, rnpb)
end

"""
    celestial_to_intermediate_frame_of_date(jd1, jd2, x, y)

Returns celestial to intermediate-frame-of-date matrix for a given 2 part Julian Date.

# Example

```jldoctest
julia> celestial_to_intermediate_frame_of_date(2.4578265e6, 0.30434616919175345, 0.2, 0.3)
3×3 Rotations.RotZYZ{Float64}(0.982794, 0.368863, -0.952794):
  0.977932    -0.0604048  0.2
 -0.00243134   0.953936   0.3
 -0.208909    -0.293866   0.932738
```
"""
function celestial_to_intermediate_frame_of_date(jd1, jd2, x, y)
    celestial_to_intermediate(x, y, s00(jd1, jd2, x, y))
end

"""
    celestial_to_intermediate_matrix(jd1, jd2, rbpn)

Returns celestial-to-intermediate matrix for given 2 part Julian date and celestial-to-true matrix.

# Example

```jldoctest
julia> celestial_to_intermediate_matrix(2.4578265e6, 0.30434616919175345, rand(3,3))
3×3 Rotations.RotZYZ{Float64}(0.921674, 1.10597, -0.729284):
  0.732866  -0.413437  0.54035
 -0.13652    0.688686  0.712092
 -0.666536  -0.595636  0.448272
```
"""
function celestial_to_intermediate_matrix(jd1, jd2, rbpn)
    x, y = rbpn[1,3], rbpn[2,3]
    celestial_to_intermediate_frame_of_date(jd1, jd2, x, y)
end

"""
    celestial_to_intermediate_matrix_a00(jd1, jd2)

Returns celestial-to-intermediate matrix for a given 2 part Julian date (TT).

# Example

```jldoctest
julia> celestial_to_intermediate_matrix_a00(2.4578265e6, 0.30434616919175345)
3×3 Rotations.RotZYZ{Float64}(-0.0240766, 0.00165637, 0.0240765):
  0.999999    6.21979e-8   0.00165589
  3.83181e-9  1.0         -3.98758e-5
 -0.00165589  3.98757e-5   0.999999
```
"""
function celestial_to_intermediate_matrix_a00(jd1, jd2)
    rbpn = precession_nutation_matrix_a00(jd1, jd2)
    celestial_to_intermediate_matrix(jd1, jd2, rbpn)
end

"""
    celestial_to_intermediate_matrix_b00(jd1, jd2)

Returns celestial-to-intermediate matrix for a given 2 part Julian date (TT).

# Example

```jldoctest
julia> celestial_to_intermediate_matrix_b00(2.4578265e6, 0.30434616919175345)
3×3 Rotations.RotZYZ{Float64}(-0.0240766, 0.00165637, 0.0240765):
0.999999    6.21979e-8   0.00165589
3.83181e-9  1.0         -3.98758e-5
-0.00165589  3.98757e-5   0.999999
```
"""
function celestial_to_intermediate_matrix_b00(jd1, jd2)
    rbpn = precession_nutation_matrix_b00(jd1, jd2)
    celestial_to_intermediate_matrix(jd1, jd2, rbpn)
end

"""
    celestial_to_terrestrial_matrix(tta, ttb, uta, utb, x, y, xp, yp)

Returns celestial to terrestrial matrix for given two, 2 part Julian dates (UT1, TT) and CIP coordinates.

# Example

```jldoctest
julia> celestial_to_terrestrial_matrix(2.4579405e6, 0.0, 2.4579405e6, -0.0007966009351851851, 0.1, 0.3, 0.5, 0.7)
3×3 Rotations.RotMatrix{3,Float64,9}:
  0.279858   0.75876   0.588186
 -0.66989   -0.284531  0.685777
  0.687697  -0.58594   0.428658
```
"""
function celestial_to_terrestrial_matrix(tta, ttb, uta, utb, x, y, xp, yp)
    rc2i = celestial_to_intermediate_frame_of_date(tta, ttb, x, y)
    era = earth_rotation_angle(uta, utb)
    sp = tio_locator(tta, ttb)

    rpom = polar_motion(xp, yp, sp)
    compose_rotation(rpom, angleaxis_to_dcm(-era, [0.0, 0.0, 1.0]), rc2i)
end

"""
    celestial_to_terrestrial_a00(tta, ttb, uta, utb, xp, yp)

Returns celestial to terrestrial matrix given two, 2 part Julian dates (UT1, TT) and coordinates of the pole.

# Example

```jldoctest
julia> celestial_to_terrestrial_a00(2.4579405e6, 0.0, 2.4579405e6, -0.0007966009351851851, 0.1, 0.3)
3×3 Rotations.RotMatrix{3,Float64,9}:
  0.235088    0.934797  0.266247
 -0.966879    0.196899  0.162409
  0.0993959  -0.295609  0.950124
```
"""
function celestial_to_terrestrial_a00(tta, ttb, uta, utb, xp, yp)
    rc2i = celestial_to_intermediate_matrix_a00(tta, ttb)
    era = earth_rotation_angle(uta, utb)
    sp = tio_locator(tta, ttb)

    rpom = polar_motion(xp, yp, sp)
    compose_rotation(rpom, angleaxis_to_dcm(-era, [0.0, 0.0, 1.0]), rc2i)
end

"""
    celestial_to_terrestrial_b00(tta, ttb, uta, utb, xp, yp)

Returns celestial to terrestrial matrix given two, 2 part Julian dates (UT1, TT) and coordinates of the pole.

# Example

```jldoctest
julia> celestial_to_terrestrial_b00(2.4579405e6, 0.0, 2.4579405e6, -0.0007966009351851851, 0.1, 0.3)
3×3 Rotations.RotMatrix{3,Float64,9}:
  0.235088    0.934797  0.266247
 -0.966879    0.196899  0.162409
  0.0993959  -0.295609  0.950124
```
"""
function celestial_to_terrestrial_b00(tta, ttb, uta, utb, xp, yp)
    rc2i = celestial_to_intermediate_matrix_b00(tta, ttb)
    era = earth_rotation_angle(uta, utb)

    rpom = polar_motion(xp, yp, 0.0)
    compose_rotation(rpom, angleaxis_to_dcm(-era, [0.0, 0.0, 1.0]), rc2i)
end

"""
    c2tpe(tta, ttb, uta, utb, dpsi, deps, xp, yp)

Returns the celestial to terrestrial matrix given two, 2 part Julian date (UT1, TT), the nutation and coordinates of the pole(radians).

# Example

```jldoctest
julia> c2tpe(2.4579405e6, 0.0, 2.4579405e6, -0.0007966009351851851, 0.1, 0.3, 0.5, 0.7)
3×3 Rotations.RotMatrix{3,Float64,9}:
  0.253143   0.777796  0.575284
 -0.676266  -0.28296   0.680146
  0.691797  -0.561219  0.454367
```
"""
function c2tpe(tta, ttb, uta, utb, dpsi, deps, xp, yp)
    epsa, rb, rp, rbp, rn, rbpn = precession_nutation00(tta, ttb, dpsi, deps)
    gmst = greenwich_mean_sidereal_time00(uta, utb, tta, ttb)
    ee = equation_of_equinoxes_00(tta, ttb, epsa, dpsi)
    sp = tio_locator(tta, ttb)

    rpom = polar_motion(xp, yp, sp)
    compose_rotation(rpom, angleaxis_to_dcm(-(gmst + ee), [0.0, 0.0, 1.0]), rbpn)
end

end


text:
using JSON
using StrTables

const VER = UInt32(1)

const inpname = "emoji_pretty.json"
const vers  = "master" # Julia used 0f0cf4ea8845eb52d26df2a48c3c31c3b8cad14e
const dpath = "https://raw.githubusercontent.com/iamcal/emoji-data/"

const fname = "emoji.dat"
const datapath = "../data"

const disp = [false]

# Get manual additions to the tables
include("../src/manual_emoji.jl")

str_to_uint32(str) = UInt32[ch%UInt32 for ch in str]

function make_tables(dpath, ver, fname)
    lname = joinpath(datapath, fname)
    if isfile(lname)
        println("Loaded: ", lname)
        src = lname
    else
        src = string(dpath, ver, '/', fname)
        download(src, lname)
        println("Saved to: ", lname)
    end
    emojidata = JSON.parsefile(lname)

    mandict = Dict(manual)
    symnam = String[ n for (n, v) in manual ]
    symval = Vector{UInt32}[ str_to_uint32(v) for (n, v) in manual ]
    ind = 0
    for emoji in emojidata
        # Make a vector of Chars out of hex data
        unified = emoji["unified"]
        unistr = UInt32[parse_hex(UInt32, str) for str in split(unified,'-')]
        strval = String(Char.(unistr))
        vecnames = emoji["short_names"]
        for name in vecnames
            manval = get(mandict, name, "")
            if manval == ""
                disp[] && println('#', ind += 1, '\t', unified, '\t', name)
                push!(symnam, name)
                push!(symval, unistr)
            else
                println(name, " => ", strval, "  overridden by manual entry: ", manval)
            end
        end
    end
    disp[] && println()

    # Get emoji names sorted
    srtnam = sortperm(symnam)
    srtval = symval[srtnam]

    # BMP characters
    l16 = Tuple{UInt16, UInt16}[]
    # non-BMP characters (in range 0x10000 - 0x1ffff)
    l32 = Tuple{UInt16, UInt16}[]
    # Vector of characters
    l2c = Tuple{String, UInt16}[]

    max2c = 1
    for i in eachindex(srtnam)
        chrs = srtval[i]
        len = length(chrs)
        if len > 1
            max2c = max(max2c, len)
            push!(l2c, (String(Char[ch for ch in chrs]), i))
        else
            ch = chrs[1]
            if ch > 0x1ffff
                error("Character $ch too large")
            elseif ch > 0x0ffff
                push!(l32, (ch%UInt16, i))
            else
                push!(l16, (ch%UInt16, i))
            end
        end
    end

    # We now have 3 vectors, for single BMP characters, for non-BMP characters, and for strings
    # each has the value and a index into the name table
    # We need to create a vector the same size as the name table, that gives the index
    # of into one of the three tables, in order to go from names to 1 or 2 output characters
    # We also need, for each of the 3 tables, a sorted vector that goes from the indices
    # in each table to the index into the name table (so that we can find multiple names for
    # each character)

    indvec = create_vector(UInt16, length(srtnam))
    vec16, ind16, base32 = sortsplit!(indvec, l16, 0)
    vec32, ind32, base2c = sortsplit!(indvec, l32, base32)
    vec2c, ind2c, basefn = sortsplit!(indvec, l2c, base2c)

    (VER, string(now()), src,
     base32%UInt32, base2c%UInt32, StrTable(symnam[srtnam]), indvec,
     vec16, ind16, vec32, ind32, StrTable(vec2c), ind2c, max2c%UInt32)
end

savfile = joinpath(datapath, fname)
if isfile(savfile)
    println("Tables already exist")
else
    tup = nothing
    println("Creating tables")
    try
        global tup
        tup = make_tables(dpath, vers, inpname)
    catch ex
        println(sprint(showerror, ex, catch_backtrace()))
    end
    println("Saving tables to ", savfile)
    StrTables.save(savfile, tup)
    println("Done")
end

ext: jl


text:
# Adjoint test for preconditioners
# Author: Ziyi Yin, ziyi.yin@gatech.edu
# Date: Dec 2020
#

### Model
model, _ , _ = setup_model()
nsrc = 4
_, _, recGeometry, _ = setup_geom(model; nsrc=nsrc)

tol = 5f-5
nrec = length(recGeometry.xloc[1])
nt = recGeometry.nt[1]

println("Test linearity of preconditioner operators")

######## Test left preconditioners

C = CumsumOp(nt,nsrc,nrec)
D = DiffOp(nt,nsrc,nrec)
FDiff = FractionalIntegrationOp(nt,nsrc,nrec,-0.8)
FInt = FractionalIntegrationOp(nt,nsrc,nrec,0.8)

left_list = [C,D,FDiff,FInt,C',D',FDiff',FInt']

data1 = Array{Array}(undef, nsrc)
data2 = Array{Array}(undef, nsrc)
for j=1:nsrc
    data1[j] = randn(Float32, nt, nrec)
    data2[j] = randn(Float32, nt, nrec)
end
d_obs1 = judiVector(recGeometry,data1)
d_obs2 = judiVector(recGeometry,data2)

for Op in left_list
    a = Op*(d_obs1+d_obs2)
    b = Op*d_obs1 + Op*d_obs2
    scalar = randn(Float32)
    c = scalar*(Op*d_obs1)
    d = Op*(scalar*d_obs1)
    @printf(" F * (a + b): %2.5e, F * a + F * b : %2.5e, relative error : %2.5e \n", norm(a), norm(b), norm(a-b)/norm(a))
    @printf(" a (F x): %2.5e, F a x : %2.5e, relative error : %2.5e \n", norm(c), norm(d), norm(c-d)/norm(c))
    @test isapprox(a, b, rtol=tol)
end

######## Test right preconditioners

M = MaskOp(model.n,randn(Float32,model.n))

right_list = [M, M']

dm1 = randn(Float32,prod(model.n))
dm2 = randn(Float32,prod(model.n))

for Op in right_list
    a = Op*(dm1+dm2)
    b = Op*dm1 + Op*dm2
    scalar = randn(Float32)
    c = scalar*(Op*dm1)
    d = Op*(scalar*dm1)
    @printf(" F * (a + b): %2.5e, F * a + F * b : %2.5e, relative error : %2.5e \n", norm(a), norm(b), norm(a-b)/norm(a))
    @printf(" a (F x): %2.5e, F a x : %2.5e, relative error : %2.5e \n", norm(c), norm(d), norm(c-d)/norm(c))
    @test isapprox(a, b, rtol=tol)
end


text:
using ATools
using Test
using Distributions
using Statistics

qx = Normal(100.0, 10.0)
qs = rand(qx, 100000)
dx = Normal(100.0, 5.0)
data = rand(dx, 100000)
xs=rand(Float64, length(qs))

include("math_test.jl")
include("util_test.jl")
include("histos_test.jl")
include("fits_test.jl")
include("raytracing_test.jl")
include("io_test.jl")
include("reco_tests.jl")

max_stars_repo_name: andLaing/ATools

max_issues_repo_name: andLaing/ATools

max_forks_repo_name: andLaing/ATools


text:
function get_stdlib_list(upstream::GitHub.Repo, branch_name::AbstractString; auth)
    return mktempdir() do temp_dir
        return cd(temp_dir) do
            list = Vector{StdlibInfo}(undef, 0)
            upstream_clone_url = "https://github.com/$(upstream.full_name).git"
            run(`git clone $(upstream_clone_url) UPSTREAM`)
            return cd("UPSTREAM") do
                run(`git checkout $(branch_name)`)
                my_regex_2 = r"^([\w]*?)_GIT_URL[\s]*?:=[\s]*?([^\n\r]*)$"
                my_regex_3 = r"^([\w]*?)\.version$"
                my_regex_4 = r"^([\w]*?)_BRANCH[\s]*?=[\s]*?([\w\-\.]*?)$"
                my_regex_5 = r"^([\w]*?)_(SHA[\w\d]*?)[\s]*?=[\s]*?([\w]*?)$"
                for filename in readdir("stdlib")
                    if isfile(joinpath("stdlib", filename))
                        m3 = match(my_regex_3, filename)
                        if !(m3 isa Nothing)
                            name = m3[1]
                            current_shas = Dict{String, String}()
                            branch_dict = Dict{Symbol, String}()
                            stdlib_version_lines = split(read(joinpath("stdlib", "$(name).version"), String), '\n')
                            for stdlib_version_line in stdlib_version_lines
                                m2 = match(my_regex_2, stdlib_version_line)
                                m4 = match(my_regex_4, stdlib_version_line)
                                m5 = match(my_regex_5, stdlib_version_line)
                                if !(m2 isa Nothing)
                                    branch_dict[:git_url] = strip(m2[2])
                                elseif !(m4 isa Nothing)
                                    branch_dict[:branch] = m4[2]
                                elseif !(m5 isa Nothing)
                                    current_shas[lowercase(m5[2])] = m5[3]
                                end
                            end
                            haskey(branch_dict, :branch) || error("$filename is missing a _BRANCH = entry")
                            haskey(branch_dict, :git_url) || error("$filename is missing a _GIT_URL := entry")
                            branch = branch_dict[:branch]
                            git_url = branch_dict[:git_url]
                            stdlib = StdlibInfo(;
                                name = name,
                                branch = branch,
                                git_url = git_url,
                                current_shas = current_shas
                            )
                            push!(list, stdlib)
                        end
                    end
                end
                return list
            end # cd("UPSTREAM")
        end # cd(temp_dir)
    end # mktempdir()
end


text:
import TransferEntropy: transferentropy

"""
    transferentropy(source::AbstractArray{<:Real, 1},
        response::AbstractArray{<:Real, 1},
        k::Int, l::Int, m::Int,
        binning_scheme::Union{RectangularBinning, Vector{RectangularBinning}};
        summary_statistic::Function = StatsBase.mean,
        η::Int = 1, τ::Int = 1,
        estimator::BinningTransferEntropyEstimator = VisitationFrequency(b = 2)) -> FLoat64

Convenience function for computing transfer entropy from a `source` time series to a
`response` time series.

Transfer entropy is computed using the provided `estimator` over
a `k + l + m`-dimensional [generalised delay reconstruction](@ref custom_delay_reconstruction)
of the input data. The embedding delay used for the reconstructions is `τ`,
and the prediction lag is `η`. Transfer entropy is computed using the following conventions
for the [generalised delay reconstructions](@ref custom_delay_reconstruction).

```math
\\begin{align}
T_f^{(k)} &= \\{(T(t+\\eta_k), \\ldots, T(t+\\eta_2), T(t+\\eta_1)) \\} \\\\
T_{pp}^{(l)} &= \\{ (T(t), T(t-\\tau_1), T(t-\\tau_2), \\ldots, T(t - \\tau_{l - 1})) \\} \\\\
S_{pp}^{(m)} &= \\{ (S(t), S(t-\\tau_1), S(t-\\tau_2), \\ldots, S(t-\\tau_{m - 1})) \\}
\\end{align}
```

where ``T_f`` denotes the `k`-dimensional set of vectors furnishing the future states of ``T``,
``T_{pp}`` denotes the `l`-dimensional set of vectors furnishing the past and present states of ``T``,
and ``S_{pp}`` denotes the `m`-dimensional set of vectors furnishing the past and present of ``S``.
``\\eta`` is the prediction lag. This convenience function uses ``\\tau_1`` = `τ`,
``\\tau_2`` = `2*τ`, ``\\tau_3`` = `3*τ`, and so on.

Combined, we get the generalised embedding ``\\mathbb{E} = (T_f^{(k)}, T_{pp}^{(l)}, S_{pp}^{(m)})``.
TE is then computed as

```math
\\begin{align}
TE_{S \\rightarrow T} = \\int_{\\mathbb{E}} P(T_f, T_{pp}, S_{pp}) \\log_{b}{\\left(\\frac{P(T_f | T_{pp}, S_{pp})}{P(T_f | T_{pp})}\\right)}
\\end{align}
```

The probabilities needed for computing ``TE_{S \\rightarrow T}`` are estimated
over a discretization of the state space reconstruction, and this discretization is
dictated by the provided `binning_scheme`(s). If there is more than one binning scheme,
then the transfer entropy is summarised using `summary_statistic` over the partitions,
and a single value for the transfer entropy is returned.

## Arguments

- **`source`**: The source data series.
- **`target`**: The target data series.
- **`k`**: The dimension of the ``T_{f}`` component of the embedding.
- **`l`**: The dimension of the ``T_{pp}`` component of the embedding.
- **`m`**: The dimension of the ``S_{pp}`` component of the embedding.
- **`binning_scheme`**: The binning scheme(s) used to construct the partitions
    over which TE is computed. Must be either one or several instances of
    `RectangularBinning`s (provided as a vector). TE is computed for each
    of the resulting partitions.

## Keyword arguments

- **`τ::Int = 1`**: The embedding lag (lags the ``T_{pp}`` component of the embedding).
- **`η::Int = 1`**: The prediction lag.
- **`estimator::BinningTransferEntropyEstimator = VisitationFrequency()`**: The
    transfer entropy estimator to use.

## Returns

A single value for the transfer entropy.
"""
function transferentropy(source::AbstractArray{<:Real, 1},
        response::AbstractArray{<:Real, 1},
        k::Int, l::Int, m::Int,
        binning_scheme::Union{RectangularBinning, Vector{RectangularBinning}};
        summary_statistic::Function = StatsBase.mean,
        η = 1, τ = 1,
        estimator = VisitationFrequency(b = 2))

    k + l + m >= 3 || throw(ArgumentError("`dim = k + l + m` must be 3 or higher for regular TE"))

    pts, vars = te_embed(source, response, k, l, m, η = η, τ = τ)

    # Compute TE over different partitions
    # ====================================
    bs = binning_scheme isa Vector{RectangularBinning} ? binning_scheme : [binning_scheme]
    tes = map(binscheme -> transferentropy(pts, vars, binscheme, estimator), bs)

    return summary_statistic(tes)
end

"""
    transferentropy(source::AbstractArray{<:Real, 1},
        response::AbstractArray{<:Real, 1},
        cond::AbstractArray{<:Real, 1},
        k::Int, l::Int, m::Int,
        binning_scheme::Union{RectangularBinning, Vector{RectangularBinning}};
        summary_statistic::Function = StatsBase.mean,
        η::Int = 1, τ::Int = 1,
        estimator::BinningTransferEntropyEstimator = VisitationFrequency(b = 2)) -> FLoat64

Convenience function for computing transfer entropy from a `source` time series to a
`response` time series, conditioned on a third time series `cond`.

Transfer entropy is computed using the provided `estimator` over
a `k + l + m + n`-dimensional [generalised delay reconstruction](@ref custom_delay_reconstruction)
of the input data. The embedding delay used for the reconstructions is `τ`,
and the prediction lag is `η`. Transfer entropy is computed using the following conventions
for the [generalised delay reconstructions](@ref custom_delay_reconstruction).

```math
\\begin{align}
T_f^{(k)} &= \\{(T(t+\\eta_k), \\ldots, T(t+\\eta_2), T(t+\\eta_1)) \\} \\\\
T_{pp}^{(l)} &= \\{ (T(t), T(t-\\tau_1), T(t-\\tau_2), \\ldots, T(t - \\tau_{l - 1})) \\} \\\\
S_{pp}^{(m)} &= \\{ (S(t), S(t-\\tau_1), S(t-\\tau_2), \\ldots, S(t-\\tau_{m - 1})) \\} \\\\
C_{pp}^{(n)} &= \\{ (C(t), C(t-\\tau_1), C(t-\\tau_2), \\ldots, C(t-\\tau_{n - 1})) \\}
\\end{align}
```

where ``T_f`` denotes the `k`-dimensional set of vectors furnishing the future states of ``T``,
``T_{pp}`` denotes the `l`-dimensional set of vectors furnishing the past and present states of ``T``,
``S_{pp}`` denotes the `m`-dimensional set of vectors furnishing the past and present of ``S``,
and ``C_{pp}`` denotes the `n`-dimensional set of vectors furnishing the past and present of ``C``.
``\\eta`` is the prediction lag. This convenience function uses ``\\tau_1`` = `τ`,
``\\tau_2`` = `2*τ`, ``\\tau_3`` = `3*τ`, and so on.

Combined, we get the generalised embedding  ``\\mathbb{E} = (T_f^{(k)}, T_{pp}^{(l)}, S_{pp}^{(m)}, C_{pp}^{(n)})``.
TE is then computed as

```math
\\begin{align}
TE_{S \\rightarrow T | C} = \\int_{\\mathbb{E}} P(T_f, T_{pp}, S_{pp}, C_{pp}) \\log_{b}{\\left(\\frac{P(T_f | T_{pp}, S_{pp}, C_{pp})}{P(T_f | T_{pp}, C_{pp})}\\right)}
\\end{align}
```

The probabilities needed for computing ``TE_{S \\rightarrow T | C}`` are estimated
over a discretization of the state space reconstruction, and this discretization is
dictated by the provided `binning_scheme`(s). If there is more than one binning scheme,
then the transfer entropy is summarised using `summary_statistic` over the partitions,
and a single value for the transfer entropy is returned.

## Arguments

- **`source`**: The source data series.
- **`target`**: The target data series.
- **`k::Int`**: The dimension of the ``T_{f}`` component of the embedding.
- **`l::Int`**: The dimension of the ``T_{pp}`` component of the embedding.
- **`m::Int`**: The dimension of the ``S_{pp}`` component of the embedding.
- **`n::Int`**: The dimension of the ``C_{pp}`` component of the embedding.
- **`binning_scheme`**: The binning scheme(s) used to construct the partitions
    over which TE is computed. Must be either one or several instances of
    `RectangularBinning`s (provided as a vector). TE is computed for each
    of the resulting partitions.

## Keyword arguments

- **`τ::Int = 1`**: The embedding lag (lags the ``T_{pp}`` component of the embedding).
- **`η::Int = 1`**: The prediction lag.
- **`estimator::BinningTransferEntropyEstimator = VisitationFrequency()`**: The
    transfer entropy estimator to use.

## Returns

A single value for the transfer entropy.
"""
function transferentropy(source::AbstractArray{<:Real, 1},
        response::AbstractArray{<:Real, 1},
        cond::AbstractArray{<:Real, 1},
        k::Int, l::Int, m::Int, n::Int,
        binning_scheme::Union{RectangularBinning, Vector{RectangularBinning}};
        summary_statistic::Function = StatsBase.mean,
        η = 1, τ = 1,
        estimator = VisitationFrequency(b = 2))

    k + l + m + n >= 4 || throw(ArgumentError("`dim = k + l + m + n` must be 4 or higher for conditional TE"))

    pts, vars = te_embed(source, response, cond, k, l, m, n, τ = τ, η = η)

    # Compute TE over different partitions
    # ====================================
    bs = binning_scheme isa Vector{RectangularBinning} ? binning_scheme : [binning_scheme]
    tes = map(binscheme -> transferentropy(pts, vars, binscheme, estimator), bs)

    return summary_statistic(tes)
end

export transferentropy


text:
using Test

# each line is "Step X must be completed before step Y can begin"
parseline(line) = (line[6], line[37])

test_input = parseline.(split(strip("""
Step C must be finished before step A can begin.
Step C must be finished before step F can begin.
Step A must be finished before step B can begin.
Step A must be finished before step D can begin.
Step B must be finished before step E can begin.
Step D must be finished before step E can begin.
Step F must be finished before step E can begin.
"""), '\n'))

input = parseline.(readlines("07.input"))

function forward(input)
    edges = Dict{Char,Vector{Char}}()
    for (before, after) in input
        push!(get!(edges, before, Char[]), after)
    end
    edges
end

backward(input) = forward(reverse(i) for i in input)

mutable struct JobQueue
    f_edges::Dict{Char,Vector{Char}}
    b_edges::Dict{Char,Set{Char}}
    queue::Vector{Char}
end

function JobQueue(input)
    f_edges = forward(input)
    b_edges = Dict(t=>Set(f) for (t,f) in backward(input))
    queue = sort!(collect(setdiff(keys(f_edges), keys(b_edges))), rev=true)
    JobQueue(f_edges, b_edges, queue)
end

start!(jq::JobQueue) = pop!(jq.queue)
function finish!(jq::JobQueue, job)
    for dep in get(jq.f_edges, job, ())
        if isempty(delete!(jq.b_edges[dep], job))
            push!(jq.queue, dep)
            # could also delete dep from back edges but not necessary
        end
    end

    sort!(jq.queue, rev=true)
    job
end
Base.isempty(jq::JobQueue) = isempty(jq.queue)

Base.iterate(jq::JobQueue) = iterate(jq, 0)
Base.iterate(jq::JobQueue, state) =
    isempty(jq) ? nothing : (finish!(jq, start!(jq)), state+1)

Base.length(jq::JobQueue) = length(union(keys(jq.f_edges), keys(jq.b_edges)))
Base.eltype(::Type{JobQueue}) = Char

star1(input) = collect(JobQueue(input))
@test String(star1(test_input)) == "CABDFE"

String(star1(input))

# star 2: two twists.  jobs have a TIME associated (61 seconds for A, 62 for B,
# etc.), and there are multiple workers who can work in parallel.  the strategy
# is very similar, except that instead of just taking jobs straight out of the
# queue, need to have an additional step where each worker takes a job and then
# reports when they'll be done.  then advance the time to the next worker to
# finish.

# what to do next depends on whether there are workers free and jobs available.
#
# while there are jobs in the master queue and the workers aren't all busy,
# start a job and hand it to the workers.
#
# if there are no jobs available, take the next job the workers will finish, and
# finish it in the main queue (and update the current time)
#
# if there are no workers free, same as above...
mutable struct MultiWorkerQueue
    jobs::JobQueue
    working_on::Vector{Char}
    t_done::Vector{Int}
    t::Int
    job_time::Function
end

function MultiWorkerQueue(input, n_worker::Int, job_time::Function)
    working_on = fill('.', n_worker)
    t_done = fill(typemax(Int), n_worker)
    jobs = JobQueue(input)
    MultiWorkerQueue(jobs, working_on, t_done, 0, job_time)
end

function Base.iterate(mwq::MultiWorkerQueue, state)
    println("t = $(mwq.t)")

    # done when no more jobs and all workers idle
    isempty(mwq.jobs) && all(isequal('.'), mwq.working_on) && return nothing

    # while there are idle workers and jobs available, fill the work queues
    for worker in findall(isequal('.'), mwq.working_on)
        isempty(mwq.jobs) && break
        job = start!(mwq.jobs)
        mwq.working_on[worker] = job
        mwq.t_done[worker] = mwq.t + mwq.job_time(job)
        println("  queueing job $job for worker $worker (done at t=$(mwq.t_done[worker]))")
    end

    # find the next job to be finished, increment the time, and return it
    mwq.t, worker = findmin(mwq.t_done)
    job = mwq.working_on[worker]
    finish!(mwq.jobs, job)
    mwq.t_done[worker] = typemax(Int)
    mwq.working_on[worker] = '.'
    return job, mwq.t
end

Base.iterate(mwq::MultiWorkerQueue) = iterate(mwq::MultiWorkerQueue, 0)
Base.length(mwq::MultiWorkerQueue) = length(mwq.jobs)

function star2(input; n_worker=5, job_t = (c) -> c-'A'+61)
    jobs = MultiWorkerQueue(input, n_worker, job_t)
    jobs_done = collect(jobs)
    jobs_done, jobs.t
end


text:
#-----------------------------------------------------------------------# Weight
@recipe function f(wt::Weight; nobs=50)
    xguide --> "n"
    yguide --> "w(n)"
    label --> name(wt)
    xlims --> (0, nobs + 1)
    ylims --> (0, 1.02)
    linewidth --> 2
    seriestype --> :scatter
    wt.(1:nobs)
end

#-----------------------------------------------------------------------# Fallback
@recipe function f(o::OnlineStat)
    legend --> false
    axis --> false
    grid --> false
    ylims --> (0, 1)
    xlims --> (0, 1)
    annotations --> [(.5, .75, name(o) * ":"), (.5, .4, string(value(o)))]
    zeros(0)
end

#-----------------------------------------------------------------------------# OrderStats
@recipe function f(o::OrderStats)
    label --> "Approximate CDF via OrderStats"
    legend --> :topleft
    xguide --> "x"
    yguide --> "P(X ≤ x)"
    a, b = value(o.ex)
    v = vcat(a, value(o), b)
    k = length(v)
    v, [1:k] ./ k
end

#-----------------------------------------------------------------------# residual plot
@recipe function f(o::OnlineStat{VectorOb}, x::AbstractMatrix, y::AbstractVector)
    ylab --> "Residual"
    xlab --> "Observation Index"
    legend --> false
    @series begin
        linetype --> :scatter
        ŷ = predict(o, x)
        eachindex(y), y - ŷ
    end
end

@recipe function f(o::OnlineStat{XY})
    ylab --> "beta_j"
    xlab --> "j"
    seriestype --> :scatter
    coef(o)
end

#-----------------------------------------------------------------------# StatCollection
@recipe function f(s::OnlineStatsBase.StatCollection)
    if :layout in keys(plotattributes)
        for stat in s.stats
            @series begin stat end
        end
    else  # hack to ensure series aren't sent to wrong subplots
        layout --> length(s.stats)
        for i in eachindex(s.stats)
            @series begin
                subplot --> i
                s.stats[i]
            end
        end
    end
end

#-----------------------------------------------------------------------# GroupBy
@recipe function f(o::GroupBy{T, <:Union{KHist, Hist}}) where {T}
    sort!(o.value)
    link --> :all
    for (k, v) in pairs(o.value)
        @series begin
            label --> k
            v
        end
    end
end
@recipe function f(o::GroupBy)
    sort!(o.value)
    link --> :all
    collect(keys(o.value)), value.(collect(values(o.value)))
end

#-----------------------------------------------------------------------# StatLag
@recipe function f(o::StatLag)
    layout --> length(o.value)
    for i in 1:length(o.value)
        @series begin
            subplot --> i
            label --> "Current - $(i-1)"
            o.value[i]
        end
    end
end

#-----------------------------------------------------------------------# AutoCov
@recipe function f(o::AutoCov)
    xguide --> "Lag"
    yguide --> "Autocorrelation"
    ylim --> (0, 1)
    seriestype --> :scatter
    autocor(o)
end

#-----------------------------------------------------------------------# CovMatrix
@recipe function f(o::CovMatrix; corr = false)
    seriestype --> :heatmap
    corr ? cor(o) : cov(o)
end

#-----------------------------------------------------------------------# Histograms
@recipe function f(o::HistogramStat; normalize=true)
    e, c = edges(o), counts(o)
    inds = findfirst(x -> x > 0, c):findlast(x -> x > 0, c)
    normalize --> normalize
    Histogram(e[vcat(inds, inds[end] + 1)], c[inds], _closed(o))
end

_closed(o::Hist) = o.left ? :left : :right
_closed(o) = :left

#-----------------------------------------------------------------------# CountMap
@recipe function f(o::CountMap, kys = keys(o); sortby = :keys)
    seriestype --> :bar
    kys = collect(kys)
    vls = [o.value[ky] for ky in kys]
    sortby in [:keys, :values] || @warn("sortby = :$sortby not recognized")
    sp = sortby == :keys ? sortperm(kys) : sortperm(vls)
    x, y = string.(kys[sp]), vls[sp]
    hover --> ["($xi, $yi)" for (xi,yi) in zip(x, y)], :quiet
    label --> "Count"
    x, y
end

#-----------------------------------------------------------------------------# Partition
_xguide(o::Partition) = "Nobs"
_xguide(o::Union{IndexedPartition, KIndexedPartition}) = "Index"

__middle(x) = x
__middle(x::Tuple) = middle(x...)

@recipe function f(o::Union{Partition, IndexedPartition, KIndexedPartition}; type=1)
    if type === 1
        @series begin
            t = o isa Partition ? "" : "Joint Distribution"
            title --> t
            xguide --> _xguide(o)
            o.parts
        end
    elseif type === 2
        @series begin
            title --> "Stat Over Entire Index"
            seriestype --> :bar
            orientation --> :h
            mergestats(o)
        end
    elseif type === 3
        @series begin
            l = o isa Partition ? "Nobs / Part" : "Index Distribution"
            label --> "Index Distribution"
            seriestype := :bar
            linewidth --> 0
            line_alpha --> 0
            data = [__middle(loc) => nobs(o) for (loc,o) in o.parts]
            first.(data), last.(data)
        end
    else
        error("type must be 1, 2, or 3.")
    end
end

@recipe function f(parts::T) where T <: AbstractVector{<:Pair{<:Any, <:OnlineStat}}
    seriestype  --> _seriestype(parts)
    seriesalpha --> _alpha(parts)
    linewidth   --> _linewidth(parts)
    linealpha   --> _linealpha(parts)
    line_z      --> _line_z(parts)
    fill_z      --> _fill_z(parts)
    ylims       --> _ylims(parts)
    xlims       --> _xlims(parts)
    yguide      --> _yguide(parts)

    hover --> string.(nobs.(last.(parts))), :quiet

    group = _group(parts)
    x, y = xy(parts)

    if isnothing(group)
        label --> _label(parts)
        return x, y
    else
        for k in unique(group)
            @series begin
                label --> _label(parts) * ": $k"
                x[group .== k], y[group .== k]
            end
        end
    end
end

@recipe function f(parts::AbstractVector{<:Pair{<:Any, <:Union{Series, Group}}})
    nstats = length(parts[1][2].stats)
    for i in 1:nstats
        @series begin
            map(x -> x[1] => x[2].stats[i], parts)
        end
    end
end

function xy(parts::AbstractVector{<:Pair{<:Any, <:OnlineStat}})
    data = xy.(parts)
    x, y = vcat(first.(data)...), vcat(last.(data)...)
end

_yguide(::Any) = ""
_yguide(::AbstractArray{<:Pair{<:Any, <:CountMap}}) = "Probability"

_seriestype(::Any) = :path
_seriestype(::AbstractVector{<:Pair{<:TwoThings, <:Union{Extrema, HistogramStat, CountMap, ProbMap, Counter}}}) = :shape
_seriestype(::AbstractVector{<:Pair{<:Number, <:Mean}}) = :scatter

_linealpha(parts) = 1
_linealpha(::AbstractVector{<:Pair{<:TwoThings, <:Union{Extrema,HistogramStat,CountMap}}}) = 0

_line_z(parts) = nothing
_line_z(parts::AbstractVector{<:Pair{<:Number, <:HistogramStat}}) = repeat(vcat(counts.(last.(parts))...), inner=3)

_alpha(parts) = 1
_alpha(::AbstractVector{<:Pair{<:Any, <:Extrema}}) = .5

_linewidth(parts) = 1
_linewidth(::AbstractVector{<:Pair{<:TwoThings, <:Union{Mean, Variance}}}) = 2
_linewidth(::AbstractVector{<:Pair{<:TwoThings, <:Union{Extrema,HistogramStat,CountMap}}}) = 0
_linewidth(::AbstractVector{<:Pair{<:Number, <:Extrema}}) = 2
_linewidth(::AbstractVector{<:Pair{<:Number, <:Union{HistogramStat,CountMap}}}) = 4

_label(parts) = name(parts[1][2], false, false)
_label(::AbstractVector{<:Pair{<:TwoThings, <:Moments}}) = ["m1" "m2" "m3" "m4"]
_label(::AbstractVector{<:Pair{<:Any, <:HistogramStat}}) = ""
_label(::AbstractVector{<:Pair{<:Any, <:Counter{T}}}) where {T} = "Count of $T"

_fill_z(parts) = nothing
function _fill_z(parts::AbstractVector{<:Pair{<:TwoThings, <:HistogramStat}})
    z = vcat(map(x -> counts(x[2]), parts)...)
    z[z .> 0]
end

_ylims(parts) = (-Inf, Inf)
_ylims(parts::AbstractVector{<:Pair{<:TwoThings, <:Union{CountMap, ProbMap, Counter}}}) = (0, Inf)

_xlims(parts) = (-Inf, Inf)
_xlims(parts::AbstractVector{<:Pair{<:TwoThings{<:Number, <:Number}, <:Any}}) = parts[1][1][1], parts[end][1][2]

_group(parts) = nothing
function _group(parts::AbstractVector{<:Pair{<:TwoThings, <:Union{CountMap, ProbMap}}})
    out = map(x -> collect(keys(sort!(value(x[2])))), parts)
    out = map(x -> repeat(x,inner=5), out)
    string.(vcat(out...))
end
function _group(parts::AbstractVector{<:Pair{<:Number, <:CountMap}})
    out = map(x -> collect(keys(sort!(value(x[2])))), parts)
    string.(repeat(vcat(out...), inner=3))
end

#-----------------------------------------------------------------------------# xy Mean and Variance
function xy(part::Pair{<:TwoThings, <:Union{Mean, Variance, Counter}})
    (a,b), o = part
    [a, b, b], [value(o), value(o), NaN]
end
xy(part::Pair{<:Number, <:Union{Mean, Variance}}) = part[1], value(part[2])
#-----------------------------------------------------------------------------# xy Counter
function xy(part::Pair{<:TwoThings, <:Counter})
    (a,b), o = part
    [a, a, b, b, b], [0, value(o), value(o), 0, NaN]
end
xy(part::Pair{<:Number, <:Counter}) = part[1], value(part[2])
#-----------------------------------------------------------------------------# xy Extrema
function xy(part::Pair{<:TwoThings, <:Extrema})
    (a,b), o = part
    low, high = extrema(o)
    [a, a, b, b, b], [low, high, high, low, NaN]
end
function xy(part::Pair{<:Number, <:Extrema})
    loc, o = part
    low, high = extrema(o)
    [loc, loc, loc], [low, high, NaN]
end
#-----------------------------------------------------------------------------# xy Moments
function xy(part::Pair{<:TwoThings, <:Moments})
    (a,b), o = part
    hcat([[a, b, b] for _ in 1:4]...), hcat([[v,v,NaN] for v in value(o)]...)
end
function xy(part::Pair{<:Number, <:Moments})
    loc, o = part
    [loc loc loc loc], hcat(value(o)...)
end
#-----------------------------------------------------------------------------# xy HistogramStat
function xy(part::Pair{<:TwoThings, <:HistogramStat})
    (a,b), o = part
    x, y = [], Float64[]
    edg = edges(o)
    cnts = counts(o)
    for i in eachindex(cnts)
        if cnts[i] > 0
            append!(x, [a, a, b, b, b])
            append!(y, [edg[i], edg[i+1], edg[i+1], edg[i], NaN])
        end
    end
    x, y
end
function xy(part::Pair{<:Number, <:HistogramStat})
    loc, o = part
    x, y = [], Float64[]
    edg = edges(o)
    cnts = counts(o)
    for i in eachindex(cnts)
        if cnts[i] > 0
            append!(x, [loc, loc, loc])
            append!(y, [edg[i], edg[i+1], NaN])
        end
    end
    x, y
end

#-----------------------------------------------------------------------------# xy CountMap
function xy(part::Pair{<:TwoThings, <:CountMap})
    (a,b), o = part
    x, y = [], Float64[]
    count = 0
    for (k,v) in sort!(o.value)
        append!(x, [a, a, b, b, b])
        append!(y, [count, count + v, count + v, count, NaN] ./ nobs(o))
        count += v
    end
    x, y
end
function xy(part::Pair{<:Number, <:CountMap})
    loc, o = part
    sort!(o)
    x, y = [], Float64[]
    count = 0
    for (k,v) in sort!(o.value)
        append!(x, [loc, loc, loc])
        append!(y, [count, count + v, NaN] ./ nobs(o))
        count += v
    end
    x, y
end

#-----------------------------------------------------------------------# NBClassifier
@recipe function f(o::NBClassifier)
    kys = collect(keys(o))
    layout --> nvars(o) + 1
    for j in 1:nvars(o)
        stats = o[j]
        for (i, s) in enumerate(stats)
            @series begin
                title --> "Var $j"
                legend --> false
                subplot --> j
                s
            end
        end
    end
    @series begin
        subplot --> nvars(o) + 1
        label --> string.(reshape(kys, 1, length(kys)))
        framestyle := :none
        zeros(0, length(kys))
    end
end


text:
using SVMLightWriter
using SparseArrays
using Test

x = [0 2; 3.1 4.2]
x_sp = sparse(x)
y = [1, -1]
io = IOBuffer()
str = "1 1:3.1\n-1 0:2 1:4.2\n"

dump_svmlight_file(x, y, io)
@test String(take!(io)) == str

dump_svmlight_file(x, y, "tmp")
@test read("tmp", String) == str
rm("tmp", force = true)

dump_svmlight_file(x_sp, y, io)
@test String(take!(io)) == str


text:
function create_spline(time, status::BitArray{1}; n_inner_knots=5)
    ## construct cubic B-spline with n_inner_knots equidistant inner knots at quantiles of observed failure times
    ev = time[status]
    outer = (0.0, maximum(ev))
    knots = quantile(ev, 0.0:(1.0/(n_inner_knots+1)):1.0)[2:(end-1)]

    ## the weights don't matter here since we are only evaluating the basis functions
    Spline.spline(outer, knots, ones(Float64, length(knots) + 4), 3)
end

 ## evaluate basis functions and integrals at each event time
eval_spline(data::Data, time::Vector{Float64}=data.time) = Spline.eval_basis(data.spline, time, time)
eval_spline(spline::Spline.Spline1D, time::Vector{Float64}) = Spline.eval_basis(spline, time, time)


text:
import DelimitedFiles
import Random

"""
    PityTimerSamples

Stochastic parameters that determine pack opening behavior.  Wildcards in packs
are assumed to have a base rate that increases by an increment each time one isn't opened
in a pack.  For the rare/mythic slot, the remaining probability not occupied by the
wildcard pity timers is divided into rare or mythic pulls by the `rare_card` parameter. This
representation is used to incorporate uncertainty in WC rates into pack estimates.
"""
struct PityTimerSamples
 common_wild_card_base::Vector{Float64}
 common_wild_card_increment::Vector{Float64}
 uncommon_wild_card_base::Vector{Float64}
 uncommon_wild_card_increment::Vector{Float64}
 rare_card::Vector{Float64}
 mythic_card::Vector{Float64}
 rare_wild_card_base::Vector{Float64}
 rare_wild_card_increment::Vector{Float64}
 mythic_wild_card_base::Vector{Float64}
 mythic_wild_card_increment::Vector{Float64}
 change_point::Vector{Float64}
 ptpower::Vector{Float64}
end

"""
    PityTimerParameters

One draw from `PityTimerSamples`
"""
struct PityTimerParameters
 common_wild_card_base::Float64
 common_wild_card_increment::Float64
 uncommon_wild_card_base::Float64
 uncommon_wild_card_increment::Float64
 rare_card::Float64
 mythic_card::Float64
 rare_wild_card_base::Float64
 rare_wild_card_increment::Float64
 mythic_wild_card_base::Float64
 mythic_wild_card_increment::Float64
 change_point::Float64
 ptpower::Float64
end

function Random.:rand(x::PityTimerSamples)
 i = rand(1:length(x.rare_card))
 PityTimerParameters(
   x.common_wild_card_base[i],
   x.common_wild_card_increment[i],
   x.uncommon_wild_card_base[i],
   x.uncommon_wild_card_increment[i],
   x.rare_card[i],
   x.mythic_card[i],
   x.rare_wild_card_base[i],
   x.rare_wild_card_increment[i],
   x.mythic_wild_card_base[i],
   x.mythic_wild_card_increment[i],
   x.change_point[i],
   x.ptpower[i]
 )
end

let
 local mat, header = DelimitedFiles.readdlm(
     joinpath(@__DIR__, "..", "data", "packmodel_samples.csv"), ',';
     header = true
 )

 @assert size(mat,2) == 12

 global const pity_timer_samps = PityTimerSamples(
     mat[:,1],
     mat[:,2],
     mat[:,3],
     mat[:,4],
     mat[:,5],
     mat[:,6],
     mat[:,7],
     mat[:,8],
     mat[:,9],
     mat[:,10],
     mat[:,11],
     mat[:,12]
 )

end

max_stars_repo_name: asieev/ArenaSim

max_issues_repo_name: asieev/ArenaSim

max_forks_repo_name: asieev/ArenaSim


text:
"""
    LinearElastic(E, ν)

Isotropic linear elasticity.
# Arguments
- `E::Float64`: Young's modulus
- `ν::Float64`: Poisson's ratio
"""
struct LinearElastic <: AbstractMaterial
    # parameters
    E::Float64 # Young's modulus
    ν::Float64 # Poisson's ratio

    # precomputed elastic stiffness tensor
    Eᵉ::SymmetricTensor{4,3,Float64,36}

    LinearElastic(E::Float64, ν::Float64) = new(E, ν, elastic_tangent_3D(E, ν))
end

# Let's always supply a constructor with keyword arguments
LinearElastic(;E::Float64, ν::Float64) = LinearElastic(E, ν)

function elastic_tangent_3D(E::T, ν::T) where T
    λ = E*ν / ((1 + ν) * (1 - 2ν))
    μ = E / (2(1 + ν))
    δ(i,j) = i == j ? 1.0 : 0.0
    f = (i,j,k,l) -> λ*δ(i,j)*δ(k,l) + μ*(δ(i,k)*δ(j,l) + δ(i,l)*δ(j,k))

    C = SymmetricTensor{4, 3, T}(f)
    return C
end

struct LinearElasticState <: AbstractMaterialState end

# define which state belongs to the material
initial_material_state(::LinearElastic) = LinearElasticState()

# for restricted stress states to create buffer
get_stress_type(::LinearElasticState) = SymmetricTensor{2,3,Float64,6}

# constitutive drivers generally operate in 3D
# (we could specialize for lower dimensions if needed for performance)
"""
    material_response(m::LinearElastic, ε::SymmetricTensor{2,3})

Return the stress tensor and the stress tangent for the given strain ε such that

```math
\\boldsymbol{\\sigma} = \\mathbf{E}^\\text{e} : \\boldsymbol{\\varepsilon} .
```
No `MaterialState` is needed for the stress computation, thus if a state is handed over to `material_response`, the same state is returned.
"""
function material_response(m::LinearElastic, ε::SymmetricTensor{2,3}, state::LinearElasticState=LinearElasticState(), Δt=nothing; cache=nothing, options=nothing)
    σ = m.Eᵉ ⊡ ε
    return σ, m.Eᵉ, state
end


text:
Pkg.init()

Pkg.add("BinDeps")

Pkg.add("Compat")
Pkg.build("Compat")

Pkg.add("PyCall")

Pkg.add("HDF5")
#Pkg.add("Silo")
Pkg.clone("https://github.com/MilosSubotic/Silo.jl")
Pkg.build("Silo")
using Silo

Pkg.add("ProfileView")
Pkg.add("FastAnonymous")
Pkg.add("Lumberjack")

Pkg.add("GZip")
Pkg.add("ZipFile")

Pkg.add("PyPlot")
Pkg.add("Winston")
Pkg.add("Gadfly")
Pkg.add("Interact")
Pkg.clone("https://github.com/jverzani/GtkInteract.jl")
Pkg.add("UnicodePlots")

Pkg.add("SymPy")

#Pkg.add("DSP")
Pkg.clone("https://github.com/JuliaDSP/DSP.jl")
Pkg.add("Wavelets")
Pkg.clone("https://github.com/JayKickliter/Radio.jl")

Pkg.clone("https://github.com/jlep/Multicombinations.jl")
Pkg.add("Combinatorics")

Pkg.add("Iterators")
Pkg.add("Cartesian")

Pkg.add("Polynomials")
Pkg.add("PowerSeries")

Pkg.add("Roots")
Pkg.clone("https://github.com/JuliaControl/Control.jl")

Pkg.add("Distributions")
Pkg.add("Regression")
Pkg.add("Optim")
Pkg.add("NLopt")
Pkg.add("JuMP")

Pkg.add("Interpolations")

Pkg.add("LRUCache")


text:
# Model hyperparameters

"""
RAPIDS Docs for Holt Winters Exponential Smoothing:
    https://docs.rapids.ai/api/cuml/stable/api.html#holtwinters

Example:
```
using RAPIDS
using MLJ

X = [1, 2, 3, 4, 5, 6,
    7, 8, 9, 10, 11, 12,
    2, 3, 4, 5, 6, 7,
    8, 9, 10, 11, 12, 13,
    3, 4, 5, 6, 7, 8, 9,
    10, 11, 12, 13, 14]

model = ExponentialSmoothing()
mach = machine(model, X)
fit!(mach)
forecast(mach, 4)
```
"""
MLJModelInterface.@mlj_model mutable struct ExponentialSmoothing <: MMI.Supervised
    handle = nothing
    endog = nothing
    seasonal::String = "additive"::(_ in ("add", "additive", "mul", "multiplicative"))
    seasonal_periods::Int = 2::(_ >= 0)
    start_periods::Int = 2::(_ >= 0)
    ts_num::Int = 1::(_ >= 0)
    eps::Float64 = 2.24e-3::(_ > 0)
    verbose::Bool = false
end

# Multiple dispatch for initializing models
model_init(X, mlj_model::ExponentialSmoothing) = cuml.ExponentialSmoothing(X; mlj_to_kwargs(mlj_model)...)

const CUML_TIME_SERIES = Union{ExponentialSmoothing, }

# add metadata
MMI.load_path(::Type{<:ExponentialSmoothing}) = "$PKG.ExponentialSmoothing"
MMI.input_scitype(::Type{<:CUML_TIME_SERIES}) = Union{AbstractMatrix, Table(Continuous)}
MMI.docstring(::Type{<:ExponentialSmoothing}) = "cuML's ExponentialSmoothing: https://docs.rapids.ai/api/cuml/stable/api.html#holtwinters"

# fit methods
function MMI.fit(mlj_model::CUML_TIME_SERIES, verbosity, X, w=nothing)
    # fit the model
    model = model_init(prepare_input(X), mlj_model)
    model.fit()
    fitresult = model

    # save result
    cache = nothing
    report = ()
    return (fitresult, cache, report)
end
# predict methods
# TODO: Figure out how to handle forecast with MMI
forecast(mach, steps) = pyconvert(Vector{Float32}, mach.fitresult.forecast(steps).to_numpy())

# Classification metadata
MMI.metadata_pkg.((ExponentialSmoothing),
    name = "cuML Time Series Methods",
    uuid = "2764e59e-7dd7-4b2d-a28d-ce06411bac13", # see your Project.toml
    url  = "https://github.com/tylerjthomas9/RAPIDS.jl",  # URL to your package repo
    julia = false,          # is it written entirely in Julia?
    license = "MIT",        # your package license
    is_wrapper = true,      # does it wrap around some other package?
)


text:
using MRIReco
using Test
using LinearAlgebra
using HTTP
using Random
using StatsBase

using Scratch
using LazyArtifacts

const datadir = joinpath(artifact"MRIRecoTestData")
@info "The test data is located at $datadir."

const tmpdir  = @get_scratch!("tmp")
@info "If you want to check the output of the tests, please head to $tmpdir."

include("testIO.jl")
include("testTrajectories.jl")
include("testSimulation.jl")
include("testOperators.jl")
include("testReconstruction.jl")
include("testBrukerFile.jl")
include("testISMRMRD.jl")
include("testCoilsens.jl")


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

JLLWrappers.@generate_main_file_header("wget")
JLLWrappers.@generate_main_file("wget", UUID("25883557-5102-5516-a11b-f84f27e871d7"))
end  # module wget_jll


text:
function dataframe(resp::Requests.Response)

    df = CSV.read(IOBuffer(Requests.text(resp)))

    # force oldest date to first row

    if issorted(df[:Date])
        return df
    else
        return sort!(df)
    end
end


text:
using Test
using Distributions
using StatsFuns
# using Random

using LRMoE
using JLD2

# Random.seed!(1234)

@load "test/X_all.JLD2" X_all
@load "test/Y_all.JLD2" Y_all

Y_all = convert(Matrix{Int}, Y_all)

# tmp = LRMoE.cmm_init(Y_all, X_all, 3, ["discrete" "discrete"]; exact_Y = true, n_random = 0)

# tmp.params_init[1][1][8]

# 7: ZI-NB
# 8: ZI-GC

model_list = []

for n_comp in 2:20
    init = LRMoE.cmm_init(Y_all, X_all, n_comp, ["discrete" "discrete"]; exact_Y = true, n_random = 0)
    # push!(model_list, (init.α_init, init.ll_best))
    tmp = vcat([hcat([init.params_init[d][i][8] for i in 1:n_comp]...) for d in 1:2]...)
    push!(model_list, (init.α_init, tmp))
end

length(model_list)

model_list[1][2]
model_list[2][2]
model_list[3][2]

@save "model_list_GC.JLD2" model_list

max_stars_repo_name: sparktseung/JuliaLRMoE

max_issues_repo_name: sparktseung/JuliaLRMoE

max_forks_repo_name: sparktseung/JuliaLRMoE


text:
using DiffEqSensitivity, Test, Random
Random.seed!(123)

function ishi(X)
    A= 7
    B= 0.1
    [sin(X[1]) + A*sin(X[2])^2+ B*sin(X[3])^4 *sin(X[1])]
end
p_range = [[0.0,1.0] for i in 1:4]
N = 60000
sobol = gsa(ishi,p_range,Sobol(N=N,order=[0,1,2]))

for (first_order, r_sens) in zip(sobol.S1,[0.0242498108,0.9801111059,0.0000346274,0.0000000000])
     @test first_order[1] ≈ r_sens atol=2e-2
end
for (second_order, r_sens) in zip(sobol.S2,[ 0.999923,0.021941,0.021867,0.974586,0.974525,-0.004465])
    @test second_order[1] ≈ r_sens atol=2e-2
end
for (total_order, r_sens) in zip(sobol.ST,[2.682048e-02,9.693216e-01,6.910431e-05,0.000000e+00])
    @test total_order[1] ≈ r_sens atol=2e-2
end

######################################################################
#Model runs: 600000
#
#First order indices:
#       original          bias   std. error     min. c.i.    max. c.i.
#X1 0.0242498108 -6.038174e-05 1.130665e-03  0.0219442397 0.0266505624
#X2 0.9801111059  1.318961e-03 9.755089e-03  0.9619890024 1.0038331890
#X3 0.0000346274  1.564015e-06 7.342856e-05 -0.0001075782 0.0001843829
#X4 0.0000000000  0.000000e+00 0.000000e+00  0.0000000000 0.0000000000
#
#Total indices:
#       original          bias   std. error     min. c.i.    max. c.i.
#X1 2.682048e-02  1.087146e-04 1.177283e-03  2.462456e-02 0.0294033066
#X2 9.693216e-01 -7.482170e-04 7.634399e-03  9.550214e-01 0.9875516033
#X3 6.910431e-05  8.949816e-06 7.376054e-05 -9.023791e-05 0.0002100689
#X4 0.000000e+00  0.000000e+00 0.000000e+00  0.000000e+00 0.0000000000

#Second order indices
#     original std. error min. c.i. max. c.i.
#X12  0.999923   0.000001  0.999922  0.999924
#X13  0.021941   0.004076  0.013953  0.029929
#X14  0.021867   0.004076  0.013878  0.029856
#X23  0.974586   0.000154  0.974284  0.974888
#X24  0.974525   0.000154  0.974223  0.974827
#X34 -0.004465   0.004083 -0.012468  0.003538
######################################################################