Skip to content
Open
Show file tree
Hide file tree
Changes from 12 commits
Commits
Show all changes
37 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
9c9d296
Added *.py[cod] to ignore list
jwoillez Jul 29, 2014
3dbc8f8
Initial commit of code analyzer
jwoillez Jul 29, 2014
7523521
Initial commit of cython_generator
jwoillez Jul 29, 2014
2a60933
Ignore erfa.pyx and erfa.c in cython_numpy_auto
jwoillez Jul 29, 2014
6e653c1
Added cython_numpy_auto extension to setup.py
jwoillez Jul 29, 2014
bcbe709
Run test against both cython_numpy and cython_numpy_auto
jwoillez Jul 29, 2014
d1dcc3b
Base all classes on `object` for attribute access.
jwoillez Jul 30, 2014
729d2e0
Turn `cython_numpy_auto` into a package
jwoillez Jul 30, 2014
b6da8b8
Use jinja2 for code generation
jwoillez Jul 30, 2014
1ec74bc
Add more ctype to dtype conversions
jwoillez Jul 31, 2014
940753a
Deal with empty input/output sections in doc
jwoillez Jul 31, 2014
41cbfe1
Argument documentations have variable leading spaces
jwoillez Jul 31, 2014
7506ff0
Deal with multidimensional C arrays
jwoillez Jul 31, 2014
7bbcabb
Pointer arguments are always of input type
jwoillez Jul 31, 2014
922cba1
Remove dead code: is_in(), is_out(), is_inout()
jwoillez Jul 31, 2014
fa3612b
Implement Return() class to use along regular arguments
jwoillez Jul 31, 2014
1ebb796
Spaces between function name and opening parenthesis
jwoillez Jul 31, 2014
3d35218
Add returned parameter at end of arguments list
jwoillez Jul 31, 2014
0de22b9
np instead of numpy collides with c code arguments
jwoillez Jul 31, 2014
763c08c
Three new jinja filters: prefix, postfix, surround
jwoillez Jul 31, 2014
d03c0c1
Process all ERFA functions
jwoillez Jul 31, 2014
c9be43e
Add eraLDBODY structure support
jwoillez Jul 31, 2014
56d80ca
Use prefix, postfix, and surround in template.
jwoillez Jul 31, 2014
e3ccefd
Re-generate extern function signature with pointers
jwoillez Jul 31, 2014
e641fd4
Add returned values
jwoillez Jul 31, 2014
dade3e0
Functions with no inputs do not need vectorisation
jwoillez Jul 31, 2014
e16ff9a
Improvements for no input functions
jwoillez Jul 31, 2014
080325d
Pointer arguments are actually not necessarily inputs
jwoillez Jul 31, 2014
9981fb0
Deal with "Given and returned" documentation sections
jwoillez Jul 31, 2014
a85fcf9
Remove unwanted print statements
jwoillez Jul 31, 2014
c43c83c
Identify sections/subsections and process only Astronomy section
jwoillez Jul 31, 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
7 changes: 6 additions & 1 deletion .gitignore
Original file line number Diff line number Diff line change
@@ -1,5 +1,10 @@
/build

.project

.pydevproject

*.py[cod]

cython_numpy_auto/erfa.c

cython_numpy_auto/erfa.pyx
170 changes: 129 additions & 41 deletions cython_numpy/erfa.pyx
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
Copy link

Choose a reason for hiding this comment

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

So, this is now the generated file? It should probably be moved so it doesn't clobber the manually written one already in the repo at this location.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Actually, no. I have just improved the cython_numpy approach independently from the auto-generation. One of the goal was to add support for aper, one of the tricky wrapper as it has an inout parameter, and also uses this eraASTROM structure.

If you want to see the auto-generated wrapper, your have to run cython_generator.py in the python_numpy_auto experiment.

Copy link

Choose a reason for hiding this comment

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

Ok. I see, that's fine.

152 changes: 152 additions & 0 deletions cython_numpy_auto/code_analyzer.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,152 @@
import os.path
import re


__all__ = ['Function']


ctype_to_dtype = {'double' : "np.double",
'double *' : "np.double",
'int' : "np.int",
'int *' : "np.int",
'int[4]' : "np.dtype([('', 'i')]*4)",
'eraASTROM *': "dt_eraASTROM",
}


class FunctionDoc:

def __init__(self, doc):
self.doc = doc.replace("**"," ").replace("/*"," ").replace("*/"," ")
self.__input = None
self.__output = None

@property
def input(self):
if self.__input is None:
self.__input = []
__input = re.search("Given:\n(.+?) \n", self.doc, re.DOTALL).group(1)
for i in __input.split("\n"):
arg_doc = ArgumentDoc(i)
if arg_doc.name is not None:
self.__input.append(arg_doc)
return self.__input

@property
def output(self):
if self.__output is None:
self.__output = []
__output = re.search("Returned:\n(.+?) \n", self.doc, re.DOTALL).group(1)
for i in __output.split("\n"):
arg_doc = ArgumentDoc(i)
if arg_doc.name is not None:
self.__output.append(arg_doc)
return self.__output

def __repr__(self):
return self.doc

class ArgumentDoc:

def __init__(self, doc):
match = re.search("^ ([^ ]+)[ ]+([^ ]+)[ ]+(.+)", doc)
if match is not None:
self.name = match.group(1)
self.type = match.group(2)
self.doc = match.group(3)
else:
self.name = None
self.type = None
self.doc = None

def __repr__(self):
return " {0:15} {1:15} {2}".format(self.name, self.type, self.doc)

class Argument:

def __init__(self, definition, doc):
self.__doc = doc
self.__inout_state = None
self.definition = definition.strip()
if self.definition[-1] == "]":
self.ctype, self.name = self.definition.split(" ",1)
self.name, arr = self.name.split("[")
self.ctype += ("["+arr)
elif "*" in self.definition:
self.ctype, self.name = self.definition.split("*", 1)
self.ctype += "*"
else:
self.ctype, self.name = self.definition.split(" ", 1)

@property
def inout_state(self):
if self.__inout_state is None:
self.__inout_state = ''
for i in self.__doc.input:
if self.name in i.name.split(','):
self.__inout_state = 'in'
for o in self.__doc.output:
if self.name in o.name.split(','):
if self.__inout_state == 'in':
self.__inout_state = 'inout'
else:
self.__inout_state = 'out'
return self.__inout_state

@property
def is_in(self):
return self.inout_state == 'in'

@property
def is_out(self):
return self.inout_state == 'out'

@property
def is_inout(self):
return self.inout_state == 'inout'

@property
def ctype_ptr(self):
if self.ctype[-1] == ']':
return self.ctype.split('[')[0]+" *"
elif self.ctype[:6] == 'const ':
return self.ctype[6:]
else:
return self.ctype

@property
def dtype(self):
return ctype_to_dtype[self.ctype]

def __repr__(self):
return "Argument('{0}', name='{1}', ctype='{2}', inout_state='{3}')".format(self.definition, self.name, self.ctype, self.inout_state)

class Function:

def __init__(self, name, source_path):
self.name = name
self.pyname = name.split('era')[-1].lower()
self.filename = name.split("era")[-1].lower()+".c"
Copy link

Choose a reason for hiding this comment

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

If I'm following along correctly, this means that it's opening the .c file to find the signature for the function? It should probably open the .h file, as a system installation of erfa will generally have the headers, but not the source code.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Correct. The reason for parsing the .c file instead is that SOFA/ERFA has properly documented the input/output status of each parameter. It would also let us generate the documentation for each auto-generated function.

Copy link

Choose a reason for hiding this comment

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

Ah -- having gone and looked at the SOFA/ERFA code, I see what you're saying now, and I see why that's the only real alternative. Not to sound like a curmudgeon, but SOFA is totally doing it wrong. The .h files should include the documentation, not the .c files.

At the end of the day, it probably matters little, because I suspect we would generate the .pyx file, and its resulting .c file through Cython, and ship that in our tarball, so the end user won't need the SOFA/ERFA .c files to build astropy. But it might be worth a note to the SOFA folks about this.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I tend to agree with you on the documentation location. We could also relocate the documentation when deriving ERFA from SOFA. This could help those that will use or develop with ERFA directly.

Will open a separate issue against ERFA.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

You already did! :)
liberfa/erfa#24

self.filepath = os.path.join(os.path.normpath(source_path), self.filename)
pattern = "\n([^\n]+{0}\([^)]+\)).+?(/\*.+?\*/)".format(name)
p = re.compile(pattern, flags=re.DOTALL|re.MULTILINE)
with open(self.filepath) as f:
search = p.search(f.read())
self.cfunc = search.group(1)
self.__doc = FunctionDoc(search.group(2))
self.args = []
for arg in re.search("\(([^)]+)\)", self.cfunc, flags=re.MULTILINE|re.DOTALL).group(1).split(','):
self.args.append(Argument(arg, self.__doc))
Copy link

Choose a reason for hiding this comment

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

I'm impressed that we're able to get away with regexs to parse the argument types here -- I suppose that's fine given the regularity of ERFA. We could take this one step further and actually find what's in the headers and generate functions for everything -- and use something like pycparser to get that done. But this is probably fine for now.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Using the .h to auto-generate a full wrapper would indeed be the way to go.


def args_by_inout(self, inout_filter, prop=None, join=None):
result = []
for arg in self.args:
if arg.inout_state in inout_filter.split('|'):
if prop is None:
result.append(arg)
else:
result.append(getattr(arg, prop))
if join is not None:
return join.join(result)
else:
return result
Loading