calculation speed worst than pyepehm?

[nachoplus]

Hello,

Thanks you very much for your work!

I am translating pyephem code to skyfield and I have noticed pyephem calculation run at least at 30x faster than skyfield one.

My program wants to calculate all the artificial satellite positions for a given lat/lon/date using the NORAD TLEs

The time consuming math is in:
sat = earth.satellite(tle)
apparent = self.here(self.date).observe(sat)

real 0m31.959s
user 0m31.550s
sys 0m0.208s

compare with the pyephem:
tle.compute(self.here)

real 0m0.766s
user 0m0.532s
sys 0m0.064s

Are you aware of that? ..Or perhaps I am doing in the wrong way ...
For my app speed is critical.. Is there any way to speed up?

Nacho

[brandon-rhodes]

The base reason for the problem is that Skyfield is written in pure-Python instead of in C. There are two possible solutions, and I will take a look at them both for you.

First, I want to try cleaning up sgp4 (the package that is doing the actual slow satellite work) enough that the Numba compiler could compile it to C code with its jit decorator.

Second, I should try Skyfield under PyPy to see whether it is able to run or not. It would be nice for users to be able to opt for PyPy and its JIT-powered machine-level efficiency, instead of having to wait for C Python to build and destroy all of the float objects required as intermediate values in the long sgp4 operations!

I will report back here with comments about my progress on both fronts.

[brandon-rhodes]

Having just re-read your issue, I see that you are getting a time of 31 seconds, not milliseconds as I had initially thought when I first read it!

I would like to be able to run this calculation on my own computer — could you provide a little script that builds the TLE and builds here so that I have an example that I can run and see it take that long? Thanks!

[nachoplus]

Hello Brando,

Thanks for your support and interest. I would like to clarify that 31 seconds is for calculate all the TLE in the NORAD data base (more than 15000 objects!)

My code is part of a huge project and has a lot of dependencies and big data files, thus I need to make a simple example from scrach to show you the problem. I hope could be ready at the end of this week. I'll keep you informed.

Thanks again!

[brandon-rhodes]

Oh! Good, I am glad that a single satellite is not causing that much slowdown. If it is many computations that are taking this long, then I can probably build a small benchmark myself to work against — I just didn't know how to make even a long list of positions for a single satellite take quite that long :)

[nachoplus]

...I am updating my system to ubuntu trusty and I have to install all the stuff.. I send you a draff of my code just in case of you can't wait for my example ...

I got NORAD full tle database from:
http://www.idb.com.au/files/TLE_DATA/ALL_TLE.ZIP

Unzip then I read and compute the position of all objects using something like:

def readTLEfile(TLEfile):
    f=open(TLEfile)
    theList=f.read().split('\n')
    f.close()
    N=3
    data = [''.join(theList[n:n+N]) for n in range(0, len(theList), N)]
    return data[:-1]

def compute(lat,lon):
    astPos=[["RA","DEC","ELEVATION"]]
    TLEs=readTLEfile(<TheFile>)
    here = earth.topos(str(lat), str(lon))
    date= datetime.datetime.strptime("2014-08-30 23:00:00", '%Y-%m-%d %H:%M:%S')
    date=date.replace(tzinfo=utc)
    date_=JulianDate(utc=date)      
    for i,tle in enumerate(TLEs):
        sat = earth.satellite(tle)
        apparent = here(date_).observe(sat)
        ra,dec,distance=apparent.radec()
        astPos.append(ra,dec,distance)
    return astPos

[nachoplus]

Hi,

Here is the benchmark script. To run first download http://www.idb.com.au/files/TLE_DATA/ALL_TLE.ZIP and unzip. Also you have to use my pyephm dirty hack: brandon-rhodes/pyephem#28 If not, pyephem calcs get stucked

#!/usr/bin/python
# -*- coding: iso-8859-15 -*-

from skyfield.api import earth,JulianDate,utc
from datetime import datetime

import ephem

import csv
import cProfile

def readTLEfile(TLEfile):
    f=open(TLEfile)
    theList=f.read().split('\n')
    f.close()
    N=3
    data = [''.join(theList[n:n+N]) for n in range(0, len(theList), N)]
    return data[:-1]

class skyFieldTest():


    def compute(self,lat,lon):
        astPos=[["ID","RA","DEC","ELEVATION"]]
        TLEs=readTLEfile("ALL_TLE.TXT")
        here = earth.topos(str(lat), str(lon))
        date= datetime.strptime("2014-08-30 23:00:00", '%Y-%m-%d %H:%M:%S')
        date=date.replace(tzinfo=utc)
        date_=JulianDate(utc=date)      
        for i,tle in enumerate(TLEs):
            sat = earth.satellite(tle)
            apparent = here(date_).observe(sat)
            ra,dec,distance=apparent.radec()
            astPos.append((i,ra,dec,distance))
        return astPos

    def do(self,lat,lon):
        astPos=self.compute(lat,lon)
        print "Total number of TLEs calculated: ",len(astPos)-1
        print "Writing CSV result on skyPOS.txt"
        with open("skyPOS.txt", 'wb') as f:
            writer = csv.writer(f)
            writer.writerows(astPos)


class ephemTest():

    def compute(self,lat,lon):
        astPos=[["ID","RA","DEC","ELEVATION"]]
        TLEs=readTLEfile("ALL_TLE.TXT")
        here = ephem.Observer()
            here.lat, here.lon, here.horizon  = str(lat), str(lon), str(10)
            here.elev = float(900)
            here.temp = 25e0
            here.compute_pressure()
        here.date=ephem.date("2014-08-30 23:00:00")
        for i,tle in enumerate(TLEs):
            lines=tle.split('\r')[:-1]
        try:    
           sat=ephem.readtle(lines[0],lines[1],lines[2])
        except:
           print "ERROR computing or reading TLE:\n",lines
           continue
                sat.compute(here)
        ra  =sat.ra
        dec =sat.dec
        distance=sat.elevation
            astPos.append((i,ra,dec,distance))
        return astPos

    def do(self,lat,lon):
        astPos=self.compute(lat,lon)
        print "Total number of TLEs calculated: ",len(astPos)-1
        print "Writing CSV result on: ephemPOS.txt"
        with open("ephemPOS.txt", 'wb') as f:
            writer = csv.writer(f)
            writer.writerows(astPos)




if __name__ == '__main__':
    print " "
    print "==== Skyfield vs pyephem benchmarck ===="
    print " "
    print "Starting.",datetime.utcnow()
    print "Creating classes.."
    sky=skyFieldTest()
    eph=ephemTest()
    print " "
    t0=datetime.utcnow()
    print "Start calculation using skyfield:",t0
    sky.do("40.340388 N","4.346194 W")
    t1=datetime.utcnow()
    print "Calculation ended",t1
    lsky=t1-t0
    print "Lasting (s)",lsky
    print " "
    t0=datetime.utcnow()
    print "start calculation using pyephem:",t0
    eph.do(lon="-4:20:46.3",lat="40:20:25.4")
    t1=datetime.utcnow()
    print "Calculation ended",t1
    leph=t1-t0
    print " "
    print "\t\tLasting (hh:mm:ss)"
    print "pyephem \t",leph
    print "slyField\t",lsky
    print " "

And this is the result in my machine (ubuntu-trusty i7 16G RAM):

==== Skyfield vs pyephem benchmarck ====

Starting. 2014-11-29 23:02:38.043484
Creating classes..

Start calculation using skyfield: 2014-11-29 23:02:38.043503
Total number of TLEs calculated:  15050
Writing CSV result on skyPOS.txt
Calculation ended 2014-11-29 23:03:08.092066
Lasting (s) 0:00:30.048563

start calculation using pyephem: 2014-11-29 23:03:08.092110
Total number of TLEs calculated:  15050
Writing CSV result on: ephemPOS.txt
Calculation ended 2014-11-29 23:03:08.349396

        Lasting (hh:mm:ss)
pyephem     0:00:00.257286
slyField    0:00:30.048563

I hope this would be usefull...

Cheers!
Nacho

[brandon-rhodes]

Even while you are waiting for me to perform some optimizations on the sgp4 module, there is one quick adjustment you can make to your code for a slight speedup: instead of running the expensive here(date_) function every single time through your loop to figure out where the Earth is pointing, do something like this before you start the loop:

here_now = here(date_)

Then you can simply say:

here_now.observe(sat)

If you make that change, are you able to measure any speedup on your machine?

[nachoplus]

Hi @brandon-rhodes ,

Ok. I follow your sugestion:

....
#Init here_now only one time outside the loop
        here_now = here(date_)    
        for i,tle in enumerate(TLEs):
            sat = earth.satellite(tle)
            apparent = here_now.observe(sat)
            ra,dec,distance=apparent.radec()
            astPos.append((i,ra,dec,distance))
.....

And this is the result:

==== Skyfield vs pyephem benchmarck (v2)====

Starting. 2014-12-10 07:47:50.232971
Creating classes..

Start calculation using skyfield: 2014-12-10 07:47:50.232989
Total number of TLEs calculated:  15050
Writing CSV result on skyPOS.txt
Calculation ended 2014-12-10 07:48:03.711814
Lasting (s) 0:00:13.478825

start calculation using pyephem: 2014-12-10 07:48:03.711860
Total number of TLEs calculated:  15050
Writing CSV result on: ephemPOS.txt
Calculation ended 2014-12-10 07:48:03.998857
Lasting (s) 0:00:00.286997

        Lasting (hh:mm:ss)
pyephem     0:00:00.286997
slyField    0:00:13.478825

Big improvement! but not enough for my needs.

It is my intention to go deeper to see what's going on ussing python line_profiler or similar tool.

[nachoplus]

The error could be due to a bad ALL_TLE.TXT file. Some time they public
with a some lines wrong. Be sure that is rigth o try different file. You
can get some from http://www.tle.info/data/ Any of this file sould work

Sorry for the bad indentation but is not mine fault, is a GitHub defect. I
just copy/paste my tested code..

Nacho

2014-12-11 21:23 GMT+01:00 wenzul notifications@github.com:

@nachoplus https://github.com/nachoplus Your code isn't working for me.
You have wrong indentations and after I fix that I get

twoline2rv() takes at least 3 arguments (2 given)

Is this really the code which you have used?

—
Reply to this email directly or view it on GitHub
#30 (comment)
.

[wenzul]

These are the top x time consuming functions.

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
    15071    9.709    0.001   20.675    0.001 nutationlib.py:208(iau2000a)
    30138    4.155    0.000    6.498    0.000 io.py:228(sscanf)
     6038    3.606    0.001    5.664    0.001 propagation.py:933(_dspace)
   180852    3.460    0.000    3.460    0.000 {method 'reshape' of 'numpy.ndarray' objects}
   120564    2.938    0.000    2.938    0.000 {method 'dot' of 'numpy.ndarray' objects}
    30138    2.401    0.000    9.043    0.000 propagation.py:1647(sgp4)
    60284    2.151    0.000    7.353    0.000 numeric.py:1106(tensordot)
  1804588    1.324    0.000    1.324    0.000 {math.sin}
   301465    1.228    0.000    1.228    0.000 {numpy.core.multiarray.array}
  1816612    1.197    0.000    1.197    0.000 {math.cos}
    15071    1.162    0.000    1.632    0.000 nutationlib.py:87(equation_of_the_equinoxes_complimentary_terms)
    15069    1.125    0.000   11.665    0.001 io.py:90(twoline2rv)
    15069    1.108    0.000    3.458    0.000 propagation.py:1276(sgp4init)
    15071    0.762    0.000    0.778    0.000 nutationlib.py:281(fundamental_arguments)

I think the numpy arrays consume a lot of time because in this scenario - creating sat objects and compute just one date - is a lot of overhead without using the advantages of the array support.

[astrojuanlu]

Hello! Are you still thinking about cleaning sgp4 up to numba-JIT it? I was considering it some days ago and stumbled upon your comment by chance.

On the other hand, Travis CI supports PyPy, so perhaps you could use it to see if sgp4 and skyfield work there.

[brandon-rhodes]

I am still very interested in seeing the SGP4 algorithm cleaned up! By pivoting special placeholder values from None to nan so that they are floating point either way, and some other adjustments, I do think it could be made Numba-friendly — but I am not likely to have time to tackle the issue for several months because of other commitments.

I have never been able to get Skyfield working under PyPy because PyPy did not, the last time I tried, support enough NumPy calls (einsum was missing, I think?) to run.

[ckuethe]

I would love to see sgp4 run faster also. It's taking just under a minute (run time: 57.44336s) to compute a single position on my i5 machine using a very simple test script:

#!/usr/bin/env python
import ephem
from datetime import datetime

from time import clock
from sgp4.earth_gravity import wgs72
from sgp4.io import twoline2rv

goes15 = ['GOES 15',
    '1 36411U 10008A   16146.69708964  .00000087  00000-0  00000+0 0  9995',
    '2 36411   0.1384 104.2770 0001334 286.4218 329.3213  1.00258510 22847']

sat_obj = twoline2rv(goes15[1], goes15[2], wgs72)
calc_at_time = ephem.Date(datetime.now())

c_start = clock()
position, velocity = sat_obj.propagate(calc_at_time)
c_elapsed = clock() - c_start

print sat_obj.error, sat_obj.error_message
print "elapsed", c_elapsed

According to perf record here's what's using up more than 1% of CPU while computing a single position

40.36%  python   python2.7           [.] PyEval_EvalFrameEx
 5.85%  python   libm-2.23.so        [.] __cos_avx
 5.59%  python   libm-2.23.so        [.] __sin_avx
 5.06%  python   libc-2.23.so        [.] __sigsetjmp
 5.04%  python   python2.7           [.] PyNumber_Multiply
 3.12%  python   python2.7           [.] PyNumber_Add
 2.82%  python   python2.7           [.] PyFloat_FromDouble
 2.26%  python   libc-2.23.so        [.] __sigjmp_save
 1.20%  python   python2.7           [.] PyFPE_dummy
 1.13%  python   python2.7           [.] PyFloat_AsDouble
 1.01%  python   python2.7           [.] _setjmp@plt

[ckuethe]

I did a quick hack to add Numba support to sgp4 and it speeds things up about 8 to 10x. I might be able to make it a little faster but it'll cause an API break.

[ckuethe]

Sadly, that's still 6 seconds to compute a position (down from 60) and you really have to be computing a lot of positions because the LLVM jit takes about 30 seconds to run and then sgp4 takes another 6 seconds to load.

[ckuethe]

More observations.. with numba support hacked in:

• NOAA 19 takes 53 microseconds to propagate (down from 174 microseconds)

1 33591U 09005A   16207.49182166  .00000083  00000-0  69882-4 0  9996
2 33591  99.0428 164.5849 0015011 109.4452 250.8342 14.12091623384542

• GOES 15 takes about 6 seconds to propagate (down from 45 seconds)

1 36411U 10008A   16206.53472402 +.00000076 +00000-0 +00000-0 0  9999
2 36411 000.2488 098.3934 0000606 080.6320 180.9684 01.00279852023449

Crud. Now I'm going to have to do an analysis of propagate() time vs orbital period.

[ckuethe]

Deep space orbits (period greater than or equal to 225 minutes) are far more intensive to compute than near space, and the deep space resonance effects are the bulk of that. Even for objects as far out as GPS satellites that only make 2 orbits per day (period ~720min) calculations are pretty fast but once you get out to geosynchronous it can take a painfully long time for the resonance to complete.

[brandon-rhodes]

I have just released a new version of sgp4 that wraps the official C++ code (at least on systems where it can be installed and compiled), and have made a routine that does a fairly fast search for satellite events — because it's measuring a relative position between two Earth locations, it's able to turn off much of Skyfield's expensive accuracy:

https://rhodesmill.org/skyfield/earth-satellites.html#finding-when-a-satellite-rises-and-sets

Is anyone on this issue interested in installing the new version of Skyfield released yesterday and confirming independently that the speed should now be closer to PyEphem’s?

[astrojuanlu]

I have just released a new version of sgp4 that wraps the official C++ code (at least on systems where it can be installed and compiled),

I tried to quickly look for which specific version is it, or some release notes or relevant commits and couldn't find them, could you please share any pointers? :)

[brandon-rhodes]

Yes, the documentation and changelog are here:

https://pypi.org/project/sgp4/

So that I can make the experience more convenient for users: where did you look that lacked a link to the documentation? I'll be happy to add more links to the docs to help folks find them!

[astrojuanlu]

Let's continue in brandon-rhodes/python-sgp4#44 :)

[brandon-rhodes]

Thanks, I moved over there!

In the meantime, I'm still interested in whether anyone can confirm here in this thread that finding satellite passes now runs more quickly for them with the new release. Thanks for any data folks can provide!

[brandon-rhodes]

Pending more activity, I am closing this on the assumption that the new C++ accelerated SGP4 library has largely fixed the high cost of satellite positions in Skyfield. Please comment if you have more to add.