From e12134fc3dbbf77c53350024ab153590721bc109 Mon Sep 17 00:00:00 2001 From: Chris Hogan Date: Tue, 10 Oct 2017 12:48:13 -0500 Subject: [PATCH 01/43] Rewrite mp_fourier in pure C --- tractor/Makefile | 8 +- tractor/c_mp_fourier.i | 261 ++++ tractor/numpy.i | 3166 ++++++++++++++++++++++++++++++++++++++++ tractor/setup-mpf.py | 7 +- 4 files changed, 3440 insertions(+), 2 deletions(-) create mode 100644 tractor/c_mp_fourier.i create mode 100644 tractor/numpy.i diff --git a/tractor/Makefile b/tractor/Makefile index 145c19bb..853bcf1c 100644 --- a/tractor/Makefile +++ b/tractor/Makefile @@ -7,12 +7,15 @@ PYTHON_SO_EXT ?= $(shell $(PYTHON) -c "from distutils import sysconfig; print(sy ceres: _ceres$(PYTHON_SO_EXT) ceres.py .PHONY: ceres -mpf: mp_fourier +mpf: mp_fourier c_mp_fourier .PHONY: mpf mp_fourier: _mp_fourier$(PYTHON_SO_EXT) mp_fourier.py .PHONY: mp_fourier +c_mp_fourier: _c_mp_fourier$(PYTHON_SO_EXT) c_mp_fourier.py +.PHONY: c_mp_fourier + mix: _mix$(PYTHON_SO_EXT) mix.py .PHONY: mix @@ -22,6 +25,9 @@ emfit: _emfit$(PYTHON_SO_EXT) emfit.py _mp_fourier$(PYTHON_SO_EXT): mp_fourier.i setup-mpf.py $(PYTHON) setup-mpf.py build_ext --inplace +_c_mp_fourier$(PYTHON_SO_EXT): c_mp_fourier.i setup-mpf.py + $(PYTHON) setup-mpf.py build_ext --inplace + mix.py _mix$(PYTHON_SO_EXT): mix.i approx3.c gauss_masked.c setup-mix.py $(PYTHON) setup-mix.py build_ext --inplace diff --git a/tractor/c_mp_fourier.i b/tractor/c_mp_fourier.i new file mode 100644 index 00000000..dd44fee8 --- /dev/null +++ b/tractor/c_mp_fourier.i @@ -0,0 +1,261 @@ +%module(package="tractor") c_mp_fourier + +%{ +#define SWIG_FILE_WITH_INIT +#define _GNU_SOURCE 1 +#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION +#include +#include +#include +%} + +%include "numpy.i" + +%init %{ + // numpy + import_array(); +%} + +%apply (double* IN_ARRAY1, int DIM1) { + (double *amps, int amps_len), + (double *v, int v_len), + (double *w, int w_len) +}; +%apply (double* IN_ARRAY2, int DIM1, int DIM2) { + (double *means, int means_dim1, int means_dim2) +}; +%apply (double* IN_ARRAY3, int DIM1, int DIM2, int DIM3) { + (double *vars, int vars_dim1, int vars_dim2, int vars_dim3) +}; +%apply (double* INPLACE_ARRAY2, int DIM1, int DIM2) { + (double *out, int out_dim1, int out_dim2) +}; + +%inline %{ + +static void mixture_profile_fourier_transform(double *amps, int amps_len, + double *means, int means_dim1, int means_dim2, + double *vars, int vars_dim1, int vars_dim2, int vars_dim3, + double *v, int v_len, + double *w, int w_len, + double *out, int out_dim1, int out_dim2) +{ + // const int D = 2; + const double twopisquare = -2. * M_PI * M_PI; + + int K = amps_len; + int NV = v_len; + int NW = w_len; + int i, j, k; + + // assert(means_dim1 == K); + // assert(means_dim2 == D); + // assert(vars_dim1 == K); + // assert(vars_dim2 == D); + // assert(vars_dim3 == D); + + // for (k = 0; k < K; k++) { + // if ((means[k*D] != means[0]) || + // (means[k*D+1] != means[1])) { + // PyErr_SetString(PyExc_ValueError, "Assume all means are equal"); + // return NULL; + // } + // } + + double mu0 = means[0]; + double mu1 = means[1]; + + int zeromean = (mu0 == 0.) && (mu1 == 0.); + + // TODO: Handle complex out + + if (zeromean) { + // out = (double*)malloc(sizeof(double) * NW * NV); + } + // else { + // out = new std::complex[NW * NV]; + // } + + double* ff = out; + + for (j = 0; j < NW; j++) { + for (i = 0; i < NV; i++) { + double s = 0; + double* V = vars; + for (k = 0; k < K; k++) { + double a, b, d; + a = *V; + V++; + b = *V; + V++; + // skip c + V++; + d = *V; + V++; + + s += amps[k] * exp(twopisquare * (a * v[i] * v[i] + + 2. * b * v[i] * w[j] + + d * w[j] * w[j])); + } + if (zeromean) { + *ff = s; + ff++; + } else { + double angle = -2. * M_PI * (mu0 * v[i] + mu1 * w[j]); + *ff = s * cos(angle); + ff++; + *ff = s * sin(angle); + ff++; + } + } + } + return; +} + +%} + +// static PyObject* mixture_profile_fourier_transform( +// PyObject* po_amps, +// PyObject* po_means, +// PyObject* po_vars, +// PyObject* po_v, +// PyObject* po_w +// ) { + + // npy_intp K, NW,NV; + // const npy_intp D = 2; + // npy_intp i,j,k; + // double* amps, *means, *vars, *vv, *ww; + // PyObject* np_F; + // double* f; + // npy_intp dims[2]; + + // PyArrayObject *np_amps, *np_means, *np_vars, *np_v, *np_w; + + // if (!PyArray_Check(po_amps) || !PyArray_Check(po_means) || + // !PyArray_Check(po_vars) || !PyArray_Check(po_v) || + // !PyArray_Check(po_w)) { + // PyErr_SetString(PyExc_ValueError, "Expected numpy arrays"); + // return NULL; + // } + + // np_amps = (PyArrayObject*)po_amps; + // np_means = (PyArrayObject*)po_means; + // np_vars = (PyArrayObject*)po_vars; + // np_v = (PyArrayObject*)po_v; + // np_w = (PyArrayObject*)po_w; + + // if ((PyArray_TYPE(np_amps) != NPY_DOUBLE) || + // (PyArray_TYPE(np_means ) != NPY_DOUBLE) || + // (PyArray_TYPE(np_vars) != NPY_DOUBLE) || + // (PyArray_TYPE(np_v) != NPY_DOUBLE) || + // (PyArray_TYPE(np_w) != NPY_DOUBLE)) { + // PyErr_SetString(PyExc_ValueError, "Expected numpy double arrays"); + // return NULL; + // } + + // if (PyArray_NDIM(np_amps) != 1) { + // PyErr_SetString(PyExc_ValueError, "Expected 'amps' to be 1-d"); + // return NULL; + // } + // K = PyArray_DIM(np_amps, 0); + // if (PyArray_NDIM(np_means) != 2) { + // PyErr_SetString(PyExc_ValueError, "Expected 'means' to be 2-d"); + // return NULL; + // } + // if ((PyArray_DIM(np_means, 0) != K) || + // (PyArray_DIM(np_means, 1) != D)) { + // PyErr_SetString(PyExc_ValueError, "Expected 'means' to be K x D"); + // return NULL; + // } + // if (PyArray_NDIM(np_vars) != 3) { + // PyErr_SetString(PyExc_ValueError, "Expected 'vars' to be 3-d"); + // return NULL; + // } + // if ((PyArray_DIM(np_vars, 0) != K) || + // (PyArray_DIM(np_vars, 1) != D) || + // (PyArray_DIM(np_vars, 2) != D)) { + // PyErr_SetString(PyExc_ValueError, "Expected 'vars' to be K x D x D"); + // return NULL; + // } + + // if (PyArray_NDIM(np_v) != 1) { + // PyErr_SetString(PyExc_ValueError, "Expected 'v' to be 1-d"); + // return NULL; + // } + // if (PyArray_NDIM(np_w) != 1) { + // PyErr_SetString(PyExc_ValueError, "Expected 'w' to be 1-d"); + // return NULL; + // } + + // means = PyArray_DATA(np_means); + + // for (k=0; k +%} + +/**********************************************************************/ + +%fragment("NumPy_Backward_Compatibility", "header") +{ +%#if NPY_API_VERSION < 0x00000007 +%#define NPY_ARRAY_DEFAULT NPY_DEFAULT +%#define NPY_ARRAY_FARRAY NPY_FARRAY +%#define NPY_FORTRANORDER NPY_FORTRAN +%#endif +} + +/**********************************************************************/ + +/* The following code originally appeared in + * enthought/kiva/agg/src/numeric.i written by Eric Jones. It was + * translated from C++ to C by John Hunter. Bill Spotz has modified + * it to fix some minor bugs, upgrade from Numeric to numpy (all + * versions), add some comments and functionality, and convert from + * direct code insertion to SWIG fragments. + */ + +%fragment("NumPy_Macros", "header") +{ +/* Macros to extract array attributes. + */ +%#if NPY_API_VERSION < 0x00000007 +%#define is_array(a) ((a) && PyArray_Check((PyArrayObject*)a)) +%#define array_type(a) (int)(PyArray_TYPE((PyArrayObject*)a)) +%#define array_numdims(a) (((PyArrayObject*)a)->nd) +%#define array_dimensions(a) (((PyArrayObject*)a)->dimensions) +%#define array_size(a,i) (((PyArrayObject*)a)->dimensions[i]) +%#define array_strides(a) (((PyArrayObject*)a)->strides) +%#define array_stride(a,i) (((PyArrayObject*)a)->strides[i]) +%#define array_data(a) (((PyArrayObject*)a)->data) +%#define array_descr(a) (((PyArrayObject*)a)->descr) +%#define array_flags(a) (((PyArrayObject*)a)->flags) +%#define array_enableflags(a,f) (((PyArrayObject*)a)->flags) = f +%#define array_is_fortran(a) (PyArray_ISFORTRAN((PyArrayObject*)a)) +%#else +%#define is_array(a) ((a) && PyArray_Check(a)) +%#define array_type(a) PyArray_TYPE((PyArrayObject*)a) +%#define array_numdims(a) PyArray_NDIM((PyArrayObject*)a) +%#define array_dimensions(a) PyArray_DIMS((PyArrayObject*)a) +%#define array_strides(a) PyArray_STRIDES((PyArrayObject*)a) +%#define array_stride(a,i) PyArray_STRIDE((PyArrayObject*)a,i) +%#define array_size(a,i) PyArray_DIM((PyArrayObject*)a,i) +%#define array_data(a) PyArray_DATA((PyArrayObject*)a) +%#define array_descr(a) PyArray_DESCR((PyArrayObject*)a) +%#define array_flags(a) PyArray_FLAGS((PyArrayObject*)a) +%#define array_enableflags(a,f) PyArray_ENABLEFLAGS((PyArrayObject*)a,f) +%#define array_is_fortran(a) (PyArray_IS_F_CONTIGUOUS((PyArrayObject*)a)) +%#endif +%#define array_is_contiguous(a) (PyArray_ISCONTIGUOUS((PyArrayObject*)a)) +%#define array_is_native(a) (PyArray_ISNOTSWAPPED((PyArrayObject*)a)) +} + +/**********************************************************************/ + +%fragment("NumPy_Utilities", + "header") +{ + /* Given a PyObject, return a string describing its type. + */ + const char* pytype_string(PyObject* py_obj) + { + if (py_obj == NULL ) return "C NULL value"; + if (py_obj == Py_None ) return "Python None" ; + if (PyCallable_Check(py_obj)) return "callable" ; + if (PyString_Check( py_obj)) return "string" ; + if (PyInt_Check( py_obj)) return "int" ; + if (PyFloat_Check( py_obj)) return "float" ; + if (PyDict_Check( py_obj)) return "dict" ; + if (PyList_Check( py_obj)) return "list" ; + if (PyTuple_Check( py_obj)) return "tuple" ; +%#if PY_MAJOR_VERSION < 3 + if (PyFile_Check( py_obj)) return "file" ; + if (PyModule_Check( py_obj)) return "module" ; + if (PyInstance_Check(py_obj)) return "instance" ; +%#endif + + return "unknown type"; + } + + /* Given a NumPy typecode, return a string describing the type. + */ + const char* typecode_string(int typecode) + { + static const char* type_names[25] = {"bool", + "byte", + "unsigned byte", + "short", + "unsigned short", + "int", + "unsigned int", + "long", + "unsigned long", + "long long", + "unsigned long long", + "float", + "double", + "long double", + "complex float", + "complex double", + "complex long double", + "object", + "string", + "unicode", + "void", + "ntypes", + "notype", + "char", + "unknown"}; + return typecode < 24 ? type_names[typecode] : type_names[24]; + } + + /* Make sure input has correct numpy type. This now just calls + PyArray_EquivTypenums(). + */ + int type_match(int actual_type, + int desired_type) + { + return PyArray_EquivTypenums(actual_type, desired_type); + } + +%#ifdef SWIGPY_USE_CAPSULE + void free_cap(PyObject * cap) + { + void* array = (void*) PyCapsule_GetPointer(cap,SWIGPY_CAPSULE_NAME); + if (array != NULL) free(array); + } +%#endif + + +} + +/**********************************************************************/ + +%fragment("NumPy_Object_to_Array", + "header", + fragment="NumPy_Backward_Compatibility", + fragment="NumPy_Macros", + fragment="NumPy_Utilities") +{ + /* Given a PyObject pointer, cast it to a PyArrayObject pointer if + * legal. If not, set the python error string appropriately and + * return NULL. + */ + PyArrayObject* obj_to_array_no_conversion(PyObject* input, + int typecode) + { + PyArrayObject* ary = NULL; + if (is_array(input) && (typecode == NPY_NOTYPE || + PyArray_EquivTypenums(array_type(input), typecode))) + { + ary = (PyArrayObject*) input; + } + else if is_array(input) + { + const char* desired_type = typecode_string(typecode); + const char* actual_type = typecode_string(array_type(input)); + PyErr_Format(PyExc_TypeError, + "Array of type '%s' required. Array of type '%s' given", + desired_type, actual_type); + ary = NULL; + } + else + { + const char* desired_type = typecode_string(typecode); + const char* actual_type = pytype_string(input); + PyErr_Format(PyExc_TypeError, + "Array of type '%s' required. A '%s' was given", + desired_type, + actual_type); + ary = NULL; + } + return ary; + } + + /* Convert the given PyObject to a NumPy array with the given + * typecode. On success, return a valid PyArrayObject* with the + * correct type. On failure, the python error string will be set and + * the routine returns NULL. + */ + PyArrayObject* obj_to_array_allow_conversion(PyObject* input, + int typecode, + int* is_new_object) + { + PyArrayObject* ary = NULL; + PyObject* py_obj; + if (is_array(input) && (typecode == NPY_NOTYPE || + PyArray_EquivTypenums(array_type(input),typecode))) + { + ary = (PyArrayObject*) input; + *is_new_object = 0; + } + else + { + py_obj = PyArray_FROMANY(input, typecode, 0, 0, NPY_ARRAY_DEFAULT); + /* If NULL, PyArray_FromObject will have set python error value.*/ + ary = (PyArrayObject*) py_obj; + *is_new_object = 1; + } + return ary; + } + + /* Given a PyArrayObject, check to see if it is contiguous. If so, + * return the input pointer and flag it as not a new object. If it is + * not contiguous, create a new PyArrayObject using the original data, + * flag it as a new object and return the pointer. + */ + PyArrayObject* make_contiguous(PyArrayObject* ary, + int* is_new_object, + int min_dims, + int max_dims) + { + PyArrayObject* result; + if (array_is_contiguous(ary)) + { + result = ary; + *is_new_object = 0; + } + else + { + result = (PyArrayObject*) PyArray_ContiguousFromObject((PyObject*)ary, + array_type(ary), + min_dims, + max_dims); + *is_new_object = 1; + } + return result; + } + + /* Given a PyArrayObject, check to see if it is Fortran-contiguous. + * If so, return the input pointer, but do not flag it as not a new + * object. If it is not Fortran-contiguous, create a new + * PyArrayObject using the original data, flag it as a new object + * and return the pointer. + */ + PyArrayObject* make_fortran(PyArrayObject* ary, + int* is_new_object) + { + PyArrayObject* result; + if (array_is_fortran(ary)) + { + result = ary; + *is_new_object = 0; + } + else + { + Py_INCREF(array_descr(ary)); + result = (PyArrayObject*) PyArray_FromArray(ary, + array_descr(ary), +%#if NPY_API_VERSION < 0x00000007 + NPY_FORTRANORDER); +%#else + NPY_ARRAY_F_CONTIGUOUS); +%#endif + *is_new_object = 1; + } + return result; + } + + /* Convert a given PyObject to a contiguous PyArrayObject of the + * specified type. If the input object is not a contiguous + * PyArrayObject, a new one will be created and the new object flag + * will be set. + */ + PyArrayObject* obj_to_array_contiguous_allow_conversion(PyObject* input, + int typecode, + int* is_new_object) + { + int is_new1 = 0; + int is_new2 = 0; + PyArrayObject* ary2; + PyArrayObject* ary1 = obj_to_array_allow_conversion(input, + typecode, + &is_new1); + if (ary1) + { + ary2 = make_contiguous(ary1, &is_new2, 0, 0); + if ( is_new1 && is_new2) + { + Py_DECREF(ary1); + } + ary1 = ary2; + } + *is_new_object = is_new1 || is_new2; + return ary1; + } + + /* Convert a given PyObject to a Fortran-ordered PyArrayObject of the + * specified type. If the input object is not a Fortran-ordered + * PyArrayObject, a new one will be created and the new object flag + * will be set. + */ + PyArrayObject* obj_to_array_fortran_allow_conversion(PyObject* input, + int typecode, + int* is_new_object) + { + int is_new1 = 0; + int is_new2 = 0; + PyArrayObject* ary2; + PyArrayObject* ary1 = obj_to_array_allow_conversion(input, + typecode, + &is_new1); + if (ary1) + { + ary2 = make_fortran(ary1, &is_new2); + if (is_new1 && is_new2) + { + Py_DECREF(ary1); + } + ary1 = ary2; + } + *is_new_object = is_new1 || is_new2; + return ary1; + } +} /* end fragment */ + +/**********************************************************************/ + +%fragment("NumPy_Array_Requirements", + "header", + fragment="NumPy_Backward_Compatibility", + fragment="NumPy_Macros") +{ + /* Test whether a python object is contiguous. If array is + * contiguous, return 1. Otherwise, set the python error string and + * return 0. + */ + int require_contiguous(PyArrayObject* ary) + { + int contiguous = 1; + if (!array_is_contiguous(ary)) + { + PyErr_SetString(PyExc_TypeError, + "Array must be contiguous. A non-contiguous array was given"); + contiguous = 0; + } + return contiguous; + } + + /* Test whether a python object is (C_ or F_) contiguous. If array is + * contiguous, return 1. Otherwise, set the python error string and + * return 0. + */ + int require_c_or_f_contiguous(PyArrayObject* ary) + { + int contiguous = 1; + if (!(array_is_contiguous(ary) || array_is_fortran(ary))) + { + PyErr_SetString(PyExc_TypeError, + "Array must be contiguous (C_ or F_). A non-contiguous array was given"); + contiguous = 0; + } + return contiguous; + } + + /* Require that a numpy array is not byte-swapped. If the array is + * not byte-swapped, return 1. Otherwise, set the python error string + * and return 0. + */ + int require_native(PyArrayObject* ary) + { + int native = 1; + if (!array_is_native(ary)) + { + PyErr_SetString(PyExc_TypeError, + "Array must have native byteorder. " + "A byte-swapped array was given"); + native = 0; + } + return native; + } + + /* Require the given PyArrayObject to have a specified number of + * dimensions. If the array has the specified number of dimensions, + * return 1. Otherwise, set the python error string and return 0. + */ + int require_dimensions(PyArrayObject* ary, + int exact_dimensions) + { + int success = 1; + if (array_numdims(ary) != exact_dimensions) + { + PyErr_Format(PyExc_TypeError, + "Array must have %d dimensions. Given array has %d dimensions", + exact_dimensions, + array_numdims(ary)); + success = 0; + } + return success; + } + + /* Require the given PyArrayObject to have one of a list of specified + * number of dimensions. If the array has one of the specified number + * of dimensions, return 1. Otherwise, set the python error string + * and return 0. + */ + int require_dimensions_n(PyArrayObject* ary, + int* exact_dimensions, + int n) + { + int success = 0; + int i; + char dims_str[255] = ""; + char s[255]; + for (i = 0; i < n && !success; i++) + { + if (array_numdims(ary) == exact_dimensions[i]) + { + success = 1; + } + } + if (!success) + { + for (i = 0; i < n-1; i++) + { + sprintf(s, "%d, ", exact_dimensions[i]); + strcat(dims_str,s); + } + sprintf(s, " or %d", exact_dimensions[n-1]); + strcat(dims_str,s); + PyErr_Format(PyExc_TypeError, + "Array must have %s dimensions. Given array has %d dimensions", + dims_str, + array_numdims(ary)); + } + return success; + } + + /* Require the given PyArrayObject to have a specified shape. If the + * array has the specified shape, return 1. Otherwise, set the python + * error string and return 0. + */ + int require_size(PyArrayObject* ary, + npy_intp* size, + int n) + { + int i; + int success = 1; + int len; + char desired_dims[255] = "["; + char s[255]; + char actual_dims[255] = "["; + for(i=0; i < n;i++) + { + if (size[i] != -1 && size[i] != array_size(ary,i)) + { + success = 0; + } + } + if (!success) + { + for (i = 0; i < n; i++) + { + if (size[i] == -1) + { + sprintf(s, "*,"); + } + else + { + sprintf(s, "%ld,", (long int)size[i]); + } + strcat(desired_dims,s); + } + len = strlen(desired_dims); + desired_dims[len-1] = ']'; + for (i = 0; i < n; i++) + { + sprintf(s, "%ld,", (long int)array_size(ary,i)); + strcat(actual_dims,s); + } + len = strlen(actual_dims); + actual_dims[len-1] = ']'; + PyErr_Format(PyExc_TypeError, + "Array must have shape of %s. Given array has shape of %s", + desired_dims, + actual_dims); + } + return success; + } + + /* Require the given PyArrayObject to to be Fortran ordered. If the + * the PyArrayObject is already Fortran ordered, do nothing. Else, + * set the Fortran ordering flag and recompute the strides. + */ + int require_fortran(PyArrayObject* ary) + { + int success = 1; + int nd = array_numdims(ary); + int i; + npy_intp * strides = array_strides(ary); + if (array_is_fortran(ary)) return success; + /* Set the Fortran ordered flag */ + array_enableflags(ary,NPY_ARRAY_FARRAY); + /* Recompute the strides */ + strides[0] = strides[nd-1]; + for (i=1; i < nd; ++i) + strides[i] = strides[i-1] * array_size(ary,i-1); + return success; + } +} + +/* Combine all NumPy fragments into one for convenience */ +%fragment("NumPy_Fragments", + "header", + fragment="NumPy_Backward_Compatibility", + fragment="NumPy_Macros", + fragment="NumPy_Utilities", + fragment="NumPy_Object_to_Array", + fragment="NumPy_Array_Requirements") +{ +} + +/* End John Hunter translation (with modifications by Bill Spotz) + */ + +/* %numpy_typemaps() macro + * + * This macro defines a family of 75 typemaps that allow C arguments + * of the form + * + * 1. (DATA_TYPE IN_ARRAY1[ANY]) + * 2. (DATA_TYPE* IN_ARRAY1, DIM_TYPE DIM1) + * 3. (DIM_TYPE DIM1, DATA_TYPE* IN_ARRAY1) + * + * 4. (DATA_TYPE IN_ARRAY2[ANY][ANY]) + * 5. (DATA_TYPE* IN_ARRAY2, DIM_TYPE DIM1, DIM_TYPE DIM2) + * 6. (DIM_TYPE DIM1, DIM_TYPE DIM2, DATA_TYPE* IN_ARRAY2) + * 7. (DATA_TYPE* IN_FARRAY2, DIM_TYPE DIM1, DIM_TYPE DIM2) + * 8. (DIM_TYPE DIM1, DIM_TYPE DIM2, DATA_TYPE* IN_FARRAY2) + * + * 9. (DATA_TYPE IN_ARRAY3[ANY][ANY][ANY]) + * 10. (DATA_TYPE* IN_ARRAY3, DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3) + * 11. (DATA_TYPE** IN_ARRAY3, DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3) + * 12. (DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, DATA_TYPE* IN_ARRAY3) + * 13. (DATA_TYPE* IN_FARRAY3, DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3) + * 14. (DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, DATA_TYPE* IN_FARRAY3) + * + * 15. (DATA_TYPE IN_ARRAY4[ANY][ANY][ANY][ANY]) + * 16. (DATA_TYPE* IN_ARRAY4, DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, DIM_TYPE DIM4) + * 17. (DATA_TYPE** IN_ARRAY4, DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, DIM_TYPE DIM4) + * 18. (DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, , DIM_TYPE DIM4, DATA_TYPE* IN_ARRAY4) + * 19. (DATA_TYPE* IN_FARRAY4, DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, DIM_TYPE DIM4) + * 20. (DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, DIM_TYPE DIM4, DATA_TYPE* IN_FARRAY4) + * + * 21. (DATA_TYPE INPLACE_ARRAY1[ANY]) + * 22. (DATA_TYPE* INPLACE_ARRAY1, DIM_TYPE DIM1) + * 23. (DIM_TYPE DIM1, DATA_TYPE* INPLACE_ARRAY1) + * + * 24. (DATA_TYPE INPLACE_ARRAY2[ANY][ANY]) + * 25. (DATA_TYPE* INPLACE_ARRAY2, DIM_TYPE DIM1, DIM_TYPE DIM2) + * 26. (DIM_TYPE DIM1, DIM_TYPE DIM2, DATA_TYPE* INPLACE_ARRAY2) + * 27. (DATA_TYPE* INPLACE_FARRAY2, DIM_TYPE DIM1, DIM_TYPE DIM2) + * 28. (DIM_TYPE DIM1, DIM_TYPE DIM2, DATA_TYPE* INPLACE_FARRAY2) + * + * 29. (DATA_TYPE INPLACE_ARRAY3[ANY][ANY][ANY]) + * 30. (DATA_TYPE* INPLACE_ARRAY3, DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3) + * 31. (DATA_TYPE** INPLACE_ARRAY3, DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3) + * 32. (DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, DATA_TYPE* INPLACE_ARRAY3) + * 33. (DATA_TYPE* INPLACE_FARRAY3, DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3) + * 34. (DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, DATA_TYPE* INPLACE_FARRAY3) + * + * 35. (DATA_TYPE INPLACE_ARRAY4[ANY][ANY][ANY][ANY]) + * 36. (DATA_TYPE* INPLACE_ARRAY4, DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, DIM_TYPE DIM4) + * 37. (DATA_TYPE** INPLACE_ARRAY4, DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, DIM_TYPE DIM4) + * 38. (DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, DIM_TYPE DIM4, DATA_TYPE* INPLACE_ARRAY4) + * 39. (DATA_TYPE* INPLACE_FARRAY4, DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, DIM_TYPE DIM4) + * 40. (DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, DIM_TYPE DIM4, DATA_TYPE* INPLACE_FARRAY4) + * + * 41. (DATA_TYPE ARGOUT_ARRAY1[ANY]) + * 42. (DATA_TYPE* ARGOUT_ARRAY1, DIM_TYPE DIM1) + * 43. (DIM_TYPE DIM1, DATA_TYPE* ARGOUT_ARRAY1) + * + * 44. (DATA_TYPE ARGOUT_ARRAY2[ANY][ANY]) + * + * 45. (DATA_TYPE ARGOUT_ARRAY3[ANY][ANY][ANY]) + * + * 46. (DATA_TYPE ARGOUT_ARRAY4[ANY][ANY][ANY][ANY]) + * + * 47. (DATA_TYPE** ARGOUTVIEW_ARRAY1, DIM_TYPE* DIM1) + * 48. (DIM_TYPE* DIM1, DATA_TYPE** ARGOUTVIEW_ARRAY1) + * + * 49. (DATA_TYPE** ARGOUTVIEW_ARRAY2, DIM_TYPE* DIM1, DIM_TYPE* DIM2) + * 50. (DIM_TYPE* DIM1, DIM_TYPE* DIM2, DATA_TYPE** ARGOUTVIEW_ARRAY2) + * 51. (DATA_TYPE** ARGOUTVIEW_FARRAY2, DIM_TYPE* DIM1, DIM_TYPE* DIM2) + * 52. (DIM_TYPE* DIM1, DIM_TYPE* DIM2, DATA_TYPE** ARGOUTVIEW_FARRAY2) + * + * 53. (DATA_TYPE** ARGOUTVIEW_ARRAY3, DIM_TYPE* DIM1, DIM_TYPE* DIM2, DIM_TYPE* DIM3) + * 54. (DIM_TYPE* DIM1, DIM_TYPE* DIM2, DIM_TYPE* DIM3, DATA_TYPE** ARGOUTVIEW_ARRAY3) + * 55. (DATA_TYPE** ARGOUTVIEW_FARRAY3, DIM_TYPE* DIM1, DIM_TYPE* DIM2, DIM_TYPE* DIM3) + * 56. (DIM_TYPE* DIM1, DIM_TYPE* DIM2, DIM_TYPE* DIM3, DATA_TYPE** ARGOUTVIEW_FARRAY3) + * + * 57. (DATA_TYPE** ARGOUTVIEW_ARRAY4, DIM_TYPE* DIM1, DIM_TYPE* DIM2, DIM_TYPE* DIM3, DIM_TYPE* DIM4) + * 58. (DIM_TYPE* DIM1, DIM_TYPE* DIM2, DIM_TYPE* DIM3, DIM_TYPE* DIM4, DATA_TYPE** ARGOUTVIEW_ARRAY4) + * 59. (DATA_TYPE** ARGOUTVIEW_FARRAY4, DIM_TYPE* DIM1, DIM_TYPE* DIM2, DIM_TYPE* DIM3, DIM_TYPE* DIM4) + * 60. (DIM_TYPE* DIM1, DIM_TYPE* DIM2, DIM_TYPE* DIM3, DIM_TYPE* DIM4, DATA_TYPE** ARGOUTVIEW_FARRAY4) + * + * 61. (DATA_TYPE** ARGOUTVIEWM_ARRAY1, DIM_TYPE* DIM1) + * 62. (DIM_TYPE* DIM1, DATA_TYPE** ARGOUTVIEWM_ARRAY1) + * + * 63. (DATA_TYPE** ARGOUTVIEWM_ARRAY2, DIM_TYPE* DIM1, DIM_TYPE* DIM2) + * 64. (DIM_TYPE* DIM1, DIM_TYPE* DIM2, DATA_TYPE** ARGOUTVIEWM_ARRAY2) + * 65. (DATA_TYPE** ARGOUTVIEWM_FARRAY2, DIM_TYPE* DIM1, DIM_TYPE* DIM2) + * 66. (DIM_TYPE* DIM1, DIM_TYPE* DIM2, DATA_TYPE** ARGOUTVIEWM_FARRAY2) + * + * 67. (DATA_TYPE** ARGOUTVIEWM_ARRAY3, DIM_TYPE* DIM1, DIM_TYPE* DIM2, DIM_TYPE* DIM3) + * 68. (DIM_TYPE* DIM1, DIM_TYPE* DIM2, DIM_TYPE* DIM3, DATA_TYPE** ARGOUTVIEWM_ARRAY3) + * 69. (DATA_TYPE** ARGOUTVIEWM_FARRAY3, DIM_TYPE* DIM1, DIM_TYPE* DIM2, DIM_TYPE* DIM3) + * 70. (DIM_TYPE* DIM1, DIM_TYPE* DIM2, DIM_TYPE* DIM3, DATA_TYPE** ARGOUTVIEWM_FARRAY3) + * + * 71. (DATA_TYPE** ARGOUTVIEWM_ARRAY4, DIM_TYPE* DIM1, DIM_TYPE* DIM2, DIM_TYPE* DIM3, DIM_TYPE* DIM4) + * 72. (DIM_TYPE* DIM1, DIM_TYPE* DIM2, DIM_TYPE* DIM3, DIM_TYPE* DIM4, DATA_TYPE** ARGOUTVIEWM_ARRAY4) + * 73. (DATA_TYPE** ARGOUTVIEWM_FARRAY4, DIM_TYPE* DIM1, DIM_TYPE* DIM2, DIM_TYPE* DIM3, DIM_TYPE* DIM4) + * 74. (DIM_TYPE* DIM1, DIM_TYPE* DIM2, DIM_TYPE* DIM3, DIM_TYPE* DIM4, DATA_TYPE** ARGOUTVIEWM_FARRAY4) + * + * 75. (DATA_TYPE* INPLACE_ARRAY_FLAT, DIM_TYPE DIM_FLAT) + * + * where "DATA_TYPE" is any type supported by the NumPy module, and + * "DIM_TYPE" is any int-like type suitable for specifying dimensions. + * The difference between "ARRAY" typemaps and "FARRAY" typemaps is + * that the "FARRAY" typemaps expect Fortran ordering of + * multidimensional arrays. In python, the dimensions will not need + * to be specified (except for the "DATA_TYPE* ARGOUT_ARRAY1" + * typemaps). The IN_ARRAYs can be a numpy array or any sequence that + * can be converted to a numpy array of the specified type. The + * INPLACE_ARRAYs must be numpy arrays of the appropriate type. The + * ARGOUT_ARRAYs will be returned as new numpy arrays of the + * appropriate type. + * + * These typemaps can be applied to existing functions using the + * %apply directive. For example: + * + * %apply (double* IN_ARRAY1, int DIM1) {(double* series, int length)}; + * double prod(double* series, int length); + * + * %apply (int DIM1, int DIM2, double* INPLACE_ARRAY2) + * {(int rows, int cols, double* matrix )}; + * void floor(int rows, int cols, double* matrix, double f); + * + * %apply (double IN_ARRAY3[ANY][ANY][ANY]) + * {(double tensor[2][2][2] )}; + * %apply (double ARGOUT_ARRAY3[ANY][ANY][ANY]) + * {(double low[2][2][2] )}; + * %apply (double ARGOUT_ARRAY3[ANY][ANY][ANY]) + * {(double upp[2][2][2] )}; + * void luSplit(double tensor[2][2][2], + * double low[2][2][2], + * double upp[2][2][2] ); + * + * or directly with + * + * double prod(double* IN_ARRAY1, int DIM1); + * + * void floor(int DIM1, int DIM2, double* INPLACE_ARRAY2, double f); + * + * void luSplit(double IN_ARRAY3[ANY][ANY][ANY], + * double ARGOUT_ARRAY3[ANY][ANY][ANY], + * double ARGOUT_ARRAY3[ANY][ANY][ANY]); + */ + +%define %numpy_typemaps(DATA_TYPE, DATA_TYPECODE, DIM_TYPE) + +/************************/ +/* Input Array Typemaps */ +/************************/ + +/* Typemap suite for (DATA_TYPE IN_ARRAY1[ANY]) + */ +%typecheck(SWIG_TYPECHECK_DOUBLE_ARRAY, + fragment="NumPy_Macros") + (DATA_TYPE IN_ARRAY1[ANY]) +{ + $1 = is_array($input) || PySequence_Check($input); +} +%typemap(in, + fragment="NumPy_Fragments") + (DATA_TYPE IN_ARRAY1[ANY]) + (PyArrayObject* array=NULL, int is_new_object=0) +{ + npy_intp size[1] = { $1_dim0 }; + array = obj_to_array_contiguous_allow_conversion($input, + DATA_TYPECODE, + &is_new_object); + if (!array || !require_dimensions(array, 1) || + !require_size(array, size, 1)) SWIG_fail; + $1 = ($1_ltype) array_data(array); +} +%typemap(freearg) + (DATA_TYPE IN_ARRAY1[ANY]) +{ + if (is_new_object$argnum && array$argnum) + { Py_DECREF(array$argnum); } +} + +/* Typemap suite for (DATA_TYPE* IN_ARRAY1, DIM_TYPE DIM1) + */ +%typecheck(SWIG_TYPECHECK_DOUBLE_ARRAY, + fragment="NumPy_Macros") + (DATA_TYPE* IN_ARRAY1, DIM_TYPE DIM1) +{ + $1 = is_array($input) || PySequence_Check($input); +} +%typemap(in, + fragment="NumPy_Fragments") + (DATA_TYPE* IN_ARRAY1, DIM_TYPE DIM1) + (PyArrayObject* array=NULL, int is_new_object=0) +{ + npy_intp size[1] = { -1 }; + array = obj_to_array_contiguous_allow_conversion($input, + DATA_TYPECODE, + &is_new_object); + if (!array || !require_dimensions(array, 1) || + !require_size(array, size, 1)) SWIG_fail; + $1 = (DATA_TYPE*) array_data(array); + $2 = (DIM_TYPE) array_size(array,0); +} +%typemap(freearg) + (DATA_TYPE* IN_ARRAY1, DIM_TYPE DIM1) +{ + if (is_new_object$argnum && array$argnum) + { Py_DECREF(array$argnum); } +} + +/* Typemap suite for (DIM_TYPE DIM1, DATA_TYPE* IN_ARRAY1) + */ +%typecheck(SWIG_TYPECHECK_DOUBLE_ARRAY, + fragment="NumPy_Macros") + (DIM_TYPE DIM1, DATA_TYPE* IN_ARRAY1) +{ + $1 = is_array($input) || PySequence_Check($input); +} +%typemap(in, + fragment="NumPy_Fragments") + (DIM_TYPE DIM1, DATA_TYPE* IN_ARRAY1) + (PyArrayObject* array=NULL, int is_new_object=0) +{ + npy_intp size[1] = {-1}; + array = obj_to_array_contiguous_allow_conversion($input, + DATA_TYPECODE, + &is_new_object); + if (!array || !require_dimensions(array, 1) || + !require_size(array, size, 1)) SWIG_fail; + $1 = (DIM_TYPE) array_size(array,0); + $2 = (DATA_TYPE*) array_data(array); +} +%typemap(freearg) + (DIM_TYPE DIM1, DATA_TYPE* IN_ARRAY1) +{ + if (is_new_object$argnum && array$argnum) + { Py_DECREF(array$argnum); } +} + +/* Typemap suite for (DATA_TYPE IN_ARRAY2[ANY][ANY]) + */ +%typecheck(SWIG_TYPECHECK_DOUBLE_ARRAY, + fragment="NumPy_Macros") + (DATA_TYPE IN_ARRAY2[ANY][ANY]) +{ + $1 = is_array($input) || PySequence_Check($input); +} +%typemap(in, + fragment="NumPy_Fragments") + (DATA_TYPE IN_ARRAY2[ANY][ANY]) + (PyArrayObject* array=NULL, int is_new_object=0) +{ + npy_intp size[2] = { $1_dim0, $1_dim1 }; + array = obj_to_array_contiguous_allow_conversion($input, + DATA_TYPECODE, + &is_new_object); + if (!array || !require_dimensions(array, 2) || + !require_size(array, size, 2)) SWIG_fail; + $1 = ($1_ltype) array_data(array); +} +%typemap(freearg) + (DATA_TYPE IN_ARRAY2[ANY][ANY]) +{ + if (is_new_object$argnum && array$argnum) + { Py_DECREF(array$argnum); } +} + +/* Typemap suite for (DATA_TYPE* IN_ARRAY2, DIM_TYPE DIM1, DIM_TYPE DIM2) + */ +%typecheck(SWIG_TYPECHECK_DOUBLE_ARRAY, + fragment="NumPy_Macros") + (DATA_TYPE* IN_ARRAY2, DIM_TYPE DIM1, DIM_TYPE DIM2) +{ + $1 = is_array($input) || PySequence_Check($input); +} +%typemap(in, + fragment="NumPy_Fragments") + (DATA_TYPE* IN_ARRAY2, DIM_TYPE DIM1, DIM_TYPE DIM2) + (PyArrayObject* array=NULL, int is_new_object=0) +{ + npy_intp size[2] = { -1, -1 }; + array = obj_to_array_contiguous_allow_conversion($input, DATA_TYPECODE, + &is_new_object); + if (!array || !require_dimensions(array, 2) || + !require_size(array, size, 2)) SWIG_fail; + $1 = (DATA_TYPE*) array_data(array); + $2 = (DIM_TYPE) array_size(array,0); + $3 = (DIM_TYPE) array_size(array,1); +} +%typemap(freearg) + (DATA_TYPE* IN_ARRAY2, DIM_TYPE DIM1, DIM_TYPE DIM2) +{ + if (is_new_object$argnum && array$argnum) + { Py_DECREF(array$argnum); } +} + +/* Typemap suite for (DIM_TYPE DIM1, DIM_TYPE DIM2, DATA_TYPE* IN_ARRAY2) + */ +%typecheck(SWIG_TYPECHECK_DOUBLE_ARRAY, + fragment="NumPy_Macros") + (DIM_TYPE DIM1, DIM_TYPE DIM2, DATA_TYPE* IN_ARRAY2) +{ + $1 = is_array($input) || PySequence_Check($input); +} +%typemap(in, + fragment="NumPy_Fragments") + (DIM_TYPE DIM1, DIM_TYPE DIM2, DATA_TYPE* IN_ARRAY2) + (PyArrayObject* array=NULL, int is_new_object=0) +{ + npy_intp size[2] = { -1, -1 }; + array = obj_to_array_contiguous_allow_conversion($input, + DATA_TYPECODE, + &is_new_object); + if (!array || !require_dimensions(array, 2) || + !require_size(array, size, 2)) SWIG_fail; + $1 = (DIM_TYPE) array_size(array,0); + $2 = (DIM_TYPE) array_size(array,1); + $3 = (DATA_TYPE*) array_data(array); +} +%typemap(freearg) + (DIM_TYPE DIM1, DIM_TYPE DIM2, DATA_TYPE* IN_ARRAY2) +{ + if (is_new_object$argnum && array$argnum) + { Py_DECREF(array$argnum); } +} + +/* Typemap suite for (DATA_TYPE* IN_FARRAY2, DIM_TYPE DIM1, DIM_TYPE DIM2) + */ +%typecheck(SWIG_TYPECHECK_DOUBLE_ARRAY, + fragment="NumPy_Macros") + (DATA_TYPE* IN_FARRAY2, DIM_TYPE DIM1, DIM_TYPE DIM2) +{ + $1 = is_array($input) || PySequence_Check($input); +} +%typemap(in, + fragment="NumPy_Fragments") + (DATA_TYPE* IN_FARRAY2, DIM_TYPE DIM1, DIM_TYPE DIM2) + (PyArrayObject* array=NULL, int is_new_object=0) +{ + npy_intp size[2] = { -1, -1 }; + array = obj_to_array_fortran_allow_conversion($input, + DATA_TYPECODE, + &is_new_object); + if (!array || !require_dimensions(array, 2) || + !require_size(array, size, 2) || !require_fortran(array)) SWIG_fail; + $1 = (DATA_TYPE*) array_data(array); + $2 = (DIM_TYPE) array_size(array,0); + $3 = (DIM_TYPE) array_size(array,1); +} +%typemap(freearg) + (DATA_TYPE* IN_FARRAY2, DIM_TYPE DIM1, DIM_TYPE DIM2) +{ + if (is_new_object$argnum && array$argnum) + { Py_DECREF(array$argnum); } +} + +/* Typemap suite for (DIM_TYPE DIM1, DIM_TYPE DIM2, DATA_TYPE* IN_FARRAY2) + */ +%typecheck(SWIG_TYPECHECK_DOUBLE_ARRAY, + fragment="NumPy_Macros") + (DIM_TYPE DIM1, DIM_TYPE DIM2, DATA_TYPE* IN_FARRAY2) +{ + $1 = is_array($input) || PySequence_Check($input); +} +%typemap(in, + fragment="NumPy_Fragments") + (DIM_TYPE DIM1, DIM_TYPE DIM2, DATA_TYPE* IN_FARRAY2) + (PyArrayObject* array=NULL, int is_new_object=0) +{ + npy_intp size[2] = { -1, -1 }; + array = obj_to_array_fortran_allow_conversion($input, + DATA_TYPECODE, + &is_new_object); + if (!array || !require_dimensions(array, 2) || + !require_size(array, size, 2) || !require_fortran(array)) SWIG_fail; + $1 = (DIM_TYPE) array_size(array,0); + $2 = (DIM_TYPE) array_size(array,1); + $3 = (DATA_TYPE*) array_data(array); +} +%typemap(freearg) + (DIM_TYPE DIM1, DIM_TYPE DIM2, DATA_TYPE* IN_FARRAY2) +{ + if (is_new_object$argnum && array$argnum) + { Py_DECREF(array$argnum); } +} + +/* Typemap suite for (DATA_TYPE IN_ARRAY3[ANY][ANY][ANY]) + */ +%typecheck(SWIG_TYPECHECK_DOUBLE_ARRAY, + fragment="NumPy_Macros") + (DATA_TYPE IN_ARRAY3[ANY][ANY][ANY]) +{ + $1 = is_array($input) || PySequence_Check($input); +} +%typemap(in, + fragment="NumPy_Fragments") + (DATA_TYPE IN_ARRAY3[ANY][ANY][ANY]) + (PyArrayObject* array=NULL, int is_new_object=0) +{ + npy_intp size[3] = { $1_dim0, $1_dim1, $1_dim2 }; + array = obj_to_array_contiguous_allow_conversion($input, + DATA_TYPECODE, + &is_new_object); + if (!array || !require_dimensions(array, 3) || + !require_size(array, size, 3)) SWIG_fail; + $1 = ($1_ltype) array_data(array); +} +%typemap(freearg) + (DATA_TYPE IN_ARRAY3[ANY][ANY][ANY]) +{ + if (is_new_object$argnum && array$argnum) + { Py_DECREF(array$argnum); } +} + +/* Typemap suite for (DATA_TYPE* IN_ARRAY3, DIM_TYPE DIM1, DIM_TYPE DIM2, + * DIM_TYPE DIM3) + */ +%typecheck(SWIG_TYPECHECK_DOUBLE_ARRAY, + fragment="NumPy_Macros") + (DATA_TYPE* IN_ARRAY3, DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3) +{ + $1 = is_array($input) || PySequence_Check($input); +} +%typemap(in, + fragment="NumPy_Fragments") + (DATA_TYPE* IN_ARRAY3, DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3) + (PyArrayObject* array=NULL, int is_new_object=0) +{ + npy_intp size[3] = { -1, -1, -1 }; + array = obj_to_array_contiguous_allow_conversion($input, DATA_TYPECODE, + &is_new_object); + if (!array || !require_dimensions(array, 3) || + !require_size(array, size, 3)) SWIG_fail; + $1 = (DATA_TYPE*) array_data(array); + $2 = (DIM_TYPE) array_size(array,0); + $3 = (DIM_TYPE) array_size(array,1); + $4 = (DIM_TYPE) array_size(array,2); +} +%typemap(freearg) + (DATA_TYPE* IN_ARRAY3, DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3) +{ + if (is_new_object$argnum && array$argnum) + { Py_DECREF(array$argnum); } +} + +/* Typemap suite for (DATA_TYPE** IN_ARRAY3, DIM_TYPE DIM1, DIM_TYPE DIM2, + * DIM_TYPE DIM3) + */ +%typecheck(SWIG_TYPECHECK_DOUBLE_ARRAY, + fragment="NumPy_Macros") + (DATA_TYPE** IN_ARRAY3, DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3) +{ + /* for now, only concerned with lists */ + $1 = PySequence_Check($input); +} +%typemap(in, + fragment="NumPy_Fragments") + (DATA_TYPE** IN_ARRAY3, DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3) + (DATA_TYPE** array=NULL, PyArrayObject** object_array=NULL, int* is_new_object_array=NULL) +{ + npy_intp size[2] = { -1, -1 }; + PyArrayObject* temp_array; + Py_ssize_t i; + int is_new_object; + + /* length of the list */ + $2 = PyList_Size($input); + + /* the arrays */ + array = (DATA_TYPE **)malloc($2*sizeof(DATA_TYPE *)); + object_array = (PyArrayObject **)calloc($2,sizeof(PyArrayObject *)); + is_new_object_array = (int *)calloc($2,sizeof(int)); + + if (array == NULL || object_array == NULL || is_new_object_array == NULL) + { + SWIG_fail; + } + + for (i=0; i<$2; i++) + { + temp_array = obj_to_array_contiguous_allow_conversion(PySequence_GetItem($input,i), DATA_TYPECODE, &is_new_object); + + /* the new array must be stored so that it can be destroyed in freearg */ + object_array[i] = temp_array; + is_new_object_array[i] = is_new_object; + + if (!temp_array || !require_dimensions(temp_array, 2)) SWIG_fail; + + /* store the size of the first array in the list, then use that for comparison. */ + if (i == 0) + { + size[0] = array_size(temp_array,0); + size[1] = array_size(temp_array,1); + } + + if (!require_size(temp_array, size, 2)) SWIG_fail; + + array[i] = (DATA_TYPE*) array_data(temp_array); + } + + $1 = (DATA_TYPE**) array; + $3 = (DIM_TYPE) size[0]; + $4 = (DIM_TYPE) size[1]; +} +%typemap(freearg) + (DATA_TYPE** IN_ARRAY3, DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3) +{ + Py_ssize_t i; + + if (array$argnum!=NULL) free(array$argnum); + + /*freeing the individual arrays if needed */ + if (object_array$argnum!=NULL) + { + if (is_new_object_array$argnum!=NULL) + { + for (i=0; i<$2; i++) + { + if (object_array$argnum[i] != NULL && is_new_object_array$argnum[i]) + { Py_DECREF(object_array$argnum[i]); } + } + free(is_new_object_array$argnum); + } + free(object_array$argnum); + } +} + +/* Typemap suite for (DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, + * DATA_TYPE* IN_ARRAY3) + */ +%typecheck(SWIG_TYPECHECK_DOUBLE_ARRAY, + fragment="NumPy_Macros") + (DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, DATA_TYPE* IN_ARRAY3) +{ + $1 = is_array($input) || PySequence_Check($input); +} +%typemap(in, + fragment="NumPy_Fragments") + (DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, DATA_TYPE* IN_ARRAY3) + (PyArrayObject* array=NULL, int is_new_object=0) +{ + npy_intp size[3] = { -1, -1, -1 }; + array = obj_to_array_contiguous_allow_conversion($input, DATA_TYPECODE, + &is_new_object); + if (!array || !require_dimensions(array, 3) || + !require_size(array, size, 3)) SWIG_fail; + $1 = (DIM_TYPE) array_size(array,0); + $2 = (DIM_TYPE) array_size(array,1); + $3 = (DIM_TYPE) array_size(array,2); + $4 = (DATA_TYPE*) array_data(array); +} +%typemap(freearg) + (DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, DATA_TYPE* IN_ARRAY3) +{ + if (is_new_object$argnum && array$argnum) + { Py_DECREF(array$argnum); } +} + +/* Typemap suite for (DATA_TYPE* IN_FARRAY3, DIM_TYPE DIM1, DIM_TYPE DIM2, + * DIM_TYPE DIM3) + */ +%typecheck(SWIG_TYPECHECK_DOUBLE_ARRAY, + fragment="NumPy_Macros") + (DATA_TYPE* IN_FARRAY3, DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3) +{ + $1 = is_array($input) || PySequence_Check($input); +} +%typemap(in, + fragment="NumPy_Fragments") + (DATA_TYPE* IN_FARRAY3, DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3) + (PyArrayObject* array=NULL, int is_new_object=0) +{ + npy_intp size[3] = { -1, -1, -1 }; + array = obj_to_array_fortran_allow_conversion($input, DATA_TYPECODE, + &is_new_object); + if (!array || !require_dimensions(array, 3) || + !require_size(array, size, 3) | !require_fortran(array)) SWIG_fail; + $1 = (DATA_TYPE*) array_data(array); + $2 = (DIM_TYPE) array_size(array,0); + $3 = (DIM_TYPE) array_size(array,1); + $4 = (DIM_TYPE) array_size(array,2); +} +%typemap(freearg) + (DATA_TYPE* IN_FARRAY3, DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3) +{ + if (is_new_object$argnum && array$argnum) + { Py_DECREF(array$argnum); } +} + +/* Typemap suite for (DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, + * DATA_TYPE* IN_FARRAY3) + */ +%typecheck(SWIG_TYPECHECK_DOUBLE_ARRAY, + fragment="NumPy_Macros") + (DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, DATA_TYPE* IN_FARRAY3) +{ + $1 = is_array($input) || PySequence_Check($input); +} +%typemap(in, + fragment="NumPy_Fragments") + (DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, DATA_TYPE* IN_FARRAY3) + (PyArrayObject* array=NULL, int is_new_object=0) +{ + npy_intp size[3] = { -1, -1, -1 }; + array = obj_to_array_fortran_allow_conversion($input, + DATA_TYPECODE, + &is_new_object); + if (!array || !require_dimensions(array, 3) || + !require_size(array, size, 3) || !require_fortran(array)) SWIG_fail; + $1 = (DIM_TYPE) array_size(array,0); + $2 = (DIM_TYPE) array_size(array,1); + $3 = (DIM_TYPE) array_size(array,2); + $4 = (DATA_TYPE*) array_data(array); +} +%typemap(freearg) + (DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, DATA_TYPE* IN_FARRAY3) +{ + if (is_new_object$argnum && array$argnum) + { Py_DECREF(array$argnum); } +} + +/* Typemap suite for (DATA_TYPE IN_ARRAY4[ANY][ANY][ANY][ANY]) + */ +%typecheck(SWIG_TYPECHECK_DOUBLE_ARRAY, + fragment="NumPy_Macros") + (DATA_TYPE IN_ARRAY4[ANY][ANY][ANY][ANY]) +{ + $1 = is_array($input) || PySequence_Check($input); +} +%typemap(in, + fragment="NumPy_Fragments") + (DATA_TYPE IN_ARRAY4[ANY][ANY][ANY][ANY]) + (PyArrayObject* array=NULL, int is_new_object=0) +{ + npy_intp size[4] = { $1_dim0, $1_dim1, $1_dim2 , $1_dim3}; + array = obj_to_array_contiguous_allow_conversion($input, DATA_TYPECODE, + &is_new_object); + if (!array || !require_dimensions(array, 4) || + !require_size(array, size, 4)) SWIG_fail; + $1 = ($1_ltype) array_data(array); +} +%typemap(freearg) + (DATA_TYPE IN_ARRAY4[ANY][ANY][ANY][ANY]) +{ + if (is_new_object$argnum && array$argnum) + { Py_DECREF(array$argnum); } +} + +/* Typemap suite for (DATA_TYPE* IN_ARRAY4, DIM_TYPE DIM1, DIM_TYPE DIM2, + * DIM_TYPE DIM3, DIM_TYPE DIM4) + */ +%typecheck(SWIG_TYPECHECK_DOUBLE_ARRAY, + fragment="NumPy_Macros") + (DATA_TYPE* IN_ARRAY4, DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, DIM_TYPE DIM4) +{ + $1 = is_array($input) || PySequence_Check($input); +} +%typemap(in, + fragment="NumPy_Fragments") + (DATA_TYPE* IN_ARRAY4, DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, DIM_TYPE DIM4) + (PyArrayObject* array=NULL, int is_new_object=0) +{ + npy_intp size[4] = { -1, -1, -1, -1 }; + array = obj_to_array_contiguous_allow_conversion($input, DATA_TYPECODE, + &is_new_object); + if (!array || !require_dimensions(array, 4) || + !require_size(array, size, 4)) SWIG_fail; + $1 = (DATA_TYPE*) array_data(array); + $2 = (DIM_TYPE) array_size(array,0); + $3 = (DIM_TYPE) array_size(array,1); + $4 = (DIM_TYPE) array_size(array,2); + $5 = (DIM_TYPE) array_size(array,3); +} +%typemap(freearg) + (DATA_TYPE* IN_ARRAY4, DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, DIM_TYPE DIM4) +{ + if (is_new_object$argnum && array$argnum) + { Py_DECREF(array$argnum); } +} + +/* Typemap suite for (DATA_TYPE** IN_ARRAY4, DIM_TYPE DIM1, DIM_TYPE DIM2, + * DIM_TYPE DIM3, DIM_TYPE DIM4) + */ +%typecheck(SWIG_TYPECHECK_DOUBLE_ARRAY, + fragment="NumPy_Macros") + (DATA_TYPE** IN_ARRAY4, DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, DIM_TYPE DIM4) +{ + /* for now, only concerned with lists */ + $1 = PySequence_Check($input); +} +%typemap(in, + fragment="NumPy_Fragments") + (DATA_TYPE** IN_ARRAY4, DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, DIM_TYPE DIM4) + (DATA_TYPE** array=NULL, PyArrayObject** object_array=NULL, int* is_new_object_array=NULL) +{ + npy_intp size[3] = { -1, -1, -1 }; + PyArrayObject* temp_array; + Py_ssize_t i; + int is_new_object; + + /* length of the list */ + $2 = PyList_Size($input); + + /* the arrays */ + array = (DATA_TYPE **)malloc($2*sizeof(DATA_TYPE *)); + object_array = (PyArrayObject **)calloc($2,sizeof(PyArrayObject *)); + is_new_object_array = (int *)calloc($2,sizeof(int)); + + if (array == NULL || object_array == NULL || is_new_object_array == NULL) + { + SWIG_fail; + } + + for (i=0; i<$2; i++) + { + temp_array = obj_to_array_contiguous_allow_conversion(PySequence_GetItem($input,i), DATA_TYPECODE, &is_new_object); + + /* the new array must be stored so that it can be destroyed in freearg */ + object_array[i] = temp_array; + is_new_object_array[i] = is_new_object; + + if (!temp_array || !require_dimensions(temp_array, 3)) SWIG_fail; + + /* store the size of the first array in the list, then use that for comparison. */ + if (i == 0) + { + size[0] = array_size(temp_array,0); + size[1] = array_size(temp_array,1); + size[2] = array_size(temp_array,2); + } + + if (!require_size(temp_array, size, 3)) SWIG_fail; + + array[i] = (DATA_TYPE*) array_data(temp_array); + } + + $1 = (DATA_TYPE**) array; + $3 = (DIM_TYPE) size[0]; + $4 = (DIM_TYPE) size[1]; + $5 = (DIM_TYPE) size[2]; +} +%typemap(freearg) + (DATA_TYPE** IN_ARRAY4, DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, DIM_TYPE DIM4) +{ + Py_ssize_t i; + + if (array$argnum!=NULL) free(array$argnum); + + /*freeing the individual arrays if needed */ + if (object_array$argnum!=NULL) + { + if (is_new_object_array$argnum!=NULL) + { + for (i=0; i<$2; i++) + { + if (object_array$argnum[i] != NULL && is_new_object_array$argnum[i]) + { Py_DECREF(object_array$argnum[i]); } + } + free(is_new_object_array$argnum); + } + free(object_array$argnum); + } +} + +/* Typemap suite for (DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, DIM_TYPE DIM4, + * DATA_TYPE* IN_ARRAY4) + */ +%typecheck(SWIG_TYPECHECK_DOUBLE_ARRAY, + fragment="NumPy_Macros") + (DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, DIM_TYPE DIM4, DATA_TYPE* IN_ARRAY4) +{ + $1 = is_array($input) || PySequence_Check($input); +} +%typemap(in, + fragment="NumPy_Fragments") + (DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, DIM_TYPE DIM4, DATA_TYPE* IN_ARRAY4) + (PyArrayObject* array=NULL, int is_new_object=0) +{ + npy_intp size[4] = { -1, -1, -1 , -1}; + array = obj_to_array_contiguous_allow_conversion($input, DATA_TYPECODE, + &is_new_object); + if (!array || !require_dimensions(array, 4) || + !require_size(array, size, 4)) SWIG_fail; + $1 = (DIM_TYPE) array_size(array,0); + $2 = (DIM_TYPE) array_size(array,1); + $3 = (DIM_TYPE) array_size(array,2); + $4 = (DIM_TYPE) array_size(array,3); + $5 = (DATA_TYPE*) array_data(array); +} +%typemap(freearg) + (DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, DIM_TYPE DIM4, DATA_TYPE* IN_ARRAY4) +{ + if (is_new_object$argnum && array$argnum) + { Py_DECREF(array$argnum); } +} + +/* Typemap suite for (DATA_TYPE* IN_FARRAY4, DIM_TYPE DIM1, DIM_TYPE DIM2, + * DIM_TYPE DIM3, DIM_TYPE DIM4) + */ +%typecheck(SWIG_TYPECHECK_DOUBLE_ARRAY, + fragment="NumPy_Macros") + (DATA_TYPE* IN_FARRAY4, DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, DIM_TYPE DIM4) +{ + $1 = is_array($input) || PySequence_Check($input); +} +%typemap(in, + fragment="NumPy_Fragments") + (DATA_TYPE* IN_FARRAY4, DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, DIM_TYPE DIM4) + (PyArrayObject* array=NULL, int is_new_object=0) +{ + npy_intp size[4] = { -1, -1, -1, -1 }; + array = obj_to_array_fortran_allow_conversion($input, DATA_TYPECODE, + &is_new_object); + if (!array || !require_dimensions(array, 4) || + !require_size(array, size, 4) | !require_fortran(array)) SWIG_fail; + $1 = (DATA_TYPE*) array_data(array); + $2 = (DIM_TYPE) array_size(array,0); + $3 = (DIM_TYPE) array_size(array,1); + $4 = (DIM_TYPE) array_size(array,2); + $5 = (DIM_TYPE) array_size(array,3); +} +%typemap(freearg) + (DATA_TYPE* IN_FARRAY4, DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, DIM_TYPE DIM4) +{ + if (is_new_object$argnum && array$argnum) + { Py_DECREF(array$argnum); } +} + +/* Typemap suite for (DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, DIM_TYPE DIM4, + * DATA_TYPE* IN_FARRAY4) + */ +%typecheck(SWIG_TYPECHECK_DOUBLE_ARRAY, + fragment="NumPy_Macros") + (DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, DIM_TYPE DIM4, DATA_TYPE* IN_FARRAY4) +{ + $1 = is_array($input) || PySequence_Check($input); +} +%typemap(in, + fragment="NumPy_Fragments") + (DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, DIM_TYPE DIM4, DATA_TYPE* IN_FARRAY4) + (PyArrayObject* array=NULL, int is_new_object=0) +{ + npy_intp size[4] = { -1, -1, -1 , -1 }; + array = obj_to_array_fortran_allow_conversion($input, DATA_TYPECODE, + &is_new_object); + if (!array || !require_dimensions(array, 4) || + !require_size(array, size, 4) || !require_fortran(array)) SWIG_fail; + $1 = (DIM_TYPE) array_size(array,0); + $2 = (DIM_TYPE) array_size(array,1); + $3 = (DIM_TYPE) array_size(array,2); + $4 = (DIM_TYPE) array_size(array,3); + $5 = (DATA_TYPE*) array_data(array); +} +%typemap(freearg) + (DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, DIM_TYPE DIM4, DATA_TYPE* IN_FARRAY4) +{ + if (is_new_object$argnum && array$argnum) + { Py_DECREF(array$argnum); } +} + +/***************************/ +/* In-Place Array Typemaps */ +/***************************/ + +/* Typemap suite for (DATA_TYPE INPLACE_ARRAY1[ANY]) + */ +%typecheck(SWIG_TYPECHECK_DOUBLE_ARRAY, + fragment="NumPy_Macros") + (DATA_TYPE INPLACE_ARRAY1[ANY]) +{ + $1 = is_array($input) && PyArray_EquivTypenums(array_type($input), + DATA_TYPECODE); +} +%typemap(in, + fragment="NumPy_Fragments") + (DATA_TYPE INPLACE_ARRAY1[ANY]) + (PyArrayObject* array=NULL) +{ + npy_intp size[1] = { $1_dim0 }; + array = obj_to_array_no_conversion($input, DATA_TYPECODE); + if (!array || !require_dimensions(array,1) || !require_size(array, size, 1) || + !require_contiguous(array) || !require_native(array)) SWIG_fail; + $1 = ($1_ltype) array_data(array); +} + +/* Typemap suite for (DATA_TYPE* INPLACE_ARRAY1, DIM_TYPE DIM1) + */ +%typecheck(SWIG_TYPECHECK_DOUBLE_ARRAY, + fragment="NumPy_Macros") + (DATA_TYPE* INPLACE_ARRAY1, DIM_TYPE DIM1) +{ + $1 = is_array($input) && PyArray_EquivTypenums(array_type($input), + DATA_TYPECODE); +} +%typemap(in, + fragment="NumPy_Fragments") + (DATA_TYPE* INPLACE_ARRAY1, DIM_TYPE DIM1) + (PyArrayObject* array=NULL, int i=1) +{ + array = obj_to_array_no_conversion($input, DATA_TYPECODE); + if (!array || !require_dimensions(array,1) || !require_contiguous(array) + || !require_native(array)) SWIG_fail; + $1 = (DATA_TYPE*) array_data(array); + $2 = 1; + for (i=0; i < array_numdims(array); ++i) $2 *= array_size(array,i); +} + +/* Typemap suite for (DIM_TYPE DIM1, DATA_TYPE* INPLACE_ARRAY1) + */ +%typecheck(SWIG_TYPECHECK_DOUBLE_ARRAY, + fragment="NumPy_Macros") + (DIM_TYPE DIM1, DATA_TYPE* INPLACE_ARRAY1) +{ + $1 = is_array($input) && PyArray_EquivTypenums(array_type($input), + DATA_TYPECODE); +} +%typemap(in, + fragment="NumPy_Fragments") + (DIM_TYPE DIM1, DATA_TYPE* INPLACE_ARRAY1) + (PyArrayObject* array=NULL, int i=0) +{ + array = obj_to_array_no_conversion($input, DATA_TYPECODE); + if (!array || !require_dimensions(array,1) || !require_contiguous(array) + || !require_native(array)) SWIG_fail; + $1 = 1; + for (i=0; i < array_numdims(array); ++i) $1 *= array_size(array,i); + $2 = (DATA_TYPE*) array_data(array); +} + +/* Typemap suite for (DATA_TYPE INPLACE_ARRAY2[ANY][ANY]) + */ +%typecheck(SWIG_TYPECHECK_DOUBLE_ARRAY, + fragment="NumPy_Macros") + (DATA_TYPE INPLACE_ARRAY2[ANY][ANY]) +{ + $1 = is_array($input) && PyArray_EquivTypenums(array_type($input), + DATA_TYPECODE); +} +%typemap(in, + fragment="NumPy_Fragments") + (DATA_TYPE INPLACE_ARRAY2[ANY][ANY]) + (PyArrayObject* array=NULL) +{ + npy_intp size[2] = { $1_dim0, $1_dim1 }; + array = obj_to_array_no_conversion($input, DATA_TYPECODE); + if (!array || !require_dimensions(array,2) || !require_size(array, size, 2) || + !require_contiguous(array) || !require_native(array)) SWIG_fail; + $1 = ($1_ltype) array_data(array); +} + +/* Typemap suite for (DATA_TYPE* INPLACE_ARRAY2, DIM_TYPE DIM1, DIM_TYPE DIM2) + */ +%typecheck(SWIG_TYPECHECK_DOUBLE_ARRAY, + fragment="NumPy_Macros") + (DATA_TYPE* INPLACE_ARRAY2, DIM_TYPE DIM1, DIM_TYPE DIM2) +{ + $1 = is_array($input) && PyArray_EquivTypenums(array_type($input), + DATA_TYPECODE); +} +%typemap(in, + fragment="NumPy_Fragments") + (DATA_TYPE* INPLACE_ARRAY2, DIM_TYPE DIM1, DIM_TYPE DIM2) + (PyArrayObject* array=NULL) +{ + array = obj_to_array_no_conversion($input, DATA_TYPECODE); + if (!array || !require_dimensions(array,2) || !require_contiguous(array) + || !require_native(array)) SWIG_fail; + $1 = (DATA_TYPE*) array_data(array); + $2 = (DIM_TYPE) array_size(array,0); + $3 = (DIM_TYPE) array_size(array,1); +} + +/* Typemap suite for (DIM_TYPE DIM1, DIM_TYPE DIM2, DATA_TYPE* INPLACE_ARRAY2) + */ +%typecheck(SWIG_TYPECHECK_DOUBLE_ARRAY, + fragment="NumPy_Macros") + (DIM_TYPE DIM1, DIM_TYPE DIM2, DATA_TYPE* INPLACE_ARRAY2) +{ + $1 = is_array($input) && PyArray_EquivTypenums(array_type($input), + DATA_TYPECODE); +} +%typemap(in, + fragment="NumPy_Fragments") + (DIM_TYPE DIM1, DIM_TYPE DIM2, DATA_TYPE* INPLACE_ARRAY2) + (PyArrayObject* array=NULL) +{ + array = obj_to_array_no_conversion($input, DATA_TYPECODE); + if (!array || !require_dimensions(array,2) || !require_contiguous(array) || + !require_native(array)) SWIG_fail; + $1 = (DIM_TYPE) array_size(array,0); + $2 = (DIM_TYPE) array_size(array,1); + $3 = (DATA_TYPE*) array_data(array); +} + +/* Typemap suite for (DATA_TYPE* INPLACE_FARRAY2, DIM_TYPE DIM1, DIM_TYPE DIM2) + */ +%typecheck(SWIG_TYPECHECK_DOUBLE_ARRAY, + fragment="NumPy_Macros") + (DATA_TYPE* INPLACE_FARRAY2, DIM_TYPE DIM1, DIM_TYPE DIM2) +{ + $1 = is_array($input) && PyArray_EquivTypenums(array_type($input), + DATA_TYPECODE); +} +%typemap(in, + fragment="NumPy_Fragments") + (DATA_TYPE* INPLACE_FARRAY2, DIM_TYPE DIM1, DIM_TYPE DIM2) + (PyArrayObject* array=NULL) +{ + array = obj_to_array_no_conversion($input, DATA_TYPECODE); + if (!array || !require_dimensions(array,2) || !require_contiguous(array) + || !require_native(array) || !require_fortran(array)) SWIG_fail; + $1 = (DATA_TYPE*) array_data(array); + $2 = (DIM_TYPE) array_size(array,0); + $3 = (DIM_TYPE) array_size(array,1); +} + +/* Typemap suite for (DIM_TYPE DIM1, DIM_TYPE DIM2, DATA_TYPE* INPLACE_FARRAY2) + */ +%typecheck(SWIG_TYPECHECK_DOUBLE_ARRAY, + fragment="NumPy_Macros") + (DIM_TYPE DIM1, DIM_TYPE DIM2, DATA_TYPE* INPLACE_FARRAY2) +{ + $1 = is_array($input) && PyArray_EquivTypenums(array_type($input), + DATA_TYPECODE); +} +%typemap(in, + fragment="NumPy_Fragments") + (DIM_TYPE DIM1, DIM_TYPE DIM2, DATA_TYPE* INPLACE_FARRAY2) + (PyArrayObject* array=NULL) +{ + array = obj_to_array_no_conversion($input, DATA_TYPECODE); + if (!array || !require_dimensions(array,2) || !require_contiguous(array) || + !require_native(array) || !require_fortran(array)) SWIG_fail; + $1 = (DIM_TYPE) array_size(array,0); + $2 = (DIM_TYPE) array_size(array,1); + $3 = (DATA_TYPE*) array_data(array); +} + +/* Typemap suite for (DATA_TYPE INPLACE_ARRAY3[ANY][ANY][ANY]) + */ +%typecheck(SWIG_TYPECHECK_DOUBLE_ARRAY, + fragment="NumPy_Macros") + (DATA_TYPE INPLACE_ARRAY3[ANY][ANY][ANY]) +{ + $1 = is_array($input) && PyArray_EquivTypenums(array_type($input), + DATA_TYPECODE); +} +%typemap(in, + fragment="NumPy_Fragments") + (DATA_TYPE INPLACE_ARRAY3[ANY][ANY][ANY]) + (PyArrayObject* array=NULL) +{ + npy_intp size[3] = { $1_dim0, $1_dim1, $1_dim2 }; + array = obj_to_array_no_conversion($input, DATA_TYPECODE); + if (!array || !require_dimensions(array,3) || !require_size(array, size, 3) || + !require_contiguous(array) || !require_native(array)) SWIG_fail; + $1 = ($1_ltype) array_data(array); +} + +/* Typemap suite for (DATA_TYPE* INPLACE_ARRAY3, DIM_TYPE DIM1, DIM_TYPE DIM2, + * DIM_TYPE DIM3) + */ +%typecheck(SWIG_TYPECHECK_DOUBLE_ARRAY, + fragment="NumPy_Macros") + (DATA_TYPE* INPLACE_ARRAY3, DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3) +{ + $1 = is_array($input) && PyArray_EquivTypenums(array_type($input), + DATA_TYPECODE); +} +%typemap(in, + fragment="NumPy_Fragments") + (DATA_TYPE* INPLACE_ARRAY3, DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3) + (PyArrayObject* array=NULL) +{ + array = obj_to_array_no_conversion($input, DATA_TYPECODE); + if (!array || !require_dimensions(array,3) || !require_contiguous(array) || + !require_native(array)) SWIG_fail; + $1 = (DATA_TYPE*) array_data(array); + $2 = (DIM_TYPE) array_size(array,0); + $3 = (DIM_TYPE) array_size(array,1); + $4 = (DIM_TYPE) array_size(array,2); +} + +/* Typemap suite for (DATA_TYPE** INPLACE_ARRAY3, DIM_TYPE DIM1, DIM_TYPE DIM2, + * DIM_TYPE DIM3) + */ +%typecheck(SWIG_TYPECHECK_DOUBLE_ARRAY, + fragment="NumPy_Macros") + (DATA_TYPE** INPLACE_ARRAY3, DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3) +{ + $1 = PySequence_Check($input); +} +%typemap(in, + fragment="NumPy_Fragments") + (DATA_TYPE** INPLACE_ARRAY3, DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3) + (DATA_TYPE** array=NULL, PyArrayObject** object_array=NULL) +{ + npy_intp size[2] = { -1, -1 }; + PyArrayObject* temp_array; + Py_ssize_t i; + + /* length of the list */ + $2 = PyList_Size($input); + + /* the arrays */ + array = (DATA_TYPE **)malloc($2*sizeof(DATA_TYPE *)); + object_array = (PyArrayObject **)calloc($2,sizeof(PyArrayObject *)); + + if (array == NULL || object_array == NULL) + { + SWIG_fail; + } + + for (i=0; i<$2; i++) + { + temp_array = obj_to_array_no_conversion(PySequence_GetItem($input,i), DATA_TYPECODE); + + /* the new array must be stored so that it can be destroyed in freearg */ + object_array[i] = temp_array; + + if ( !temp_array || !require_dimensions(temp_array, 2) || + !require_contiguous(temp_array) || + !require_native(temp_array) || + !PyArray_EquivTypenums(array_type(temp_array), DATA_TYPECODE) + ) SWIG_fail; + + /* store the size of the first array in the list, then use that for comparison. */ + if (i == 0) + { + size[0] = array_size(temp_array,0); + size[1] = array_size(temp_array,1); + } + + if (!require_size(temp_array, size, 2)) SWIG_fail; + + array[i] = (DATA_TYPE*) array_data(temp_array); + } + + $1 = (DATA_TYPE**) array; + $3 = (DIM_TYPE) size[0]; + $4 = (DIM_TYPE) size[1]; +} +%typemap(freearg) + (DATA_TYPE** INPLACE_ARRAY3, DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3) +{ + if (array$argnum!=NULL) free(array$argnum); + if (object_array$argnum!=NULL) free(object_array$argnum); +} + +/* Typemap suite for (DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, + * DATA_TYPE* INPLACE_ARRAY3) + */ +%typecheck(SWIG_TYPECHECK_DOUBLE_ARRAY, + fragment="NumPy_Macros") + (DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, DATA_TYPE* INPLACE_ARRAY3) +{ + $1 = is_array($input) && PyArray_EquivTypenums(array_type($input), + DATA_TYPECODE); +} +%typemap(in, + fragment="NumPy_Fragments") + (DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, DATA_TYPE* INPLACE_ARRAY3) + (PyArrayObject* array=NULL) +{ + array = obj_to_array_no_conversion($input, DATA_TYPECODE); + if (!array || !require_dimensions(array,3) || !require_contiguous(array) + || !require_native(array)) SWIG_fail; + $1 = (DIM_TYPE) array_size(array,0); + $2 = (DIM_TYPE) array_size(array,1); + $3 = (DIM_TYPE) array_size(array,2); + $4 = (DATA_TYPE*) array_data(array); +} + +/* Typemap suite for (DATA_TYPE* INPLACE_FARRAY3, DIM_TYPE DIM1, DIM_TYPE DIM2, + * DIM_TYPE DIM3) + */ +%typecheck(SWIG_TYPECHECK_DOUBLE_ARRAY, + fragment="NumPy_Macros") + (DATA_TYPE* INPLACE_FARRAY3, DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3) +{ + $1 = is_array($input) && PyArray_EquivTypenums(array_type($input), + DATA_TYPECODE); +} +%typemap(in, + fragment="NumPy_Fragments") + (DATA_TYPE* INPLACE_FARRAY3, DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3) + (PyArrayObject* array=NULL) +{ + array = obj_to_array_no_conversion($input, DATA_TYPECODE); + if (!array || !require_dimensions(array,3) || !require_contiguous(array) || + !require_native(array) || !require_fortran(array)) SWIG_fail; + $1 = (DATA_TYPE*) array_data(array); + $2 = (DIM_TYPE) array_size(array,0); + $3 = (DIM_TYPE) array_size(array,1); + $4 = (DIM_TYPE) array_size(array,2); +} + +/* Typemap suite for (DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, + * DATA_TYPE* INPLACE_FARRAY3) + */ +%typecheck(SWIG_TYPECHECK_DOUBLE_ARRAY, + fragment="NumPy_Macros") + (DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, DATA_TYPE* INPLACE_FARRAY3) +{ + $1 = is_array($input) && PyArray_EquivTypenums(array_type($input), + DATA_TYPECODE); +} +%typemap(in, + fragment="NumPy_Fragments") + (DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, DATA_TYPE* INPLACE_FARRAY3) + (PyArrayObject* array=NULL) +{ + array = obj_to_array_no_conversion($input, DATA_TYPECODE); + if (!array || !require_dimensions(array,3) || !require_contiguous(array) + || !require_native(array) || !require_fortran(array)) SWIG_fail; + $1 = (DIM_TYPE) array_size(array,0); + $2 = (DIM_TYPE) array_size(array,1); + $3 = (DIM_TYPE) array_size(array,2); + $4 = (DATA_TYPE*) array_data(array); +} + +/* Typemap suite for (DATA_TYPE INPLACE_ARRAY4[ANY][ANY][ANY][ANY]) + */ +%typecheck(SWIG_TYPECHECK_DOUBLE_ARRAY, + fragment="NumPy_Macros") + (DATA_TYPE INPLACE_ARRAY4[ANY][ANY][ANY][ANY]) +{ + $1 = is_array($input) && PyArray_EquivTypenums(array_type($input), + DATA_TYPECODE); +} +%typemap(in, + fragment="NumPy_Fragments") + (DATA_TYPE INPLACE_ARRAY4[ANY][ANY][ANY][ANY]) + (PyArrayObject* array=NULL) +{ + npy_intp size[4] = { $1_dim0, $1_dim1, $1_dim2 , $1_dim3 }; + array = obj_to_array_no_conversion($input, DATA_TYPECODE); + if (!array || !require_dimensions(array,4) || !require_size(array, size, 4) || + !require_contiguous(array) || !require_native(array)) SWIG_fail; + $1 = ($1_ltype) array_data(array); +} + +/* Typemap suite for (DATA_TYPE* INPLACE_ARRAY4, DIM_TYPE DIM1, DIM_TYPE DIM2, + * DIM_TYPE DIM3, DIM_TYPE DIM4) + */ +%typecheck(SWIG_TYPECHECK_DOUBLE_ARRAY, + fragment="NumPy_Macros") + (DATA_TYPE* INPLACE_ARRAY4, DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, DIM_TYPE DIM4) +{ + $1 = is_array($input) && PyArray_EquivTypenums(array_type($input), + DATA_TYPECODE); +} +%typemap(in, + fragment="NumPy_Fragments") + (DATA_TYPE* INPLACE_ARRAY4, DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, DIM_TYPE DIM4) + (PyArrayObject* array=NULL) +{ + array = obj_to_array_no_conversion($input, DATA_TYPECODE); + if (!array || !require_dimensions(array,4) || !require_contiguous(array) || + !require_native(array)) SWIG_fail; + $1 = (DATA_TYPE*) array_data(array); + $2 = (DIM_TYPE) array_size(array,0); + $3 = (DIM_TYPE) array_size(array,1); + $4 = (DIM_TYPE) array_size(array,2); + $5 = (DIM_TYPE) array_size(array,3); +} + +/* Typemap suite for (DATA_TYPE** INPLACE_ARRAY4, DIM_TYPE DIM1, DIM_TYPE DIM2, + * DIM_TYPE DIM3, DIM_TYPE DIM4) + */ +%typecheck(SWIG_TYPECHECK_DOUBLE_ARRAY, + fragment="NumPy_Macros") + (DATA_TYPE** INPLACE_ARRAY4, DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, DIM_TYPE DIM4) +{ + $1 = PySequence_Check($input); +} +%typemap(in, + fragment="NumPy_Fragments") + (DATA_TYPE** INPLACE_ARRAY4, DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, DIM_TYPE DIM4) + (DATA_TYPE** array=NULL, PyArrayObject** object_array=NULL) +{ + npy_intp size[3] = { -1, -1, -1 }; + PyArrayObject* temp_array; + Py_ssize_t i; + + /* length of the list */ + $2 = PyList_Size($input); + + /* the arrays */ + array = (DATA_TYPE **)malloc($2*sizeof(DATA_TYPE *)); + object_array = (PyArrayObject **)calloc($2,sizeof(PyArrayObject *)); + + if (array == NULL || object_array == NULL) + { + SWIG_fail; + } + + for (i=0; i<$2; i++) + { + temp_array = obj_to_array_no_conversion(PySequence_GetItem($input,i), DATA_TYPECODE); + + /* the new array must be stored so that it can be destroyed in freearg */ + object_array[i] = temp_array; + + if ( !temp_array || !require_dimensions(temp_array, 3) || + !require_contiguous(temp_array) || + !require_native(temp_array) || + !PyArray_EquivTypenums(array_type(temp_array), DATA_TYPECODE) + ) SWIG_fail; + + /* store the size of the first array in the list, then use that for comparison. */ + if (i == 0) + { + size[0] = array_size(temp_array,0); + size[1] = array_size(temp_array,1); + size[2] = array_size(temp_array,2); + } + + if (!require_size(temp_array, size, 3)) SWIG_fail; + + array[i] = (DATA_TYPE*) array_data(temp_array); + } + + $1 = (DATA_TYPE**) array; + $3 = (DIM_TYPE) size[0]; + $4 = (DIM_TYPE) size[1]; + $5 = (DIM_TYPE) size[2]; +} +%typemap(freearg) + (DATA_TYPE** INPLACE_ARRAY4, DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, DIM_TYPE DIM4) +{ + if (array$argnum!=NULL) free(array$argnum); + if (object_array$argnum!=NULL) free(object_array$argnum); +} + +/* Typemap suite for (DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, DIM_TYPE DIM4, + * DATA_TYPE* INPLACE_ARRAY4) + */ +%typecheck(SWIG_TYPECHECK_DOUBLE_ARRAY, + fragment="NumPy_Macros") + (DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, DIM_TYPE DIM4, DATA_TYPE* INPLACE_ARRAY4) +{ + $1 = is_array($input) && PyArray_EquivTypenums(array_type($input), + DATA_TYPECODE); +} +%typemap(in, + fragment="NumPy_Fragments") + (DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, DIM_TYPE DIM4, DATA_TYPE* INPLACE_ARRAY4) + (PyArrayObject* array=NULL) +{ + array = obj_to_array_no_conversion($input, DATA_TYPECODE); + if (!array || !require_dimensions(array,4) || !require_contiguous(array) + || !require_native(array)) SWIG_fail; + $1 = (DIM_TYPE) array_size(array,0); + $2 = (DIM_TYPE) array_size(array,1); + $3 = (DIM_TYPE) array_size(array,2); + $4 = (DIM_TYPE) array_size(array,3); + $5 = (DATA_TYPE*) array_data(array); +} + +/* Typemap suite for (DATA_TYPE* INPLACE_FARRAY4, DIM_TYPE DIM1, DIM_TYPE DIM2, + * DIM_TYPE DIM3, DIM_TYPE DIM4) + */ +%typecheck(SWIG_TYPECHECK_DOUBLE_ARRAY, + fragment="NumPy_Macros") + (DATA_TYPE* INPLACE_FARRAY4, DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, DIM_TYPE DIM4) +{ + $1 = is_array($input) && PyArray_EquivTypenums(array_type($input), + DATA_TYPECODE); +} +%typemap(in, + fragment="NumPy_Fragments") + (DATA_TYPE* INPLACE_FARRAY4, DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, DIM_TYPE DIM4) + (PyArrayObject* array=NULL) +{ + array = obj_to_array_no_conversion($input, DATA_TYPECODE); + if (!array || !require_dimensions(array,4) || !require_contiguous(array) || + !require_native(array) || !require_fortran(array)) SWIG_fail; + $1 = (DATA_TYPE*) array_data(array); + $2 = (DIM_TYPE) array_size(array,0); + $3 = (DIM_TYPE) array_size(array,1); + $4 = (DIM_TYPE) array_size(array,2); + $5 = (DIM_TYPE) array_size(array,3); +} + +/* Typemap suite for (DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, + * DATA_TYPE* INPLACE_FARRAY4) + */ +%typecheck(SWIG_TYPECHECK_DOUBLE_ARRAY, + fragment="NumPy_Macros") + (DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, DIM_TYPE DIM4, DATA_TYPE* INPLACE_FARRAY4) +{ + $1 = is_array($input) && PyArray_EquivTypenums(array_type($input), + DATA_TYPECODE); +} +%typemap(in, + fragment="NumPy_Fragments") + (DIM_TYPE DIM1, DIM_TYPE DIM2, DIM_TYPE DIM3, DIM_TYPE DIM4, DATA_TYPE* INPLACE_FARRAY4) + (PyArrayObject* array=NULL) +{ + array = obj_to_array_no_conversion($input, DATA_TYPECODE); + if (!array || !require_dimensions(array,4) || !require_contiguous(array) + || !require_native(array) || !require_fortran(array)) SWIG_fail; + $1 = (DIM_TYPE) array_size(array,0); + $2 = (DIM_TYPE) array_size(array,1); + $3 = (DIM_TYPE) array_size(array,2); + $4 = (DIM_TYPE) array_size(array,3); + $5 = (DATA_TYPE*) array_data(array); +} + +/*************************/ +/* Argout Array Typemaps */ +/*************************/ + +/* Typemap suite for (DATA_TYPE ARGOUT_ARRAY1[ANY]) + */ +%typemap(in,numinputs=0, + fragment="NumPy_Backward_Compatibility,NumPy_Macros") + (DATA_TYPE ARGOUT_ARRAY1[ANY]) + (PyObject* array = NULL) +{ + npy_intp dims[1] = { $1_dim0 }; + array = PyArray_SimpleNew(1, dims, DATA_TYPECODE); + if (!array) SWIG_fail; + $1 = ($1_ltype) array_data(array); +} +%typemap(argout) + (DATA_TYPE ARGOUT_ARRAY1[ANY]) +{ + $result = SWIG_Python_AppendOutput($result,(PyObject*)array$argnum); +} + +/* Typemap suite for (DATA_TYPE* ARGOUT_ARRAY1, DIM_TYPE DIM1) + */ +%typemap(in,numinputs=1, + fragment="NumPy_Fragments") + (DATA_TYPE* ARGOUT_ARRAY1, DIM_TYPE DIM1) + (PyObject* array = NULL) +{ + npy_intp dims[1]; + if (!PyInt_Check($input)) + { + const char* typestring = pytype_string($input); + PyErr_Format(PyExc_TypeError, + "Int dimension expected. '%s' given.", + typestring); + SWIG_fail; + } + $2 = (DIM_TYPE) PyInt_AsLong($input); + dims[0] = (npy_intp) $2; + array = PyArray_SimpleNew(1, dims, DATA_TYPECODE); + if (!array) SWIG_fail; + $1 = (DATA_TYPE*) array_data(array); +} +%typemap(argout) + (DATA_TYPE* ARGOUT_ARRAY1, DIM_TYPE DIM1) +{ + $result = SWIG_Python_AppendOutput($result,(PyObject*)array$argnum); +} + +/* Typemap suite for (DIM_TYPE DIM1, DATA_TYPE* ARGOUT_ARRAY1) + */ +%typemap(in,numinputs=1, + fragment="NumPy_Fragments") + (DIM_TYPE DIM1, DATA_TYPE* ARGOUT_ARRAY1) + (PyObject* array = NULL) +{ + npy_intp dims[1]; + if (!PyInt_Check($input)) + { + const char* typestring = pytype_string($input); + PyErr_Format(PyExc_TypeError, + "Int dimension expected. '%s' given.", + typestring); + SWIG_fail; + } + $1 = (DIM_TYPE) PyInt_AsLong($input); + dims[0] = (npy_intp) $1; + array = PyArray_SimpleNew(1, dims, DATA_TYPECODE); + if (!array) SWIG_fail; + $2 = (DATA_TYPE*) array_data(array); +} +%typemap(argout) + (DIM_TYPE DIM1, DATA_TYPE* ARGOUT_ARRAY1) +{ + $result = SWIG_Python_AppendOutput($result,(PyObject*)array$argnum); +} + +/* Typemap suite for (DATA_TYPE ARGOUT_ARRAY2[ANY][ANY]) + */ +%typemap(in,numinputs=0, + fragment="NumPy_Backward_Compatibility,NumPy_Macros") + (DATA_TYPE ARGOUT_ARRAY2[ANY][ANY]) + (PyObject* array = NULL) +{ + npy_intp dims[2] = { $1_dim0, $1_dim1 }; + array = PyArray_SimpleNew(2, dims, DATA_TYPECODE); + if (!array) SWIG_fail; + $1 = ($1_ltype) array_data(array); +} +%typemap(argout) + (DATA_TYPE ARGOUT_ARRAY2[ANY][ANY]) +{ + $result = SWIG_Python_AppendOutput($result,(PyObject*)array$argnum); +} + +/* Typemap suite for (DATA_TYPE ARGOUT_ARRAY3[ANY][ANY][ANY]) + */ +%typemap(in,numinputs=0, + fragment="NumPy_Backward_Compatibility,NumPy_Macros") + (DATA_TYPE ARGOUT_ARRAY3[ANY][ANY][ANY]) + (PyObject* array = NULL) +{ + npy_intp dims[3] = { $1_dim0, $1_dim1, $1_dim2 }; + array = PyArray_SimpleNew(3, dims, DATA_TYPECODE); + if (!array) SWIG_fail; + $1 = ($1_ltype) array_data(array); +} +%typemap(argout) + (DATA_TYPE ARGOUT_ARRAY3[ANY][ANY][ANY]) +{ + $result = SWIG_Python_AppendOutput($result,(PyObject*)array$argnum); +} + +/* Typemap suite for (DATA_TYPE ARGOUT_ARRAY4[ANY][ANY][ANY][ANY]) + */ +%typemap(in,numinputs=0, + fragment="NumPy_Backward_Compatibility,NumPy_Macros") + (DATA_TYPE ARGOUT_ARRAY4[ANY][ANY][ANY][ANY]) + (PyObject* array = NULL) +{ + npy_intp dims[4] = { $1_dim0, $1_dim1, $1_dim2, $1_dim3 }; + array = PyArray_SimpleNew(4, dims, DATA_TYPECODE); + if (!array) SWIG_fail; + $1 = ($1_ltype) array_data(array); +} +%typemap(argout) + (DATA_TYPE ARGOUT_ARRAY4[ANY][ANY][ANY][ANY]) +{ + $result = SWIG_Python_AppendOutput($result,(PyObject*)array$argnum); +} + +/*****************************/ +/* Argoutview Array Typemaps */ +/*****************************/ + +/* Typemap suite for (DATA_TYPE** ARGOUTVIEW_ARRAY1, DIM_TYPE* DIM1) + */ +%typemap(in,numinputs=0) + (DATA_TYPE** ARGOUTVIEW_ARRAY1, DIM_TYPE* DIM1 ) + (DATA_TYPE* data_temp = NULL , DIM_TYPE dim_temp) +{ + $1 = &data_temp; + $2 = &dim_temp; +} +%typemap(argout, + fragment="NumPy_Backward_Compatibility") + (DATA_TYPE** ARGOUTVIEW_ARRAY1, DIM_TYPE* DIM1) +{ + npy_intp dims[1] = { *$2 }; + PyObject* obj = PyArray_SimpleNewFromData(1, dims, DATA_TYPECODE, (void*)(*$1)); + PyArrayObject* array = (PyArrayObject*) obj; + + if (!array) SWIG_fail; + $result = SWIG_Python_AppendOutput($result,obj); +} + +/* Typemap suite for (DIM_TYPE* DIM1, DATA_TYPE** ARGOUTVIEW_ARRAY1) + */ +%typemap(in,numinputs=0) + (DIM_TYPE* DIM1 , DATA_TYPE** ARGOUTVIEW_ARRAY1) + (DIM_TYPE dim_temp, DATA_TYPE* data_temp = NULL ) +{ + $1 = &dim_temp; + $2 = &data_temp; +} +%typemap(argout, + fragment="NumPy_Backward_Compatibility") + (DIM_TYPE* DIM1, DATA_TYPE** ARGOUTVIEW_ARRAY1) +{ + npy_intp dims[1] = { *$1 }; + PyObject* obj = PyArray_SimpleNewFromData(1, dims, DATA_TYPECODE, (void*)(*$2)); + PyArrayObject* array = (PyArrayObject*) obj; + + if (!array) SWIG_fail; + $result = SWIG_Python_AppendOutput($result,obj); +} + +/* Typemap suite for (DATA_TYPE** ARGOUTVIEW_ARRAY2, DIM_TYPE* DIM1, DIM_TYPE* DIM2) + */ +%typemap(in,numinputs=0) + (DATA_TYPE** ARGOUTVIEW_ARRAY2, DIM_TYPE* DIM1 , DIM_TYPE* DIM2 ) + (DATA_TYPE* data_temp = NULL , DIM_TYPE dim1_temp, DIM_TYPE dim2_temp) +{ + $1 = &data_temp; + $2 = &dim1_temp; + $3 = &dim2_temp; +} +%typemap(argout, + fragment="NumPy_Backward_Compatibility") + (DATA_TYPE** ARGOUTVIEW_ARRAY2, DIM_TYPE* DIM1, DIM_TYPE* DIM2) +{ + npy_intp dims[2] = { *$2, *$3 }; + PyObject* obj = PyArray_SimpleNewFromData(2, dims, DATA_TYPECODE, (void*)(*$1)); + PyArrayObject* array = (PyArrayObject*) obj; + + if (!array) SWIG_fail; + $result = SWIG_Python_AppendOutput($result,obj); +} + +/* Typemap suite for (DIM_TYPE* DIM1, DIM_TYPE* DIM2, DATA_TYPE** ARGOUTVIEW_ARRAY2) + */ +%typemap(in,numinputs=0) + (DIM_TYPE* DIM1 , DIM_TYPE* DIM2 , DATA_TYPE** ARGOUTVIEW_ARRAY2) + (DIM_TYPE dim1_temp, DIM_TYPE dim2_temp, DATA_TYPE* data_temp = NULL ) +{ + $1 = &dim1_temp; + $2 = &dim2_temp; + $3 = &data_temp; +} +%typemap(argout, + fragment="NumPy_Backward_Compatibility") + (DIM_TYPE* DIM1, DIM_TYPE* DIM2, DATA_TYPE** ARGOUTVIEW_ARRAY2) +{ + npy_intp dims[2] = { *$1, *$2 }; + PyObject* obj = PyArray_SimpleNewFromData(2, dims, DATA_TYPECODE, (void*)(*$3)); + PyArrayObject* array = (PyArrayObject*) obj; + + if (!array) SWIG_fail; + $result = SWIG_Python_AppendOutput($result,obj); +} + +/* Typemap suite for (DATA_TYPE** ARGOUTVIEW_FARRAY2, DIM_TYPE* DIM1, DIM_TYPE* DIM2) + */ +%typemap(in,numinputs=0) + (DATA_TYPE** ARGOUTVIEW_FARRAY2, DIM_TYPE* DIM1 , DIM_TYPE* DIM2 ) + (DATA_TYPE* data_temp = NULL , DIM_TYPE dim1_temp, DIM_TYPE dim2_temp) +{ + $1 = &data_temp; + $2 = &dim1_temp; + $3 = &dim2_temp; +} +%typemap(argout, + fragment="NumPy_Backward_Compatibility,NumPy_Array_Requirements") + (DATA_TYPE** ARGOUTVIEW_FARRAY2, DIM_TYPE* DIM1, DIM_TYPE* DIM2) +{ + npy_intp dims[2] = { *$2, *$3 }; + PyObject* obj = PyArray_SimpleNewFromData(2, dims, DATA_TYPECODE, (void*)(*$1)); + PyArrayObject* array = (PyArrayObject*) obj; + + if (!array || !require_fortran(array)) SWIG_fail; + $result = SWIG_Python_AppendOutput($result,obj); +} + +/* Typemap suite for (DIM_TYPE* DIM1, DIM_TYPE* DIM2, DATA_TYPE** ARGOUTVIEW_FARRAY2) + */ +%typemap(in,numinputs=0) + (DIM_TYPE* DIM1 , DIM_TYPE* DIM2 , DATA_TYPE** ARGOUTVIEW_FARRAY2) + (DIM_TYPE dim1_temp, DIM_TYPE dim2_temp, DATA_TYPE* data_temp = NULL ) +{ + $1 = &dim1_temp; + $2 = &dim2_temp; + $3 = &data_temp; +} +%typemap(argout, + fragment="NumPy_Backward_Compatibility,NumPy_Array_Requirements") + (DIM_TYPE* DIM1, DIM_TYPE* DIM2, DATA_TYPE** ARGOUTVIEW_FARRAY2) +{ + npy_intp dims[2] = { *$1, *$2 }; + PyObject* obj = PyArray_SimpleNewFromData(2, dims, DATA_TYPECODE, (void*)(*$3)); + PyArrayObject* array = (PyArrayObject*) obj; + + if (!array || !require_fortran(array)) SWIG_fail; + $result = SWIG_Python_AppendOutput($result,obj); +} + +/* Typemap suite for (DATA_TYPE** ARGOUTVIEW_ARRAY3, DIM_TYPE* DIM1, DIM_TYPE* DIM2, + DIM_TYPE* DIM3) + */ +%typemap(in,numinputs=0) + (DATA_TYPE** ARGOUTVIEW_ARRAY3, DIM_TYPE* DIM1 , DIM_TYPE* DIM2 , DIM_TYPE* DIM3 ) + (DATA_TYPE* data_temp = NULL , DIM_TYPE dim1_temp, DIM_TYPE dim2_temp, DIM_TYPE dim3_temp) +{ + $1 = &data_temp; + $2 = &dim1_temp; + $3 = &dim2_temp; + $4 = &dim3_temp; +} +%typemap(argout, + fragment="NumPy_Backward_Compatibility") + (DATA_TYPE** ARGOUTVIEW_ARRAY3, DIM_TYPE* DIM1, DIM_TYPE* DIM2, DIM_TYPE* DIM3) +{ + npy_intp dims[3] = { *$2, *$3, *$4 }; + PyObject* obj = PyArray_SimpleNewFromData(3, dims, DATA_TYPECODE, (void*)(*$1)); + PyArrayObject* array = (PyArrayObject*) obj; + + if (!array) SWIG_fail; + $result = SWIG_Python_AppendOutput($result,obj); +} + +/* Typemap suite for (DIM_TYPE* DIM1, DIM_TYPE* DIM2, DIM_TYPE* DIM3, + DATA_TYPE** ARGOUTVIEW_ARRAY3) + */ +%typemap(in,numinputs=0) + (DIM_TYPE* DIM1, DIM_TYPE* DIM2, DIM_TYPE* DIM3, DATA_TYPE** ARGOUTVIEW_ARRAY3) + (DIM_TYPE dim1_temp, DIM_TYPE dim2_temp, DIM_TYPE dim3_temp, DATA_TYPE* data_temp = NULL) +{ + $1 = &dim1_temp; + $2 = &dim2_temp; + $3 = &dim3_temp; + $4 = &data_temp; +} +%typemap(argout, + fragment="NumPy_Backward_Compatibility") + (DIM_TYPE* DIM1, DIM_TYPE* DIM2, DIM_TYPE* DIM3, DATA_TYPE** ARGOUTVIEW_ARRAY3) +{ + npy_intp dims[3] = { *$1, *$2, *$3 }; + PyObject* obj = PyArray_SimpleNewFromData(3, dims, DATA_TYPECODE, (void*)(*$4)); + PyArrayObject* array = (PyArrayObject*) obj; + + if (!array) SWIG_fail; + $result = SWIG_Python_AppendOutput($result,obj); +} + +/* Typemap suite for (DATA_TYPE** ARGOUTVIEW_FARRAY3, DIM_TYPE* DIM1, DIM_TYPE* DIM2, + DIM_TYPE* DIM3) + */ +%typemap(in,numinputs=0) + (DATA_TYPE** ARGOUTVIEW_FARRAY3, DIM_TYPE* DIM1 , DIM_TYPE* DIM2 , DIM_TYPE* DIM3 ) + (DATA_TYPE* data_temp = NULL , DIM_TYPE dim1_temp, DIM_TYPE dim2_temp, DIM_TYPE dim3_temp) +{ + $1 = &data_temp; + $2 = &dim1_temp; + $3 = &dim2_temp; + $4 = &dim3_temp; +} +%typemap(argout, + fragment="NumPy_Backward_Compatibility,NumPy_Array_Requirements") + (DATA_TYPE** ARGOUTVIEW_FARRAY3, DIM_TYPE* DIM1, DIM_TYPE* DIM2, DIM_TYPE* DIM3) +{ + npy_intp dims[3] = { *$2, *$3, *$4 }; + PyObject* obj = PyArray_SimpleNewFromData(3, dims, DATA_TYPECODE, (void*)(*$1)); + PyArrayObject* array = (PyArrayObject*) obj; + + if (!array || !require_fortran(array)) SWIG_fail; + $result = SWIG_Python_AppendOutput($result,obj); +} + +/* Typemap suite for (DIM_TYPE* DIM1, DIM_TYPE* DIM2, DIM_TYPE* DIM3, + DATA_TYPE** ARGOUTVIEW_FARRAY3) + */ +%typemap(in,numinputs=0) + (DIM_TYPE* DIM1 , DIM_TYPE* DIM2 , DIM_TYPE* DIM3 , DATA_TYPE** ARGOUTVIEW_FARRAY3) + (DIM_TYPE dim1_temp, DIM_TYPE dim2_temp, DIM_TYPE dim3_temp, DATA_TYPE* data_temp = NULL ) +{ + $1 = &dim1_temp; + $2 = &dim2_temp; + $3 = &dim3_temp; + $4 = &data_temp; +} +%typemap(argout, + fragment="NumPy_Backward_Compatibility,NumPy_Array_Requirements") + (DIM_TYPE* DIM1, DIM_TYPE* DIM2, DIM_TYPE* DIM3, DATA_TYPE** ARGOUTVIEW_FARRAY3) +{ + npy_intp dims[3] = { *$1, *$2, *$3 }; + PyObject* obj = PyArray_SimpleNewFromData(3, dims, DATA_TYPECODE, (void*)(*$4)); + PyArrayObject* array = (PyArrayObject*) obj; + + if (!array || !require_fortran(array)) SWIG_fail; + $result = SWIG_Python_AppendOutput($result,obj); +} + +/* Typemap suite for (DATA_TYPE** ARGOUTVIEW_ARRAY4, DIM_TYPE* DIM1, DIM_TYPE* DIM2, + DIM_TYPE* DIM3, DIM_TYPE* DIM4) + */ +%typemap(in,numinputs=0) + (DATA_TYPE** ARGOUTVIEW_ARRAY4, DIM_TYPE* DIM1 , DIM_TYPE* DIM2 , DIM_TYPE* DIM3 , DIM_TYPE* DIM4 ) + (DATA_TYPE* data_temp = NULL , DIM_TYPE dim1_temp, DIM_TYPE dim2_temp, DIM_TYPE dim3_temp, DIM_TYPE dim4_temp) +{ + $1 = &data_temp; + $2 = &dim1_temp; + $3 = &dim2_temp; + $4 = &dim3_temp; + $5 = &dim4_temp; +} +%typemap(argout, + fragment="NumPy_Backward_Compatibility") + (DATA_TYPE** ARGOUTVIEW_ARRAY4, DIM_TYPE* DIM1, DIM_TYPE* DIM2, DIM_TYPE* DIM3, DIM_TYPE* DIM4) +{ + npy_intp dims[4] = { *$2, *$3, *$4 , *$5 }; + PyObject* obj = PyArray_SimpleNewFromData(4, dims, DATA_TYPECODE, (void*)(*$1)); + PyArrayObject* array = (PyArrayObject*) obj; + + if (!array) SWIG_fail; + $result = SWIG_Python_AppendOutput($result,obj); +} + +/* Typemap suite for (DIM_TYPE* DIM1, DIM_TYPE* DIM2, DIM_TYPE* DIM3, DIM_TYPE* DIM4, + DATA_TYPE** ARGOUTVIEW_ARRAY4) + */ +%typemap(in,numinputs=0) + (DIM_TYPE* DIM1 , DIM_TYPE* DIM2 , DIM_TYPE* DIM3 , DIM_TYPE* DIM4 , DATA_TYPE** ARGOUTVIEW_ARRAY4) + (DIM_TYPE dim1_temp, DIM_TYPE dim2_temp, DIM_TYPE dim3_temp, DIM_TYPE dim4_temp, DATA_TYPE* data_temp = NULL ) +{ + $1 = &dim1_temp; + $2 = &dim2_temp; + $3 = &dim3_temp; + $4 = &dim4_temp; + $5 = &data_temp; +} +%typemap(argout, + fragment="NumPy_Backward_Compatibility") + (DIM_TYPE* DIM1, DIM_TYPE* DIM2, DIM_TYPE* DIM3, DIM_TYPE* DIM4, DATA_TYPE** ARGOUTVIEW_ARRAY4) +{ + npy_intp dims[4] = { *$1, *$2, *$3 , *$4 }; + PyObject* obj = PyArray_SimpleNewFromData(4, dims, DATA_TYPECODE, (void*)(*$5)); + PyArrayObject* array = (PyArrayObject*) obj; + + if (!array) SWIG_fail; + $result = SWIG_Python_AppendOutput($result,obj); +} + +/* Typemap suite for (DATA_TYPE** ARGOUTVIEW_FARRAY4, DIM_TYPE* DIM1, DIM_TYPE* DIM2, + DIM_TYPE* DIM3, DIM_TYPE* DIM4) + */ +%typemap(in,numinputs=0) + (DATA_TYPE** ARGOUTVIEW_FARRAY4, DIM_TYPE* DIM1 , DIM_TYPE* DIM2 , DIM_TYPE* DIM3 , DIM_TYPE* DIM4 ) + (DATA_TYPE* data_temp = NULL , DIM_TYPE dim1_temp, DIM_TYPE dim2_temp, DIM_TYPE dim3_temp, DIM_TYPE dim4_temp) +{ + $1 = &data_temp; + $2 = &dim1_temp; + $3 = &dim2_temp; + $4 = &dim3_temp; + $5 = &dim4_temp; +} +%typemap(argout, + fragment="NumPy_Backward_Compatibility,NumPy_Array_Requirements") + (DATA_TYPE** ARGOUTVIEW_FARRAY4, DIM_TYPE* DIM1, DIM_TYPE* DIM2, DIM_TYPE* DIM3, DIM_TYPE* DIM4) +{ + npy_intp dims[4] = { *$2, *$3, *$4 , *$5 }; + PyObject* obj = PyArray_SimpleNewFromData(4, dims, DATA_TYPECODE, (void*)(*$1)); + PyArrayObject* array = (PyArrayObject*) obj; + + if (!array || !require_fortran(array)) SWIG_fail; + $result = SWIG_Python_AppendOutput($result,obj); +} + +/* Typemap suite for (DIM_TYPE* DIM1, DIM_TYPE* DIM2, DIM_TYPE* DIM3, DIM_TYPE* DIM4, + DATA_TYPE** ARGOUTVIEW_FARRAY4) + */ +%typemap(in,numinputs=0) + (DIM_TYPE* DIM1 , DIM_TYPE* DIM2 , DIM_TYPE* DIM3 , DIM_TYPE* DIM4 , DATA_TYPE** ARGOUTVIEW_FARRAY4) + (DIM_TYPE dim1_temp, DIM_TYPE dim2_temp, DIM_TYPE dim3_temp, DIM_TYPE dim4_temp, DATA_TYPE* data_temp = NULL ) +{ + $1 = &dim1_temp; + $2 = &dim2_temp; + $3 = &dim3_temp; + $4 = &dim4_temp; + $5 = &data_temp; +} +%typemap(argout, + fragment="NumPy_Backward_Compatibility,NumPy_Array_Requirements") + (DIM_TYPE* DIM1, DIM_TYPE* DIM2, DIM_TYPE* DIM3, DIM_TYPE* DIM4, DATA_TYPE** ARGOUTVIEW_FARRAY4) +{ + npy_intp dims[4] = { *$1, *$2, *$3 , *$4 }; + PyObject* obj = PyArray_SimpleNewFromData(4, dims, DATA_TYPECODE, (void*)(*$5)); + PyArrayObject* array = (PyArrayObject*) obj; + + if (!array || !require_fortran(array)) SWIG_fail; + $result = SWIG_Python_AppendOutput($result,obj); +} + +/*************************************/ +/* Managed Argoutview Array Typemaps */ +/*************************************/ + +/* Typemap suite for (DATA_TYPE** ARGOUTVIEWM_ARRAY1, DIM_TYPE* DIM1) + */ +%typemap(in,numinputs=0) + (DATA_TYPE** ARGOUTVIEWM_ARRAY1, DIM_TYPE* DIM1 ) + (DATA_TYPE* data_temp = NULL , DIM_TYPE dim_temp) +{ + $1 = &data_temp; + $2 = &dim_temp; +} +%typemap(argout, + fragment="NumPy_Backward_Compatibility,NumPy_Utilities") + (DATA_TYPE** ARGOUTVIEWM_ARRAY1, DIM_TYPE* DIM1) +{ + npy_intp dims[1] = { *$2 }; + PyObject* obj = PyArray_SimpleNewFromData(1, dims, DATA_TYPECODE, (void*)(*$1)); + PyArrayObject* array = (PyArrayObject*) obj; + + if (!array) SWIG_fail; + +%#ifdef SWIGPY_USE_CAPSULE + PyObject* cap = PyCapsule_New((void*)(*$1), SWIGPY_CAPSULE_NAME, free_cap); +%#else + PyObject* cap = PyCObject_FromVoidPtr((void*)(*$1), free); +%#endif + +%#if NPY_API_VERSION < 0x00000007 + PyArray_BASE(array) = cap; +%#else + PyArray_SetBaseObject(array,cap); +%#endif + + $result = SWIG_Python_AppendOutput($result,obj); +} + +/* Typemap suite for (DIM_TYPE* DIM1, DATA_TYPE** ARGOUTVIEWM_ARRAY1) + */ +%typemap(in,numinputs=0) + (DIM_TYPE* DIM1 , DATA_TYPE** ARGOUTVIEWM_ARRAY1) + (DIM_TYPE dim_temp, DATA_TYPE* data_temp = NULL ) +{ + $1 = &dim_temp; + $2 = &data_temp; +} +%typemap(argout, + fragment="NumPy_Backward_Compatibility,NumPy_Utilities") + (DIM_TYPE* DIM1, DATA_TYPE** ARGOUTVIEWM_ARRAY1) +{ + npy_intp dims[1] = { *$1 }; + PyObject* obj = PyArray_SimpleNewFromData(1, dims, DATA_TYPECODE, (void*)(*$2)); + PyArrayObject* array = (PyArrayObject*) obj; + + if (!array) SWIG_fail; + +%#ifdef SWIGPY_USE_CAPSULE + PyObject* cap = PyCapsule_New((void*)(*$1), SWIGPY_CAPSULE_NAME, free_cap); +%#else + PyObject* cap = PyCObject_FromVoidPtr((void*)(*$1), free); +%#endif + +%#if NPY_API_VERSION < 0x00000007 + PyArray_BASE(array) = cap; +%#else + PyArray_SetBaseObject(array,cap); +%#endif + + $result = SWIG_Python_AppendOutput($result,obj); +} + +/* Typemap suite for (DATA_TYPE** ARGOUTVIEWM_ARRAY2, DIM_TYPE* DIM1, DIM_TYPE* DIM2) + */ +%typemap(in,numinputs=0) + (DATA_TYPE** ARGOUTVIEWM_ARRAY2, DIM_TYPE* DIM1 , DIM_TYPE* DIM2 ) + (DATA_TYPE* data_temp = NULL , DIM_TYPE dim1_temp, DIM_TYPE dim2_temp) +{ + $1 = &data_temp; + $2 = &dim1_temp; + $3 = &dim2_temp; +} +%typemap(argout, + fragment="NumPy_Backward_Compatibility,NumPy_Utilities") + (DATA_TYPE** ARGOUTVIEWM_ARRAY2, DIM_TYPE* DIM1, DIM_TYPE* DIM2) +{ + npy_intp dims[2] = { *$2, *$3 }; + PyObject* obj = PyArray_SimpleNewFromData(2, dims, DATA_TYPECODE, (void*)(*$1)); + PyArrayObject* array = (PyArrayObject*) obj; + + if (!array) SWIG_fail; + +%#ifdef SWIGPY_USE_CAPSULE + PyObject* cap = PyCapsule_New((void*)(*$1), SWIGPY_CAPSULE_NAME, free_cap); +%#else + PyObject* cap = PyCObject_FromVoidPtr((void*)(*$1), free); +%#endif + +%#if NPY_API_VERSION < 0x00000007 + PyArray_BASE(array) = cap; +%#else + PyArray_SetBaseObject(array,cap); +%#endif + + $result = SWIG_Python_AppendOutput($result,obj); +} + +/* Typemap suite for (DIM_TYPE* DIM1, DIM_TYPE* DIM2, DATA_TYPE** ARGOUTVIEWM_ARRAY2) + */ +%typemap(in,numinputs=0) + (DIM_TYPE* DIM1 , DIM_TYPE* DIM2 , DATA_TYPE** ARGOUTVIEWM_ARRAY2) + (DIM_TYPE dim1_temp, DIM_TYPE dim2_temp, DATA_TYPE* data_temp = NULL ) +{ + $1 = &dim1_temp; + $2 = &dim2_temp; + $3 = &data_temp; +} +%typemap(argout, + fragment="NumPy_Backward_Compatibility,NumPy_Utilities") + (DIM_TYPE* DIM1, DIM_TYPE* DIM2, DATA_TYPE** ARGOUTVIEWM_ARRAY2) +{ + npy_intp dims[2] = { *$1, *$2 }; + PyObject* obj = PyArray_SimpleNewFromData(2, dims, DATA_TYPECODE, (void*)(*$3)); + PyArrayObject* array = (PyArrayObject*) obj; + + if (!array) SWIG_fail; + +%#ifdef SWIGPY_USE_CAPSULE + PyObject* cap = PyCapsule_New((void*)(*$1), SWIGPY_CAPSULE_NAME, free_cap); +%#else + PyObject* cap = PyCObject_FromVoidPtr((void*)(*$1), free); +%#endif + +%#if NPY_API_VERSION < 0x00000007 + PyArray_BASE(array) = cap; +%#else + PyArray_SetBaseObject(array,cap); +%#endif + + $result = SWIG_Python_AppendOutput($result,obj); +} + +/* Typemap suite for (DATA_TYPE** ARGOUTVIEWM_FARRAY2, DIM_TYPE* DIM1, DIM_TYPE* DIM2) + */ +%typemap(in,numinputs=0) + (DATA_TYPE** ARGOUTVIEWM_FARRAY2, DIM_TYPE* DIM1 , DIM_TYPE* DIM2 ) + (DATA_TYPE* data_temp = NULL , DIM_TYPE dim1_temp, DIM_TYPE dim2_temp) +{ + $1 = &data_temp; + $2 = &dim1_temp; + $3 = &dim2_temp; +} +%typemap(argout, + fragment="NumPy_Backward_Compatibility,NumPy_Array_Requirements,NumPy_Utilities") + (DATA_TYPE** ARGOUTVIEWM_FARRAY2, DIM_TYPE* DIM1, DIM_TYPE* DIM2) +{ + npy_intp dims[2] = { *$2, *$3 }; + PyObject* obj = PyArray_SimpleNewFromData(2, dims, DATA_TYPECODE, (void*)(*$1)); + PyArrayObject* array = (PyArrayObject*) obj; + + if (!array || !require_fortran(array)) SWIG_fail; + +%#ifdef SWIGPY_USE_CAPSULE + PyObject* cap = PyCapsule_New((void*)(*$1), SWIGPY_CAPSULE_NAME, free_cap); +%#else + PyObject* cap = PyCObject_FromVoidPtr((void*)(*$1), free); +%#endif + +%#if NPY_API_VERSION < 0x00000007 + PyArray_BASE(array) = cap; +%#else + PyArray_SetBaseObject(array,cap); +%#endif + + $result = SWIG_Python_AppendOutput($result,obj); +} + +/* Typemap suite for (DIM_TYPE* DIM1, DIM_TYPE* DIM2, DATA_TYPE** ARGOUTVIEWM_FARRAY2) + */ +%typemap(in,numinputs=0) + (DIM_TYPE* DIM1 , DIM_TYPE* DIM2 , DATA_TYPE** ARGOUTVIEWM_FARRAY2) + (DIM_TYPE dim1_temp, DIM_TYPE dim2_temp, DATA_TYPE* data_temp = NULL ) +{ + $1 = &dim1_temp; + $2 = &dim2_temp; + $3 = &data_temp; +} +%typemap(argout, + fragment="NumPy_Backward_Compatibility,NumPy_Array_Requirements,NumPy_Utilities") + (DIM_TYPE* DIM1, DIM_TYPE* DIM2, DATA_TYPE** ARGOUTVIEWM_FARRAY2) +{ + npy_intp dims[2] = { *$1, *$2 }; + PyObject* obj = PyArray_SimpleNewFromData(2, dims, DATA_TYPECODE, (void*)(*$3)); + PyArrayObject* array = (PyArrayObject*) obj; + + if (!array || !require_fortran(array)) SWIG_fail; + +%#ifdef SWIGPY_USE_CAPSULE + PyObject* cap = PyCapsule_New((void*)(*$1), SWIGPY_CAPSULE_NAME, free_cap); +%#else + PyObject* cap = PyCObject_FromVoidPtr((void*)(*$1), free); +%#endif + +%#if NPY_API_VERSION < 0x00000007 + PyArray_BASE(array) = cap; +%#else + PyArray_SetBaseObject(array,cap); +%#endif + + $result = SWIG_Python_AppendOutput($result,obj); +} + +/* Typemap suite for (DATA_TYPE** ARGOUTVIEWM_ARRAY3, DIM_TYPE* DIM1, DIM_TYPE* DIM2, + DIM_TYPE* DIM3) + */ +%typemap(in,numinputs=0) + (DATA_TYPE** ARGOUTVIEWM_ARRAY3, DIM_TYPE* DIM1 , DIM_TYPE* DIM2 , DIM_TYPE* DIM3 ) + (DATA_TYPE* data_temp = NULL , DIM_TYPE dim1_temp, DIM_TYPE dim2_temp, DIM_TYPE dim3_temp) +{ + $1 = &data_temp; + $2 = &dim1_temp; + $3 = &dim2_temp; + $4 = &dim3_temp; +} +%typemap(argout, + fragment="NumPy_Backward_Compatibility,NumPy_Utilities") + (DATA_TYPE** ARGOUTVIEWM_ARRAY3, DIM_TYPE* DIM1, DIM_TYPE* DIM2, DIM_TYPE* DIM3) +{ + npy_intp dims[3] = { *$2, *$3, *$4 }; + PyObject* obj = PyArray_SimpleNewFromData(3, dims, DATA_TYPECODE, (void*)(*$1)); + PyArrayObject* array = (PyArrayObject*) obj; + + if (!array) SWIG_fail; + +%#ifdef SWIGPY_USE_CAPSULE + PyObject* cap = PyCapsule_New((void*)(*$1), SWIGPY_CAPSULE_NAME, free_cap); +%#else + PyObject* cap = PyCObject_FromVoidPtr((void*)(*$1), free); +%#endif + +%#if NPY_API_VERSION < 0x00000007 + PyArray_BASE(array) = cap; +%#else + PyArray_SetBaseObject(array,cap); +%#endif + + $result = SWIG_Python_AppendOutput($result,obj); +} + +/* Typemap suite for (DIM_TYPE* DIM1, DIM_TYPE* DIM2, DIM_TYPE* DIM3, + DATA_TYPE** ARGOUTVIEWM_ARRAY3) + */ +%typemap(in,numinputs=0) + (DIM_TYPE* DIM1 , DIM_TYPE* DIM2 , DIM_TYPE* DIM3 , DATA_TYPE** ARGOUTVIEWM_ARRAY3) + (DIM_TYPE dim1_temp, DIM_TYPE dim2_temp, DIM_TYPE dim3_temp, DATA_TYPE* data_temp = NULL ) +{ + $1 = &dim1_temp; + $2 = &dim2_temp; + $3 = &dim3_temp; + $4 = &data_temp; +} +%typemap(argout, + fragment="NumPy_Backward_Compatibility,NumPy_Utilities") + (DIM_TYPE* DIM1, DIM_TYPE* DIM2, DIM_TYPE* DIM3, DATA_TYPE** ARGOUTVIEWM_ARRAY3) +{ + npy_intp dims[3] = { *$1, *$2, *$3 }; + PyObject* obj= PyArray_SimpleNewFromData(3, dims, DATA_TYPECODE, (void*)(*$4)); + PyArrayObject* array = (PyArrayObject*) obj; + + if (!array) SWIG_fail; + +%#ifdef SWIGPY_USE_CAPSULE + PyObject* cap = PyCapsule_New((void*)(*$1), SWIGPY_CAPSULE_NAME, free_cap); +%#else + PyObject* cap = PyCObject_FromVoidPtr((void*)(*$1), free); +%#endif + +%#if NPY_API_VERSION < 0x00000007 + PyArray_BASE(array) = cap; +%#else + PyArray_SetBaseObject(array,cap); +%#endif + + $result = SWIG_Python_AppendOutput($result,obj); +} + +/* Typemap suite for (DATA_TYPE** ARGOUTVIEWM_FARRAY3, DIM_TYPE* DIM1, DIM_TYPE* DIM2, + DIM_TYPE* DIM3) + */ +%typemap(in,numinputs=0) + (DATA_TYPE** ARGOUTVIEWM_FARRAY3, DIM_TYPE* DIM1 , DIM_TYPE* DIM2 , DIM_TYPE* DIM3 ) + (DATA_TYPE* data_temp = NULL , DIM_TYPE dim1_temp, DIM_TYPE dim2_temp, DIM_TYPE dim3_temp) +{ + $1 = &data_temp; + $2 = &dim1_temp; + $3 = &dim2_temp; + $4 = &dim3_temp; +} +%typemap(argout, + fragment="NumPy_Backward_Compatibility,NumPy_Array_Requirements,NumPy_Utilities") + (DATA_TYPE** ARGOUTVIEWM_FARRAY3, DIM_TYPE* DIM1, DIM_TYPE* DIM2, DIM_TYPE* DIM3) +{ + npy_intp dims[3] = { *$2, *$3, *$4 }; + PyObject* obj = PyArray_SimpleNewFromData(3, dims, DATA_TYPECODE, (void*)(*$1)); + PyArrayObject* array = (PyArrayObject*) obj; + + if (!array || !require_fortran(array)) SWIG_fail; + +%#ifdef SWIGPY_USE_CAPSULE + PyObject* cap = PyCapsule_New((void*)(*$1), SWIGPY_CAPSULE_NAME, free_cap); +%#else + PyObject* cap = PyCObject_FromVoidPtr((void*)(*$1), free); +%#endif + +%#if NPY_API_VERSION < 0x00000007 + PyArray_BASE(array) = cap; +%#else + PyArray_SetBaseObject(array,cap); +%#endif + + $result = SWIG_Python_AppendOutput($result,obj); +} + +/* Typemap suite for (DIM_TYPE* DIM1, DIM_TYPE* DIM2, DIM_TYPE* DIM3, + DATA_TYPE** ARGOUTVIEWM_FARRAY3) + */ +%typemap(in,numinputs=0) + (DIM_TYPE* DIM1 , DIM_TYPE* DIM2 , DIM_TYPE* DIM3 , DATA_TYPE** ARGOUTVIEWM_FARRAY3) + (DIM_TYPE dim1_temp, DIM_TYPE dim2_temp, DIM_TYPE dim3_temp, DATA_TYPE* data_temp = NULL ) +{ + $1 = &dim1_temp; + $2 = &dim2_temp; + $3 = &dim3_temp; + $4 = &data_temp; +} +%typemap(argout, + fragment="NumPy_Backward_Compatibility,NumPy_Array_Requirements,NumPy_Utilities") + (DIM_TYPE* DIM1, DIM_TYPE* DIM2, DIM_TYPE* DIM3, DATA_TYPE** ARGOUTVIEWM_FARRAY3) +{ + npy_intp dims[3] = { *$1, *$2, *$3 }; + PyObject* obj = PyArray_SimpleNewFromData(3, dims, DATA_TYPECODE, (void*)(*$4)); + PyArrayObject* array = (PyArrayObject*) obj; + + if (!array || !require_fortran(array)) SWIG_fail; + +%#ifdef SWIGPY_USE_CAPSULE + PyObject* cap = PyCapsule_New((void*)(*$1), SWIGPY_CAPSULE_NAME, free_cap); +%#else + PyObject* cap = PyCObject_FromVoidPtr((void*)(*$1), free); +%#endif + +%#if NPY_API_VERSION < 0x00000007 + PyArray_BASE(array) = cap; +%#else + PyArray_SetBaseObject(array,cap); +%#endif + + $result = SWIG_Python_AppendOutput($result,obj); +} + +/* Typemap suite for (DATA_TYPE** ARGOUTVIEWM_ARRAY4, DIM_TYPE* DIM1, DIM_TYPE* DIM2, + DIM_TYPE* DIM3, DIM_TYPE* DIM4) + */ +%typemap(in,numinputs=0) + (DATA_TYPE** ARGOUTVIEWM_ARRAY4, DIM_TYPE* DIM1 , DIM_TYPE* DIM2 , DIM_TYPE* DIM3 , DIM_TYPE* DIM4 ) + (DATA_TYPE* data_temp = NULL , DIM_TYPE dim1_temp, DIM_TYPE dim2_temp, DIM_TYPE dim3_temp, DIM_TYPE dim4_temp) +{ + $1 = &data_temp; + $2 = &dim1_temp; + $3 = &dim2_temp; + $4 = &dim3_temp; + $5 = &dim4_temp; +} +%typemap(argout, + fragment="NumPy_Backward_Compatibility,NumPy_Utilities") + (DATA_TYPE** ARGOUTVIEWM_ARRAY4, DIM_TYPE* DIM1, DIM_TYPE* DIM2, DIM_TYPE* DIM3, DIM_TYPE* DIM4) +{ + npy_intp dims[4] = { *$2, *$3, *$4 , *$5 }; + PyObject* obj = PyArray_SimpleNewFromData(4, dims, DATA_TYPECODE, (void*)(*$1)); + PyArrayObject* array = (PyArrayObject*) obj; + + if (!array) SWIG_fail; + +%#ifdef SWIGPY_USE_CAPSULE + PyObject* cap = PyCapsule_New((void*)(*$1), SWIGPY_CAPSULE_NAME, free_cap); +%#else + PyObject* cap = PyCObject_FromVoidPtr((void*)(*$1), free); +%#endif + +%#if NPY_API_VERSION < 0x00000007 + PyArray_BASE(array) = cap; +%#else + PyArray_SetBaseObject(array,cap); +%#endif + + $result = SWIG_Python_AppendOutput($result,obj); +} + +/* Typemap suite for (DIM_TYPE* DIM1, DIM_TYPE* DIM2, DIM_TYPE* DIM3, DIM_TYPE* DIM4, + DATA_TYPE** ARGOUTVIEWM_ARRAY4) + */ +%typemap(in,numinputs=0) + (DIM_TYPE* DIM1 , DIM_TYPE* DIM2 , DIM_TYPE* DIM3 , DIM_TYPE* DIM4 , DATA_TYPE** ARGOUTVIEWM_ARRAY4) + (DIM_TYPE dim1_temp, DIM_TYPE dim2_temp, DIM_TYPE dim3_temp, DIM_TYPE dim4_temp, DATA_TYPE* data_temp = NULL ) +{ + $1 = &dim1_temp; + $2 = &dim2_temp; + $3 = &dim3_temp; + $4 = &dim4_temp; + $5 = &data_temp; +} +%typemap(argout, + fragment="NumPy_Backward_Compatibility,NumPy_Utilities") + (DIM_TYPE* DIM1, DIM_TYPE* DIM2, DIM_TYPE* DIM3, DIM_TYPE* DIM4, DATA_TYPE** ARGOUTVIEWM_ARRAY4) +{ + npy_intp dims[4] = { *$1, *$2, *$3 , *$4 }; + PyObject* obj = PyArray_SimpleNewFromData(4, dims, DATA_TYPECODE, (void*)(*$5)); + PyArrayObject* array = (PyArrayObject*) obj; + + if (!array) SWIG_fail; + +%#ifdef SWIGPY_USE_CAPSULE + PyObject* cap = PyCapsule_New((void*)(*$1), SWIGPY_CAPSULE_NAME, free_cap); +%#else + PyObject* cap = PyCObject_FromVoidPtr((void*)(*$1), free); +%#endif + +%#if NPY_API_VERSION < 0x00000007 + PyArray_BASE(array) = cap; +%#else + PyArray_SetBaseObject(array,cap); +%#endif + + $result = SWIG_Python_AppendOutput($result,obj); +} + +/* Typemap suite for (DATA_TYPE** ARGOUTVIEWM_FARRAY4, DIM_TYPE* DIM1, DIM_TYPE* DIM2, + DIM_TYPE* DIM3, DIM_TYPE* DIM4) + */ +%typemap(in,numinputs=0) + (DATA_TYPE** ARGOUTVIEWM_FARRAY4, DIM_TYPE* DIM1 , DIM_TYPE* DIM2 , DIM_TYPE* DIM3 , DIM_TYPE* DIM4 ) + (DATA_TYPE* data_temp = NULL , DIM_TYPE dim1_temp, DIM_TYPE dim2_temp, DIM_TYPE dim3_temp, DIM_TYPE dim4_temp) +{ + $1 = &data_temp; + $2 = &dim1_temp; + $3 = &dim2_temp; + $4 = &dim3_temp; + $5 = &dim4_temp; +} +%typemap(argout, + fragment="NumPy_Backward_Compatibility,NumPy_Array_Requirements,NumPy_Utilities") + (DATA_TYPE** ARGOUTVIEWM_FARRAY4, DIM_TYPE* DIM1, DIM_TYPE* DIM2, DIM_TYPE* DIM3) +{ + npy_intp dims[4] = { *$2, *$3, *$4 , *$5 }; + PyObject* obj = PyArray_SimpleNewFromData(4, dims, DATA_TYPECODE, (void*)(*$1)); + PyArrayObject* array = (PyArrayObject*) obj; + + if (!array || !require_fortran(array)) SWIG_fail; + +%#ifdef SWIGPY_USE_CAPSULE + PyObject* cap = PyCapsule_New((void*)(*$1), SWIGPY_CAPSULE_NAME, free_cap); +%#else + PyObject* cap = PyCObject_FromVoidPtr((void*)(*$1), free); +%#endif + +%#if NPY_API_VERSION < 0x00000007 + PyArray_BASE(array) = cap; +%#else + PyArray_SetBaseObject(array,cap); +%#endif + + $result = SWIG_Python_AppendOutput($result,obj); +} + +/* Typemap suite for (DIM_TYPE* DIM1, DIM_TYPE* DIM2, DIM_TYPE* DIM3, DIM_TYPE* DIM4, + DATA_TYPE** ARGOUTVIEWM_FARRAY4) + */ +%typemap(in,numinputs=0) + (DIM_TYPE* DIM1 , DIM_TYPE* DIM2 , DIM_TYPE* DIM3 , DIM_TYPE* DIM4 , DATA_TYPE** ARGOUTVIEWM_FARRAY4) + (DIM_TYPE dim1_temp, DIM_TYPE dim2_temp, DIM_TYPE dim3_temp, DIM_TYPE dim4_temp, DATA_TYPE* data_temp = NULL ) +{ + $1 = &dim1_temp; + $2 = &dim2_temp; + $3 = &dim3_temp; + $4 = &dim4_temp; + $5 = &data_temp; +} +%typemap(argout, + fragment="NumPy_Backward_Compatibility,NumPy_Array_Requirements,NumPy_Utilities") + (DIM_TYPE* DIM1, DIM_TYPE* DIM2, DIM_TYPE* DIM3, DIM_TYPE* DIM4, DATA_TYPE** ARGOUTVIEWM_FARRAY4) +{ + npy_intp dims[4] = { *$1, *$2, *$3 , *$4 }; + PyObject* obj = PyArray_SimpleNewFromData(4, dims, DATA_TYPECODE, (void*)(*$5)); + PyArrayObject* array = (PyArrayObject*) obj; + + if (!array || !require_fortran(array)) SWIG_fail; + +%#ifdef SWIGPY_USE_CAPSULE + PyObject* cap = PyCapsule_New((void*)(*$1), SWIGPY_CAPSULE_NAME, free_cap); +%#else + PyObject* cap = PyCObject_FromVoidPtr((void*)(*$1), free); +%#endif + +%#if NPY_API_VERSION < 0x00000007 + PyArray_BASE(array) = cap; +%#else + PyArray_SetBaseObject(array,cap); +%#endif + + $result = SWIG_Python_AppendOutput($result,obj); +} + +/* Typemap suite for (DATA_TYPE** ARGOUTVIEWM_ARRAY4, DIM_TYPE* DIM1, DIM_TYPE* DIM2, + DIM_TYPE* DIM3, DIM_TYPE* DIM4) + */ +%typemap(in,numinputs=0) + (DATA_TYPE** ARGOUTVIEWM_ARRAY4, DIM_TYPE* DIM1 , DIM_TYPE* DIM2 , DIM_TYPE* DIM3 , DIM_TYPE* DIM4 ) + (DATA_TYPE* data_temp = NULL , DIM_TYPE dim1_temp, DIM_TYPE dim2_temp, DIM_TYPE dim3_temp, DIM_TYPE dim4_temp) +{ + $1 = &data_temp; + $2 = &dim1_temp; + $3 = &dim2_temp; + $4 = &dim3_temp; + $5 = &dim4_temp; +} +%typemap(argout, + fragment="NumPy_Backward_Compatibility,NumPy_Utilities") + (DATA_TYPE** ARGOUTVIEWM_ARRAY4, DIM_TYPE* DIM1, DIM_TYPE* DIM2, DIM_TYPE* DIM3, DIM_TYPE* DIM4) +{ + npy_intp dims[4] = { *$2, *$3, *$4 , *$5 }; + PyObject* obj = PyArray_SimpleNewFromData(4, dims, DATA_TYPECODE, (void*)(*$1)); + PyArrayObject* array = (PyArrayObject*) obj; + + if (!array) SWIG_fail; + +%#ifdef SWIGPY_USE_CAPSULE + PyObject* cap = PyCapsule_New((void*)(*$1), SWIGPY_CAPSULE_NAME, free_cap); +%#else + PyObject* cap = PyCObject_FromVoidPtr((void*)(*$1), free); +%#endif + +%#if NPY_API_VERSION < 0x00000007 + PyArray_BASE(array) = cap; +%#else + PyArray_SetBaseObject(array,cap); +%#endif + + $result = SWIG_Python_AppendOutput($result,obj); +} + +/* Typemap suite for (DIM_TYPE* DIM1, DIM_TYPE* DIM2, DIM_TYPE* DIM3, DIM_TYPE* DIM4, + DATA_TYPE** ARGOUTVIEWM_ARRAY4) + */ +%typemap(in,numinputs=0) + (DIM_TYPE* DIM1 , DIM_TYPE* DIM2 , DIM_TYPE* DIM3 , DIM_TYPE* DIM4 , DATA_TYPE** ARGOUTVIEWM_ARRAY4) + (DIM_TYPE dim1_temp, DIM_TYPE dim2_temp, DIM_TYPE dim3_temp, DIM_TYPE dim4_temp, DATA_TYPE* data_temp = NULL ) +{ + $1 = &dim1_temp; + $2 = &dim2_temp; + $3 = &dim3_temp; + $4 = &dim4_temp; + $5 = &data_temp; +} +%typemap(argout, + fragment="NumPy_Backward_Compatibility,NumPy_Utilities") + (DIM_TYPE* DIM1, DIM_TYPE* DIM2, DIM_TYPE* DIM3, DIM_TYPE* DIM4, DATA_TYPE** ARGOUTVIEWM_ARRAY4) +{ + npy_intp dims[4] = { *$1, *$2, *$3 , *$4 }; + PyObject* obj = PyArray_SimpleNewFromData(4, dims, DATA_TYPECODE, (void*)(*$5)); + PyArrayObject* array = (PyArrayObject*) obj; + + if (!array) SWIG_fail; + +%#ifdef SWIGPY_USE_CAPSULE + PyObject* cap = PyCapsule_New((void*)(*$1), SWIGPY_CAPSULE_NAME, free_cap); +%#else + PyObject* cap = PyCObject_FromVoidPtr((void*)(*$1), free); +%#endif + +%#if NPY_API_VERSION < 0x00000007 + PyArray_BASE(array) = cap; +%#else + PyArray_SetBaseObject(array,cap); +%#endif + + $result = SWIG_Python_AppendOutput($result,obj); +} + +/* Typemap suite for (DATA_TYPE** ARGOUTVIEWM_FARRAY4, DIM_TYPE* DIM1, DIM_TYPE* DIM2, + DIM_TYPE* DIM3, DIM_TYPE* DIM4) + */ +%typemap(in,numinputs=0) + (DATA_TYPE** ARGOUTVIEWM_FARRAY4, DIM_TYPE* DIM1 , DIM_TYPE* DIM2 , DIM_TYPE* DIM3 , DIM_TYPE* DIM4 ) + (DATA_TYPE* data_temp = NULL , DIM_TYPE dim1_temp, DIM_TYPE dim2_temp, DIM_TYPE dim3_temp, DIM_TYPE dim4_temp) +{ + $1 = &data_temp; + $2 = &dim1_temp; + $3 = &dim2_temp; + $4 = &dim3_temp; + $5 = &dim4_temp; +} +%typemap(argout, + fragment="NumPy_Backward_Compatibility,NumPy_Array_Requirements,NumPy_Utilities") + (DATA_TYPE** ARGOUTVIEWM_FARRAY4, DIM_TYPE* DIM1, DIM_TYPE* DIM2, DIM_TYPE* DIM3, DIM_TYPE* DIM4) +{ + npy_intp dims[4] = { *$2, *$3, *$4 , *$5 }; + PyObject* obj = PyArray_SimpleNewFromData(4, dims, DATA_TYPECODE, (void*)(*$1)); + PyArrayObject* array = (PyArrayObject*) obj; + + if (!array || !require_fortran(array)) SWIG_fail; + +%#ifdef SWIGPY_USE_CAPSULE + PyObject* cap = PyCapsule_New((void*)(*$1), SWIGPY_CAPSULE_NAME, free_cap); +%#else + PyObject* cap = PyCObject_FromVoidPtr((void*)(*$1), free); +%#endif + +%#if NPY_API_VERSION < 0x00000007 + PyArray_BASE(array) = cap; +%#else + PyArray_SetBaseObject(array,cap); +%#endif + + $result = SWIG_Python_AppendOutput($result,obj); +} + +/* Typemap suite for (DIM_TYPE* DIM1, DIM_TYPE* DIM2, DIM_TYPE* DIM3, DIM_TYPE* DIM4, + DATA_TYPE** ARGOUTVIEWM_FARRAY4) + */ +%typemap(in,numinputs=0) + (DIM_TYPE* DIM1 , DIM_TYPE* DIM2 , DIM_TYPE* DIM3 , DIM_TYPE* DIM4 , DATA_TYPE** ARGOUTVIEWM_FARRAY4) + (DIM_TYPE dim1_temp, DIM_TYPE dim2_temp, DIM_TYPE dim3_temp, DIM_TYPE dim4_temp, DATA_TYPE* data_temp = NULL ) +{ + $1 = &dim1_temp; + $2 = &dim2_temp; + $3 = &dim3_temp; + $4 = &dim4_temp; + $5 = &data_temp; +} +%typemap(argout, + fragment="NumPy_Backward_Compatibility,NumPy_Array_Requirements,NumPy_Utilities") + (DIM_TYPE* DIM1, DIM_TYPE* DIM2, DIM_TYPE* DIM3, DIM_TYPE* DIM4, DATA_TYPE** ARGOUTVIEWM_FARRAY4) +{ + npy_intp dims[4] = { *$1, *$2, *$3 , *$4 }; + PyObject* obj = PyArray_SimpleNewFromData(4, dims, DATA_TYPECODE, (void*)(*$5)); + PyArrayObject* array = (PyArrayObject*) obj; + + if (!array || !require_fortran(array)) SWIG_fail; + +%#ifdef SWIGPY_USE_CAPSULE + PyObject* cap = PyCapsule_New((void*)(*$1), SWIGPY_CAPSULE_NAME, free_cap); +%#else + PyObject* cap = PyCObject_FromVoidPtr((void*)(*$1), free); +%#endif + +%#if NPY_API_VERSION < 0x00000007 + PyArray_BASE(array) = cap; +%#else + PyArray_SetBaseObject(array,cap); +%#endif + + $result = SWIG_Python_AppendOutput($result,obj); +} + +/**************************************/ +/* In-Place Array Typemap - flattened */ +/**************************************/ + +/* Typemap suite for (DATA_TYPE* INPLACE_ARRAY_FLAT, DIM_TYPE DIM_FLAT) + */ +%typecheck(SWIG_TYPECHECK_DOUBLE_ARRAY, + fragment="NumPy_Macros") + (DATA_TYPE* INPLACE_ARRAY_FLAT, DIM_TYPE DIM_FLAT) +{ + $1 = is_array($input) && PyArray_EquivTypenums(array_type($input), + DATA_TYPECODE); +} +%typemap(in, + fragment="NumPy_Fragments") + (DATA_TYPE* INPLACE_ARRAY_FLAT, DIM_TYPE DIM_FLAT) + (PyArrayObject* array=NULL, int i=1) +{ + array = obj_to_array_no_conversion($input, DATA_TYPECODE); + if (!array || !require_c_or_f_contiguous(array) + || !require_native(array)) SWIG_fail; + $1 = (DATA_TYPE*) array_data(array); + $2 = 1; + for (i=0; i < array_numdims(array); ++i) $2 *= array_size(array,i); +} + +%enddef /* %numpy_typemaps() macro */ +/* *************************************************************** */ + +/* Concrete instances of the %numpy_typemaps() macro: Each invocation + * below applies all of the typemaps above to the specified data type. + */ +%numpy_typemaps(signed char , NPY_BYTE , int) +%numpy_typemaps(unsigned char , NPY_UBYTE , int) +%numpy_typemaps(short , NPY_SHORT , int) +%numpy_typemaps(unsigned short , NPY_USHORT , int) +%numpy_typemaps(int , NPY_INT , int) +%numpy_typemaps(unsigned int , NPY_UINT , int) +%numpy_typemaps(long , NPY_LONG , int) +%numpy_typemaps(unsigned long , NPY_ULONG , int) +%numpy_typemaps(long long , NPY_LONGLONG , int) +%numpy_typemaps(unsigned long long, NPY_ULONGLONG, int) +%numpy_typemaps(float , NPY_FLOAT , int) +%numpy_typemaps(double , NPY_DOUBLE , int) + +/* *************************************************************** + * The follow macro expansion does not work, because C++ bool is 4 + * bytes and NPY_BOOL is 1 byte + * + * %numpy_typemaps(bool, NPY_BOOL, int) + */ + +/* *************************************************************** + * On my Mac, I get the following warning for this macro expansion: + * 'swig/python detected a memory leak of type 'long double *', no destructor found.' + * + * %numpy_typemaps(long double, NPY_LONGDOUBLE, int) + */ + +#ifdef __cplusplus + +%include + +%numpy_typemaps(std::complex, NPY_CFLOAT , int) +%numpy_typemaps(std::complex, NPY_CDOUBLE, int) + +#endif + +#endif /* SWIGPYTHON */ diff --git a/tractor/setup-mpf.py b/tractor/setup-mpf.py index a6b98d4c..82c6f3b2 100644 --- a/tractor/setup-mpf.py +++ b/tractor/setup-mpf.py @@ -9,6 +9,11 @@ extra_compile_args=['-g'], extra_link_args=['-g'], ) +my_swig_mod = Extension('_c_mp_fourier', + sources=['c_mp_fourier.i'], + include_dirs=numpy_inc, + extra_compile_args=['-g', '-O0'], + extra_link_args=['-g']) setup(name = 'Gaussian mixtures -- Fourier transform', version = '1.0', @@ -16,4 +21,4 @@ author = 'Lang & Hogg', author_email = 'dstndstn@gmail.com', url = 'http://astrometry.net', - ext_modules = [c_swig_module]) + ext_modules = [c_swig_module, my_swig_mod]) From cb57eef98d451576b5811ee9f441f046728b5a51 Mon Sep 17 00:00:00 2001 From: Chris Hogan Date: Tue, 10 Oct 2017 16:26:19 -0500 Subject: [PATCH 02/43] Better vectorization --- tractor/c_mp_fourier.i | 61 +++++++----------------------------------- tractor/setup-mpf.py | 2 +- 2 files changed, 11 insertions(+), 52 deletions(-) diff --git a/tractor/c_mp_fourier.i b/tractor/c_mp_fourier.i index dd44fee8..d2feae94 100644 --- a/tractor/c_mp_fourier.i +++ b/tractor/c_mp_fourier.i @@ -40,7 +40,6 @@ static void mixture_profile_fourier_transform(double *amps, int amps_len, double *w, int w_len, double *out, int out_dim1, int out_dim2) { - // const int D = 2; const double twopisquare = -2. * M_PI * M_PI; int K = amps_len; @@ -48,67 +47,27 @@ static void mixture_profile_fourier_transform(double *amps, int amps_len, int NW = w_len; int i, j, k; - // assert(means_dim1 == K); - // assert(means_dim2 == D); - // assert(vars_dim1 == K); - // assert(vars_dim2 == D); - // assert(vars_dim3 == D); - - // for (k = 0; k < K; k++) { - // if ((means[k*D] != means[0]) || - // (means[k*D+1] != means[1])) { - // PyErr_SetString(PyExc_ValueError, "Assume all means are equal"); - // return NULL; - // } - // } - - double mu0 = means[0]; - double mu1 = means[1]; - - int zeromean = (mu0 == 0.) && (mu1 == 0.); - - // TODO: Handle complex out - - if (zeromean) { - // out = (double*)malloc(sizeof(double) * NW * NV); - } - // else { - // out = new std::complex[NW * NV]; - // } - - double* ff = out; + double *s = (double*)malloc(sizeof(double) * NV); for (j = 0; j < NW; j++) { for (i = 0; i < NV; i++) { - double s = 0; - double* V = vars; + s[i] = 0; for (k = 0; k < K; k++) { double a, b, d; - a = *V; - V++; - b = *V; - V++; - // skip c - V++; - d = *V; - V++; + int offset = k * 4; - s += amps[k] * exp(twopisquare * (a * v[i] * v[i] + + a = vars[offset]; + b = vars[offset + 1]; + d = vars[offset + 3]; + + s[i] += amps[k] * exp(twopisquare * (a * v[i] * v[i] + 2. * b * v[i] * w[j] + d * w[j] * w[j])); } - if (zeromean) { - *ff = s; - ff++; - } else { - double angle = -2. * M_PI * (mu0 * v[i] + mu1 * w[j]); - *ff = s * cos(angle); - ff++; - *ff = s * sin(angle); - ff++; - } + out[NV * j + i] = s[i]; } } + free(s); return; } diff --git a/tractor/setup-mpf.py b/tractor/setup-mpf.py index 82c6f3b2..c8ad10a3 100644 --- a/tractor/setup-mpf.py +++ b/tractor/setup-mpf.py @@ -12,7 +12,7 @@ my_swig_mod = Extension('_c_mp_fourier', sources=['c_mp_fourier.i'], include_dirs=numpy_inc, - extra_compile_args=['-g', '-O0'], + extra_compile_args=['-g'], extra_link_args=['-g']) setup(name = 'Gaussian mixtures -- Fourier transform', From fd2201621490500483217eb1d8d485c99e696743 Mon Sep 17 00:00:00 2001 From: Chris Hogan Date: Tue, 10 Oct 2017 16:28:45 -0500 Subject: [PATCH 03/43] Icc options --- tractor/setup-mpf.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tractor/setup-mpf.py b/tractor/setup-mpf.py index c8ad10a3..8ceacead 100644 --- a/tractor/setup-mpf.py +++ b/tractor/setup-mpf.py @@ -12,8 +12,8 @@ my_swig_mod = Extension('_c_mp_fourier', sources=['c_mp_fourier.i'], include_dirs=numpy_inc, - extra_compile_args=['-g'], - extra_link_args=['-g']) + extra_compile_args=['-g', '-xhost'], + extra_link_args=['-g', '-limf', '-lsvml']) setup(name = 'Gaussian mixtures -- Fourier transform', version = '1.0', From ef20660429b272530397ae5c947a5d2e158441f9 Mon Sep 17 00:00:00 2001 From: Chris Hogan Date: Wed, 11 Oct 2017 09:12:23 -0500 Subject: [PATCH 04/43] Final vectorized version --- tractor/c_mp_fourier.i | 162 ++--------------------------------------- 1 file changed, 7 insertions(+), 155 deletions(-) diff --git a/tractor/c_mp_fourier.i b/tractor/c_mp_fourier.i index d2feae94..348f0c76 100644 --- a/tractor/c_mp_fourier.i +++ b/tractor/c_mp_fourier.i @@ -47,24 +47,22 @@ static void mixture_profile_fourier_transform(double *amps, int amps_len, int NW = w_len; int i, j, k; - double *s = (double*)malloc(sizeof(double) * NV); + double *s = (double*)malloc(sizeof(double) * NV * NW); + memset(s, 0, sizeof(double) * NV * NW); for (j = 0; j < NW; j++) { for (i = 0; i < NV; i++) { - s[i] = 0; for (k = 0; k < K; k++) { - double a, b, d; int offset = k * 4; + double a = vars[offset]; + double b = vars[offset + 1]; + double d = vars[offset + 3]; - a = vars[offset]; - b = vars[offset + 1]; - d = vars[offset + 3]; - - s[i] += amps[k] * exp(twopisquare * (a * v[i] * v[i] + + s[NV * j + i] += amps[k] * exp(twopisquare * (a * v[i] * v[i] + 2. * b * v[i] * w[j] + d * w[j] * w[j])); } - out[NV * j + i] = s[i]; + out[NV * j + i] = s[NV * j + i]; } } free(s); @@ -72,149 +70,3 @@ static void mixture_profile_fourier_transform(double *amps, int amps_len, } %} - -// static PyObject* mixture_profile_fourier_transform( -// PyObject* po_amps, -// PyObject* po_means, -// PyObject* po_vars, -// PyObject* po_v, -// PyObject* po_w -// ) { - - // npy_intp K, NW,NV; - // const npy_intp D = 2; - // npy_intp i,j,k; - // double* amps, *means, *vars, *vv, *ww; - // PyObject* np_F; - // double* f; - // npy_intp dims[2]; - - // PyArrayObject *np_amps, *np_means, *np_vars, *np_v, *np_w; - - // if (!PyArray_Check(po_amps) || !PyArray_Check(po_means) || - // !PyArray_Check(po_vars) || !PyArray_Check(po_v) || - // !PyArray_Check(po_w)) { - // PyErr_SetString(PyExc_ValueError, "Expected numpy arrays"); - // return NULL; - // } - - // np_amps = (PyArrayObject*)po_amps; - // np_means = (PyArrayObject*)po_means; - // np_vars = (PyArrayObject*)po_vars; - // np_v = (PyArrayObject*)po_v; - // np_w = (PyArrayObject*)po_w; - - // if ((PyArray_TYPE(np_amps) != NPY_DOUBLE) || - // (PyArray_TYPE(np_means ) != NPY_DOUBLE) || - // (PyArray_TYPE(np_vars) != NPY_DOUBLE) || - // (PyArray_TYPE(np_v) != NPY_DOUBLE) || - // (PyArray_TYPE(np_w) != NPY_DOUBLE)) { - // PyErr_SetString(PyExc_ValueError, "Expected numpy double arrays"); - // return NULL; - // } - - // if (PyArray_NDIM(np_amps) != 1) { - // PyErr_SetString(PyExc_ValueError, "Expected 'amps' to be 1-d"); - // return NULL; - // } - // K = PyArray_DIM(np_amps, 0); - // if (PyArray_NDIM(np_means) != 2) { - // PyErr_SetString(PyExc_ValueError, "Expected 'means' to be 2-d"); - // return NULL; - // } - // if ((PyArray_DIM(np_means, 0) != K) || - // (PyArray_DIM(np_means, 1) != D)) { - // PyErr_SetString(PyExc_ValueError, "Expected 'means' to be K x D"); - // return NULL; - // } - // if (PyArray_NDIM(np_vars) != 3) { - // PyErr_SetString(PyExc_ValueError, "Expected 'vars' to be 3-d"); - // return NULL; - // } - // if ((PyArray_DIM(np_vars, 0) != K) || - // (PyArray_DIM(np_vars, 1) != D) || - // (PyArray_DIM(np_vars, 2) != D)) { - // PyErr_SetString(PyExc_ValueError, "Expected 'vars' to be K x D x D"); - // return NULL; - // } - - // if (PyArray_NDIM(np_v) != 1) { - // PyErr_SetString(PyExc_ValueError, "Expected 'v' to be 1-d"); - // return NULL; - // } - // if (PyArray_NDIM(np_w) != 1) { - // PyErr_SetString(PyExc_ValueError, "Expected 'w' to be 1-d"); - // return NULL; - // } - - // means = PyArray_DATA(np_means); - - // for (k=0; k Date: Tue, 17 Oct 2017 14:05:15 -0500 Subject: [PATCH 05/43] Call my c_mp_fourier --- tractor/mixture_profiles.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tractor/mixture_profiles.py b/tractor/mixture_profiles.py index 4b20ee94..ba4750ba 100644 --- a/tractor/mixture_profiles.py +++ b/tractor/mixture_profiles.py @@ -9,7 +9,7 @@ #import scipy.spatial.distance as scp try: - from tractor import mp_fourier + from tractor import c_mp_fourier as mp_fourier except: mp_fourier = None From 88e76ea1feea7bbbffc78e48374a0e5d0a6315ff Mon Sep 17 00:00:00 2001 From: Chris Hogan Date: Tue, 17 Oct 2017 14:12:42 -0500 Subject: [PATCH 06/43] Adjust fourier_transform call --- tractor/mixture_profiles.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/tractor/mixture_profiles.py b/tractor/mixture_profiles.py index ba4750ba..0383ff02 100644 --- a/tractor/mixture_profiles.py +++ b/tractor/mixture_profiles.py @@ -159,8 +159,9 @@ def getFourierTransform(self, v, w, use_mp_fourier=True): ''' if mp_fourier and use_mp_fourier: - f = mp_fourier.mixture_profile_fourier_transform( - self.amp, self.mean, self.var, v, w) + f = np.zeros((self.w, self.v)) + mp_fourier.mixture_profile_fourier_transform( + self.amp, self.mean, self.var, v, w, f) return f Fsum = None From 42232672e00a43dc24b9cc9e1c62d1c28dd4e649 Mon Sep 17 00:00:00 2001 From: Chris Hogan Date: Wed, 18 Oct 2017 16:29:18 -0500 Subject: [PATCH 07/43] Enable auto parallelization --- tractor/c_mp_fourier.i | 10 ++++++---- tractor/mixture_profiles.py | 2 +- tractor/setup-mpf.py | 4 ++-- 3 files changed, 9 insertions(+), 7 deletions(-) diff --git a/tractor/c_mp_fourier.i b/tractor/c_mp_fourier.i index 348f0c76..2d80e80f 100644 --- a/tractor/c_mp_fourier.i +++ b/tractor/c_mp_fourier.i @@ -50,19 +50,21 @@ static void mixture_profile_fourier_transform(double *amps, int amps_len, double *s = (double*)malloc(sizeof(double) * NV * NW); memset(s, 0, sizeof(double) * NV * NW); + #pragma parallel for (j = 0; j < NW; j++) { for (i = 0; i < NV; i++) { + int index = NV * j + i; for (k = 0; k < K; k++) { int offset = k * 4; double a = vars[offset]; double b = vars[offset + 1]; double d = vars[offset + 3]; - s[NV * j + i] += amps[k] * exp(twopisquare * (a * v[i] * v[i] + - 2. * b * v[i] * w[j] + - d * w[j] * w[j])); + s[index] += amps[k] * exp(twopisquare * (a * v[i] * v[i] + + 2. * b * v[i] * w[j] + + d * w[j] * w[j])); } - out[NV * j + i] = s[NV * j + i]; + out[index] = s[index]; } } free(s); diff --git a/tractor/mixture_profiles.py b/tractor/mixture_profiles.py index 0383ff02..675dbd55 100644 --- a/tractor/mixture_profiles.py +++ b/tractor/mixture_profiles.py @@ -159,7 +159,7 @@ def getFourierTransform(self, v, w, use_mp_fourier=True): ''' if mp_fourier and use_mp_fourier: - f = np.zeros((self.w, self.v)) + f = np.zeros((len(w), len(v))) mp_fourier.mixture_profile_fourier_transform( self.amp, self.mean, self.var, v, w, f) return f diff --git a/tractor/setup-mpf.py b/tractor/setup-mpf.py index 8ceacead..5d19e4cd 100644 --- a/tractor/setup-mpf.py +++ b/tractor/setup-mpf.py @@ -12,8 +12,8 @@ my_swig_mod = Extension('_c_mp_fourier', sources=['c_mp_fourier.i'], include_dirs=numpy_inc, - extra_compile_args=['-g', '-xhost'], - extra_link_args=['-g', '-limf', '-lsvml']) + extra_compile_args=['-g', '-xhost', '-parallel', '-qopt-report=5'], + extra_link_args=['-g', '-lsvml']) setup(name = 'Gaussian mixtures -- Fourier transform', version = '1.0', From 9b7137b4e93d3eefe8743715174ea0f8ca2304ce Mon Sep 17 00:00:00 2001 From: Chris Hogan Date: Thu, 19 Oct 2017 15:28:36 -0500 Subject: [PATCH 08/43] Final c_mp_fourier.i --- tractor/c_mp_fourier.i | 10 +++++----- tractor/mixture_profiles.py | 0 tractor/setup-mpf.py | 2 +- 3 files changed, 6 insertions(+), 6 deletions(-) mode change 100644 => 100755 tractor/c_mp_fourier.i mode change 100644 => 100755 tractor/mixture_profiles.py mode change 100644 => 100755 tractor/setup-mpf.py diff --git a/tractor/c_mp_fourier.i b/tractor/c_mp_fourier.i old mode 100644 new mode 100755 index 2d80e80f..c39d6698 --- a/tractor/c_mp_fourier.i +++ b/tractor/c_mp_fourier.i @@ -6,7 +6,6 @@ #define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION #include #include -#include %} %include "numpy.i" @@ -50,19 +49,20 @@ static void mixture_profile_fourier_transform(double *amps, int amps_len, double *s = (double*)malloc(sizeof(double) * NV * NW); memset(s, 0, sizeof(double) * NV * NW); - #pragma parallel for (j = 0; j < NW; j++) { + double w_j = w[j]; + double w_j_sqr = w_j * w_j; for (i = 0; i < NV; i++) { int index = NV * j + i; + double v_i = v[i]; + double v_i_sqr = v_i * v_i; for (k = 0; k < K; k++) { int offset = k * 4; double a = vars[offset]; double b = vars[offset + 1]; double d = vars[offset + 3]; - s[index] += amps[k] * exp(twopisquare * (a * v[i] * v[i] + - 2. * b * v[i] * w[j] + - d * w[j] * w[j])); + s[index] += amps[k] * exp(twopisquare * (a * v_i_sqr + 2. * b * v_i * w_j + d * w_j_sqr)); } out[index] = s[index]; } diff --git a/tractor/mixture_profiles.py b/tractor/mixture_profiles.py old mode 100644 new mode 100755 diff --git a/tractor/setup-mpf.py b/tractor/setup-mpf.py old mode 100644 new mode 100755 index 5d19e4cd..8df6a558 --- a/tractor/setup-mpf.py +++ b/tractor/setup-mpf.py @@ -12,7 +12,7 @@ my_swig_mod = Extension('_c_mp_fourier', sources=['c_mp_fourier.i'], include_dirs=numpy_inc, - extra_compile_args=['-g', '-xhost', '-parallel', '-qopt-report=5'], + extra_compile_args=['-g', '-xhost', '-qopt-report=5'], extra_link_args=['-g', '-lsvml']) setup(name = 'Gaussian mixtures -- Fourier transform', From 61ff9a0e3db5f419e7ce4caa4db881a7fd296e72 Mon Sep 17 00:00:00 2001 From: Chris Hogan Date: Mon, 23 Oct 2017 09:13:51 -0500 Subject: [PATCH 09/43] Intel fft patch --- tractor/galaxy.py | 16 ++++++++++++++++ tractor/psf.py | 17 +++++++++++++++++ tractor/psfex.py | 16 ++++++++++++++++ 3 files changed, 49 insertions(+) diff --git a/tractor/galaxy.py b/tractor/galaxy.py index 6dbe9690..a01611ef 100644 --- a/tractor/galaxy.py +++ b/tractor/galaxy.py @@ -16,6 +16,22 @@ from __future__ import division import numpy as np +import numpy.fft_intel.libifft as m +def irfftn_numpy(x, s=None, axes=None): + a = np.asarray(x) + no_trim = (s is None) and (axes is None) + s, axes = m._cook_nd_args(a, s, axes, invreal=True) + la = axes[-1] + ovr_x = False + if len(s) > 1: + if not no_trim: + a = m._fix_dimensions(a, s, axes) + for ii in range(len(axes)-1): + a = m.ifft(a, s[ii], axes[ii], overwrite_x=ovr_x) + ovr_x = True + a = m.irfft_numpy(a, n = s[-1], axis=la) + return a +m.irfftn_numpy = irfftn_numpy from astrometry.util.miscutils import get_overlapping_region diff --git a/tractor/psf.py b/tractor/psf.py index c0675165..238e53ee 100644 --- a/tractor/psf.py +++ b/tractor/psf.py @@ -2,6 +2,23 @@ from __future__ import division import numpy as np +import numpy.fft_intel.libifft as m +def irfftn_numpy(x, s=None, axes=None): + a = np.asarray(x) + no_trim = (s is None) and (axes is None) + s, axes = m._cook_nd_args(a, s, axes, invreal=True) + la = axes[-1] + ovr_x = False + if len(s) > 1: + if not no_trim: + a = m._fix_dimensions(a, s, axes) + for ii in range(len(axes)-1): + a = m.ifft(a, s[ii], axes[ii], overwrite_x=ovr_x) + ovr_x = True + a = m.irfft_numpy(a, n = s[-1], axis=la) + return a +m.irfftn_numpy = irfftn_numpy + from astrometry.util.miscutils import lanczos_filter from tractor.image import Image diff --git a/tractor/psfex.py b/tractor/psfex.py index 4d53c400..96e2c262 100644 --- a/tractor/psfex.py +++ b/tractor/psfex.py @@ -1,6 +1,22 @@ from __future__ import print_function import numpy as np +import numpy.fft_intel.libifft as m +def irfftn_numpy(x, s=None, axes=None): + a = np.asarray(x) + no_trim = (s is None) and (axes is None) + s, axes = m._cook_nd_args(a, s, axes, invreal=True) + la = axes[-1] + ovr_x = False + if len(s) > 1: + if not no_trim: + a = m._fix_dimensions(a, s, axes) + for ii in range(len(axes)-1): + a = m.ifft(a, s[ii], axes[ii], overwrite_x=ovr_x) + ovr_x = True + a = m.irfft_numpy(a, n = s[-1], axis=la) + return a +m.irfftn_numpy = irfftn_numpy from tractor.utils import MultiParams, getClassName from tractor.psf import GaussianMixturePSF, PixelizedPSF From 65018341ae24950c4c41ada1daa51b31192a29ae Mon Sep 17 00:00:00 2001 From: Dustin Lang Date: Wed, 25 Oct 2017 11:48:13 -0700 Subject: [PATCH 10/43] add swigged correlate and correlate7 functions --- tractor/c_mp_fourier.i | 262 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 262 insertions(+) diff --git a/tractor/c_mp_fourier.i b/tractor/c_mp_fourier.i index c39d6698..74297af0 100755 --- a/tractor/c_mp_fourier.i +++ b/tractor/c_mp_fourier.i @@ -6,6 +6,9 @@ #define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION #include #include + + //#include + //#include %} %include "numpy.i" @@ -30,8 +33,267 @@ (double *out, int out_dim1, int out_dim2) }; + +%apply (double* INPLACE_ARRAY2, int DIM1, int DIM2) { + (double *work, int work_dim1, int work_dim2) +}; +%apply (double* INPLACE_ARRAY2, int DIM1, int DIM2) { + (double *img, int img_dim1, int img_dim2) +}; +%apply (double* INPLACE_ARRAY1, int DIM1) { + (double *filtx, int filtx_dim) +}; +%apply (double* INPLACE_ARRAY1, int DIM1) { + (double *filty, int filty_dim) +}; + + %inline %{ +#if 0 + } +#endif + +static void correlate(double* img, int img_dim1, int img_dim2, + double* filtx, int filtx_dim, + double* filty, int filty_dim, + double* work, int work_dim1, int work_dim2, + double* out, int out_dim1, int out_dim2) { + + //printf("img %i x %i, filter x %i, filter y %i\n", img_dim1, img_dim2, filtx_dim, filty_dim); + + __assume_aligned(img, 64); + __assume_aligned(work, 64); + __assume_aligned(out, 64); + + assert(filtx_dim <= 8); + assert(filty_dim <= 8); + + assert(img_dim1 == out_dim1); + assert(img_dim2 == out_dim2); + assert(img_dim1 == work_dim1); + assert(img_dim2 == work_dim2); + + double filter[8]; + int i, j, k; + int W = img_dim2; + int H = img_dim1; + + assert(W > 8); + assert(H > 8); + + /* + for (i=0; i 8); + assert(H > 8); + + memcpy(filter, filtx, 7 * sizeof(double)); + + // first run the filtx over image rows + for (j=0; j Date: Wed, 25 Oct 2017 12:41:30 -0700 Subject: [PATCH 11/43] use correlate7 --- tractor/c_mp_fourier.i | 8 ++++---- tractor/galaxy.py | 26 +++++++++++++++++++++----- 2 files changed, 25 insertions(+), 9 deletions(-) diff --git a/tractor/c_mp_fourier.i b/tractor/c_mp_fourier.i index 74297af0..8688e7ec 100755 --- a/tractor/c_mp_fourier.i +++ b/tractor/c_mp_fourier.i @@ -71,8 +71,8 @@ static void correlate(double* img, int img_dim1, int img_dim2, assert(img_dim1 == out_dim1); assert(img_dim2 == out_dim2); - assert(img_dim1 == work_dim1); - assert(img_dim2 == work_dim2); + assert(work_dim1 >= img_dim1); + assert(work_dim2 >= img_dim2); double filter[8]; int i, j, k; @@ -199,8 +199,8 @@ static void correlate7(double* img, int img_dim1, int img_dim2, assert(img_dim1 == out_dim1); assert(img_dim2 == out_dim2); - assert(img_dim1 == work_dim1); - assert(img_dim2 == work_dim2); + assert(work_dim1 >= img_dim1); + assert(work_dim2 >= img_dim2); double filter[7]; int i, j, k; diff --git a/tractor/galaxy.py b/tractor/galaxy.py index a01611ef..a9000925 100644 --- a/tractor/galaxy.py +++ b/tractor/galaxy.py @@ -485,7 +485,7 @@ def run_mog(amix=None, mm=None): # Lanczos-3 interpolation in ~ the same way we do for # pixelized PSFs. from astrometry.util.miscutils import lanczos_filter - from scipy.ndimage.filters import correlate1d + #from scipy.ndimage.filters import correlate1d #L = 3 L = fft_lanczos_order Lx = lanczos_filter(L, np.arange(-L, L+1) + mux) @@ -496,14 +496,25 @@ def run_mog(amix=None, mm=None): #print('Lx centroid', np.sum(Lx * (np.arange(-L,L+1)))) #print('Ly centroid', np.sum(Ly * (np.arange(-L,L+1)))) - #print('kernels:', Lx, Ly) + assert(len(Lx) == 7) + assert(len(Ly) == 7) + #cx = correlate1d(G, Lx, axis=1, mode='constant') + #G = correlate1d(cx, Ly, axis=0, mode='constant') + #del cx - cx = correlate1d(G, Lx, axis=1, mode='constant') - G = correlate1d(cx, Ly, axis=0, mode='constant') - del cx + G = np.require(G, requirements=['A']) + + np.save('G.npy', G) + np.save('Lx.npy', Lx) + np.save('Ly.npy', Ly) + sys.exit(0) + + from tractor.c_mp_fourier import correlate7 + correlate7(G, Lx, Ly, work_corr7, G) else: G = np.zeros((pH,pW), np.float32) + assert(False) if modelMask is not None: gh,gw = G.shape @@ -553,6 +564,11 @@ def run_mog(amix=None, mm=None): return Patch(ix0, iy0, G) + +work_corr7 = np.zeros((1024,1024), np.float64) +work_corr7 = np.require(work_corr7, requirements=['A']) + + fft_lanczos_order = 3 def _fourier_galaxy_debug_plots(G, shG, xi,yi,xo,yo, P, Fsum, From 7fdf71b5ae5901057b17a636cca57f99da722f38 Mon Sep 17 00:00:00 2001 From: Dustin Lang Date: Wed, 25 Oct 2017 12:42:03 -0700 Subject: [PATCH 12/43] compiler args --- tractor/galaxy.py | 4 ---- tractor/psf.py | 32 ++++++++++++++++---------------- tractor/setup-mpf.py | 2 +- 3 files changed, 17 insertions(+), 21 deletions(-) diff --git a/tractor/galaxy.py b/tractor/galaxy.py index a9000925..5513b621 100644 --- a/tractor/galaxy.py +++ b/tractor/galaxy.py @@ -410,10 +410,6 @@ def run_mog(amix=None, mm=None): # print('Tim shape:', img.shape) P,(cx,cy),(pH,pW),(v,w) = psf.getFourierTransform(px, py, halfsize) - preal = np.fft.irfft2(P, s=(pH,pW)) - # print('Sum of inverse-Fourier-transformed PSF model:', preal.sum()) - - dx = px - cx dy = py - cy if modelMask is not None: diff --git a/tractor/psf.py b/tractor/psf.py index 238e53ee..5398a1d4 100644 --- a/tractor/psf.py +++ b/tractor/psf.py @@ -2,22 +2,22 @@ from __future__ import division import numpy as np -import numpy.fft_intel.libifft as m -def irfftn_numpy(x, s=None, axes=None): - a = np.asarray(x) - no_trim = (s is None) and (axes is None) - s, axes = m._cook_nd_args(a, s, axes, invreal=True) - la = axes[-1] - ovr_x = False - if len(s) > 1: - if not no_trim: - a = m._fix_dimensions(a, s, axes) - for ii in range(len(axes)-1): - a = m.ifft(a, s[ii], axes[ii], overwrite_x=ovr_x) - ovr_x = True - a = m.irfft_numpy(a, n = s[-1], axis=la) - return a -m.irfftn_numpy = irfftn_numpy +# import numpy.fft_intel.libifft as m +# def irfftn_numpy(x, s=None, axes=None): +# a = np.asarray(x) +# no_trim = (s is None) and (axes is None) +# s, axes = m._cook_nd_args(a, s, axes, invreal=True) +# la = axes[-1] +# ovr_x = False +# if len(s) > 1: +# if not no_trim: +# a = m._fix_dimensions(a, s, axes) +# for ii in range(len(axes)-1): +# a = m.ifft(a, s[ii], axes[ii], overwrite_x=ovr_x) +# ovr_x = True +# a = m.irfft_numpy(a, n = s[-1], axis=la) +# return a +# m.irfftn_numpy = irfftn_numpy from astrometry.util.miscutils import lanczos_filter diff --git a/tractor/setup-mpf.py b/tractor/setup-mpf.py index 8df6a558..59212254 100755 --- a/tractor/setup-mpf.py +++ b/tractor/setup-mpf.py @@ -12,7 +12,7 @@ my_swig_mod = Extension('_c_mp_fourier', sources=['c_mp_fourier.i'], include_dirs=numpy_inc, - extra_compile_args=['-g', '-xhost', '-qopt-report=5'], + extra_compile_args=['-g', '-xhost', '-qopt-report=5', '-axMIC-AVX512'], #,AVX'], extra_link_args=['-g', '-lsvml']) setup(name = 'Gaussian mixtures -- Fourier transform', From 1a707efac6224f71da71239bca0331ed1981a4aa Mon Sep 17 00:00:00 2001 From: Dustin Lang Date: Wed, 25 Oct 2017 12:45:10 -0700 Subject: [PATCH 13/43] remove fix for intel np fft --- tractor/galaxy.py | 32 ++++++++++++++++---------------- tractor/psfex.py | 32 ++++++++++++++++---------------- 2 files changed, 32 insertions(+), 32 deletions(-) diff --git a/tractor/galaxy.py b/tractor/galaxy.py index 5513b621..c349d06d 100644 --- a/tractor/galaxy.py +++ b/tractor/galaxy.py @@ -16,22 +16,22 @@ from __future__ import division import numpy as np -import numpy.fft_intel.libifft as m -def irfftn_numpy(x, s=None, axes=None): - a = np.asarray(x) - no_trim = (s is None) and (axes is None) - s, axes = m._cook_nd_args(a, s, axes, invreal=True) - la = axes[-1] - ovr_x = False - if len(s) > 1: - if not no_trim: - a = m._fix_dimensions(a, s, axes) - for ii in range(len(axes)-1): - a = m.ifft(a, s[ii], axes[ii], overwrite_x=ovr_x) - ovr_x = True - a = m.irfft_numpy(a, n = s[-1], axis=la) - return a -m.irfftn_numpy = irfftn_numpy +# import numpy.fft_intel.libifft as m +# def irfftn_numpy(x, s=None, axes=None): +# a = np.asarray(x) +# no_trim = (s is None) and (axes is None) +# s, axes = m._cook_nd_args(a, s, axes, invreal=True) +# la = axes[-1] +# ovr_x = False +# if len(s) > 1: +# if not no_trim: +# a = m._fix_dimensions(a, s, axes) +# for ii in range(len(axes)-1): +# a = m.ifft(a, s[ii], axes[ii], overwrite_x=ovr_x) +# ovr_x = True +# a = m.irfft_numpy(a, n = s[-1], axis=la) +# return a +# m.irfftn_numpy = irfftn_numpy from astrometry.util.miscutils import get_overlapping_region diff --git a/tractor/psfex.py b/tractor/psfex.py index 96e2c262..a9d97b6c 100644 --- a/tractor/psfex.py +++ b/tractor/psfex.py @@ -1,22 +1,22 @@ from __future__ import print_function import numpy as np -import numpy.fft_intel.libifft as m -def irfftn_numpy(x, s=None, axes=None): - a = np.asarray(x) - no_trim = (s is None) and (axes is None) - s, axes = m._cook_nd_args(a, s, axes, invreal=True) - la = axes[-1] - ovr_x = False - if len(s) > 1: - if not no_trim: - a = m._fix_dimensions(a, s, axes) - for ii in range(len(axes)-1): - a = m.ifft(a, s[ii], axes[ii], overwrite_x=ovr_x) - ovr_x = True - a = m.irfft_numpy(a, n = s[-1], axis=la) - return a -m.irfftn_numpy = irfftn_numpy +# import numpy.fft_intel.libifft as m +# def irfftn_numpy(x, s=None, axes=None): +# a = np.asarray(x) +# no_trim = (s is None) and (axes is None) +# s, axes = m._cook_nd_args(a, s, axes, invreal=True) +# la = axes[-1] +# ovr_x = False +# if len(s) > 1: +# if not no_trim: +# a = m._fix_dimensions(a, s, axes) +# for ii in range(len(axes)-1): +# a = m.ifft(a, s[ii], axes[ii], overwrite_x=ovr_x) +# ovr_x = True +# a = m.irfft_numpy(a, n = s[-1], axis=la) +# return a +# m.irfftn_numpy = irfftn_numpy from tractor.utils import MultiParams, getClassName from tractor.psf import GaussianMixturePSF, PixelizedPSF From 68617084356aec7dfdbf24f3637df0529eaadcaf Mon Sep 17 00:00:00 2001 From: Dustin Lang Date: Wed, 1 Nov 2017 13:06:09 -0700 Subject: [PATCH 14/43] make correlate7 use the 'img' as input-output array; make PSF use correlate7 also --- tractor/c_mp_fourier.i | 17 +++++++++-------- tractor/galaxy.py | 4 ++-- tractor/psf.py | 30 ++++++++++++++++++++++++------ 3 files changed, 35 insertions(+), 16 deletions(-) diff --git a/tractor/c_mp_fourier.i b/tractor/c_mp_fourier.i index 8688e7ec..d3c0bda7 100755 --- a/tractor/c_mp_fourier.i +++ b/tractor/c_mp_fourier.i @@ -40,10 +40,10 @@ %apply (double* INPLACE_ARRAY2, int DIM1, int DIM2) { (double *img, int img_dim1, int img_dim2) }; -%apply (double* INPLACE_ARRAY1, int DIM1) { +%apply (double* IN_ARRAY1, int DIM1) { (double *filtx, int filtx_dim) }; -%apply (double* INPLACE_ARRAY1, int DIM1) { +%apply (double* IN_ARRAY1, int DIM1) { (double *filty, int filty_dim) }; @@ -187,12 +187,11 @@ static void correlate(double* img, int img_dim1, int img_dim2, static void correlate7(double* img, int img_dim1, int img_dim2, double* filtx, int filtx_dim, double* filty, int filty_dim, - double* work, int work_dim1, int work_dim2, - double* out, int out_dim1, int out_dim2) { + double* work, int work_dim1, int work_dim2) { + // Output goes back into "img"! __assume_aligned(img, 64); __assume_aligned(work, 64); - __assume_aligned(out, 64); assert(filtx_dim == 7); assert(filty_dim == 7); @@ -259,6 +258,8 @@ static void correlate7(double* img, int img_dim1, int img_dim2, int workH = W; int workW = H; + // Output goes back into "img"! + // Now run filty over rows of the 'work' array for (j=0; j Date: Wed, 1 Nov 2017 13:06:36 -0700 Subject: [PATCH 15/43] remove debugging --- tractor/galaxy.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/tractor/galaxy.py b/tractor/galaxy.py index 849e0e73..6bd3df29 100644 --- a/tractor/galaxy.py +++ b/tractor/galaxy.py @@ -502,10 +502,10 @@ def run_mog(amix=None, mm=None): G = np.require(G, requirements=['A']) correlate7(G, Lx, Ly, work_corr7) - np.save('G.npy', G) - np.save('Lx.npy', Lx) - np.save('Ly.npy', Ly) - sys.exit(0) + # np.save('G.npy', G) + # np.save('Lx.npy', Lx) + # np.save('Ly.npy', Ly) + # sys.exit(0) else: From 8eca064083f4bd7dea10e423e55e2521e77f6a4a Mon Sep 17 00:00:00 2001 From: Dustin Lang Date: Thu, 2 Nov 2017 08:24:40 -0700 Subject: [PATCH 16/43] add float version of correlate7; compute integer offsets and use "restrict" keyword to make compiler optimizer happy --- tractor/c_mp_fourier.i | 159 ++++++++++++++++++++++++++++++++++------- tractor/psf.py | 23 ++++-- 2 files changed, 153 insertions(+), 29 deletions(-) diff --git a/tractor/c_mp_fourier.i b/tractor/c_mp_fourier.i index d3c0bda7..4374685b 100755 --- a/tractor/c_mp_fourier.i +++ b/tractor/c_mp_fourier.i @@ -48,6 +48,21 @@ }; +%apply (float* INPLACE_ARRAY2, int DIM1, int DIM2) { + (float *work, int work_dim1, int work_dim2) +}; +%apply (float* INPLACE_ARRAY2, int DIM1, int DIM2) { + (float *img, int img_dim1, int img_dim2) +}; +/* +%apply (float* IN_ARRAY1, int DIM1) { + (float *filtx, int filtx_dim) +}; +%apply (float* IN_ARRAY1, int DIM1) { + (float *filty, int filty_dim) +}; +*/ + %inline %{ #if 0 @@ -184,10 +199,10 @@ static void correlate(double* img, int img_dim1, int img_dim2, -static void correlate7(double* img, int img_dim1, int img_dim2, - double* filtx, int filtx_dim, - double* filty, int filty_dim, - double* work, int work_dim1, int work_dim2) { +static void correlate7(double* restrict img, int img_dim1, int img_dim2, + double* restrict filtx, int filtx_dim, + double* restrict filty, int filty_dim, + double* restrict work, int work_dim1, int work_dim2) { // Output goes back into "img"! __assume_aligned(img, 64); @@ -214,46 +229,143 @@ static void correlate7(double* img, int img_dim1, int img_dim2, // first run the filtx over image rows for (j=0; j= img_dim1); + assert(work_dim2 >= img_dim2); + + float filter[7]; + int i, j, k; + int W = img_dim2; + int H = img_dim1; + + assert(W > 8); + assert(H > 8); + + //memcpy(filter, filtx, 7 * sizeof(float)); + // double -> float + for (i=0; i<7; i++) + filter[i] = filtx[i]; + + // first run the filtx over image rows for (j=0; j float + for (i=0; i<7; i++) + filter[i] = filty[i]; int workH = W; int workW = H; @@ -263,16 +375,16 @@ static void correlate7(double* img, int img_dim1, int img_dim2, // Now run filty over rows of the 'work' array for (j=0; j Date: Thu, 2 Nov 2017 09:03:39 -0700 Subject: [PATCH 17/43] correlate7: change loop iteration from i=3... to i=0... --- tractor/c_mp_fourier.i | 79 +++++++++++++++++++++++++----------------- 1 file changed, 47 insertions(+), 32 deletions(-) diff --git a/tractor/c_mp_fourier.i b/tractor/c_mp_fourier.i index 4374685b..7a9ab8c6 100755 --- a/tractor/c_mp_fourier.i +++ b/tractor/c_mp_fourier.i @@ -239,17 +239,19 @@ static void correlate7(double* restrict img, int img_dim1, int img_dim2, work[i*H + j] = sum; } // middle section - offset = j*W - 3; - for (i=3; i<=(W-4); i++) { + //offset = j*W - 3; + //for (i=3; i<=(W-4); i++) { + for (i=0; i<=(W-7); i++) { double sum = 0.0; for (k=0; k<7; k++) //sum += filter[k] * img_row[i-3+k]; sum += filter[k] * img[offset + i+k]; - work[i*H + j] = sum; + //work[i*H + j] = sum; + work[(i+3)*H + j] = sum; } // special handling of right edge // i=0 is the rightmost pixel - offset = j*W; + //offset = j*W; for (i=0; i<3; i++) { double sum = 0.0; for (k=0; k<4+i; k++) @@ -279,17 +281,19 @@ static void correlate7(double* restrict img, int img_dim1, int img_dim2, img[i*W + j] = sum; } // middle section - offset = j * workW - 3; - for (i=3; i<=(workW-4); i++) { + //offset = j * workW - 3; + //for (i=3; i<=(workW-4); i++) { + for (i=0; i<=(workW-7); i++) { double sum = 0.0; for (k=0; k<7; k++) //sum += filter[k] * work_row[i-3+k]; sum += filter[k] * work[offset + i+k]; - img[i*W + j] = sum; + //img[i*W + j] = sum; + img[(i+3)*W + j] = sum; } // special handling of right edge // i=0 is the rightmost pixel - offset = j * workW; + //offset = j * workW; for (i=0; i<3; i++) { double sum = 0.0; for (k=0; k<4+i; k++) @@ -303,12 +307,11 @@ static void correlate7(double* restrict img, int img_dim1, int img_dim2, -static void correlate7f(float* img, int img_dim1, int img_dim2, - //float* filtx, int filtx_dim, - //float* filty, int filty_dim, - double* filtx, int filtx_dim, - double* filty, int filty_dim, - float* work, int work_dim1, int work_dim2) { + +static void correlate7f(float* restrict img, int img_dim1, int img_dim2, + double* restrict filtx, int filtx_dim, + double* restrict filty, int filty_dim, + float* restrict work, int work_dim1, int work_dim2) { // Output goes back into "img"! __assume_aligned(img, 64); @@ -330,40 +333,46 @@ static void correlate7f(float* img, int img_dim1, int img_dim2, assert(W > 8); assert(H > 8); - //memcpy(filter, filtx, 7 * sizeof(float)); - // double -> float + //memcpy(filter, filtx, 7 * sizeof(double)); for (i=0; i<7; i++) filter[i] = filtx[i]; // first run the filtx over image rows for (j=0; j float + //memcpy(filter, filty, 7 * sizeof(double)); for (i=0; i<7; i++) filter[i] = filty[i]; @@ -375,26 +384,34 @@ static void correlate7f(float* img, int img_dim1, int img_dim2, // Now run filty over rows of the 'work' array for (j=0; j Date: Thu, 2 Nov 2017 09:04:29 -0700 Subject: [PATCH 18/43] convert PsfEx values to float32 upon read; align PixelizedPsf image --- tractor/psf.py | 16 ++++++++++------ tractor/psfex.py | 2 ++ 2 files changed, 12 insertions(+), 6 deletions(-) diff --git a/tractor/psf.py b/tractor/psf.py index 5468712a..18a03beb 100644 --- a/tractor/psf.py +++ b/tractor/psf.py @@ -53,7 +53,8 @@ def __init__(self, img, Lorder=3): - *Lorder* is the order of the Lanczos interpolant used for shifting the image to subpixel positions. ''' - self.img = img + # align + self.img = np.require(img, requirements=['A']) H,W = img.shape assert((H % 2) == 1) assert((W % 2) == 1) @@ -132,20 +133,21 @@ def getPointSourcePatch(self, px, py, minval=0., modelMask=None, **kwargs): Lx /= Lx.sum() Ly /= Ly.sum() - print('img', img.dtype) - print('Lx', Lx.dtype) - print('Ly', Ly.dtype) + # print('img', img.dtype) + # print('Lx', Lx.dtype) + # print('Ly', Ly.dtype) if modelMask is None: assert(len(Lx) == 7) assert(len(Ly) == 7) - img = np.require(img, requirements=['A']) from tractor.c_mp_fourier import correlate7, correlate7f + #print('modelMask is None, dtype', img.dtype) if img.dtype == np.float32: correlate7f(img, Lx, Ly, work_corr7f) else: correlate7(img, Lx, Ly, work_corr7) + #print('modelMask is None, dtype', img.dtype, 'done') # from tractor.c_mp_fourier import correlate7 # correlate7(img, Lx, Ly, work_corr7) return Patch(x0, y0, img) @@ -165,12 +167,14 @@ def getPointSourcePatch(self, px, py, minval=0., modelMask=None, **kwargs): assert(len(Lx) == 7) assert(len(Ly) == 7) - mm = np.require(mm, requirements=['A']) + #mm = np.require(mm, requirements=['A']) from tractor.c_mp_fourier import correlate7, correlate7f + #print('mm, dtype', mm.dtype) if mm.dtype == np.float32: correlate7f(mm, Lx, Ly, work_corr7f) else: correlate7(mm, Lx, Ly, work_corr7) + #print('mm, dtype', mm.dtype, 'done') #sx = correlate1d(mm, Lx, axis=1, mode='constant') #mm = correlate1d(sx, Ly, axis=0, mode='constant') diff --git a/tractor/psfex.py b/tractor/psfex.py index a9d97b6c..aa8994e2 100644 --- a/tractor/psfex.py +++ b/tractor/psfex.py @@ -185,6 +185,7 @@ def __init__(self, fn=None, ext=1, Ti=None): from astrometry.util.fits import fits_table T = fits_table(fn, ext=ext) ims = T.psf_mask[0] + ims = ims.astype(np.float32) #print('Got', ims.shape, 'PSF images') hdr = T.get_header() @@ -240,6 +241,7 @@ def __init__(self, fn=None, ext=1, Ti=None): ne = (degree + 1) * (degree + 2) / 2 assert(Ti.psfaxis3 == ne) ims = Ti.psf_mask + ims = ims.astype(np.float32) assert(len(ims.shape) == 3) assert(ims.shape[0] == ne) self.psfbases = ims From 8aa4c48e6bcf1ef3e92063bd32b7428af509979b Mon Sep 17 00:00:00 2001 From: Dustin Lang Date: Thu, 2 Nov 2017 09:06:23 -0700 Subject: [PATCH 19/43] tidy --- tractor/c_mp_fourier.i | 48 ------------------------------------------ 1 file changed, 48 deletions(-) diff --git a/tractor/c_mp_fourier.i b/tractor/c_mp_fourier.i index 7a9ab8c6..1fc11604 100755 --- a/tractor/c_mp_fourier.i +++ b/tractor/c_mp_fourier.i @@ -6,9 +6,6 @@ #define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION #include #include - - //#include - //#include %} %include "numpy.i" @@ -54,14 +51,6 @@ %apply (float* INPLACE_ARRAY2, int DIM1, int DIM2) { (float *img, int img_dim1, int img_dim2) }; -/* -%apply (float* IN_ARRAY1, int DIM1) { - (float *filtx, int filtx_dim) -}; -%apply (float* IN_ARRAY1, int DIM1) { - (float *filty, int filty_dim) -}; -*/ %inline %{ @@ -229,33 +218,25 @@ static void correlate7(double* restrict img, int img_dim1, int img_dim2, // first run the filtx over image rows for (j=0; j Date: Thu, 2 Nov 2017 09:33:04 -0700 Subject: [PATCH 20/43] remove assertion --- tractor/galaxy.py | 11 ----------- 1 file changed, 11 deletions(-) diff --git a/tractor/galaxy.py b/tractor/galaxy.py index 6bd3df29..5866850e 100644 --- a/tractor/galaxy.py +++ b/tractor/galaxy.py @@ -502,15 +502,8 @@ def run_mog(amix=None, mm=None): G = np.require(G, requirements=['A']) correlate7(G, Lx, Ly, work_corr7) - # np.save('G.npy', G) - # np.save('Lx.npy', Lx) - # np.save('Ly.npy', Ly) - # sys.exit(0) - - else: G = np.zeros((pH,pW), np.float32) - assert(False) if modelMask is not None: gh,gw = G.shape @@ -522,10 +515,6 @@ def run_mog(amix=None, mm=None): shG = np.zeros((mh,mw), G.dtype) shG[yo,xo] = G[yi,xi] - # print('shift:', (sx,sy), 'mm', (mw,mh), 'g', (gw,gh)) - # print('yi,xi,', yi,xi) - # print('yo,xo,', yo,xo) - if debug_ps is not None: _fourier_galaxy_debug_plots(G, shG, xi,yi,xo,yo, P, Fsum, pW,pH, psf) From d3fe462babce0889457f88259b6d04d8ceabdfaa Mon Sep 17 00:00:00 2001 From: Dustin Lang Date: Fri, 3 Nov 2017 05:59:25 -0700 Subject: [PATCH 21/43] correlate7: create separate output array; ensure float32 --- tractor/psf.py | 26 ++++++++++++++++---------- 1 file changed, 16 insertions(+), 10 deletions(-) diff --git a/tractor/psf.py b/tractor/psf.py index 18a03beb..6afddf05 100644 --- a/tractor/psf.py +++ b/tractor/psf.py @@ -143,14 +143,19 @@ def getPointSourcePatch(self, px, py, minval=0., modelMask=None, **kwargs): assert(len(Ly) == 7) from tractor.c_mp_fourier import correlate7, correlate7f #print('modelMask is None, dtype', img.dtype) - if img.dtype == np.float32: - correlate7f(img, Lx, Ly, work_corr7f) - else: - correlate7(img, Lx, Ly, work_corr7) + + outimg = np.empty(img.shape, np.float32) + outimg[:,:] = img + correlate7f(outimg, Lx, Ly, work_corr7f) + return Patch(x0, y0, outimg) + # if img.dtype == np.float32: + # correlate7f(img, Lx, Ly, work_corr7f) + # else: + # correlate7(img, Lx, Ly, work_corr7) + #return Patch(x0, y0, img) #print('modelMask is None, dtype', img.dtype, 'done') # from tractor.c_mp_fourier import correlate7 # correlate7(img, Lx, Ly, work_corr7) - return Patch(x0, y0, img) # sx = correlate1d(img, Lx, axis=1, mode='constant') # shifted = correlate1d(sx, Ly, axis=0, mode='constant') # assert(np.all(np.isfinite(shifted))) @@ -159,7 +164,8 @@ def getPointSourcePatch(self, px, py, minval=0., modelMask=None, **kwargs): # padding = L # Create a modelMask + padding sized stamp and insert PSF image into it - mm = np.zeros((mh+2*padding, mw+2*padding), img.dtype) + #mm = np.zeros((mh+2*padding, mw+2*padding), img.dtype) + mm = np.zeros((mh+2*padding, mw+2*padding), np.float32) yi,yo = get_overlapping_region(my0-y0-padding, my0-y0+mh-1+padding, 0, H-1) xi,xo = get_overlapping_region(mx0-x0-padding, mx0-x0+mw-1+padding, 0, W-1) @@ -170,10 +176,10 @@ def getPointSourcePatch(self, px, py, minval=0., modelMask=None, **kwargs): #mm = np.require(mm, requirements=['A']) from tractor.c_mp_fourier import correlate7, correlate7f #print('mm, dtype', mm.dtype) - if mm.dtype == np.float32: - correlate7f(mm, Lx, Ly, work_corr7f) - else: - correlate7(mm, Lx, Ly, work_corr7) + #if mm.dtype == np.float32: + correlate7f(mm, Lx, Ly, work_corr7f) + #else: + # correlate7(mm, Lx, Ly, work_corr7) #print('mm, dtype', mm.dtype, 'done') #sx = correlate1d(mm, Lx, axis=1, mode='constant') From bf2d45bd4e836f61fd174587877a685a3e7d66b9 Mon Sep 17 00:00:00 2001 From: Dustin Lang Date: Fri, 3 Nov 2017 14:58:35 -0400 Subject: [PATCH 22/43] ceres: update 'gflags' pkg-config command; and cast a type for PyArray_TYPE --- tractor/Makefile | 4 ++-- tractor/ceres.i | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/tractor/Makefile b/tractor/Makefile index 145c19bb..df7ab815 100644 --- a/tractor/Makefile +++ b/tractor/Makefile @@ -58,8 +58,8 @@ BLAS_LIB ?= -lblas GLOG_LIB ?= $(shell pkg-config --libs libglog) #-lglog GLOG_INC ?= $(shell pkg-config --cflags libglog) -GFLAGS_LIB ?= $(shell pkg-config --libs libgflags) -GFLAGS_INC ?= $(shell pkg-config --cflags libgflags) +GFLAGS_LIB ?= $(shell pkg-config --libs gflags) +GFLAGS_INC ?= $(shell pkg-config --cflags gflags) CERES_LIB_DIR ?= /usr/local/lib diff --git a/tractor/ceres.i b/tractor/ceres.i index d377bfad..9019b8f9 100644 --- a/tractor/ceres.i +++ b/tractor/ceres.i @@ -57,7 +57,7 @@ static PyObject* real_ceres_forced_phot(PyObject* blocks, Nblocks = PyList_Size(blocks); //printf("N blocks: %i\n", (int)Nblocks); assert(PyArray_Check(np_fluxes)); - assert(PyArray_TYPE(np_fluxes) == NPY_DOUBLE); + assert(PyArray_TYPE((PyArrayObject*)np_fluxes) == NPY_DOUBLE); int Nfluxes = (int)PyArray_Size(np_fluxes); double* realfluxes = (double*)PyArray_DATA((PyArrayObject*)np_fluxes); T* mod0data; From b421de4468ac17c1d0fc1e40e2fc717556795b8c Mon Sep 17 00:00:00 2001 From: Dustin Lang Date: Mon, 6 Nov 2017 13:18:04 -0800 Subject: [PATCH 23/43] merge c_mp_fourier into mp_fourier --- tractor/mp_fourier.i | 386 ++++++++++++++++++++++++++++--------------- 1 file changed, 252 insertions(+), 134 deletions(-) diff --git a/tractor/mp_fourier.i b/tractor/mp_fourier.i index a496e719..94233f90 100644 --- a/tractor/mp_fourier.i +++ b/tractor/mp_fourier.i @@ -1,17 +1,58 @@ %module(package="tractor") mp_fourier %{ +#define SWIG_FILE_WITH_INIT #define _GNU_SOURCE 1 #define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION #include #include #include - %} +%} + +%include "numpy.i" %init %{ // numpy import_array(); - %} +%} + +%apply (double* IN_ARRAY1, int DIM1) { + (double *amps, int amps_len), + (double *v, int v_len), + (double *w, int w_len) +}; +%apply (double* IN_ARRAY2, int DIM1, int DIM2) { + (double *means, int means_dim1, int means_dim2) +}; +%apply (double* IN_ARRAY3, int DIM1, int DIM2, int DIM3) { + (double *vars, int vars_dim1, int vars_dim2, int vars_dim3) +}; +%apply (double* INPLACE_ARRAY2, int DIM1, int DIM2) { + (double *out, int out_dim1, int out_dim2) +}; + + +%apply (double* INPLACE_ARRAY2, int DIM1, int DIM2) { + (double *work, int work_dim1, int work_dim2) +}; +%apply (double* INPLACE_ARRAY2, int DIM1, int DIM2) { + (double *img, int img_dim1, int img_dim2) +}; +%apply (double* IN_ARRAY1, int DIM1) { + (double *filtx, int filtx_dim) +}; +%apply (double* IN_ARRAY1, int DIM1) { + (double *filty, int filty_dim) +}; + + +%apply (float* INPLACE_ARRAY2, int DIM1, int DIM2) { + (float *work, int work_dim1, int work_dim2) +}; +%apply (float* INPLACE_ARRAY2, int DIM1, int DIM2) { + (float *img, int img_dim1, int img_dim2) +}; + %inline %{ @@ -19,150 +60,227 @@ } // fool emacs indenter #endif -static PyObject* mixture_profile_fourier_transform( - PyObject* po_amps, - PyObject* po_means, - PyObject* po_vars, - PyObject* po_v, - PyObject* po_w - ) { - - npy_intp K, NW,NV; - const npy_intp D = 2; - npy_intp i,j,k; - double* amps, *means, *vars, *vv, *ww; - PyObject* np_F; - double* f; - npy_intp dims[2]; - - PyArrayObject *np_amps, *np_means, *np_vars, *np_v, *np_w; - - if (!PyArray_Check(po_amps) || !PyArray_Check(po_means) || - !PyArray_Check(po_vars) || !PyArray_Check(po_v) || - !PyArray_Check(po_w)) { - PyErr_SetString(PyExc_ValueError, "Expected numpy arrays"); - return NULL; - } +static void correlate7(double* restrict img, int img_dim1, int img_dim2, + double* restrict filtx, int filtx_dim, + double* restrict filty, int filty_dim, + double* restrict work, int work_dim1, int work_dim2) { + // Output goes back into "img"! - np_amps = (PyArrayObject*)po_amps; - np_means = (PyArrayObject*)po_means; - np_vars = (PyArrayObject*)po_vars; - np_v = (PyArrayObject*)po_v; - np_w = (PyArrayObject*)po_w; - - if ((PyArray_TYPE(np_amps) != NPY_DOUBLE) || - (PyArray_TYPE(np_means ) != NPY_DOUBLE) || - (PyArray_TYPE(np_vars) != NPY_DOUBLE) || - (PyArray_TYPE(np_v) != NPY_DOUBLE) || - (PyArray_TYPE(np_w) != NPY_DOUBLE)) { - PyErr_SetString(PyExc_ValueError, "Expected numpy double arrays"); - return NULL; - } + __assume_aligned(img, 64); + __assume_aligned(work, 64); - if (PyArray_NDIM(np_amps) != 1) { - PyErr_SetString(PyExc_ValueError, "Expected 'amps' to be 1-d"); - return NULL; - } - K = PyArray_DIM(np_amps, 0); - if (PyArray_NDIM(np_means) != 2) { - PyErr_SetString(PyExc_ValueError, "Expected 'means' to be 2-d"); - return NULL; - } - if ((PyArray_DIM(np_means, 0) != K) || - (PyArray_DIM(np_means, 1) != D)) { - PyErr_SetString(PyExc_ValueError, "Expected 'means' to be K x D"); - return NULL; - } - if (PyArray_NDIM(np_vars) != 3) { - PyErr_SetString(PyExc_ValueError, "Expected 'vars' to be 3-d"); - return NULL; - } - if ((PyArray_DIM(np_vars, 0) != K) || - (PyArray_DIM(np_vars, 1) != D) || - (PyArray_DIM(np_vars, 2) != D)) { - PyErr_SetString(PyExc_ValueError, "Expected 'vars' to be K x D x D"); - return NULL; + assert(filtx_dim == 7); + assert(filty_dim == 7); + + assert(img_dim1 == out_dim1); + assert(img_dim2 == out_dim2); + assert(work_dim1 >= img_dim1); + assert(work_dim2 >= img_dim2); + + double filter[7]; + int i, j, k; + int W = img_dim2; + int H = img_dim1; + + assert(W > 8); + assert(H > 8); + + memcpy(filter, filtx, 7 * sizeof(double)); + + // first run the filtx over image rows + for (j=0; j= img_dim1); + assert(work_dim2 >= img_dim2); + + float filter[7]; + int i, j, k; + int W = img_dim2; + int H = img_dim1; + + assert(W > 8); + assert(H > 8); + + //memcpy(filter, filtx, 7 * sizeof(double)); + for (i=0; i<7; i++) + filter[i] = filtx[i]; + + // first run the filtx over image rows + for (j=0; j Date: Mon, 6 Nov 2017 13:26:03 -0800 Subject: [PATCH 24/43] modify setup-mpf.py to produce 'generic' and 'intel' versions of mp_fourier module --- tractor/mp_fourier.i | 4 ++++ tractor/setup-mpf.py | 26 +++++++++++++++----------- 2 files changed, 19 insertions(+), 11 deletions(-) diff --git a/tractor/mp_fourier.i b/tractor/mp_fourier.i index 94233f90..3cf1eb50 100644 --- a/tractor/mp_fourier.i +++ b/tractor/mp_fourier.i @@ -66,8 +66,10 @@ static void correlate7(double* restrict img, int img_dim1, int img_dim2, double* restrict work, int work_dim1, int work_dim2) { // Output goes back into "img"! +#ifdef __INTEL_COMPILER __assume_aligned(img, 64); __assume_aligned(work, 64); +#endif assert(filtx_dim == 7); assert(filty_dim == 7); @@ -156,8 +158,10 @@ static void correlate7f(float* restrict img, int img_dim1, int img_dim2, float* restrict work, int work_dim1, int work_dim2) { // Output goes back into "img"! +#ifdef __INTEL_COMPILER __assume_aligned(img, 64); __assume_aligned(work, 64); +#endif assert(filtx_dim == 7); assert(filty_dim == 7); diff --git a/tractor/setup-mpf.py b/tractor/setup-mpf.py index 59212254..7edd865d 100755 --- a/tractor/setup-mpf.py +++ b/tractor/setup-mpf.py @@ -1,19 +1,23 @@ from distutils.core import setup, Extension from numpy.distutils.misc_util import get_numpy_include_dirs +import os numpy_inc = get_numpy_include_dirs() -c_swig_module = Extension('_mp_fourier', - sources = ['mp_fourier.i' ], - include_dirs = numpy_inc, - extra_compile_args=['-g'], - extra_link_args=['-g'], +if os.environ.get('CC') == 'icc': + mpf_module = Extension('_intel_mp_fourier', + sources=['mp_fourier.i'], + include_dirs=numpy_inc, + extra_compile_args=['-g', '-xhost', '-qopt-report=5', '-axMIC-AVX512'], + extra_link_args=['-g', '-lsvml'] + ) +else: + mpf_module = Extension('_mp_fourier', + sources = ['mp_fourier.i' ], + include_dirs = numpy_inc, + extra_compile_args=['-g'], + extra_link_args=['-g'], ) -my_swig_mod = Extension('_c_mp_fourier', - sources=['c_mp_fourier.i'], - include_dirs=numpy_inc, - extra_compile_args=['-g', '-xhost', '-qopt-report=5', '-axMIC-AVX512'], #,AVX'], - extra_link_args=['-g', '-lsvml']) setup(name = 'Gaussian mixtures -- Fourier transform', version = '1.0', @@ -21,4 +25,4 @@ author = 'Lang & Hogg', author_email = 'dstndstn@gmail.com', url = 'http://astrometry.net', - ext_modules = [c_swig_module, my_swig_mod]) + ext_modules = [mpf_module]) From 7678892ad645bbcb8c05439fed71b38cf20cc82d Mon Sep 17 00:00:00 2001 From: Dustin Lang Date: Mon, 6 Nov 2017 13:58:21 -0800 Subject: [PATCH 25/43] makefile: build intel and generic versions of mp_fourier --- tractor/Makefile | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/tractor/Makefile b/tractor/Makefile index 853bcf1c..5d17f6b0 100644 --- a/tractor/Makefile +++ b/tractor/Makefile @@ -7,14 +7,14 @@ PYTHON_SO_EXT ?= $(shell $(PYTHON) -c "from distutils import sysconfig; print(sy ceres: _ceres$(PYTHON_SO_EXT) ceres.py .PHONY: ceres -mpf: mp_fourier c_mp_fourier +mpf: mp_fourier intel_mp_fourier .PHONY: mpf mp_fourier: _mp_fourier$(PYTHON_SO_EXT) mp_fourier.py .PHONY: mp_fourier -c_mp_fourier: _c_mp_fourier$(PYTHON_SO_EXT) c_mp_fourier.py -.PHONY: c_mp_fourier +intel_mp_fourier: _intel_mp_fourier$(PYTHON_SO_EXT) mp_fourier.py +.PHONY: intel_mp_fourier mix: _mix$(PYTHON_SO_EXT) mix.py .PHONY: mix @@ -25,8 +25,8 @@ emfit: _emfit$(PYTHON_SO_EXT) emfit.py _mp_fourier$(PYTHON_SO_EXT): mp_fourier.i setup-mpf.py $(PYTHON) setup-mpf.py build_ext --inplace -_c_mp_fourier$(PYTHON_SO_EXT): c_mp_fourier.i setup-mpf.py - $(PYTHON) setup-mpf.py build_ext --inplace +_intel_mp_fourier$(PYTHON_SO_EXT): mp_fourier.i setup-mpf.py + CC=icc $(PYTHON) setup-mpf.py build_ext --inplace mix.py _mix$(PYTHON_SO_EXT): mix.i approx3.c gauss_masked.c setup-mix.py $(PYTHON) setup-mix.py build_ext --inplace From 072261c864ddd4aa1cb6a26a13ec23b8296c20f5 Mon Sep 17 00:00:00 2001 From: Dustin Lang Date: Mon, 6 Nov 2017 13:58:55 -0800 Subject: [PATCH 26/43] start plugging in generic (non-zero-mean) gaussian_fourier, but then backpedal --- tractor/galaxy.py | 2 +- tractor/mixture_profiles.py | 12 +++--- tractor/mp_fourier.i | 81 +++++++++++++++++++++++++++++++------ 3 files changed, 76 insertions(+), 19 deletions(-) diff --git a/tractor/galaxy.py b/tractor/galaxy.py index 5866850e..66f7133b 100644 --- a/tractor/galaxy.py +++ b/tractor/galaxy.py @@ -470,7 +470,7 @@ def run_mog(amix=None, mm=None): if fftmix is not None: #print('fftmix; mux,muy=', mux,muy) - Fsum = fftmix.getFourierTransform(v, w) + Fsum = fftmix.getFourierTransform(v, w, zero_mean=True) # print('inverse Fourier-transforming into result size:', pH,pW) G = np.fft.irfft2(Fsum * P, s=(pH,pW)) diff --git a/tractor/mixture_profiles.py b/tractor/mixture_profiles.py index 675dbd55..38fc8190 100755 --- a/tractor/mixture_profiles.py +++ b/tractor/mixture_profiles.py @@ -144,7 +144,7 @@ def convolve(self, other): newk = nextnewk return MixtureOfGaussians(newamp, newmean, newvar) - def getFourierTransform(self, v, w, use_mp_fourier=True): + def getFourierTransform(self, v, w, use_mp_fourier=True, zero_mean=False): ''' v: FFT frequencies in the x direction w: FFT frequencies in the y direction @@ -156,12 +156,14 @@ def getFourierTransform(self, v, w, use_mp_fourier=True): F = np.fft.rfft2(img) v = np.fft.rfftfreq(W) w = np.fft.fftfreq(H) + + If zero_mean is *True*, ignore the *mean* of this mixture of Gaussians. ''' - if mp_fourier and use_mp_fourier: + if mp_fourier and use_mp_fourier and zero_mean: f = np.zeros((len(w), len(v))) - mp_fourier.mixture_profile_fourier_transform( - self.amp, self.mean, self.var, v, w, f) + mp_fourier.gaussian_fourier_transform_zero_mean( + self.amp, self.var, v, w, f) return f Fsum = None @@ -187,7 +189,7 @@ def getFourierTransform(self, v, w, use_mp_fourier=True): Fsum += amp * F return Fsum - + # ideally pos is a numpy array shape (N, self.D) # returns a numpy array shape (N) # may fail for self.D == 1 diff --git a/tractor/mp_fourier.i b/tractor/mp_fourier.i index 3cf1eb50..fe21bf34 100644 --- a/tractor/mp_fourier.i +++ b/tractor/mp_fourier.i @@ -247,23 +247,20 @@ static void correlate7f(float* restrict img, int img_dim1, int img_dim2, } -static void mixture_profile_fourier_transform(double *amps, int amps_len, - double *means, int means_dim1, int means_dim2, - double *vars, int vars_dim1, int vars_dim2, int vars_dim3, - double *v, int v_len, - double *w, int w_len, - double *out, int out_dim1, int out_dim2) +static void gaussian_fourier_transform_zero_mean( + double * restrict amps, int amps_len, + double * restrict vars, int vars_dim1, int vars_dim2, int vars_dim3, + double * restrict v, int v_len, + double * restrict w, int w_len, + double * restrict out, int out_dim1, int out_dim2) { - const double twopisquare = -2. * M_PI * M_PI; + const double negtwopisquare = -2. * M_PI * M_PI; int K = amps_len; int NV = v_len; int NW = w_len; int i, j, k; - double *s = (double*)malloc(sizeof(double) * NV * NW); - memset(s, 0, sizeof(double) * NV * NW); - for (j = 0; j < NW; j++) { double w_j = w[j]; double w_j_sqr = w_j * w_j; @@ -271,20 +268,78 @@ static void mixture_profile_fourier_transform(double *amps, int amps_len, int index = NV * j + i; double v_i = v[i]; double v_i_sqr = v_i * v_i; + double sum = 0.0; for (k = 0; k < K; k++) { int offset = k * 4; double a = vars[offset]; double b = vars[offset + 1]; double d = vars[offset + 3]; - s[index] += amps[k] * exp(twopisquare * (a * v_i_sqr + 2. * b * v_i * w_j + d * w_j_sqr)); + sum += amps[k] * exp(negtwopisquare * + (a * v_i_sqr + 2. * b * v_i * w_j + d * w_j_sqr)); } - out[index] = s[index]; + out[index] = sum; } } - free(s); return; } +static void gaussian_fourier_transform( + double * restrict amps, int amps_len, + double * restrict means, int means_dim1, int means_dim2, + double * restrict vars, int vars_dim1, int vars_dim2, int vars_dim3, + double * restrict v, int v_len, + double * restrict w, int w_len, + double * restrict out, int out_dim1, int out_dim2) +{ + // General case with non-zero means; output must be complex128. + // ASSUMES 2-d data, and that all means are equal. + + const double negtwopisquare = -2. * M_PI * M_PI; + + int K = amps_len; + int NV = v_len; + int NW = w_len; + int i, j, k; + + assert( + + double mu0 = means[0]; + double mu1 = means[1]; + + // + + for (j = 0; j < NW; j++) { + double w_j = w[j]; + double w_j_sqr = w_j * w_j; + for (i = 0; i < NV; i++) { + int index = 2 * (NV * j + i); + double v_i = v[i]; + double v_i_sqr = v_i * v_i; + double sum = 0.0; + for (k = 0; k < K; k++) { + int offset = k * 4; + double a = vars[offset]; + double b = vars[offset + 1]; + double d = vars[offset + 3]; + + sum += amps[k] * exp(negtwopisquare * + (a * v_i_sqr + 2. * b * v_i * w_j + d * w_j_sqr)); + } + + double angle = -2. * M_PI * (mu0 * vv[i] + mu1 * ww[j]); + out[index+0] = sum * cos(angle); + out[index+1] = sum * sin(angle); + } + } + return; +} + + + + + + + %} From 19eaed4967215244be3631220eada4dea90a9fd4 Mon Sep 17 00:00:00 2001 From: Dustin Lang Date: Mon, 6 Nov 2017 13:59:39 -0800 Subject: [PATCH 27/43] continue backpedal: remove non-zero-mean gaussian fourier transform --- tractor/mp_fourier.i | 57 -------------------------------------------- 1 file changed, 57 deletions(-) diff --git a/tractor/mp_fourier.i b/tractor/mp_fourier.i index fe21bf34..193f42dc 100644 --- a/tractor/mp_fourier.i +++ b/tractor/mp_fourier.i @@ -285,61 +285,4 @@ static void gaussian_fourier_transform_zero_mean( } -static void gaussian_fourier_transform( - double * restrict amps, int amps_len, - double * restrict means, int means_dim1, int means_dim2, - double * restrict vars, int vars_dim1, int vars_dim2, int vars_dim3, - double * restrict v, int v_len, - double * restrict w, int w_len, - double * restrict out, int out_dim1, int out_dim2) -{ - // General case with non-zero means; output must be complex128. - // ASSUMES 2-d data, and that all means are equal. - - const double negtwopisquare = -2. * M_PI * M_PI; - - int K = amps_len; - int NV = v_len; - int NW = w_len; - int i, j, k; - - assert( - - double mu0 = means[0]; - double mu1 = means[1]; - - // - - for (j = 0; j < NW; j++) { - double w_j = w[j]; - double w_j_sqr = w_j * w_j; - for (i = 0; i < NV; i++) { - int index = 2 * (NV * j + i); - double v_i = v[i]; - double v_i_sqr = v_i * v_i; - double sum = 0.0; - for (k = 0; k < K; k++) { - int offset = k * 4; - double a = vars[offset]; - double b = vars[offset + 1]; - double d = vars[offset + 3]; - - sum += amps[k] * exp(negtwopisquare * - (a * v_i_sqr + 2. * b * v_i * w_j + d * w_j_sqr)); - } - - double angle = -2. * M_PI * (mu0 * vv[i] + mu1 * ww[j]); - out[index+0] = sum * cos(angle); - out[index+1] = sum * sin(angle); - } - } - return; -} - - - - - - - %} From 26b465e84a8760822166397dff565f6ff174b8aa Mon Sep 17 00:00:00 2001 From: Dustin Lang Date: Mon, 6 Nov 2017 14:15:17 -0800 Subject: [PATCH 28/43] swig mp_fourier.i into intel_mp_fourier.py also --- tractor/Makefile | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tractor/Makefile b/tractor/Makefile index 5d17f6b0..11099b27 100644 --- a/tractor/Makefile +++ b/tractor/Makefile @@ -13,7 +13,7 @@ mpf: mp_fourier intel_mp_fourier mp_fourier: _mp_fourier$(PYTHON_SO_EXT) mp_fourier.py .PHONY: mp_fourier -intel_mp_fourier: _intel_mp_fourier$(PYTHON_SO_EXT) mp_fourier.py +intel_mp_fourier: _intel_mp_fourier$(PYTHON_SO_EXT) #intel_mp_fourier.py .PHONY: intel_mp_fourier mix: _mix$(PYTHON_SO_EXT) mix.py @@ -26,6 +26,7 @@ _mp_fourier$(PYTHON_SO_EXT): mp_fourier.i setup-mpf.py $(PYTHON) setup-mpf.py build_ext --inplace _intel_mp_fourier$(PYTHON_SO_EXT): mp_fourier.i setup-mpf.py + cp mp_fourier.i intel_mp_fourier.i CC=icc $(PYTHON) setup-mpf.py build_ext --inplace mix.py _mix$(PYTHON_SO_EXT): mix.i approx3.c gauss_masked.c setup-mix.py From ab5dfe556f5422fb73274822f9abd3d0f5c7a25d Mon Sep 17 00:00:00 2001 From: Dustin Lang Date: Mon, 6 Nov 2017 14:15:38 -0800 Subject: [PATCH 29/43] be explicit about gaussian_fourier_transform output type; assert array sizes --- tractor/mixture_profiles.py | 2 +- tractor/mp_fourier.i | 6 ++++++ 2 files changed, 7 insertions(+), 1 deletion(-) diff --git a/tractor/mixture_profiles.py b/tractor/mixture_profiles.py index 38fc8190..622e64e6 100755 --- a/tractor/mixture_profiles.py +++ b/tractor/mixture_profiles.py @@ -161,7 +161,7 @@ def getFourierTransform(self, v, w, use_mp_fourier=True, zero_mean=False): ''' if mp_fourier and use_mp_fourier and zero_mean: - f = np.zeros((len(w), len(v))) + f = np.zeros((len(w), len(v)), np.float64) mp_fourier.gaussian_fourier_transform_zero_mean( self.amp, self.var, v, w, f) return f diff --git a/tractor/mp_fourier.i b/tractor/mp_fourier.i index 193f42dc..4c1ada27 100644 --- a/tractor/mp_fourier.i +++ b/tractor/mp_fourier.i @@ -261,6 +261,12 @@ static void gaussian_fourier_transform_zero_mean( int NW = w_len; int i, j, k; + assert(vars_dim1 == K); + assert(vars_dim2 == 2); + assert(vars_dim3 == 2); + assert(out_dim1 == NV); + assert(out_dim2 == NW); + for (j = 0; j < NW; j++) { double w_j = w[j]; double w_j_sqr = w_j * w_j; From 0ae6fe10fcfade500fa2a6176f6903fbee03846a Mon Sep 17 00:00:00 2001 From: Dustin Lang Date: Mon, 6 Nov 2017 14:27:36 -0800 Subject: [PATCH 30/43] try pulling out lanczos_shift_image function --- tractor/mixture_profiles.py | 8 ++-- tractor/psf.py | 95 +++++++++++++------------------------ 2 files changed, 38 insertions(+), 65 deletions(-) diff --git a/tractor/mixture_profiles.py b/tractor/mixture_profiles.py index 622e64e6..7cde470b 100755 --- a/tractor/mixture_profiles.py +++ b/tractor/mixture_profiles.py @@ -6,12 +6,14 @@ import pylab as plt import matplotlib.cm as cm import numpy as np -#import scipy.spatial.distance as scp try: - from tractor import c_mp_fourier as mp_fourier + from tractor import intel_mp_fourier as mp_fourier except: - mp_fourier = None + try: + from tractor import mp_fourier + except: + mp_fourier = None from tractor.patch import Patch diff --git a/tractor/psf.py b/tractor/psf.py index 6afddf05..4ae509e5 100644 --- a/tractor/psf.py +++ b/tractor/psf.py @@ -2,23 +2,6 @@ from __future__ import division import numpy as np -# import numpy.fft_intel.libifft as m -# def irfftn_numpy(x, s=None, axes=None): -# a = np.asarray(x) -# no_trim = (s is None) and (axes is None) -# s, axes = m._cook_nd_args(a, s, axes, invreal=True) -# la = axes[-1] -# ovr_x = False -# if len(s) > 1: -# if not no_trim: -# a = m._fix_dimensions(a, s, axes) -# for ii in range(len(axes)-1): -# a = m.ifft(a, s[ii], axes[ii], overwrite_x=ovr_x) -# ovr_x = True -# a = m.irfft_numpy(a, n = s[-1], axis=la) -# return a -# m.irfftn_numpy = irfftn_numpy - from astrometry.util.miscutils import lanczos_filter from tractor.image import Image @@ -31,6 +14,37 @@ from tractor import mixture_profiles as mp from tractor import ducks +def lanczos_shift_image(img, dx, dy): + L = 3 + Lx = lanczos_filter(L, np.arange(-L, L+1) + dx) + Ly = lanczos_filter(L, np.arange(-L, L+1) + dy) + # Normalize the Lanczos interpolants (preserve flux) + Lx /= Lx.sum() + Ly /= Ly.sum() + + sx = correlate1d(img, Lx, axis=1, mode='constant') + shifted = correlate1d(sx, Ly, axis=0, mode='constant') + assert(np.all(np.isfinite(shifted))) + return shifted + +# assert(len(Lx) == 7) +# assert(len(Ly) == 7) +# from tractor.c_mp_fourier import correlate7, correlate7f +# outimg = np.empty(img.shape, np.float32) +# outimg[:,:] = img +# correlate7f(outimg, Lx, Ly, work_corr7f) +# +# +# try: +# from tractor import intel_mp_fourier as mp_fourier +# except: +# try: +# from tractor import mp_fourier +# except: +# mp_fourier = None +# from tractor.c_mp_fourier import correlate7, correlate7f + + class HybridPSF(object): pass @@ -126,40 +140,9 @@ def getPointSourcePatch(self, px, py, minval=0., modelMask=None, **kwargs): # image as usual, and then copy it into the modelMask # space. - L = self.Lorder - Lx = lanczos_filter(L, np.arange(-L, L+1) + dx) - Ly = lanczos_filter(L, np.arange(-L, L+1) + dy) - # Normalize the Lanczos interpolants (preserve flux) - Lx /= Lx.sum() - Ly /= Ly.sum() - - # print('img', img.dtype) - # print('Lx', Lx.dtype) - # print('Ly', Ly.dtype) - if modelMask is None: - - assert(len(Lx) == 7) - assert(len(Ly) == 7) - from tractor.c_mp_fourier import correlate7, correlate7f - #print('modelMask is None, dtype', img.dtype) - - outimg = np.empty(img.shape, np.float32) - outimg[:,:] = img - correlate7f(outimg, Lx, Ly, work_corr7f) + outimg = lanczos_shift_image(img, dx, dy) return Patch(x0, y0, outimg) - # if img.dtype == np.float32: - # correlate7f(img, Lx, Ly, work_corr7f) - # else: - # correlate7(img, Lx, Ly, work_corr7) - #return Patch(x0, y0, img) - #print('modelMask is None, dtype', img.dtype, 'done') - # from tractor.c_mp_fourier import correlate7 - # correlate7(img, Lx, Ly, work_corr7) - # sx = correlate1d(img, Lx, axis=1, mode='constant') - # shifted = correlate1d(sx, Ly, axis=0, mode='constant') - # assert(np.all(np.isfinite(shifted))) - # return Patch(x0, y0, shifted) # padding = L @@ -171,19 +154,7 @@ def getPointSourcePatch(self, px, py, minval=0., modelMask=None, **kwargs): xi,xo = get_overlapping_region(mx0-x0-padding, mx0-x0+mw-1+padding, 0, W-1) mm[yo,xo] = img[yi,xi] - assert(len(Lx) == 7) - assert(len(Ly) == 7) - #mm = np.require(mm, requirements=['A']) - from tractor.c_mp_fourier import correlate7, correlate7f - #print('mm, dtype', mm.dtype) - #if mm.dtype == np.float32: - correlate7f(mm, Lx, Ly, work_corr7f) - #else: - # correlate7(mm, Lx, Ly, work_corr7) - #print('mm, dtype', mm.dtype, 'done') - - #sx = correlate1d(mm, Lx, axis=1, mode='constant') - #mm = correlate1d(sx, Ly, axis=0, mode='constant') + mm = lanczos_shift_image(mm, dx, dy) mm = mm[padding:-padding, padding:-padding] assert(np.all(np.isfinite(mm))) From 75ac41f5fd20a084a405420c1ef4239a3d3f5f8d Mon Sep 17 00:00:00 2001 From: Dustin Lang Date: Tue, 7 Nov 2017 11:45:00 -0800 Subject: [PATCH 31/43] more intel vs generic mp_fourier work; factor out lanczos_shift_image() function --- tractor/Makefile | 2 +- tractor/galaxy.py | 67 +++++++++++++++++--------------------------- tractor/psf.py | 51 +++++++++++++++++---------------- tractor/setup-mpf.py | 2 +- 4 files changed, 55 insertions(+), 67 deletions(-) diff --git a/tractor/Makefile b/tractor/Makefile index 11099b27..7956c1a8 100644 --- a/tractor/Makefile +++ b/tractor/Makefile @@ -26,7 +26,7 @@ _mp_fourier$(PYTHON_SO_EXT): mp_fourier.i setup-mpf.py $(PYTHON) setup-mpf.py build_ext --inplace _intel_mp_fourier$(PYTHON_SO_EXT): mp_fourier.i setup-mpf.py - cp mp_fourier.i intel_mp_fourier.i + cat mp_fourier.i | sed s/mp_fourier/intel_mp_fourier/g > intel_mp_fourier.i CC=icc $(PYTHON) setup-mpf.py build_ext --inplace mix.py _mix$(PYTHON_SO_EXT): mix.i approx3.c gauss_masked.c setup-mix.py diff --git a/tractor/galaxy.py b/tractor/galaxy.py index 66f7133b..af1fcbf2 100644 --- a/tractor/galaxy.py +++ b/tractor/galaxy.py @@ -16,22 +16,6 @@ from __future__ import division import numpy as np -# import numpy.fft_intel.libifft as m -# def irfftn_numpy(x, s=None, axes=None): -# a = np.asarray(x) -# no_trim = (s is None) and (axes is None) -# s, axes = m._cook_nd_args(a, s, axes, invreal=True) -# la = axes[-1] -# ovr_x = False -# if len(s) > 1: -# if not no_trim: -# a = m._fix_dimensions(a, s, axes) -# for ii in range(len(axes)-1): -# a = m.ifft(a, s[ii], axes[ii], overwrite_x=ovr_x) -# ovr_x = True -# a = m.irfft_numpy(a, n = s[-1], axis=la) -# return a -# m.irfftn_numpy = irfftn_numpy from astrometry.util.miscutils import get_overlapping_region @@ -40,8 +24,6 @@ from tractor.patch import Patch, add_patches, ModelMask from tractor.basics import SingleProfileSource, BasicSource -#from .cache import Cache - debug_ps = None def get_galaxy_cache(): @@ -478,29 +460,32 @@ def run_mog(amix=None, mm=None): # after cutting G down to nearly its final size... tricky # tho - # Lanczos-3 interpolation in ~ the same way we do for - # pixelized PSFs. - from astrometry.util.miscutils import lanczos_filter - #from scipy.ndimage.filters import correlate1d - #L = 3 - L = fft_lanczos_order - Lx = lanczos_filter(L, np.arange(-L, L+1) + mux) - Ly = lanczos_filter(L, np.arange(-L, L+1) + muy) - # Normalize the Lanczos interpolants (preserve flux) - Lx /= Lx.sum() - Ly /= Ly.sum() - #print('Lx centroid', np.sum(Lx * (np.arange(-L,L+1)))) - #print('Ly centroid', np.sum(Ly * (np.arange(-L,L+1)))) - - assert(len(Lx) == 7) - assert(len(Ly) == 7) - #cx = correlate1d(G, Lx, axis=1, mode='constant') - #G = correlate1d(cx, Ly, axis=0, mode='constant') - #del cx - - from tractor.c_mp_fourier import correlate7 - G = np.require(G, requirements=['A']) - correlate7(G, Lx, Ly, work_corr7) + from tractor.psf import lanczos_shift_image + G = lanczos_shift_image(G, mux, muy) + + # # Lanczos-3 interpolation in ~ the same way we do for + # # pixelized PSFs. + # from astrometry.util.miscutils import lanczos_filter + # #from scipy.ndimage.filters import correlate1d + # #L = 3 + # L = fft_lanczos_order + # Lx = lanczos_filter(L, np.arange(-L, L+1) + mux) + # Ly = lanczos_filter(L, np.arange(-L, L+1) + muy) + # # Normalize the Lanczos interpolants (preserve flux) + # Lx /= Lx.sum() + # Ly /= Ly.sum() + # #print('Lx centroid', np.sum(Lx * (np.arange(-L,L+1)))) + # #print('Ly centroid', np.sum(Ly * (np.arange(-L,L+1)))) + # + # assert(len(Lx) == 7) + # assert(len(Ly) == 7) + # #cx = correlate1d(G, Lx, axis=1, mode='constant') + # #G = correlate1d(cx, Ly, axis=0, mode='constant') + # #del cx + # + # from tractor.c_mp_fourier import correlate7 + # G = np.require(G, requirements=['A']) + # correlate7(G, Lx, Ly, work_corr7) else: G = np.zeros((pH,pW), np.float32) diff --git a/tractor/psf.py b/tractor/psf.py index 4ae509e5..47c25403 100644 --- a/tractor/psf.py +++ b/tractor/psf.py @@ -14,7 +14,18 @@ from tractor import mixture_profiles as mp from tractor import ducks + +try: + from tractor import intel_mp_fourier as mp_fourier +except: + try: + from tractor import mp_fourier as mp_fourier + except: + mp_fourier = None +# from tractor.c_mp_fourier import correlate7, correlate7f + def lanczos_shift_image(img, dx, dy): + from scipy.ndimage import correlate1d L = 3 Lx = lanczos_filter(L, np.arange(-L, L+1) + dx) Ly = lanczos_filter(L, np.arange(-L, L+1) + dy) @@ -22,27 +33,20 @@ def lanczos_shift_image(img, dx, dy): Lx /= Lx.sum() Ly /= Ly.sum() - sx = correlate1d(img, Lx, axis=1, mode='constant') - shifted = correlate1d(sx, Ly, axis=0, mode='constant') - assert(np.all(np.isfinite(shifted))) - return shifted - -# assert(len(Lx) == 7) -# assert(len(Ly) == 7) -# from tractor.c_mp_fourier import correlate7, correlate7f -# outimg = np.empty(img.shape, np.float32) -# outimg[:,:] = img -# correlate7f(outimg, Lx, Ly, work_corr7f) -# -# -# try: -# from tractor import intel_mp_fourier as mp_fourier -# except: -# try: -# from tractor import mp_fourier -# except: -# mp_fourier = None -# from tractor.c_mp_fourier import correlate7, correlate7f + #print('mp_fourier:', mp_fourier) + if mp_fourier is None: + sx = correlate1d(img, Lx, axis=1, mode='constant') + outimg = correlate1d(sx, Ly, axis=0, mode='constant') + else: + assert(len(Lx) == 7) + assert(len(Ly) == 7) + outimg = np.empty(img.shape, np.float32) + outimg[:,:] = img + mp_fourier.correlate7f(outimg, Lx, Ly, work_corr7f) + + assert(np.all(np.isfinite(outimg))) + return outimg + class HybridPSF(object): @@ -145,19 +149,18 @@ def getPointSourcePatch(self, px, py, minval=0., modelMask=None, **kwargs): return Patch(x0, y0, outimg) # + L = 3 padding = L # Create a modelMask + padding sized stamp and insert PSF image into it #mm = np.zeros((mh+2*padding, mw+2*padding), img.dtype) mm = np.zeros((mh+2*padding, mw+2*padding), np.float32) - yi,yo = get_overlapping_region(my0-y0-padding, my0-y0+mh-1+padding, 0, H-1) xi,xo = get_overlapping_region(mx0-x0-padding, mx0-x0+mw-1+padding, 0, W-1) mm[yo,xo] = img[yi,xi] - mm = lanczos_shift_image(mm, dx, dy) - mm = mm[padding:-padding, padding:-padding] assert(np.all(np.isfinite(mm))) + return Patch(mx0, my0, mm) def getFourierTransformSize(self, radius): diff --git a/tractor/setup-mpf.py b/tractor/setup-mpf.py index 7edd865d..0817117b 100755 --- a/tractor/setup-mpf.py +++ b/tractor/setup-mpf.py @@ -6,7 +6,7 @@ if os.environ.get('CC') == 'icc': mpf_module = Extension('_intel_mp_fourier', - sources=['mp_fourier.i'], + sources=['intel_mp_fourier.i'], include_dirs=numpy_inc, extra_compile_args=['-g', '-xhost', '-qopt-report=5', '-axMIC-AVX512'], extra_link_args=['-g', '-lsvml'] From 9eb0b02c21ed362bc318b47a2fa15d568a7a335d Mon Sep 17 00:00:00 2001 From: Dustin Lang Date: Tue, 7 Nov 2017 12:00:06 -0800 Subject: [PATCH 32/43] enable in-place lanczos shifting --- tractor/galaxy.py | 29 ++++------------------------- tractor/psf.py | 21 ++++++++++++--------- 2 files changed, 16 insertions(+), 34 deletions(-) diff --git a/tractor/galaxy.py b/tractor/galaxy.py index af1fcbf2..7124979d 100644 --- a/tractor/galaxy.py +++ b/tractor/galaxy.py @@ -460,32 +460,11 @@ def run_mog(amix=None, mm=None): # after cutting G down to nearly its final size... tricky # tho + # Lanczos-3 interpolation in ~the same way we do for + # pixelized PSFs. from tractor.psf import lanczos_shift_image - G = lanczos_shift_image(G, mux, muy) - - # # Lanczos-3 interpolation in ~ the same way we do for - # # pixelized PSFs. - # from astrometry.util.miscutils import lanczos_filter - # #from scipy.ndimage.filters import correlate1d - # #L = 3 - # L = fft_lanczos_order - # Lx = lanczos_filter(L, np.arange(-L, L+1) + mux) - # Ly = lanczos_filter(L, np.arange(-L, L+1) + muy) - # # Normalize the Lanczos interpolants (preserve flux) - # Lx /= Lx.sum() - # Ly /= Ly.sum() - # #print('Lx centroid', np.sum(Lx * (np.arange(-L,L+1)))) - # #print('Ly centroid', np.sum(Ly * (np.arange(-L,L+1)))) - # - # assert(len(Lx) == 7) - # assert(len(Ly) == 7) - # #cx = correlate1d(G, Lx, axis=1, mode='constant') - # #G = correlate1d(cx, Ly, axis=0, mode='constant') - # #del cx - # - # from tractor.c_mp_fourier import correlate7 - # G = np.require(G, requirements=['A']) - # correlate7(G, Lx, Ly, work_corr7) + G = G.astype(np.float32) + lanczos_shift_image(G, mux, muy, inplace=True) else: G = np.zeros((pH,pW), np.float32) diff --git a/tractor/psf.py b/tractor/psf.py index 47c25403..41b7f644 100644 --- a/tractor/psf.py +++ b/tractor/psf.py @@ -24,7 +24,7 @@ mp_fourier = None # from tractor.c_mp_fourier import correlate7, correlate7f -def lanczos_shift_image(img, dx, dy): +def lanczos_shift_image(img, dx, dy, inplace=False): from scipy.ndimage import correlate1d L = 3 Lx = lanczos_filter(L, np.arange(-L, L+1) + dx) @@ -40,13 +40,22 @@ def lanczos_shift_image(img, dx, dy): else: assert(len(Lx) == 7) assert(len(Ly) == 7) - outimg = np.empty(img.shape, np.float32) - outimg[:,:] = img + if inplace: + assert(img.dtype == np.float32) + outimg = img + else: + outimg = np.empty(img.shape, np.float32) + outimg[:,:] = img mp_fourier.correlate7f(outimg, Lx, Ly, work_corr7f) assert(np.all(np.isfinite(outimg))) return outimg +#work_corr7 = np.zeros((1024,1024), np.float64) +#work_corr7 = np.require(work_corr7, requirements=['A']) + +work_corr7f = np.zeros((1024,1024), np.float32) +work_corr7f = np.require(work_corr7f, requirements=['A']) class HybridPSF(object): @@ -233,12 +242,6 @@ def getFourierTransform(self, px, py, radius): return rtn -work_corr7 = np.zeros((1024,1024), np.float64) -work_corr7 = np.require(work_corr7, requirements=['A']) - -work_corr7f = np.zeros((1024,1024), np.float32) -work_corr7f = np.require(work_corr7f, requirements=['A']) - class GaussianMixturePSF(MogParams, ducks.ImageCalibration): ''' A PSF model that is a mixture of general 2-D Gaussians From d1d559815f01f939bbd869464c32d8a46b7b59d7 Mon Sep 17 00:00:00 2001 From: Dustin Lang Date: Tue, 7 Nov 2017 12:01:51 -0800 Subject: [PATCH 33/43] tidy unused code --- tractor/galaxy.py | 7 ------- tractor/psf.py | 2 +- 2 files changed, 1 insertion(+), 8 deletions(-) diff --git a/tractor/galaxy.py b/tractor/galaxy.py index 7124979d..053ce27d 100644 --- a/tractor/galaxy.py +++ b/tractor/galaxy.py @@ -513,13 +513,6 @@ def run_mog(amix=None, mm=None): return Patch(ix0, iy0, G) - -work_corr7 = np.zeros((1024,1024), np.float64) -work_corr7 = np.require(work_corr7, requirements=['A']) - - -fft_lanczos_order = 3 - def _fourier_galaxy_debug_plots(G, shG, xi,yi,xo,yo, P, Fsum, pW,pH, psf): import pylab as plt diff --git a/tractor/psf.py b/tractor/psf.py index 41b7f644..267a16bb 100644 --- a/tractor/psf.py +++ b/tractor/psf.py @@ -19,7 +19,7 @@ from tractor import intel_mp_fourier as mp_fourier except: try: - from tractor import mp_fourier as mp_fourier + from tractor import mp_fourier except: mp_fourier = None # from tractor.c_mp_fourier import correlate7, correlate7f From 43571a775cad737cabd29e817297e9b1ddeeb570 Mon Sep 17 00:00:00 2001 From: Dustin Lang Date: Tue, 7 Nov 2017 12:27:00 -0800 Subject: [PATCH 34/43] add force_python option to lanczos_shift_image; debugging on import of intel and generic C versions --- tractor/psf.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/tractor/psf.py b/tractor/psf.py index 267a16bb..22f36f4b 100644 --- a/tractor/psf.py +++ b/tractor/psf.py @@ -2,7 +2,6 @@ from __future__ import division import numpy as np -from astrometry.util.miscutils import lanczos_filter from tractor.image import Image from tractor.pointsource import PointSource @@ -18,14 +17,17 @@ try: from tractor import intel_mp_fourier as mp_fourier except: + print('tractor.psf: failed import intel version of mp_fourier library. Falling back to generic version.') try: from tractor import mp_fourier except: + print('tractor.psf: failed import C version of mp_fourier library. Falling back to python version.') mp_fourier = None -# from tractor.c_mp_fourier import correlate7, correlate7f -def lanczos_shift_image(img, dx, dy, inplace=False): +def lanczos_shift_image(img, dx, dy, inplace=False, force_python=False): from scipy.ndimage import correlate1d + from astrometry.util.miscutils import lanczos_filter + L = 3 Lx = lanczos_filter(L, np.arange(-L, L+1) + dx) Ly = lanczos_filter(L, np.arange(-L, L+1) + dy) @@ -34,7 +36,7 @@ def lanczos_shift_image(img, dx, dy, inplace=False): Ly /= Ly.sum() #print('mp_fourier:', mp_fourier) - if mp_fourier is None: + if mp_fourier is None or force_python: sx = correlate1d(img, Lx, axis=1, mode='constant') outimg = correlate1d(sx, Ly, axis=0, mode='constant') else: From 3a14bb504712d77de94ab5f45a452719526904cf Mon Sep 17 00:00:00 2001 From: Dustin Lang Date: Tue, 7 Nov 2017 13:00:20 -0800 Subject: [PATCH 35/43] why does travis have trouble w/ 'restrict' keyword? --- .travis.yml | 2 ++ Makefile | 1 + tractor/Makefile | 2 +- tractor/mp_fourier.i | 11 +++++++++++ 4 files changed, 15 insertions(+), 1 deletion(-) diff --git a/.travis.yml b/.travis.yml index 49245ead..591ac8c4 100644 --- a/.travis.yml +++ b/.travis.yml @@ -41,6 +41,8 @@ install: # - which fitsverify script: + - echo CC is $CC + - $CC --version - python --version - python setup.py install - python setup.py build_ext --inplace diff --git a/Makefile b/Makefile index 3bf0686e..6fe837fe 100644 --- a/Makefile +++ b/Makefile @@ -10,6 +10,7 @@ emfit: FORCE $(MAKE) -C tractor emfit mpf: FORCE $(MAKE) -C tractor mpf + -$(MAKE) -C tractor intel_mp_fourier cov: coverage erase diff --git a/tractor/Makefile b/tractor/Makefile index 7956c1a8..315e2f14 100644 --- a/tractor/Makefile +++ b/tractor/Makefile @@ -7,7 +7,7 @@ PYTHON_SO_EXT ?= $(shell $(PYTHON) -c "from distutils import sysconfig; print(sy ceres: _ceres$(PYTHON_SO_EXT) ceres.py .PHONY: ceres -mpf: mp_fourier intel_mp_fourier +mpf: mp_fourier .PHONY: mpf mp_fourier: _mp_fourier$(PYTHON_SO_EXT) mp_fourier.py diff --git a/tractor/mp_fourier.i b/tractor/mp_fourier.i index 4c1ada27..08dd8b03 100644 --- a/tractor/mp_fourier.i +++ b/tractor/mp_fourier.i @@ -7,6 +7,7 @@ #include #include #include + %} %include "numpy.i" @@ -56,6 +57,16 @@ %inline %{ +// #ifdef __INTEL_COMPILER +// #define RESTRICT restrict +// #else +// #define RESTRICT __restrict +// #endif + +// #ifndef __INTEL_COMPILER +// #define restrict __restrict +// #endif + #if 0 } // fool emacs indenter #endif From 95122f822ae4db4e6158e4065f085188c95f5dcd Mon Sep 17 00:00:00 2001 From: Dustin Lang Date: Tue, 7 Nov 2017 13:08:05 -0800 Subject: [PATCH 36/43] travis: gcc version? --- .travis.yml | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/.travis.yml b/.travis.yml index 591ac8c4..ac3f37ba 100644 --- a/.travis.yml +++ b/.travis.yml @@ -41,8 +41,7 @@ install: # - which fitsverify script: - - echo CC is $CC - - $CC --version + - gcc --version - python --version - python setup.py install - python setup.py build_ext --inplace From d9d6495c2e12e8207566be1ecf59327d9d70d7d7 Mon Sep 17 00:00:00 2001 From: Dustin Lang Date: Tue, 7 Nov 2017 13:17:59 -0800 Subject: [PATCH 37/43] 'restrict' is in c99 --- tractor/setup-mpf.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tractor/setup-mpf.py b/tractor/setup-mpf.py index 0817117b..76f1736b 100755 --- a/tractor/setup-mpf.py +++ b/tractor/setup-mpf.py @@ -15,7 +15,7 @@ mpf_module = Extension('_mp_fourier', sources = ['mp_fourier.i' ], include_dirs = numpy_inc, - extra_compile_args=['-g'], + extra_compile_args=['-g', '-std=c99'], extra_link_args=['-g'], ) From feb840b800be34acdf1f24adaac89d63eb9df7cb Mon Sep 17 00:00:00 2001 From: Dustin Lang Date: Tue, 7 Nov 2017 13:25:48 -0800 Subject: [PATCH 38/43] 'restrict' is in c99 --- setup.py | 1 + 1 file changed, 1 insertion(+) diff --git a/setup.py b/setup.py index 76896c83..22a591ad 100644 --- a/setup.py +++ b/setup.py @@ -61,6 +61,7 @@ class CustomBuild(build): module_fourier = Extension('tractor._mp_fourier', sources = ['tractor/mp_fourier.i'], include_dirs = numpy_inc, + extra_compile_args=['-std=c99'], extra_objects = [], undef_macros=['NDEBUG'], ) From 674dd7189a7ca6340fd253390e2f21ddfc945e36 Mon Sep 17 00:00:00 2001 From: Dustin Lang Date: Tue, 7 Nov 2017 13:36:11 -0800 Subject: [PATCH 39/43] oops, remove asserts for non-existent variables --- tractor/mp_fourier.i | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/tractor/mp_fourier.i b/tractor/mp_fourier.i index 08dd8b03..e89448b0 100644 --- a/tractor/mp_fourier.i +++ b/tractor/mp_fourier.i @@ -85,8 +85,6 @@ static void correlate7(double* restrict img, int img_dim1, int img_dim2, assert(filtx_dim == 7); assert(filty_dim == 7); - assert(img_dim1 == out_dim1); - assert(img_dim2 == out_dim2); assert(work_dim1 >= img_dim1); assert(work_dim2 >= img_dim2); @@ -163,10 +161,10 @@ static void correlate7(double* restrict img, int img_dim1, int img_dim2, } -static void correlate7f(float* restrict img, int img_dim1, int img_dim2, +static void correlate7f(float* restrict img, int img_dim1, int img_dim2, double* restrict filtx, int filtx_dim, double* restrict filty, int filty_dim, - float* restrict work, int work_dim1, int work_dim2) { + float* restrict work, int work_dim1, int work_dim2) { // Output goes back into "img"! #ifdef __INTEL_COMPILER @@ -177,8 +175,6 @@ static void correlate7f(float* restrict img, int img_dim1, int img_dim2, assert(filtx_dim == 7); assert(filty_dim == 7); - assert(img_dim1 == out_dim1); - assert(img_dim2 == out_dim2); assert(work_dim1 >= img_dim1); assert(work_dim2 >= img_dim2); From 2e0087b1da37fe7df7c23192e8a5ae5a7c09347e Mon Sep 17 00:00:00 2001 From: Dustin Lang Date: Tue, 7 Nov 2017 13:44:45 -0800 Subject: [PATCH 40/43] NV/NW dimensions reversed in assert --- tractor/mp_fourier.i | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tractor/mp_fourier.i b/tractor/mp_fourier.i index e89448b0..680a7958 100644 --- a/tractor/mp_fourier.i +++ b/tractor/mp_fourier.i @@ -271,8 +271,8 @@ static void gaussian_fourier_transform_zero_mean( assert(vars_dim1 == K); assert(vars_dim2 == 2); assert(vars_dim3 == 2); - assert(out_dim1 == NV); - assert(out_dim2 == NW); + assert(out_dim1 == NW); + assert(out_dim2 == NV); for (j = 0; j < NW; j++) { double w_j = w[j]; From 9b11c75d2c70be986fd2966af25e56dc43e3b002 Mon Sep 17 00:00:00 2001 From: Dustin Lang Date: Fri, 17 Nov 2017 15:09:48 -0500 Subject: [PATCH 41/43] remove c_mp_fourier.i --- tractor/c_mp_fourier.i | 413 ----------------------------------------- 1 file changed, 413 deletions(-) delete mode 100755 tractor/c_mp_fourier.i diff --git a/tractor/c_mp_fourier.i b/tractor/c_mp_fourier.i deleted file mode 100755 index 1fc11604..00000000 --- a/tractor/c_mp_fourier.i +++ /dev/null @@ -1,413 +0,0 @@ -%module(package="tractor") c_mp_fourier - -%{ -#define SWIG_FILE_WITH_INIT -#define _GNU_SOURCE 1 -#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION -#include -#include -%} - -%include "numpy.i" - -%init %{ - // numpy - import_array(); -%} - -%apply (double* IN_ARRAY1, int DIM1) { - (double *amps, int amps_len), - (double *v, int v_len), - (double *w, int w_len) -}; -%apply (double* IN_ARRAY2, int DIM1, int DIM2) { - (double *means, int means_dim1, int means_dim2) -}; -%apply (double* IN_ARRAY3, int DIM1, int DIM2, int DIM3) { - (double *vars, int vars_dim1, int vars_dim2, int vars_dim3) -}; -%apply (double* INPLACE_ARRAY2, int DIM1, int DIM2) { - (double *out, int out_dim1, int out_dim2) -}; - - -%apply (double* INPLACE_ARRAY2, int DIM1, int DIM2) { - (double *work, int work_dim1, int work_dim2) -}; -%apply (double* INPLACE_ARRAY2, int DIM1, int DIM2) { - (double *img, int img_dim1, int img_dim2) -}; -%apply (double* IN_ARRAY1, int DIM1) { - (double *filtx, int filtx_dim) -}; -%apply (double* IN_ARRAY1, int DIM1) { - (double *filty, int filty_dim) -}; - - -%apply (float* INPLACE_ARRAY2, int DIM1, int DIM2) { - (float *work, int work_dim1, int work_dim2) -}; -%apply (float* INPLACE_ARRAY2, int DIM1, int DIM2) { - (float *img, int img_dim1, int img_dim2) -}; - -%inline %{ - -#if 0 - } -#endif - -static void correlate(double* img, int img_dim1, int img_dim2, - double* filtx, int filtx_dim, - double* filty, int filty_dim, - double* work, int work_dim1, int work_dim2, - double* out, int out_dim1, int out_dim2) { - - //printf("img %i x %i, filter x %i, filter y %i\n", img_dim1, img_dim2, filtx_dim, filty_dim); - - __assume_aligned(img, 64); - __assume_aligned(work, 64); - __assume_aligned(out, 64); - - assert(filtx_dim <= 8); - assert(filty_dim <= 8); - - assert(img_dim1 == out_dim1); - assert(img_dim2 == out_dim2); - assert(work_dim1 >= img_dim1); - assert(work_dim2 >= img_dim2); - - double filter[8]; - int i, j, k; - int W = img_dim2; - int H = img_dim1; - - assert(W > 8); - assert(H > 8); - - /* - for (i=0; i= img_dim1); - assert(work_dim2 >= img_dim2); - - double filter[7]; - int i, j, k; - int W = img_dim2; - int H = img_dim1; - - assert(W > 8); - assert(H > 8); - - memcpy(filter, filtx, 7 * sizeof(double)); - - // first run the filtx over image rows - for (j=0; j= img_dim1); - assert(work_dim2 >= img_dim2); - - float filter[7]; - int i, j, k; - int W = img_dim2; - int H = img_dim1; - - assert(W > 8); - assert(H > 8); - - //memcpy(filter, filtx, 7 * sizeof(double)); - for (i=0; i<7; i++) - filter[i] = filtx[i]; - - // first run the filtx over image rows - for (j=0; j Date: Fri, 17 Nov 2017 15:30:03 -0500 Subject: [PATCH 42/43] remove special intel_mp_fourier library; add intel compile instructions to setup.py --- setup.py | 30 +++++++++++++++--------------- tractor/Makefile | 8 -------- tractor/mixture_profiles.py | 7 ++----- tractor/mp_fourier.i | 10 ---------- tractor/psf.py | 10 +++------- tractor/setup-mpf.py | 22 ++++++++++------------ 6 files changed, 30 insertions(+), 57 deletions(-) diff --git a/setup.py b/setup.py index 22a591ad..be3475dd 100644 --- a/setup.py +++ b/setup.py @@ -1,9 +1,7 @@ from distutils.core import setup, Extension from distutils.command.build_ext import * from distutils.dist import Distribution - -# import sys -# py3 = (sys.version_info[0] >= 3) +import os # from http://stackoverflow.com/questions/12491328/python-distutils-not-include-the-swig-generated-module from distutils.command.build import build @@ -58,22 +56,26 @@ class CustomBuild(build): #extra_compile_args=['-O0','-g'], #extra_link_args=['-O0', '-g'], -module_fourier = Extension('tractor._mp_fourier', - sources = ['tractor/mp_fourier.i'], - include_dirs = numpy_inc, - extra_compile_args=['-std=c99'], - extra_objects = [], - undef_macros=['NDEBUG'], - ) - module_em = Extension('tractor._emfit', sources = ['tractor/emfit.i' ], include_dirs = numpy_inc, extra_objects = [], undef_macros=['NDEBUG'], ) -#extra_compile_args=['-O0','-g'], -#extra_link_args=['-O0', '-g'], + +kwargs = {} +if os.environ.get('CC') == 'icc': + kwargs.update(extra_compile_args=['-g', '-xhost', '-axMIC-AVX512'], + extra_link_args=['-g', '-lsvml']) +else: + kwargs.update(extra_compile_args=['-g', '-std=c99'], + extra_link_args=['-g']) + +module_fourier = Extension('tractor._mp_fourier', + sources = ['tractor/mp_fourier.i'], + include_dirs = numpy_inc, + undef_macros=['NDEBUG'], + **kwargs) class MyDistribution(Distribution): display_options = Distribution.display_options + [ @@ -101,8 +103,6 @@ class MyDistribution(Distribution): author_email="dstndstn@gmail.com", packages=['tractor', 'wise'], ext_modules = mods, -# py_modules = pymods, -# data_files=[('lib/python/wise', ['wise/wise-psf-avg.fits'])], package_data={'wise':['wise-psf-avg.fits', 'allsky-atlas.fits']}, package_dir={'wise':'wise', 'tractor':'tractor'}, url="http://theTractor.org/", diff --git a/tractor/Makefile b/tractor/Makefile index bf74fb20..20775e6e 100644 --- a/tractor/Makefile +++ b/tractor/Makefile @@ -13,9 +13,6 @@ mpf: mp_fourier mp_fourier: _mp_fourier$(PYTHON_SO_EXT) mp_fourier.py .PHONY: mp_fourier -intel_mp_fourier: _intel_mp_fourier$(PYTHON_SO_EXT) #intel_mp_fourier.py -.PHONY: intel_mp_fourier - mix: _mix$(PYTHON_SO_EXT) mix.py .PHONY: mix @@ -25,17 +22,12 @@ emfit: _emfit$(PYTHON_SO_EXT) emfit.py _mp_fourier$(PYTHON_SO_EXT): mp_fourier.i setup-mpf.py $(PYTHON) setup-mpf.py build_ext --inplace -_intel_mp_fourier$(PYTHON_SO_EXT): mp_fourier.i setup-mpf.py - cat mp_fourier.i | sed s/mp_fourier/intel_mp_fourier/g > intel_mp_fourier.i - CC=icc $(PYTHON) setup-mpf.py build_ext --inplace - mix.py _mix$(PYTHON_SO_EXT): mix.i approx3.c gauss_masked.c setup-mix.py $(PYTHON) setup-mix.py build_ext --inplace _emfit$(PYTHON_SO_EXT): emfit.i emfit2.c setup-emfit.py $(PYTHON) setup-emfit.py build_ext --inplace - NUMPY_INC := $(shell $(PYTHON) -c "from __future__ import print_function; from numpy.distutils.misc_util import get_numpy_include_dirs as d; print(' '.join('-I'+x for x in d()))") PYMOD_LIB ?= -L$(shell $(PYTHON_CONFIG) --prefix)/lib $(shell $(PYTHON_CONFIG) --libs) diff --git a/tractor/mixture_profiles.py b/tractor/mixture_profiles.py index e90a28bc..4e1e6a49 100755 --- a/tractor/mixture_profiles.py +++ b/tractor/mixture_profiles.py @@ -8,12 +8,9 @@ import numpy as np try: - from tractor import intel_mp_fourier as mp_fourier + from tractor import mp_fourier except: - try: - from tractor import mp_fourier - except: - mp_fourier = None + mp_fourier = None from tractor.patch import Patch diff --git a/tractor/mp_fourier.i b/tractor/mp_fourier.i index 680a7958..b8df11c2 100644 --- a/tractor/mp_fourier.i +++ b/tractor/mp_fourier.i @@ -57,16 +57,6 @@ %inline %{ -// #ifdef __INTEL_COMPILER -// #define RESTRICT restrict -// #else -// #define RESTRICT __restrict -// #endif - -// #ifndef __INTEL_COMPILER -// #define restrict __restrict -// #endif - #if 0 } // fool emacs indenter #endif diff --git a/tractor/psf.py b/tractor/psf.py index caedd553..bb7643f3 100644 --- a/tractor/psf.py +++ b/tractor/psf.py @@ -15,14 +15,10 @@ try: - from tractor import intel_mp_fourier as mp_fourier + from tractor import mp_fourier except: - print('tractor.psf: failed import intel version of mp_fourier library. Falling back to generic version.') - try: - from tractor import mp_fourier - except: - print('tractor.psf: failed import C version of mp_fourier library. Falling back to python version.') - mp_fourier = None + print('tractor.psf: failed to import C version of mp_fourier library. Falling back to python version.') + mp_fourier = None def lanczos_shift_image(img, dx, dy, inplace=False, force_python=False): from scipy.ndimage import correlate1d diff --git a/tractor/setup-mpf.py b/tractor/setup-mpf.py index 76f1736b..696c87f3 100755 --- a/tractor/setup-mpf.py +++ b/tractor/setup-mpf.py @@ -4,20 +4,18 @@ numpy_inc = get_numpy_include_dirs() +kwargs = {} if os.environ.get('CC') == 'icc': - mpf_module = Extension('_intel_mp_fourier', - sources=['intel_mp_fourier.i'], - include_dirs=numpy_inc, - extra_compile_args=['-g', '-xhost', '-qopt-report=5', '-axMIC-AVX512'], - extra_link_args=['-g', '-lsvml'] - ) + kwargs.update(extra_compile_args=['-g', '-xhost', '-axMIC-AVX512'], + extra_link_args=['-g', '-lsvml']) else: - mpf_module = Extension('_mp_fourier', - sources = ['mp_fourier.i' ], - include_dirs = numpy_inc, - extra_compile_args=['-g', '-std=c99'], - extra_link_args=['-g'], - ) + kwargs.update(extra_compile_args=['-g', '-std=c99'], + extra_link_args=['-g']) + +mpf_module = Extension('_mp_fourier', + sources = ['mp_fourier.i' ], + include_dirs = numpy_inc, + **kwargs) setup(name = 'Gaussian mixtures -- Fourier transform', version = '1.0', From 8dd319074007a837ccf6b3a162abc08b74546847 Mon Sep 17 00:00:00 2001 From: Dustin Lang Date: Fri, 17 Nov 2017 15:31:29 -0500 Subject: [PATCH 43/43] remove more intel_mp_fourier --- Makefile | 1 - 1 file changed, 1 deletion(-) diff --git a/Makefile b/Makefile index 6fe837fe..3bf0686e 100644 --- a/Makefile +++ b/Makefile @@ -10,7 +10,6 @@ emfit: FORCE $(MAKE) -C tractor emfit mpf: FORCE $(MAKE) -C tractor mpf - -$(MAKE) -C tractor intel_mp_fourier cov: coverage erase