Skip to content

More work on the python_numpy experiment #7

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
jwoillez wants to merge 6 commits into from
Closed

More work on the python_numpy experiment #7

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
627719c
Vectorize all d2dtf() imputs.
jwoillez Jul 28, 2014
b815ba7
Added aper()
jwoillez Jul 28, 2014
f1b7ce9
Added a test for aper()
jwoillez Jul 28, 2014
3209705
Postfixed all output arguments with '_out'
jwoillez Jul 28, 2014
64e3068
cdef and assign all C arguments before use
jwoillez Jul 28, 2014
0cdfa62
Use copy constructor array() and broadcast_arrays()...
jwoillez Jul 28, 2014
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
170 changes: 129 additions & 41 deletions cython_numpy/erfa.pyx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,8 @@ cimport numpy as np

np.import_array()

__all__ = ['atco13', 'd2dtf']
__all__ = ['atco13', 'd2dtf', 'aper']


cdef extern from "erfa.h":
int eraAtco13(double rc, double dc,
Expand All @@ -15,28 +16,86 @@ cdef extern from "erfa.h":
double *dob, double *rob, double *eo)
int eraD2dtf(const char *scale, int ndp, double d1, double d2,
int *iy, int *im, int *id, int ihmsf[4])
void eraAper(double theta, eraASTROM *astrom)

cdef struct eraASTROM:
double pmt
double eb[3]
double eh[3]
double em
double v[3]
double bm1
double bpn[3][3]
double along
double phi
double xpl
double ypl
double sphi
double cphi
double diurab
double eral
double refa
double refb

dt_eraASTROM = np.dtype([('pmt','d'),
('eb','d',(3,)),
('eh','d',(3,)),
('em','d'),
('v','d',(3,)),
('bm1 ','d'),
('bpn','d',(3,3)),
('along','d'),
('phi','d'),
('xpl','d'),
('ypl','d'),
('sphi','d'),
('cphi','d'),
('diurab','d'),
('eral','d'),
('refa','d'),
('refb','d')], align=True)


#Note: the pattern used here follows https://github.com/cython/cython/wiki/tutorials-numpy#dimensionally-simple-functions



def atco13(rc, dc, pr, pd, px, rv, utc1, utc2, dut1, elong, phi, hm, xp, yp, phpa, tk, rh, wl):

shape = np.broadcast(rc, dc, pr, pd, px, rv, utc1, utc2, dut1, elong, phi, hm, xp, yp, phpa, tk, rh, wl).shape
aob = np.empty(shape, dtype=np.double)
zob = np.empty(shape, dtype=np.double)
hob = np.empty(shape, dtype=np.double)
dob = np.empty(shape, dtype=np.double)
rob = np.empty(shape, dtype=np.double)
eo = np.empty(shape, dtype=np.double)
aob_out = np.empty(shape, dtype=np.double)
zob_out = np.empty(shape, dtype=np.double)
hob_out = np.empty(shape, dtype=np.double)
dob_out = np.empty(shape, dtype=np.double)
rob_out = np.empty(shape, dtype=np.double)
eo_out = np.empty(shape, dtype=np.double)

cdef np.broadcast it = np.broadcast(rc, dc, pr, pd, px, rv, utc1, utc2, dut1, elong, phi, hm, xp, yp, phpa, tk, rh, wl, aob, zob, hob, dob, rob, eo)
cdef np.broadcast it = np.broadcast(rc, dc, pr, pd, px, rv, utc1, utc2, dut1, elong, phi, hm, xp, yp, phpa, tk, rh, wl, aob_out, zob_out, hob_out, dob_out, rob_out, eo_out)

cdef double _aob
cdef double _zob
cdef double _hob
cdef double _dob
cdef double _rob
cdef double _eo
cdef double _rc
cdef double _dc
cdef double _pr
cdef double _pd
cdef double _px
cdef double _rv
cdef double _utc1
cdef double _utc2
cdef double _dut1
cdef double _elong
cdef double _phi
cdef double _hm
cdef double _xp
cdef double _yp
cdef double _phpa
cdef double _tk
cdef double _rh
cdef double _wl
cdef double *_aob
cdef double *_zob
cdef double *_hob
cdef double *_dob
cdef double *_rob
cdef double *_eo

while np.PyArray_MultiIter_NOTDONE(it):

Expand All @@ -58,47 +117,76 @@ def atco13(rc, dc, pr, pd, px, rv, utc1, utc2, dut1, elong, phi, hm, xp, yp, php
_tk = (<double*>np.PyArray_MultiIter_DATA(it, 15))[0]
_rh = (<double*>np.PyArray_MultiIter_DATA(it, 16))[0]
_wl = (<double*>np.PyArray_MultiIter_DATA(it, 17))[0]
_aob = (<double*>np.PyArray_MultiIter_DATA(it, 18))
_zob = (<double*>np.PyArray_MultiIter_DATA(it, 19))
_hob = (<double*>np.PyArray_MultiIter_DATA(it, 20))
_dob = (<double*>np.PyArray_MultiIter_DATA(it, 21))
_rob = (<double*>np.PyArray_MultiIter_DATA(it, 22))
_eo = (<double*>np.PyArray_MultiIter_DATA(it, 23))

ret = eraAtco13(_rc, _dc, _pr, _pd, _px, _rv, _utc1, _utc2, _dut1, _elong, _phi, _hm, _xp, _yp, _phpa, _tk, _rh, _wl, &_aob, &_zob, &_hob, &_dob, &_rob, &_eo)

(<double*>np.PyArray_MultiIter_DATA(it, 18))[0] = _aob
(<double*>np.PyArray_MultiIter_DATA(it, 19))[0] = _zob
(<double*>np.PyArray_MultiIter_DATA(it, 20))[0] = _hob
(<double*>np.PyArray_MultiIter_DATA(it, 21))[0] = _dob
(<double*>np.PyArray_MultiIter_DATA(it, 22))[0] = _rob
(<double*>np.PyArray_MultiIter_DATA(it, 23))[0] = _eo
ret = eraAtco13(_rc, _dc, _pr, _pd, _px, _rv, _utc1, _utc2, _dut1, _elong, _phi, _hm, _xp, _yp, _phpa, _tk, _rh, _wl, _aob, _zob, _hob, _dob, _rob, _eo)

np.PyArray_MultiIter_NEXT(it)

return aob, zob, hob, dob, rob, eo
return aob_out, zob_out, hob_out, dob_out, rob_out, eo_out



def d2dtf(scale, ndp, d1, d2):

shape = np.broadcast(d1, d2).shape
iy = np.empty(shape, dtype=np.int)
im = np.empty(shape, dtype=np.int)
id = np.empty(shape, dtype=np.int)
ihmsf = np.empty(shape, dtype=[('h','i'),('m','i'),('s','i'),('f','i')])
shape = np.broadcast(scale, ndp, d1, d2).shape
iy_out = np.empty(shape, dtype=np.int)
im_out = np.empty(shape, dtype=np.int)
id_out = np.empty(shape, dtype=np.int)
ihmsf_out = np.empty(shape, dtype=[('h','i'),('m','i'),('s','i'),('f','i')])

cdef np.broadcast it = np.broadcast(d1, d2, iy, im, id, ihmsf)
cdef np.broadcast it = np.broadcast(scale, ndp, d1, d2, iy_out, im_out, id_out, ihmsf_out)

cdef int _iy
cdef int _im
cdef int _id
cdef int _ihmsf[4]
cdef char *_scale
cdef int _ndp
cdef double _d1
cdef double _d2
cdef int *_iy
cdef int *_im
cdef int *_id
cdef int *_ihmsf

while np.PyArray_MultiIter_NOTDONE(it):

_d1 = (<double*>np.PyArray_MultiIter_DATA(it, 0))[0]
_d2 = (<double*>np.PyArray_MultiIter_DATA(it, 1))[0]
_scale = ( <char*>np.PyArray_MultiIter_DATA(it, 0))
_ndp = ( <int*>np.PyArray_MultiIter_DATA(it, 1))[0]
_d1 = (<double*>np.PyArray_MultiIter_DATA(it, 2))[0]
_d2 = (<double*>np.PyArray_MultiIter_DATA(it, 3))[0]
_iy = ( <int*>np.PyArray_MultiIter_DATA(it, 4))
_im = ( <int*>np.PyArray_MultiIter_DATA(it, 5))
_id = ( <int*>np.PyArray_MultiIter_DATA(it, 6))
_ihmsf = ( <int*>np.PyArray_MultiIter_DATA(it, 7))

ret = eraD2dtf(_scale, _ndp, _d1, _d2, _iy, _im, _id, _ihmsf)

np.PyArray_MultiIter_NEXT(it)

return iy_out, im_out, id_out, ihmsf_out



def aper(theta, astrom):

shape = np.broadcast(theta, astrom).shape
astrom_out = np.array(np.broadcast_arrays(theta, astrom)[1], dtype=dt_eraASTROM)

cdef np.broadcast it = np.broadcast(theta, astrom_out)

cdef double _theta
cdef eraASTROM *_astrom

while np.PyArray_MultiIter_NOTDONE(it):

ret = eraD2dtf(scale, ndp, _d1, _d2, &_iy, &_im, &_id, _ihmsf)
_theta = ( <double *>np.PyArray_MultiIter_DATA(it, 0))[0]
_astrom = (<eraASTROM *>np.PyArray_MultiIter_DATA(it, 1))

(<int*>np.PyArray_MultiIter_DATA(it, 2))[0] = _iy
(<int*>np.PyArray_MultiIter_DATA(it, 3))[0] = _im
(<int*>np.PyArray_MultiIter_DATA(it, 4))[0] = _id
(<int*>np.PyArray_MultiIter_DATA(it, 5))[0:4] = _ihmsf
eraAper(_theta, _astrom)

np.PyArray_MultiIter_NEXT(it)

return iy, im, id, ihmsf
return astrom_out
6 changes: 6 additions & 0 deletions tests/test_erfa.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,3 +22,9 @@
print(iy.shape, ihmsf.shape, ihmsf.dtype)
print(ihmsf)

astrom = np.zeros([2],dtype=erfa.dt_eraASTROM)
theta = np.arange(0,10.0)
print(theta.shape)
print(astrom.shape)
astrom = erfa.aper(theta[:,None], astrom[None,:])
print(astrom)