search

phn / pytpm

Python interface to the Telescope Pointing Machine C library.
http://phn.github.io/pytpm/
Other
11 stars 4 forks source link

performance issue after few calls to convertv6 #2

Open zonca opened 12 years ago

zonca commented 12 years ago

I noticed that after few calls, convertv6 slows down considerably, see a test case in this gist:

https://gist.github.com/1657277

I am performing the same call few times, at the beginning it is fast, but then slows down significantly. Can you please take a look at this issue?

thanks!

zonca commented 12 years ago

the first few calls give also a different result with respect to the others.

phn commented 12 years ago

Oh wow! This is really bad.

Speed

This is much worse on my system.

TIME[0]:4.9 s
[ 138.3027046   -6.8823031]
TIME[1]:0 s
[ 139.11112245   -0.73624182]
TIME[2]:0 s
[ 139.10879366   -0.75403965]
TIME[3]:0 s
[ 139.10874357   -0.75442249]
TIME[4]:0.01 s
[ 139.1087425    -0.75443064]
TIME[5]:0 s
[ 139.10874248   -0.75443081]
TIME[6]:0 s
[ 139.10874248   -0.75443082]
TIME[7]:0 s
[ 139.10874248   -0.75443082]
TIME[8]:0.01 s
[ 139.10874248   -0.75443082]
TIME[9]:18 s
[ 139.10874248   -0.75443082]
TIME[10]:36 s
[ 139.10874248   -0.75443082]
TIME[11]:36 s
[ 139.10874248   -0.75443082]
TIME[12]:36 s
[ 139.10874248   -0.75443082]
TIME[13]:36 s
[ 139.10874248   -0.75443082]
TIME[14]:36 s
[ 139.10874248   -0.75443082]

The speed issue can be solved by using a list of V6C vectors. This will run the loop at C level and not at Python level.

# The following line is ridiculous; cat2v6 should just accept a list.
# It is supposed to, but doesn't. It does return a list though. Sorry about this.
# I will fix this.
v6_list = [convert.cat2v6(alpha=az, delta=el)[0] for i in range(0, 15)]

start_clock = time.clock()
v6c_list = convert.convertv6(v6=v6_list,
        utc=ut,
        s1=tpm.TPM_S19, s2=tpm.TPM_S07,
        epoch=tpm.J2000, equinox=tpm.J2000,
        lon=lon, lat=lat, alt=alt,
        xpole=0.0, ypole=0.0,
        T=273.15, P=1013.25, H=0.0, wavelength=19986.16386)
print("TIME:%.2g s" % (time.clock() - start_clock))
cat = convert.v62cat(v6c_list)
print(np.degrees([(c['alpha'], c['delta']) for c in cat]))

# run-2 to see if output changes
v6_list = [convert.cat2v6(alpha=az, delta=el)[0] for i in range(0, 15)]
start_clock = time.clock()
v6c_list = convert.convertv6(v6=v6_list,
        utc=ut,
        s1=tpm.TPM_S19, s2=tpm.TPM_S07,
        epoch=tpm.J2000, equinox=tpm.J2000,
        lon=lon, lat=lat, alt=alt,
        xpole=0.0, ypole=0.0,
        T=273.15, P=1013.25, H=0.0, wavelength=19986.16386)
print("TIME:%.2g s" % (time.clock() - start_clock))
cat = convert.v62cat(v6c_list)
print(np.degrees([(c['alpha'], c['delta']) for c in cat]))

Output is:

TIME:5 s
[[ 138.3027046   -6.8823031]
 [ 138.3027046   -6.8823031]
 [ 138.3027046   -6.8823031]
 [ 138.3027046   -6.8823031]
 [ 138.3027046   -6.8823031]
 [ 138.3027046   -6.8823031]
 [ 138.3027046   -6.8823031]
 [ 138.3027046   -6.8823031]
 [ 138.3027046   -6.8823031]
 [ 138.3027046   -6.8823031]
 [ 138.3027046   -6.8823031]
 [ 138.3027046   -6.8823031]
 [ 138.3027046   -6.8823031]
 [ 138.3027046   -6.8823031]
 [ 138.3027046   -6.8823031]]
TIME:0 s
[[ 139.11112245   -0.73624182]
 [ 139.11112245   -0.73624182]
 [ 139.11112245   -0.73624182]
 [ 139.11112245   -0.73624182]
 [ 139.11112245   -0.73624182]
 [ 139.11112245   -0.73624182]
 [ 139.11112245   -0.73624182]
 [ 139.11112245   -0.73624182]
 [ 139.11112245   -0.73624182]
 [ 139.11112245   -0.73624182]
 [ 139.11112245   -0.73624182]
 [ 139.11112245   -0.73624182]
 [ 139.11112245   -0.73624182]
 [ 139.11112245   -0.73624182]
 [ 139.11112245   -0.73624182]]

Repeatability

This is just unacceptable. My initial guess is that this has something to do with how I use Python lists. But the problem doesn't go away even if I delete the lists, and start new.

Even the following code doesn't solve the problem:

def f():
    v6_list = [convert.cat2v6(alpha=az, delta=el)[0] for i in range(0, 15)]
    start_clock = time.clock()
    v6c_list = convert.convertv6(v6=v6_list,
            utc=ut,
            s1=tpm.TPM_S19, s2=tpm.TPM_S07,
            epoch=tpm.J2000, equinox=tpm.J2000,
            lon=lon, lat=lat, alt=alt,
            xpole=0.0, ypole=0.0,
            T=273.15, P=1013.25, H=0.0, wavelength=19986.16386)
    print("TIME:%.2g s" % (time.clock() - start_clock))
    cat = convert.v62cat(v6c_list)
    print(np.degrees([(c['alpha'], c['delta']) for c in cat]))

f()
f()

Output is:

TIME:5 s
[[ 138.3027046   -6.8823031]
 [ 138.3027046   -6.8823031]
 [ 138.3027046   -6.8823031]
 [ 138.3027046   -6.8823031]
 [ 138.3027046   -6.8823031]
 [ 138.3027046   -6.8823031]
 [ 138.3027046   -6.8823031]
 [ 138.3027046   -6.8823031]
 [ 138.3027046   -6.8823031]
 [ 138.3027046   -6.8823031]
 [ 138.3027046   -6.8823031]
 [ 138.3027046   -6.8823031]
 [ 138.3027046   -6.8823031]
 [ 138.3027046   -6.8823031]
 [ 138.3027046   -6.8823031]]
TIME:0 s
[[ 139.11112245   -0.73624182]
 [ 139.11112245   -0.73624182]
 [ 139.11112245   -0.73624182]
 [ 139.11112245   -0.73624182]
 [ 139.11112245   -0.73624182]
 [ 139.11112245   -0.73624182]
 [ 139.11112245   -0.73624182]
 [ 139.11112245   -0.73624182]
 [ 139.11112245   -0.73624182]
 [ 139.11112245   -0.73624182]
 [ 139.11112245   -0.73624182]
 [ 139.11112245   -0.73624182]
 [ 139.11112245   -0.73624182]
 [ 139.11112245   -0.73624182]
 [ 139.11112245   -0.73624182]]

Plan

For some reason, probably "backward compatibility" and/or laziness, I have been reluctant to use Numpy ndarray. Perhaps using Numpy arrays, and accessing the underlying C data from within Cython, will speed up the code. And it might also solve the repeatability issue.

I will start working on these issues within the next few days.

Thanks very much for the feedback. I don't have any opportunities for using the code, even though I have written it. This makes it hard to see if there is anything wrong with the code, or to even find out if it is user friendly. This is the first feedback I have received (for the usable/useful version of the code); so thanks again.

Prasanth

phn commented 12 years ago

Hello,

I was wondering which observatory and telescope you are referring to in the calculations i.e., altitude of 36 Km and wavelength of 2cm?

I tried the same calculation with smaller altitudes, a few 1000 meters, and there is no problem; at-least no apparent problem.

Code:

import datetime
import time
import numpy as np

from pytpm import convert, tpm

az = 3.30084818
el = 0.94610742
lat = 34.64
lon = -103.7
alt = 2500.0  # 35800.26
ut = 2455822.20000367
wavelength = 19986.16386

for i in range(0, 15):
    v6 = convert.cat2v6(alpha=az, delta=el)
    start_clock = time.clock()
    v6c = convert.convertv6(v6=v6,
        utc=ut,
        s1=tpm.TPM_S19, s2=tpm.TPM_S07,
        epoch=tpm.J2000, equinox=tpm.J2000,
        lon=lon, lat=lat, alt=alt,
        xpole=0.0, ypole=0.0,
        T=273.15, P=1013.25, H=0.0, wavelength=wavelength)
    print("TIME[%d]:%.2g s" % (i, time.clock() - start_clock))
    cat = convert.v62cat(v6c)
    print(np.degrees([cat['alpha'], cat['delta']]))

Output:

TIME[0]:0 s
[ 139.10720416   -0.76619028]
TIME[1]:0 s
[ 139.10720416   -0.76619028]
TIME[2]:0 s
[ 139.10720416   -0.76619028]
TIME[3]:0 s
[ 139.10720416   -0.76619028]
TIME[4]:0 s
[ 139.10720416   -0.76619028]
TIME[5]:0 s
[ 139.10720416   -0.76619028]
TIME[6]:0 s
[ 139.10720416   -0.76619028]
TIME[7]:0 s
[ 139.10720416   -0.76619028]
TIME[8]:0 s
[ 139.10720416   -0.76619028]
TIME[9]:0 s
[ 139.10720416   -0.76619028]
TIME[10]:0 s
[ 139.10720416   -0.76619028]
TIME[11]:0 s
[ 139.10720416   -0.76619028]
TIME[12]:0 s
[ 139.10720416   -0.76619028]
TIME[13]:0 s
[ 139.10720416   -0.76619028]
TIME[14]:0 s
[ 139.10720416   -0.76619028]

This is of-course not a solution, but is definitely something to keep in mind while trying to see what caused the problem.

Thanks, Prasanth

zonca commented 12 years ago

It is a balloon-borne microwave radiometer On Jan 26, 2012 11:04 PM, "Prasanth" < reply@reply.github.com> wrote:

Hello,

I was wondering which observatory and telescope you are referring to in the calculations i.e., altitude of 36 Km and wavelength of 2cm?

I tried the same calculation with smaller altitudes, a few 1000 meters, and there is no problem; at-least no apparent problem.

Code:

import datetime import time import numpy as np

from pytpm import convert, tpm

az = 3.30084818 el = 0.94610742 lat = 34.64 lon = -103.7 alt = 2500.0 # 35800.26 ut = 2455822.20000367 wavelength = 19986.16386

for i in range(0, 15): v6 = convert.cat2v6(alpha=az, delta=el) start_clock = time.clock() v6c = convert.convertv6(v6=v6, utc=ut, s1=tpm.TPM_S19, s2=tpm.TPM_S07, epoch=tpm.J2000, equinox=tpm.J2000, lon=lon, lat=lat, alt=alt, xpole=0.0, ypole=0.0, T=273.15, P=1013.25, H=0.0, wavelength=wavelength) print("TIME[%d]:%.2g s" % (i, time.clock() - start_clock)) cat = convert.v62cat(v6c) print(np.degrees([cat['alpha'], cat['delta']]))

Output:

TIME[0]:0 s [ 139.10720416 -0.76619028] TIME[1]:0 s [ 139.10720416 -0.76619028] TIME[2]:0 s [ 139.10720416 -0.76619028] TIME[3]:0 s [ 139.10720416 -0.76619028] TIME[4]:0 s [ 139.10720416 -0.76619028] TIME[5]:0 s [ 139.10720416 -0.76619028] TIME[6]:0 s [ 139.10720416 -0.76619028] TIME[7]:0 s [ 139.10720416 -0.76619028] TIME[8]:0 s [ 139.10720416 -0.76619028] TIME[9]:0 s [ 139.10720416 -0.76619028] TIME[10]:0 s [ 139.10720416 -0.76619028] TIME[11]:0 s [ 139.10720416 -0.76619028] TIME[12]:0 s [ 139.10720416 -0.76619028] TIME[13]:0 s [ 139.10720416 -0.76619028] TIME[14]:0 s [ 139.10720416 -0.76619028]

This is of-course not a solution, but is definitely something to keep in mind while trying to see what caused the problem.

Thanks, Prasanth

Reply to this email directly or view it on GitHub: https://github.com/phn/pytpm/issues/2#issuecomment-3682870

phn commented 12 years ago

Hello,

Of the coordinates you have, are there any that correspond to known objects. It would be great if so could you send me the observed coordinates for a few of these. This will let me test if the conversion, even if slow, is given valid results.

Thanks. Prasanth

zonca commented 12 years ago

Hi, I have a test available here:

https://gist.github.com/1672906

Comparing PyEphem to JPL Horizons, you can use this as a test.