Ported Tdbtt #14

prakharcode · 2018-03-18T16:00:17Z

Ported eraTdbtt to pure Julia.
Ported dependency eraDtdb to pure Julia.

Used approximation due to precision error.

julia> ref = TDBEpoch(2013, 3, 18, 12)
2013-03-18T12:00:00.000 TDB

julia> ifu = AstronomicalTime.Epochs.dtdb(ref.jd1, ref.jd2, 0.0, 0.0, 0.0, 0.0)
0.0015985185850907002

julia> ou = ERFA.dtdb(ref.jd1, ref.jd2, 0.0, 0.0, 0.0, 0.0)
0.0015774019690913835

giordano · 2018-03-18T16:05:44Z

src/conversions.jl

+    # Combine time argument (millennia) with deg/arcsec factor.
+     w = t / 3600.0
+    # Sun Mean Longitude.
+     elsun = mod(280.46645683 + 1296027711.03429 * w, 360.0) * DEG_2R


Use deg2rad to convert degrees to radians and remove the DEG_2R constant. Here and below

giordano · 2018-03-18T16:09:16Z

test/runtests.jl

-        @test ref == TDBEpoch(TTEpoch("2013-03-18T11:59:59.998"))
-        @test ref == TDBEpoch(TCBEpoch("2013-03-18T12:00:17.718"))
-        @test ref == TDBEpoch(TCGEpoch("2013-03-18T12:00:00.795"))
+        @test UT1Epoch(ref) ≈ UT1Epoch("2013-03-18T11:58:52.994")


Why approximated equality?

See below ;-)

coveralls · 2018-03-18T16:09:33Z

Coverage increased (+0.5%) to 97.288% when pulling 8d7a9bf on prakharcode:tbdtt into 11eb598 on JuliaAstro:master.

helgee · 2018-03-18T16:09:58Z

The discrepancy seems too high to me.

helgee · 2018-03-18T16:12:58Z

src/conversions.jl

+    end
+    # T**1
+     w1 = 0
+     for j in 679:-1:1


The lower bound of this loop and of the following ones are wrong.

The loops all work on non-overlapping segements of the coefficients (or would if they were correct 😉). Thus, the first step to simplify this would be to split coefficients into five separate arrays.

giordano

As a side note, the fairhd array is quite big and is just a container of numbers. @helgee would be a StaticArray an option? I don't know how it'd perform in this case.

codecov-io · 2018-03-18T16:14:47Z

Codecov Report

Merging #14 into master will increase coverage by 0.46%.
The diff coverage is 100%.

@@            Coverage Diff             @@
##           master      #14      +/-   ##
==========================================
+ Coverage   96.82%   97.28%   +0.46%     
==========================================
  Files           6        6              
  Lines         252      295      +43     
==========================================
+ Hits          244      287      +43     
  Misses          8        8

Impacted Files	Coverage Δ
src/AstronomicalTime.jl	`100% <ø> (ø)`	⬆️
src/conversions.jl	`99.43% <100%> (+0.17%)`	⬆️

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 11eb598...8d7a9bf. Read the comment docs.

helgee · 2018-03-18T16:19:18Z

@giordano If I am not mistaken, a mutable array is better in this case since it will be heap-allocated and there are no copies when it is passed around.

helgee · 2018-03-18T16:22:11Z

src/conversions.jl

+     w4 = 0
+     for j in 787:-1:1
+        w4 += fairhd[j][1] * sin(fairhd[j][2] * t + fairhd[j][3])
+    end


I wonder if there might be a more straight-forward way to compute this. These countdown-loops are bit confusing too me 😜

helgee · 2018-03-18T16:23:16Z

src/conversions.jl

+        w4 += fairhd[j][1] * sin(fairhd[j][2] * t + fairhd[j][3])
+    end
+    # Multiply by powers of T and combine.
+     wf = t * (t * (t * (t * w4 + w3) + w2) + w1) + w0


Same here. This is awful! Not your fault, of course @prakharcode .

@evalpoly t w0 w1 w2 w3 w4

helgee · 2018-03-18T16:24:37Z

src/conversions.jl

+function dtdb(jd1, jd2, ut, elong, u, v)
+
+    t = ((jd1 - DJ00) + jd2) / DJM
+    # Convert UT to local solar time in radians.


Indentation is off.

helgee · 2018-03-18T16:25:20Z

src/constants.jl

+
+const DEG_2R = 1.745329251994329576923691e-2
+
+const fairhd = [


Maybe move this to a separate file?

helgee · 2018-03-18T16:26:41Z

The array should probably not be accessed in row-major order.

helgee · 2018-03-18T16:37:22Z

test/runtests.jl

@@ -170,6 +172,9 @@ AstronomicalTime.update()
        @test AstronomicalTime.Epochs.tttcg(julian2(tt), julian1(tt)) == ERFA.tttcg(julian2(tt), julian1(tt))
        @test AstronomicalTime.Epochs.taiutc(julian1(tai), julian2(tai)) == ERFA.taiutc(julian1(tai), julian2(tai))
        @test AstronomicalTime.Epochs.taiutc(julian2(tai), julian1(tai)) == ERFA.taiutc(julian2(tai), julian1(tai))
+        @test AstronomicalTime.Epochs.dtdb(julian1(tdb), julian2(tdb), 0.0, 0.0, 0.0, 0.0) ≈ ERFA.dtdb(julian1(tdb), julian2(tdb), 0.0, 0.0, 0.0, 0.0)


helgee · 2018-03-18T16:38:29Z

src/constants.jl

@@ -4,7 +4,7 @@ export MJD, J2000, J1950,
    HOURS_PER_DAY, HOURS_PER_YEAR, HOURS_PER_CENTURY,
    DAYS_PER_YEAR, DAYS_PER_CENTURY,
    YEARS_PER_CENTURY,
-    OFFSET_TT_TAI, MOD_JD_77, ELG
+    OFFSET_TT_TAI, MOD_JD_77, ELG, fairhd, DJM, DJ00


The coefficients should not be exported.

helgee · 2018-03-18T16:47:20Z

I really want to get rid of these funky backwards counting loops. Does anyone have a theory why the original code is written this way?
The coefficients could be split into five separate arrays in any case. The loops do not overlap.

giordano · 2018-03-18T17:01:22Z

src/constants.jl

@@ -34,3 +34,8 @@ const OFFSET_TT_TAI = 32.184

 const MOD_JD_77 = 43144.0
 const ELG = 6.969290134e-10
+const DJM = 365250.0


Should this be defined in terms of DAYS_PER_YEAR? I guess these are days per millennium, maybe the constant can be given a more meaningful name.

benelsen · 2018-03-18T17:13:31Z

I assume the backward counting loops are for numerical stability in the summation, the coefficients are (mostly) decreasing in magnitude. The alternative is rewriting it to use a kbn summation or something alike.

helgee · 2018-03-18T17:21:38Z

Since Julia uses pairwise summation this might not be needed (see JuliaLang/julia#4039).

giordano · 2018-03-18T17:37:06Z

But sum is not (yet?) used here

helgee · 2018-03-18T17:46:55Z

w0 = sum(fairhd0[:,1] .* sin(fairhd0[:,2] .* t .+ fairhd0[:,3]))

giordano · 2018-03-18T20:27:07Z

Before merging this, we may want to check its performance. Currently it's slower than the ERFA function:

julia> tdb = TDBEpoch(2000, 1, 1, 12, 0, 0.0)
2000-01-01T12:00:00.000 TDB

julia> a = julian1(tdb)
2.4515445e6

julia> b = julian2(tdb)
0.5

julia> @benchmark ERFA.dtdb($a, $b, 0.0, 0.0, 0.0, 0.0)
BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     16.025 μs (0.00% GC)
  median time:      16.182 μs (0.00% GC)
  mean time:        16.353 μs (0.00% GC)
  maximum time:     46.679 μs (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     1

julia> @benchmark AstronomicalTime.Epochs.dtdb($a, $b, 0.0, 0.0, 0.0, 0.0)
BenchmarkTools.Trial: 
  memory estimate:  37.05 KiB
  allocs estimate:  2371
  --------------
  minimum time:     36.573 μs (0.00% GC)
  median time:      37.647 μs (0.00% GC)
  mean time:        40.081 μs (3.31% GC)
  maximum time:     974.220 μs (89.31% GC)
  --------------
  samples:          10000
  evals/sample:     1

benelsen · 2018-03-19T07:14:19Z

src/conversions.jl

+    # =====================
+
+    # T**0
+     w0 = 0


w0 … w4 are first assigned an Int64 (w0 = 0) and then reassigned a Float64 (w0 += ) two lines later.
Using w0 = 0.0 (same for w1 … w4) one can prevent a costly type-instability.

Good catch! 10% of improvement is very nice

benelsen · 2018-03-19T07:20:37Z

This is due to a type-instability in the variables used for summation, once those are initialized as floats the allocations are gone and the performance is better.

julia> b1 = @benchmark AstronomicalTime.Epochs.dtdb($a, $b, 0.123, 0.234, 0.456, 0.789)
BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     10.916 μs (0.00% GC)
  median time:      11.054 μs (0.00% GC)
  mean time:        11.063 μs (0.00% GC)
  maximum time:     27.996 μs (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     1

julia> b2 = @benchmark ERFA.dtdb($a, $b, 0.123, 0.234, 0.456, 0.789)
BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     12.242 μs (0.00% GC)
  median time:      12.295 μs (0.00% GC)
  mean time:        12.302 μs (0.00% GC)
  maximum time:     21.806 μs (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     1

julia> judge(median(b1), median(b2))
BenchmarkTools.TrialJudgement: 
  time:   -10.09% => improvement (5.00% tolerance)
  memory: +0.00% => invariant (1.00% tolerance)

giordano · 2018-03-19T10:30:56Z

src/conversions.jl

+        w4 += fairhd[j][1] * sin(fairhd[j][2] * t + fairhd[j][3])
+    end
+    # Multiply by powers of T and combine.
+     wf = t * (t * (t * (t * w4 + w3) + w2) + w1) + w0


I already wrote this before but got lost: replace this line with

wf = @evalpoly t w0 w1 w2 w3 w4

helgee · 2018-03-31T06:31:00Z

I did a few experiments. The fastest version I was able to come up with is very close to ERFA in terms of performance (consistently within 10%). It uses StaticArrays.jl after all. Originally I did not realize that the algorithm does not work on a matrix of coefficients but actually iterates over an array of small vectors.
With the coefficients partioned in five tuples of SVectors the loops can be written like this:

    w0 = 0.0
    for j in reverse(eachindex(fairhd0))
        w0 += fairhd0[j][1] * sin(fairhd0[j][2] * t + fairhd0[j][3])
    end

Which is much nicer IMHO.

giordano · 2018-03-31T06:34:47Z

The vector themselves can be written in the reverse order (or reversed at compile time), rather than using reverse at runtime.

helgee · 2018-03-31T15:48:47Z

Does not make a significant difference.

helgee · 2018-03-31T16:33:02Z

Final version, with tuples of tuples instead of StaticArrays.jl:

function dtdb_tuples2(jd1, jd2, ut, elong, u, v)
    t = ((jd1 - DJ00) + jd2) / DJM
    # Convert UT to local solar time in radians.
    tsol = mod(ut, 1.0) * 2π  + elong

    # FUNDAMENTAL ARGUMENTS:  Simon et al. 1994.
    # Combine time argument (millennia) with deg/arcsec factor.
    w = t / 3600.0
    # Sun Mean Longitude.
    elsun = deg2rad(mod(280.46645683 + 1296027711.03429 * w, 360.0))
    # Sun Mean Anomaly.
    emsun = deg2rad(mod(357.52910918 + 1295965810.481 * w, 360.0))
    # Mean Elongation of Moon from Sun.
    d = deg2rad(mod(297.85019547 + 16029616012.090 * w, 360.0))
    # Mean Longitude of Jupiter.
    elj = deg2rad(mod(34.35151874 + 109306899.89453 * w, 360.0))
    # Mean Longitude of Saturn.
    els = deg2rad(mod(50.07744430 + 44046398.47038 * w, 360.0))
    # TOPOCENTRIC TERMS:  Moyer 1981 and Murray 1983.
    wt =   +  0.00029e-10 * u * sin(tsol + elsun - els)
    +  0.00100e-10 * u * sin(tsol - 2.0 * emsun)
    +  0.00133e-10 * u * sin(tsol - d)
    +  0.00133e-10 * u * sin(tsol + elsun - elj)
    -  0.00229e-10 * u * sin(tsol + 2.0 * elsun + emsun)
    -  0.02200e-10 * v * cos(elsun + emsun)
    +  0.05312e-10 * u * sin(tsol - emsun)
    -  0.13677e-10 * u * sin(tsol + 2.0 * elsun)
    -  1.31840e-10 * v * cos(elsun)
    +  3.17679e-10 * u * sin(tsol)
    # =====================
    # Fairhead et al. model
    # =====================

    # T**0
    w0 = 0.0
    for j in eachindex(fairhd0_4)
        @muladd w0 += fairhd0_4[j][1] * sin(fairhd0_4[j][2] * t + fairhd0_4[j][3])
    end
    # T**1
    w1 = 0.0
    for j in eachindex(fairhd1_4)
        @muladd w1 += fairhd1_4[j][1] * sin(fairhd1_4[j][2] * t + fairhd1_4[j][3])
    end
    # T**2
    w2 = 0.0
    for j in eachindex(fairhd2_4)
        @muladd w2 += fairhd2_4[j][1] * sin(fairhd2_4[j][2] * t + fairhd2_4[j][3])
    end
    # T**3
    w3 = 0.0
    for j in eachindex(fairhd3_4)
        @muladd w3 += fairhd3_4[j][1] * sin(fairhd3_4[j][2] * t + fairhd3_4[j][3])
    end
    # T**4
    w4 = 0.0
    for j in eachindex(fairhd4_4)
        @muladd w4 += fairhd4_4[j][1] * sin(fairhd4_4[j][2] * t + fairhd4_4[j][3])
    end
    # Multiply by powers of T and combine.
    #= wf = t * (t * (t * (t * w4 + w3) + w2) + w1) + w0 =#
    wf = @evalpoly t w0 w1 w2 w3 w4
    # Adjustments to use JPL planetary masses instead of IAU.
    wj = 0.00065e-6 * sin(6069.776754 * t + 4.021194) +
        0.00033e-6 * sin( 213.299095 * t + 5.543132) +
        (-0.00196e-6 * sin(6208.294251 * t + 5.696701)) +
        (-0.00173e-6 * sin(  74.781599 * t + 2.435900)) +
        0.03638e-6 * t * t
    # ============
    # Final result
    # ============
    # TDB-TT in seconds.
    w = wt + wf + wj
end

0.6: 10.504 μs
0.7: 7.166 μs 🎉

giordano · 2018-03-31T17:25:50Z

Good! I'm surprised long tuples (more than 16 elements) don't have some overhead. What's performance of the function from ERFA on your machine?

helgee · 2018-03-31T18:29:00Z

~9 μs

…

Am 31.03.2018 um 19:25 schrieb Mosè Giordano ***@***.***>: Good! I'm surprised long tuples (more than 16 elements) don't have some overhead. What's performance of the function from ERFA on your machine? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or mute the thread.

benelsen · 2018-03-31T20:29:02Z

AFAIR StaticArrays are basically NTuples with added sugar.
You can maybe try sprinkling in some @inbounds for the last few ‰ as it's not really possible to get out of bounds there.

prakharcode · 2018-04-03T14:04:46Z

@helgee is this error tolerable: -9.930719894379461e-5 == -9.930719894379447e-5, when comparing AstronomicalTime.Epochs.dtdb(julian1(tdb), julian2(tdb), 0.0, 0.0, 0.0, 0.0) == ERFA.dtdb(julian1(tdb), julian2(tdb), 0.0, 0.0, 0.0, 0.0)
?

prakharcode · 2018-04-03T14:29:44Z

@giordano the edits suggested by you, @helgee and @benelsen helped a lot in making the function efficient. Now the function benchmarks as follows:

julia> b1 = @benchmark ERFA.dtdb($a, $b, 0.0, 0.0, 0.0, 0.0)
BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     16.120 μs (0.00% GC)
  median time:      16.244 μs (0.00% GC)
  mean time:        16.277 μs (0.00% GC)
  maximum time:     36.475 μs (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     1

julia> b2 = @benchmark AstronomicalTime.Epochs.dtdb($a, $b, 0.0, 0.0, 0.0, 0.0)
BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     11.470 μs (0.00% GC)
  median time:      11.766 μs (0.00% GC)
  mean time:        11.866 μs (0.00% GC)
  maximum time:     23.690 μs (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     1

julia> judge(median(b2), median(b1))
BenchmarkTools.TrialJudgement: 
  time:   -27.57% => improvement (5.00% tolerance)
  memory: +0.00% => invariant (1.00% tolerance)

helgee · 2018-04-03T14:51:43Z

There is no error on my machine. Did you change the direction of the loops?

prakharcode · 2018-04-03T15:25:02Z

Yes @helgee, I resolved it!

helgee · 2018-04-03T15:47:31Z

src/constants.jl

@@ -34,3 +34,8 @@ const OFFSET_TT_TAI = 32.184

 const MOD_JD_77 = 43144.0
 const ELG = 6.969290134e-10
+const DAYS_PER_MILLENNIUM = 365250.0
+const DJ00 = 2451545.0


Same as J2000

helgee · 2018-04-03T15:48:21Z

You need to add the new dependency to REQUIRE.

helgee · 2018-04-03T17:15:55Z

Nice work everybody 🎉

) * Use legal generated functions for conversions * Enable precompilation * Minor refactor of Period type and whitespace cleanup * Revert "Minor refactor of Period type and whitespace cleanup" This reverts commit b375073. * Make parameter for transform macro explicit * Replacing the calls to ERFA by pure Julia functions for eraTaiutc (#13) * Replacing the calls to ERFA by pure Julia functions for eraUtctai (#15) * Ported Tdbtt (#14) * Use legal generated functions for conversions * Allow additional arguments in transform macro * Fix transform and various cleanups * Cleanup

prakharcode added 5 commits March 18, 2018 19:50

added dtdb

983380d

Added tdbtt and dtdb

800527e

resolved

74c9f58

tdbtt done

045e04a

tdbtt docs

8624c61

giordano reviewed Mar 18, 2018

View reviewed changes

helgee reviewed Mar 18, 2018

View reviewed changes

giordano reviewed Mar 18, 2018

View reviewed changes

helgee reviewed Mar 18, 2018

View reviewed changes

src/constants.jl Outdated

const DEG_2R = 1.745329251994329576923691e-2

const fairhd = [

Copy link

Member

helgee Mar 18, 2018

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe move this to a separate file?

fix

7a9d837

helgee reviewed Mar 18, 2018

View reviewed changes

fix_1

addb1c9

helgee mentioned this pull request Mar 18, 2018

Excise ERFA #1

Closed

20 tasks

giordano reviewed Mar 18, 2018

View reviewed changes

benelsen reviewed Mar 19, 2018

View reviewed changes

giordano reviewed Mar 19, 2018

View reviewed changes

helgee and others added 2 commits April 3, 2018 08:01

Merge branch 'master' into tbdtt

8521136

merged master

6b4bb36

prakharcode added 3 commits April 3, 2018 20:29

updated

35f4755

res

5b57e78

reversed

9bbb5f2

required changed

1685828

helgee reviewed Apr 3, 2018

View reviewed changes

require

8d7a9bf

helgee merged commit fe85130 into JuliaAstro:master Apr 3, 2018

prakharcode deleted the tbdtt branch April 3, 2018 17:33

helgee pushed a commit that referenced this pull request Apr 4, 2018

Ported Tdbtt (#14)

a116172


		const DEG_2R = 1.745329251994329576923691e-2

		const fairhd = [

Ported Tdbtt #14

Ported Tdbtt #14

Conversation

prakharcode commented Mar 18, 2018

giordano Mar 18, 2018 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

coveralls commented Mar 18, 2018 • edited Loading

helgee commented Mar 18, 2018

Choose a reason for hiding this comment

Choose a reason for hiding this comment

giordano left a comment

Choose a reason for hiding this comment

codecov-io commented Mar 18, 2018 • edited Loading

Codecov Report

helgee commented Mar 18, 2018

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

helgee commented Mar 18, 2018

Choose a reason for hiding this comment

Choose a reason for hiding this comment

helgee commented Mar 18, 2018

giordano Mar 18, 2018 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

benelsen commented Mar 18, 2018

helgee commented Mar 18, 2018

giordano commented Mar 18, 2018

helgee commented Mar 18, 2018

giordano commented Mar 18, 2018

Choose a reason for hiding this comment

Choose a reason for hiding this comment

benelsen commented Mar 19, 2018

Choose a reason for hiding this comment

helgee commented Mar 31, 2018 • edited Loading

giordano commented Mar 31, 2018

helgee commented Mar 31, 2018

helgee commented Mar 31, 2018

giordano commented Mar 31, 2018

helgee commented Mar 31, 2018 via email

benelsen commented Mar 31, 2018

prakharcode commented Apr 3, 2018

prakharcode commented Apr 3, 2018

helgee commented Apr 3, 2018

prakharcode commented Apr 3, 2018

Choose a reason for hiding this comment

helgee commented Apr 3, 2018

helgee commented Apr 3, 2018

giordano Mar 18, 2018 •

edited

Loading

coveralls commented Mar 18, 2018 •

edited

Loading

codecov-io commented Mar 18, 2018 •

edited

Loading

giordano Mar 18, 2018 •

edited

Loading

helgee commented Mar 31, 2018 •

edited

Loading