New Python iterator approach: vectorise without dtype

.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

cython_npyiter/.gitignore

@@ -0,0 +1 @@
+/erfa.c

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 (<NewNpyArrayIterObject*>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 = (<double *>(dataptrarray[0]))
+        _v = (<double *>(dataptrarray[1]))
+        _s = (<double *>(dataptrarray[2]))[0]
+        _bm1 = (<double *>(dataptrarray[3]))[0]
+        _ppr = (<double *>(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 = (<double *>(dataptrarray[0]))[0]
+        _astrom = (<eraASTROM *>(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 = (<double *>(dataptrarray[0]))[0]
+        _dc = (<double *>(dataptrarray[1]))[0]
+        _pr = (<double *>(dataptrarray[2]))[0]
+        _pd = (<double *>(dataptrarray[3]))[0]
+        _px = (<double *>(dataptrarray[4]))[0]
+        _rv = (<double *>(dataptrarray[5]))[0]
+        _utc1 = (<double *>(dataptrarray[6]))[0]
+        _utc2 = (<double *>(dataptrarray[7]))[0]
+        _dut1 = (<double *>(dataptrarray[8]))[0]
+        _elong = (<double *>(dataptrarray[9]))[0]
+        _phi = (<double *>(dataptrarray[10]))[0]
+        _hm = (<double *>(dataptrarray[11]))[0]
+        _xp = (<double *>(dataptrarray[12]))[0]
+        _yp = (<double *>(dataptrarray[13]))[0]
+        _phpa = (<double *>(dataptrarray[14]))[0]
+        _tc = (<double *>(dataptrarray[15]))[0]
+        _rh = (<double *>(dataptrarray[16]))[0]
+        _wl = (<double *>(dataptrarray[17]))[0]
+        _aob = (<double *>(dataptrarray[18]))
+        _zob = (<double *>(dataptrarray[19]))
+        _hob = (<double *>(dataptrarray[20]))
+        _dob = (<double *>(dataptrarray[21]))
+        _rob = (<double *>(dataptrarray[22]))
+        _eo = (<double *>(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 = (<int *>(dataptrarray[0]))[0]
+        _d1 = (<double *>(dataptrarray[1]))[0]
+        _d2 = (<double *>(dataptrarray[2]))[0]
+        _iy = (<int *>(dataptrarray[3]))
+        _im = (<int *>(dataptrarray[4]))
+        _id = (<int *>(dataptrarray[5]))
+        _ihmsf = (<int *>(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 = (<double *>(dataptrarray[0]))[0]
+        _r = (<double *>(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 = (<double *>(dataptrarray[0]))
+        eraIr(_r)
+        status = iternext(GetNpyIter(it))
+
+    return r_out