diff --git a/.gitignore b/.gitignore index 1bb95c0..0041ebd 100644 --- a/.gitignore +++ b/.gitignore @@ -1,5 +1,14 @@ /build .project - .pydevproject + +*.py[cod] + +cython_numpy_auto/erfa.c + +cython_numpy_auto/erfa.pyx + +cython_npyiter_auto/erfa.c + +cython_npyiter_auto/erfa.pyx diff --git a/cython_npyiter/.gitignore b/cython_npyiter/.gitignore new file mode 100644 index 0000000..a72b5f9 --- /dev/null +++ b/cython_npyiter/.gitignore @@ -0,0 +1 @@ +/erfa.c diff --git a/cython_npyiter/erfa.pyx b/cython_npyiter/erfa.pyx new file mode 100644 index 0000000..55606ca --- /dev/null +++ b/cython_npyiter/erfa.pyx @@ -0,0 +1,325 @@ +import numpy as np +cimport numpy as np +from cpython.ref cimport PyObject + +np.import_array() + +__all__ = ['ab', 'aper', 'atco13', 'd2dtf', 'rx'] + + +cdef extern from "erfam.h": + cdef struct eraASTROM: + pass + +cdef extern from "erfa.h": + void eraAb(double pnat[3], double v[3], double s, double bm1, double ppr[3]) + void eraAper(double theta, eraASTROM *astrom) + int eraAtco13(double rc, double dc, + double pr, double pd, double px, double rv, + double utc1, double utc2, double dut1, + double elong, double phi, double hm, double xp, double yp, + double phpa, double tc, double rh, double wl, + double *aob, double *zob, double *hob, + 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 eraRx(double phi, double r[9]) + void eraIr(double*) + +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) + +ctypedef void NpyIter +ctypedef int (*IterNextFunc)(NpyIter * iter) nogil + +cdef extern from "numpy/arrayobject.h": + IterNextFunc GetIterNext "NpyIter_GetIterNext" (NpyIter *iter, char **) + char** GetDataPtrArray "NpyIter_GetDataPtrArray" (NpyIter* iter) + +ctypedef struct NewNpyArrayIterObject: + PyObject base + NpyIter *iter + +cdef inline NpyIter* GetNpyIter(object iter): + return (iter).iter + + +def ab(pnat, v, s, bm1): + + #Turn all inputs into arrays + pnat_in = np.array(pnat, dtype=np.double, copy=False, subok=True) + v_in = np.array(v, dtype=np.double, copy=False, subok=True) + s_in = np.array(s, dtype=np.double, copy=False, subok=True) + bm1_in = np.array(bm1, dtype=np.double, copy=False, subok=True) + + #Create the output array, based on the broadcasted shape, adding the generated dimensions if needed + broadcast = np.broadcast(pnat_in[...,0], v_in[...,0], s_in, bm1_in) + ppr_out = np.empty(broadcast.shape+(3,), dtype=np.double) + + #Create the iterator, broadcasting on all but the consumed dimensions + op_axes = [[-1]*(broadcast.nd-arr.ndim)+range(arr.ndim) for arr in [pnat_in[...,0], v_in[...,0], s_in, bm1_in, ppr_out[...,0]]] + op_flags = [['readonly']]*4+[['readwrite']]*1 + it = np.nditer([pnat_in, v_in, s_in, bm1_in, ppr_out], op_axes=op_axes, op_flags=op_flags) + + #Iterate + cdef double* _pnat + cdef double* _v + cdef double _s + cdef double _bm1 + cdef double* _ppr + cdef char** dataptrarray = GetDataPtrArray(GetNpyIter(it)) + cdef IterNextFunc iternext = GetIterNext(GetNpyIter(it), NULL) + cdef int status = 1 + while status: + _pnat = ((dataptrarray[0])) + _v = ((dataptrarray[1])) + _s = ((dataptrarray[2]))[0] + _bm1 = ((dataptrarray[3]))[0] + _ppr = ((dataptrarray[4])) + eraAb(_pnat, _v, _s, _bm1, _ppr) + status = iternext(GetNpyIter(it)) + + return ppr_out + + +def aper(theta, astrom): + + #Turn all inputs into arrays + theta_in = np.array(theta, dtype=np.double, copy=False, subok=True) + astrom_in = np.array(astrom, dtype=dt_eraASTROM, copy=False, subok=True) + + #Create the output array, based on the broadcasted shape, adding the generated dimensions if needed + broadcast = np.broadcast(theta_in, astrom_in) + astrom_out = np.empty(broadcast.shape, dtype=dt_eraASTROM) + np.copyto(astrom_out, astrom_in) + + #Create the iterator, broadcasting on all but the consumed dimensions + op_axes = [[-1]*(broadcast.nd-arr.ndim)+range(arr.ndim) for arr in [theta_in, astrom_out]] + op_flags = [['readonly']]*1+[['readwrite']]*1 + it = np.nditer([theta_in, astrom_out], op_axes=op_axes, op_flags=op_flags) + + #Iterate + cdef double _theta + cdef eraASTROM* _astrom + cdef char** dataptrarray = GetDataPtrArray(GetNpyIter(it)) + cdef IterNextFunc iternext = GetIterNext(GetNpyIter(it), NULL) + cdef int status = 1 + while status: + _theta = ((dataptrarray[0]))[0] + _astrom = ((dataptrarray[1])) + eraAper(_theta, _astrom) + status = iternext(GetNpyIter(it)) + + return astrom_out + + +def atco13(rc, dc, pr, pd, px, rv, utc1, utc2, dut1, elong, phi, hm, xp, yp, phpa, tc, rh, wl): + + #Turn all inputs into arrays + rc_in = np.array(rc, dtype=np.double, copy=False, subok=True) + dc_in = np.array(dc, dtype=np.double, copy=False, subok=True) + pr_in = np.array(pr, dtype=np.double, copy=False, subok=True) + pd_in = np.array(pd, dtype=np.double, copy=False, subok=True) + px_in = np.array(px, dtype=np.double, copy=False, subok=True) + rv_in = np.array(rv, dtype=np.double, copy=False, subok=True) + utc1_in = np.array(utc1, dtype=np.double, copy=False, subok=True) + utc2_in = np.array(utc2, dtype=np.double, copy=False, subok=True) + dut1_in = np.array(dut1, dtype=np.double, copy=False, subok=True) + elong_in = np.array(elong, dtype=np.double, copy=False, subok=True) + phi_in = np.array(phi, dtype=np.double, copy=False, subok=True) + hm_in = np.array(hm, dtype=np.double, copy=False, subok=True) + xp_in = np.array(xp, dtype=np.double, copy=False, subok=True) + yp_in = np.array(yp, dtype=np.double, copy=False, subok=True) + phpa_in = np.array(phpa, dtype=np.double, copy=False, subok=True) + tc_in = np.array(tc, dtype=np.double, copy=False, subok=True) + rh_in = np.array(rh, dtype=np.double, copy=False, subok=True) + wl_in = np.array(wl, dtype=np.double, copy=False, subok=True) + + #Create the output array, based on the broadcasted shape, adding the generated dimensions if needed + broadcast = np.broadcast(rc_in, dc_in, pr_in, pd_in, px_in, rv_in, utc1_in, utc2_in, dut1_in, elong_in, phi_in, hm_in, xp_in, yp_in, phpa_in, tc_in, rh_in, wl_in) + aob_out = np.empty(broadcast.shape, dtype=np.double) + zob_out = np.empty(broadcast.shape, dtype=np.double) + hob_out = np.empty(broadcast.shape, dtype=np.double) + dob_out = np.empty(broadcast.shape, dtype=np.double) + rob_out = np.empty(broadcast.shape, dtype=np.double) + eo_out = np.empty(broadcast.shape, dtype=np.double) + + #Create the iterator, broadcasting on all but the consumed dimensions + op_axes = [[-1]*(broadcast.nd-arr.ndim)+range(arr.ndim) for arr in [rc_in, dc_in, pr_in, pd_in, px_in, rv_in, utc1_in, utc2_in, dut1_in, elong_in, phi_in, hm_in, xp_in, yp_in, phpa_in, tc_in, rh_in, wl_in, aob_out, zob_out, hob_out, dob_out, rob_out, eo_out]] + op_flags = [['readonly']]*18+[['readwrite']]*6 + it = np.nditer([rc_in, dc_in, pr_in, pd_in, px_in, rv_in, utc1_in, utc2_in, dut1_in, elong_in, phi_in, hm_in, xp_in, yp_in, phpa_in, tc_in, rh_in, wl_in, aob_out, zob_out, hob_out, dob_out, rob_out, eo_out], op_axes=op_axes, op_flags=op_flags) + + #Iterate + 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 _tc + cdef double _rh + cdef double _wl + cdef double* _aob + cdef double* _zob + cdef double* _hob + cdef double* _dob + cdef double* _rob + cdef double* _eo + cdef char** dataptrarray = GetDataPtrArray(GetNpyIter(it)) + cdef IterNextFunc iternext = GetIterNext(GetNpyIter(it), NULL) + cdef int status = 1 + while status: + _rc = ((dataptrarray[0]))[0] + _dc = ((dataptrarray[1]))[0] + _pr = ((dataptrarray[2]))[0] + _pd = ((dataptrarray[3]))[0] + _px = ((dataptrarray[4]))[0] + _rv = ((dataptrarray[5]))[0] + _utc1 = ((dataptrarray[6]))[0] + _utc2 = ((dataptrarray[7]))[0] + _dut1 = ((dataptrarray[8]))[0] + _elong = ((dataptrarray[9]))[0] + _phi = ((dataptrarray[10]))[0] + _hm = ((dataptrarray[11]))[0] + _xp = ((dataptrarray[12]))[0] + _yp = ((dataptrarray[13]))[0] + _phpa = ((dataptrarray[14]))[0] + _tc = ((dataptrarray[15]))[0] + _rh = ((dataptrarray[16]))[0] + _wl = ((dataptrarray[17]))[0] + _aob = ((dataptrarray[18])) + _zob = ((dataptrarray[19])) + _hob = ((dataptrarray[20])) + _dob = ((dataptrarray[21])) + _rob = ((dataptrarray[22])) + _eo = ((dataptrarray[23])) + ret = eraAtco13(_rc, _dc, _pr, _pd, _px, _rv, _utc1, _utc2, _dut1, _elong, _phi, _hm, _xp, _yp, _phpa, _tc, _rh, _wl, _aob, _zob, _hob, _dob, _rob, _eo) + status = iternext(GetNpyIter(it)) + + return aob_out, zob_out, hob_out, dob_out, rob_out, eo_out + + +def d2dtf(scale, ndp, d1, d2): + + #Turn all inputs into arrays + ndp_in = np.array(ndp, dtype=np.int32, copy=False, subok=True) + d1_in = np.array(d1, dtype=np.double, copy=False, subok=True) + d2_in = np.array(d2, dtype=np.double, copy=False, subok=True) + + #Create the output array, based on the broadcasted shape, adding the generated dimensions if needed + broadcast = np.broadcast(ndp_in, d1_in, d2_in) + iy_out = np.empty(broadcast.shape, dtype=np.int32) + im_out = np.empty(broadcast.shape, dtype=np.int32) + id_out = np.empty(broadcast.shape, dtype=np.int32) + ihmsf_out = np.empty(broadcast.shape+(4,), dtype=np.int32) + + #Create the iterator, broadcasting on all but the consumed dimensions + op_axes = [[-1]*(broadcast.nd-arr.ndim)+range(arr.ndim) for arr in [ndp_in, d1_in, d2_in, iy_out, im_out, id_out, ihmsf_out[...,0]]] + op_flags = [['readonly']]*3+[['readwrite']]*4 + it = np.nditer([ndp_in, d1_in, d2_in, iy_out, im_out, id_out, ihmsf_out], op_axes=op_axes, op_flags=op_flags) + + #Iterate + cdef char* _scale = scale + cdef int _ndp + cdef double _d1 + cdef double _d2 + cdef int* _iy + cdef int* _im + cdef int* _id + cdef int* _ihmsf + cdef char** dataptrarray = GetDataPtrArray(GetNpyIter(it)) + cdef IterNextFunc iternext = GetIterNext(GetNpyIter(it), NULL) + cdef int status = 1 + while status: + _ndp = ((dataptrarray[0]))[0] + _d1 = ((dataptrarray[1]))[0] + _d2 = ((dataptrarray[2]))[0] + _iy = ((dataptrarray[3])) + _im = ((dataptrarray[4])) + _id = ((dataptrarray[5])) + _ihmsf = ((dataptrarray[6])) + ret = eraD2dtf(_scale, _ndp, _d1, _d2, _iy, _im, _id, _ihmsf) + status = iternext(GetNpyIter(it)) + + return iy_out, im_out, id_out, ihmsf_out + + +def rx(phi, r): + + #Turn all inputs into arrays + phi_in = np.array(phi, dtype=np.double, copy=False, subok=True) + r_in = np.array(r, dtype=np.double, copy=False, subok=True) + + #Create the output array, based on the broadcasted shape, adding the generated dimensions if needed + broadcast = np.broadcast(phi_in, r_in[...,0,0]) + r_out = np.empty(broadcast.shape+(3, 3), dtype=np.double) + np.copyto(r_out, r_in) + + #Create the iterator, broadcasting on all but the consumed dimensions + op_axes = [[-1]*(broadcast.nd-arr.ndim)+range(arr.ndim) for arr in [phi_in, r_out[...,0,0]]] + op_flags = [['readonly']]*1+[['readwrite']]*1 + it = np.nditer([phi_in, r_out], op_axes=op_axes, op_flags=op_flags) + + #Iterate + cdef double _phi + cdef double* _r + cdef char** dataptrarray = GetDataPtrArray(GetNpyIter(it)) + cdef IterNextFunc iternext = GetIterNext(GetNpyIter(it), NULL) + cdef int status = 1 + while status: + _phi = ((dataptrarray[0]))[0] + _r = ((dataptrarray[1])) + eraRx(_phi, _r) + status = iternext(GetNpyIter(it)) + + return r_out + + +def ir(): + + #Turn all inputs into arrays + + #Create the output array, based on the broadcasted shape, adding the generated dimensions if needed + broadcast = np.broadcast(np.int32(0.0), np.int32(0.0)) + r_out = np.empty(broadcast.shape+(3, 3), dtype=np.double) + + #Create the iterator, broadcasting on all but the consumed dimensions + op_axes = [[-1]*(broadcast.nd-arr.ndim)+range(arr.ndim) for arr in []+[r_out[...,0,0]]] + op_flags = [['readonly']]*0+[['readwrite']]*1 + it = np.nditer([]+[r_out], op_axes=op_axes, op_flags=op_flags) + + #Iterate + cdef double* _r + cdef char** dataptrarray = GetDataPtrArray(GetNpyIter(it)) + cdef IterNextFunc iternext = GetIterNext(GetNpyIter(it), NULL) + cdef int status = 1 + while status: + _r = ((dataptrarray[0])) + eraIr(_r) + status = iternext(GetNpyIter(it)) + + return r_out diff --git a/cython_npyiter_auto/code_analyzer.py b/cython_npyiter_auto/code_analyzer.py new file mode 100644 index 0000000..a53f65e --- /dev/null +++ b/cython_npyiter_auto/code_analyzer.py @@ -0,0 +1,132 @@ +import os.path +import re + + +__all__ = ['Function'] + + +ctype_to_dtype = {'double' : "np.double", + 'double *' : "np.double", + 'int' : "np.int32", + 'int *' : "np.int32", + 'eraASTROM' : "dt_eraASTROM", + 'char' : "np.dtype('S1')", + } + + +class Argument(object): + + def __init__(self, definition, inout_state): + self.inout_state = inout_state + self.definition = definition.strip() + if self.definition[:5] == "const": + self.definition = self.definition[6:].strip() + if "*" in self.definition: + self.ctype, self.name = self.definition.split("*", 1) + self.ptr = "*" + self.shape = () + elif "[]" in self.definition: + self.ctype, self.name = self.definition.split("[]", 1) + self.ptr = "*" + self.shape = () + else: + self.ctype, self.name = self.definition.rsplit(" ", 1) + if "[" in self.name: + self.name, arr = self.name.split("[", 1) + self.ptr = "*" + self.shape = tuple([int(s) for s in arr[:-1].split("][")]) + else: + self.ptr = "" + self.shape = () + self.ctype = self.ctype.strip() + + @property + def name_in(self): + return self.name+"_in" + + @property + def name_out(self): + return self.name+"_out" + + @property + def name_in_broadcast(self): + if len(self.shape)>0: + return self.name_in+"[..."+",0"*len(self.shape)+"]" + else: + return self.name_in + + @property + def name_out_broadcast(self): + if len(self.shape)>0: + return self.name_out+"[..."+",0"*len(self.shape)+"]" + else: + return self.name_out + + @property + def ctype_ptr(self): + return self.ctype+self.ptr + + @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 Return(object): + + def __init__(self, ctype): + self.name = 'ret' + self.ctype = ctype + self.inout_state = '=' + self.ctype_ptr = ctype + + @property + def dtype(self): + return ctype_to_dtype[self.ctype] + + +class Function(object): + """ + The following attributes are derived from the supplied ERFA function + name: C name of the function + pyname: Python name of the function + signature: Signature of the function, as extracted from the source file + documentation: Documentation block for the function, as extracted from the source file + args: List of arguments of function + ret: Return parameter of function + """ + + def __init__(self, name, inout_state, source_path): + self.name = name + self.pyname = name.split('era')[-1].lower() + #From source get: signature and documentation + pattern = "\n([^\n]+{0} ?\([^)]+\)).+?(/\*.+?\*/)".format(name) + p = re.compile(pattern, flags=re.DOTALL|re.MULTILINE) + filepath = os.path.join(os.path.normpath(source_path), self.pyname+".c") + with open(filepath) as f: + search = p.search(f.read()) + self.signature = " ".join(search.group(1).split()) #Single line with single spaces + self.documentation = search.group(2).replace("**"," ").replace("/*"," ").replace("*/"," ") + #Extract arguments from signature + self.args = [] + arg_list = self.signature.split(')')[0].split('(')[1].split(', ') + for arg_index in range(len(arg_list)): + self.args.append(Argument(arg_list[arg_index], inout_state[arg_index])) + self.ret = re.search("^(.*){0}".format(name), self.signature).group(1).strip() + if self.ret == 'double': + self.args.append(Return(self.ret)) + + 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 diff --git a/cython_npyiter_auto/config.py b/cython_npyiter_auto/config.py new file mode 100644 index 0000000..0bccb8c --- /dev/null +++ b/cython_npyiter_auto/config.py @@ -0,0 +1,26 @@ +#Convention is the following: +# i: input +# o: output +# m: modify +# r: return + +wrappers = {} +#int eraAtco13(double rc, double dc, +# double pr, double pd, double px, double rv, +# double utc1, double utc2, double dut1, +# double elong, double phi, double hm, double xp, double yp, +# double phpa, double tk, double rh, double wl, +# double *aob, double *zob, double *hob, +# double *dob, double *rob, double *eo) +wrappers["eraAtco13"] = "iiiiiiiiiiiiiiiiiioooooo" +#void eraRx(double phi, double r[3][3]) +wrappers["eraRx"] = "im" +#int eraD2dtf(const char *scale, int ndp, double d1, double d2, int *iy, int *im, int *id, int ihmsf[4]) +wrappers["eraD2dtf"] = "iiiioooo" +#void eraAper(double theta, eraASTROM *astrom) +wrappers["eraAper"] = "im" +#void eraAb(double pnat[3], double v[3], double s, double bm1, +# double ppr[3]) +wrappers["eraAb"] = "iiiio" +#void eraIr(double r[3][3]) +wrappers["eraIr"] = "o" \ No newline at end of file diff --git a/cython_npyiter_auto/cython_generator.py b/cython_npyiter_auto/cython_generator.py new file mode 100644 index 0000000..cf1d51e --- /dev/null +++ b/cython_npyiter_auto/cython_generator.py @@ -0,0 +1,54 @@ +import re +from jinja2 import Template, Environment, PackageLoader +from code_analyzer import Function + +import config + +ERFA_SOURCES = "../../erfa/src" + + +#Prepare the jinja2 templating environment + +def prefix(a_list, pre): + return [pre+'{0}'.format(an_element) for an_element in a_list] + +def postfix(a_list, post): + return ['{0}'.format(an_element)+post for an_element in a_list] + +def surround(a_list, pre, post): + return [pre+'{0}'.format(an_element)+post for an_element in a_list] + +env = Environment() +env.filters['prefix'] = prefix +env.filters['postfix'] = postfix +env.filters['surround'] = surround + + +#Extract all the ERFA function names from erfa.h + +with open(ERFA_SOURCES+"/erfa.h", "r") as f: + erfa_h = f.read() +print("Parse ERFA code:") +section_subsection_functions = re.findall('/\* (\w*)/(\w*) \*/\n(.*?)\n\n', erfa_h, flags=re.DOTALL|re.MULTILINE) +funcs = [] +for section, subsection, functions in section_subsection_functions: + print(" {0}.{1}".format(section, subsection)) + func_names = re.findall(' (\w+)\(.*?\);', functions, flags=re.DOTALL) + for name in func_names: + if name in config.wrappers.keys(): + function = Function(name, config.wrappers[name], ERFA_SOURCES) + funcs.append(function) + print(" {0}.{1}.{2} >{3}<...".format(section, subsection, name, function.signature)) +print("Done!") + + +#Render the template + +with open('erfa.pyx.jinja2', 'r') as template_file: + erfa_pyx = env.from_string(template_file.read()).render(funcs=funcs) + + +#Save the result + +with open("erfa.pyx", "w") as f: + f.write(erfa_pyx) diff --git a/cython_npyiter_auto/erfa.pyx.jinja2 b/cython_npyiter_auto/erfa.pyx.jinja2 new file mode 100644 index 0000000..3029490 --- /dev/null +++ b/cython_npyiter_auto/erfa.pyx.jinja2 @@ -0,0 +1,90 @@ +import numpy as np +cimport numpy as np +from cpython.ref cimport PyObject + +np.import_array() + +__all__ = [{{ funcs|map(attribute='pyname')|surround("'","'")|join(", ") }}] + + +cdef extern from "erfam.h": + cdef struct eraASTROM: + pass + +cdef extern from "erfa.h": +{%- for func in funcs %} + {{ func.ret }} {{ func.name }}({{ func.args_by_inout('i|m|o')|map(attribute='ctype_ptr')|join(', ') }}) +{%- endfor %} + +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) + +ctypedef void NpyIter +ctypedef int (*IterNextFunc)(NpyIter * iter) nogil + +cdef extern from "numpy/arrayobject.h": + IterNextFunc GetIterNext "NpyIter_GetIterNext" (NpyIter *iter, char **) + char** GetDataPtrArray "NpyIter_GetDataPtrArray" (NpyIter* iter) + +ctypedef struct NewNpyArrayIterObject: + PyObject base + NpyIter *iter + +cdef inline NpyIter* GetNpyIter(object iter): + return (iter).iter + + +{% for func in funcs %} +def {{ func.pyname }}({{ func.args_by_inout('i|m')|map(attribute='name')|join(', ') }}): + + #Turn all inputs into arrays + {%- for arg in func.args_by_inout('i|m') %} + {{ arg.name_in }} = np.array({{ arg.name }}, dtype={{ arg.dtype }}, copy=False, subok=True) + {%- endfor %} + + #Create the output array, based on the broadcasted shape, adding the generated dimensions if needed + broadcast = np.broadcast(np.int32(0.0), np.int32(0.0), {{ func.args_by_inout('i|m')|map(attribute='name_in_broadcast')|join(', ') }}) + {%- for arg in func.args_by_inout('m|o|r') %} + {{ arg.name_out }} = np.empty(broadcast.shape+{{ arg.shape }}, dtype={{ arg.dtype }}) + {%- endfor %} + {%- for arg in func.args_by_inout('m') %} + np.copyto({{ arg.name_out }}, {{ arg.name_in }}) + {%- endfor %} + + #Create the iterator, broadcasting on all but the consumed dimensions + op_axes = [[-1]*(broadcast.nd-arr.ndim)+range(arr.ndim) for arr in [{{ func.args_by_inout('i')|map(attribute='name_in_broadcast')|join(', ') }}]+[{{ func.args_by_inout('m|o')|map(attribute='name_out_broadcast')|join(', ') }}]] + op_flags = [['readonly']]*{{ func.args_by_inout('i')|count }}+[['readwrite']]*{{ func.args_by_inout('m|o')|count }} + it = np.nditer([{{ func.args_by_inout('i')|map(attribute='name_in')|join(', ') }}]+[{{ func.args_by_inout('m|o')|map(attribute='name_out')|join(', ') }}], op_axes=op_axes, op_flags=op_flags) + + #Iterate + {%- for arg in func.args_by_inout('i|m|o') %} + cdef {{ arg.ctype_ptr }} _{{ arg.name }} + {%- endfor %} + cdef char** dataptrarray = GetDataPtrArray(GetNpyIter(it)) + cdef IterNextFunc iternext = GetIterNext(GetNpyIter(it), NULL) + cdef int status = 1 + while status: + {%- for arg in func.args_by_inout('i|m|o') %} + _{{ arg.name }} = (<{{ arg.ctype }} *>(dataptrarray[{{ func.args.index(arg) }}])){%- if arg.ctype_ptr[-1] != '*' %}[0]{% endif %} + {%- endfor %} + {{ func.args_by_inout('ret')|map(attribute='name')|surround('_',' = ')|join }}{{ func.name }}({{ func.args_by_inout('i|m|o')|map(attribute='name')|prefix('_')|join(', ') }}) + status = iternext(GetNpyIter(it)) + + return {{ func.args_by_inout('m|o|r')|map(attribute='name_out')|join(', ') }} + +{% endfor %} diff --git a/cython_numpy/erfa.pyx b/cython_numpy/erfa.pyx index 6504f38..87a3493 100644 --- a/cython_numpy/erfa.pyx +++ b/cython_numpy/erfa.pyx @@ -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, @@ -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): @@ -58,47 +117,76 @@ def atco13(rc, dc, pr, pd, px, rv, utc1, utc2, dut1, elong, phi, hm, xp, yp, php _tk = (np.PyArray_MultiIter_DATA(it, 15))[0] _rh = (np.PyArray_MultiIter_DATA(it, 16))[0] _wl = (np.PyArray_MultiIter_DATA(it, 17))[0] + _aob = (np.PyArray_MultiIter_DATA(it, 18)) + _zob = (np.PyArray_MultiIter_DATA(it, 19)) + _hob = (np.PyArray_MultiIter_DATA(it, 20)) + _dob = (np.PyArray_MultiIter_DATA(it, 21)) + _rob = (np.PyArray_MultiIter_DATA(it, 22)) + _eo = (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) - - (np.PyArray_MultiIter_DATA(it, 18))[0] = _aob - (np.PyArray_MultiIter_DATA(it, 19))[0] = _zob - (np.PyArray_MultiIter_DATA(it, 20))[0] = _hob - (np.PyArray_MultiIter_DATA(it, 21))[0] = _dob - (np.PyArray_MultiIter_DATA(it, 22))[0] = _rob - (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.int32) + im_out = np.empty(shape, dtype=np.int32) + id_out = np.empty(shape, dtype=np.int32) + ihmsf_out = np.empty(shape, dtype=[('h','i4'),('m','i4'),('s','i4'),('f','i4')]) - 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 = (np.PyArray_MultiIter_DATA(it, 0))[0] - _d2 = (np.PyArray_MultiIter_DATA(it, 1))[0] + _scale = ( np.PyArray_MultiIter_DATA(it, 0)) + _ndp = ( np.PyArray_MultiIter_DATA(it, 1))[0] + _d1 = (np.PyArray_MultiIter_DATA(it, 2))[0] + _d2 = (np.PyArray_MultiIter_DATA(it, 3))[0] + _iy = ( np.PyArray_MultiIter_DATA(it, 4)) + _im = ( np.PyArray_MultiIter_DATA(it, 5)) + _id = ( np.PyArray_MultiIter_DATA(it, 6)) + _ihmsf = ( 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 = ( np.PyArray_MultiIter_DATA(it, 0))[0] + _astrom = (np.PyArray_MultiIter_DATA(it, 1)) - (np.PyArray_MultiIter_DATA(it, 2))[0] = _iy - (np.PyArray_MultiIter_DATA(it, 3))[0] = _im - (np.PyArray_MultiIter_DATA(it, 4))[0] = _id - (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 \ No newline at end of file diff --git a/cython_numpy_auto/__init__.py b/cython_numpy_auto/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/cython_numpy_auto/code_analyzer.py b/cython_numpy_auto/code_analyzer.py new file mode 100644 index 0000000..bde0348 --- /dev/null +++ b/cython_numpy_auto/code_analyzer.py @@ -0,0 +1,177 @@ +import os.path +import re + + +__all__ = ['Function'] + + +ctype_to_dtype = {'double' : "numpy.double", + 'double *' : "numpy.double", + 'int' : "numpy.int", + 'int *' : "numpy.int", + 'int[4]' : "numpy.dtype([('', 'i', (4,))])", + 'double[2]' : "numpy.dtype([('', 'd', (2,))])", + 'double[3]' : "numpy.dtype([('p', 'd', (3,))])", + 'double[2][3]' : "numpy.dtype([('pv', 'd', (2,3))])", + 'double[3][3]' : "numpy.dtype([('r', 'd', (3,3))])", + 'eraASTROM *' : "dt_eraASTROM", + 'char *' : "numpy.dtype('S1')", + } + + +class FunctionDoc(object): + + 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 = [] + result = re.search("Given([^\n]*):\n(.+?) \n", self.doc, re.DOTALL) + if result is not None: + __input = result.group(2) + for i in __input.split("\n"): + arg_doc = ArgumentDoc(i) + if arg_doc.name is not None: + self.__input.append(arg_doc) + result = re.search("Given and returned([^\n]*):\n(.+?) \n", self.doc, re.DOTALL) + if result is not None: + __input = result.group(2) + 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 = [] + result = re.search("Returned([^\n]*):\n(.+?) \n", self.doc, re.DOTALL) + if result is not None: + __output = result.group(2) + for i in __output.split("\n"): + arg_doc = ArgumentDoc(i) + if arg_doc.name is not None: + self.__output.append(arg_doc) + result = re.search("Given and returned([^\n]*):\n(.+?) \n", self.doc, re.DOTALL) + if result is not None: + __output = result.group(2) + 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(object): + + 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(object): + + def __init__(self, definition, doc): + self.__doc = doc + self.__inout_state = None + self.definition = definition.strip() + if "*" in self.definition: + self.ctype, self.name = self.definition.split("*", 1) + self.ctype += "*" + else: + self.ctype, self.name = self.definition.rsplit(" ", 1) + if "[" in self.name: + self.name, arr = self.name.split("[", 1) + self.ctype += ("["+arr) + + @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 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 Return(object): + + def __init__(self, ctype, doc): + self.name = 'ret' + self.ctype = ctype + self.inout_state = 'ret' + self.ctype_ptr = ctype + + @property + def dtype(self): + return ctype_to_dtype[self.ctype] + +class Function(object): + + def __init__(self, name, source_path): + self.name = name + self.pyname = name.split('era')[-1].lower() + self.filename = name.split("era")[-1].lower()+".c" + 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)) + self.ret = re.search("^(.*){0}".format(name), self.cfunc).group(1).strip() + if self.ret == 'double': + self.args.append(Return(self.ret, self.__doc)) + + 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 diff --git a/cython_numpy_auto/cython_generator.py b/cython_numpy_auto/cython_generator.py new file mode 100644 index 0000000..3b227fe --- /dev/null +++ b/cython_numpy_auto/cython_generator.py @@ -0,0 +1,41 @@ +import re +from jinja2 import Environment, PackageLoader +from cython_numpy_auto.code_analyzer import Function + +ERFA_SOURCES = "../../erfa/src" + +#Prepare the jinja2 templating environment +env = Environment(loader=PackageLoader('cython_numpy_auto', '.')) + +def prefix(a_list, pre): + return [pre+'{0}'.format(an_element) for an_element in a_list] +def postfix(a_list, post): + return ['{0}'.format(an_element)+post for an_element in a_list] +def surround(a_list, pre, post): + return [pre+'{0}'.format(an_element)+post for an_element in a_list] +env.filters['prefix'] = prefix +env.filters['postfix'] = postfix +env.filters['surround'] = surround + +erfa_pyx_in = env.get_template('erfa.pyx.in') + + +#Extract all the ERFA function names from erfa.h +with open(ERFA_SOURCES+"/erfa.h", "r") as f: + + erfa_h = f.read() + + funcs = [] + section_subsection_functions = re.findall('/\* (\w*)/(\w*) \*/\n(.*?)\n\n', erfa_h, flags=re.DOTALL|re.MULTILINE) + for section, subsection, functions in section_subsection_functions: + print("{0}.{1}".format(section, subsection)) + if section == "Astronomy": + func_names = re.findall(' (\w+)\(.*?\);', functions, flags=re.DOTALL) + for name in func_names: + print("{0}.{1}.{2}...".format(section, subsection, name)) + funcs.append(Function(name, ERFA_SOURCES)) + print("Done!") + #Render the template and save + erfa_pyx = erfa_pyx_in.render(funcs=funcs) + with open("erfa.pyx", "w") as f: + f.write(erfa_pyx) diff --git a/cython_numpy_auto/erfa.pyx.in b/cython_numpy_auto/erfa.pyx.in new file mode 100644 index 0000000..3b65992 --- /dev/null +++ b/cython_numpy_auto/erfa.pyx.in @@ -0,0 +1,99 @@ +import numpy +cimport numpy + +numpy.import_array() + +__all__ = [{{ funcs|map(attribute='name')|surround("'","'")|join(", ") }}] + + +cdef extern from "erfa.h": + struct eraASTROM: + pass + struct eraLDBODY: + pass +{%- for func in funcs %} + {{ func.ret }} {{ func.name }}({{ func.args_by_inout('in|inout|out')|map(attribute='ctype_ptr')|join(', ') }}) +{%- endfor %} + +dt_eraASTROM = numpy.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) + + +{% for func in funcs %} +def {{ func.pyname }}({{ func.args_by_inout('in|inout')|map(attribute='name')|join(', ') }}): + + {%- if func.args_by_inout('in|inout') %} + shape = numpy.broadcast({{ func.args_by_inout('in|inout')|map(attribute='name')|join(', ') }}).shape + {%- for arg in func.args_by_inout('out|ret') %} + {{ arg.name }}_out = numpy.empty(shape, dtype={{ arg.dtype }}) + {%- endfor %} + {%- for arg in func.args_by_inout('inout') %} + {{ arg.name }}_out = numpy.array(numpy.broadcast_arrays({{ func.args_by_inout('in|inout')|map(attribute='name')|join(', ') }})[{{ func.args.index(arg) }}], dtype={{ arg.dtype }}) + {%- endfor %} + + cdef numpy.broadcast it = numpy.broadcast({{ (func.args_by_inout('in')|map(attribute='name')|list + func.args_by_inout('inout|out|ret')|map(attribute='name')|postfix('_out')|list) |join(', ') }}) + {%- for arg in func.args_by_inout('in|inout|out|ret') %} + cdef {{ arg.ctype_ptr }} _{{ arg.name }} + {%- endfor %} + + while numpy.PyArray_MultiIter_NOTDONE(it): + + {%- for arg in func.args_by_inout('in|inout|out') %} + {%- if arg.ctype_ptr[-1] == '*' %} + _{{ arg.name }} = (<{{ arg.ctype_ptr }}>numpy.PyArray_MultiIter_DATA(it, {{ func.args.index(arg) }})) + {%- else %} + _{{ arg.name }} = (<{{ arg.ctype_ptr }}*>numpy.PyArray_MultiIter_DATA(it, {{ func.args.index(arg) }}))[0] + {%- endif %} + {%- endfor %} + + {{ func.args_by_inout('ret')|map(attribute='name')|surround('_',' = ')|join }}{{ func.name }}({{ func.args_by_inout('in|out|inout')|map(attribute='name')|prefix('_')|join(', ') }}) + + {%- for arg in func.args_by_inout('ret') %} + (<{{ arg.ctype_ptr }}*>numpy.PyArray_MultiIter_DATA(it, {{ func.args.index(arg) }}))[0] = _{{ arg.name }} + {%- endfor %} + + numpy.PyArray_MultiIter_NEXT(it) + {%- else %} + shape = [] + {%- for arg in func.args_by_inout('out') %} + {{ arg.name }}_out = numpy.empty(shape, dtype={{ arg.dtype }}) + {%- endfor %} + + cdef numpy.broadcast it = numpy.broadcast({{ (func.args_by_inout('in')|map(attribute='name')|list + func.args_by_inout('inout|out|ret')|map(attribute='name')|postfix('_out')|list) |join(', ') }}) + {%- for arg in func.args_by_inout('out') %} + cdef {{ arg.ctype_ptr }} _{{ arg.name }} + {%- endfor %} + + while numpy.PyArray_MultiIter_NOTDONE(it): + + {%- for arg in func.args_by_inout('in|inout|out') %} + {%- if arg.ctype_ptr[-1] == '*' %} + _{{ arg.name }} = (<{{ arg.ctype_ptr }}>numpy.PyArray_MultiIter_DATA(it, {{ func.args.index(arg) }})) + {%- else %} + _{{ arg.name }} = (<{{ arg.ctype_ptr }}*>numpy.PyArray_MultiIter_DATA(it, {{ func.args.index(arg) }}))[0] + {%- endif %} + {%- endfor %} + + {{ func.name }}({{ func.args_by_inout('in|out|inout')|map(attribute='name')|prefix('_')|join(', ') }}) + + numpy.PyArray_MultiIter_NEXT(it) + {%- endif %} + + return {{ func.args_by_inout('inout|out|ret')|map(attribute='name')|postfix('_out')|join(', ') }} + +{% endfor %} diff --git a/setup.py b/setup.py index acfd3b0..8f298bf 100644 --- a/setup.py +++ b/setup.py @@ -10,8 +10,26 @@ include_dirs = [np.get_include(), '/opt/local/include'], library_dirs = ['/opt/local/lib']) +cython_numpy_auto_ext = Extension("erfa.cython_numpy_auto", + ["cython_numpy_auto/erfa.pyx"], + libraries=['erfa'], + include_dirs = [np.get_include(), '/opt/local/include'], + library_dirs = ['/opt/local/lib']) + +cython_npyiter_ext = Extension("erfa.cython_npyiter", + ["cython_npyiter/erfa.pyx"], + libraries=['erfa'], + include_dirs = [np.get_include(), '/opt/local/include'], + library_dirs = ['/opt/local/lib']) + +cython_npyiter_auto_ext = Extension("erfa.cython_npyiter_auto", + ["cython_npyiter_auto/erfa.pyx"], + libraries=['erfa'], + include_dirs = [np.get_include(), '/opt/local/include'], + library_dirs = ['/opt/local/lib']) + setup(name = "erfa", cmdclass = { "build_ext": build_ext }, packages = ["erfa"], package_dir = {"erfa":"."}, - ext_modules = [cython_numpy_ext]) + ext_modules = [cython_numpy_ext, cython_numpy_auto_ext, cython_npyiter_ext, cython_npyiter_auto_ext]) diff --git a/tests/test_erfa.py b/tests/test_erfa.py index 80653bc..cd4b52a 100644 --- a/tests/test_erfa.py +++ b/tests/test_erfa.py @@ -1,24 +1,37 @@ +import importlib import numpy as np -import erfa.cython_numpy as erfa - -jd = np.linspace(2456855.5, 2456855.5+1.0/24.0/60.0, 60*2+1) -ra = np.linspace(0.0,np.pi*2.0,5) -dec = np.linspace(-np.pi/2.0,+np.pi/2.0,4) - -aob, zob, hob, dob, rob, eo = erfa.atco13(0.0,0.0,0.0,0.0,0.0,0.0,jd,0.0,0.0,0.0,np.pi/4.0,0.0,0.0,0.0,1014.0,0.0,0.0,0.5) -print(aob.shape) - -aob, zob, hob, dob, rob, eo = erfa.atco13(0.0,0.0,0.0,0.0,0.0,0.0,jd[0],0.0,0.0,0.0,np.pi/4.0,0.0,0.0,0.0,1014.0,0.0,0.0,0.5) -print(aob.shape) - -aob, zob, hob, dob, rob, eo = erfa.atco13(ra[:,None,None],dec[None,:,None],0.0,0.0,0.0,0.0,jd[None,None,:],0.0,0.0,0.0,np.pi/4.0,0.0,0.0,0.0,1014.0,0.0,0.0,0.5) -print(aob.shape) - -iy, im, id, ihmsf = erfa.d2dtf("UTC", 3, jd, 0.0) -print(iy.shape, ihmsf.shape, ihmsf.dtype) -print(ihmsf) - -iy, im, id, ihmsf = erfa.d2dtf("UTC", 3, jd[0], 0.0) -print(iy.shape, ihmsf.shape, ihmsf.dtype) -print(ihmsf) +for erfa_wrapper_name in ['cython_numpy', 'cython_numpy_auto', 'cython_npyiter']: + + erfa_module_name = '.'.join(['erfa', erfa_wrapper_name]) + print("Testing with {0}...".format(erfa_module_name)) + + erfa = importlib.import_module(erfa_module_name) + + jd = np.linspace(2456855.5, 2456855.5+1.0/24.0/60.0, 60*2+1) + ra = np.linspace(0.0,np.pi*2.0,5) + dec = np.linspace(-np.pi/2.0,+np.pi/2.0,4) + + aob, zob, hob, dob, rob, eo = erfa.atco13(0.0,0.0,0.0,0.0,0.0,0.0,jd,0.0,0.0,0.0,np.pi/4.0,0.0,0.0,0.0,1014.0,0.0,0.0,0.5) + print(aob.shape) + + aob, zob, hob, dob, rob, eo = erfa.atco13(0.0,0.0,0.0,0.0,0.0,0.0,jd[0],0.0,0.0,0.0,np.pi/4.0,0.0,0.0,0.0,1014.0,0.0,0.0,0.5) + print(aob.shape) + + aob, zob, hob, dob, rob, eo = erfa.atco13(ra[:,None,None],dec[None,:,None],0.0,0.0,0.0,0.0,jd[None,None,:],0.0,0.0,0.0,np.pi/4.0,0.0,0.0,0.0,1014.0,0.0,0.0,0.5) + print(aob.shape) + + iy, im, id, ihmsf = erfa.d2dtf("UTC", 3, jd, 0.0) + print(iy.shape, ihmsf.shape, ihmsf.dtype) + print(ihmsf) + + iy, im, id, ihmsf = erfa.d2dtf("UTC", 3, jd[0], 0.0) + 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)