NPY_NO_EXPORT int
_IsAligned(PyArrayObject *ap)
{
    unsigned int i, aligned = 1;
    const unsigned int alignment = PyArray_DESCR(ap)->alignment;

    /* The special casing for STRING and VOID types was removed
     * in accordance with http://projects.scipy.org/numpy/ticket/1227
     * It used to be that IsAligned always returned True for these
     * types, which is indeed the case when they are created using
     * PyArray_DescrConverter(), but not necessarily when using
     * PyArray_DescrAlignConverter(). */

    if (alignment == 1) {
        return 1;
    }
    aligned = npy_is_aligned(PyArray_DATA(ap), alignment);

    for (i = 0; i < PyArray_NDIM(ap); i++) {
#if NPY_RELAXED_STRIDES_CHECKING
        if (PyArray_DIM(ap, i) > 1) {
            /* if shape[i] == 1, the stride is never used */
            aligned &= npy_is_aligned((void*)PyArray_STRIDES(ap)[i],
                                      alignment);
        }
        else if (PyArray_DIM(ap, i) == 0) {
            /* an array with zero elements is always aligned */
            return 1;
        }
#else /* not NPY_RELAXED_STRIDES_CHECKING */
        aligned &= npy_is_aligned((void*)PyArray_STRIDES(ap)[i], alignment);
#endif /* not NPY_RELAXED_STRIDES_CHECKING */
    }
    return aligned != 0;
}

/* 0-strided arrays are not contiguous (even if dimension == 1) */
static int
_IsFortranContiguous(PyArrayObject *ap)
{
    npy_intp sd;
    npy_intp dim;
    int i;

    if (PyArray_NDIM(ap) == 0) {
        return 1;
    }
    sd = PyArray_DESCR(ap)->elsize;
    if (PyArray_NDIM(ap) == 1) {
        return PyArray_DIMS(ap)[0] == 1 || sd == PyArray_STRIDES(ap)[0];
    }
    for (i = 0; i < PyArray_NDIM(ap); ++i) {
        dim = PyArray_DIMS(ap)[i];
        /* fortran contiguous by definition */
        if (dim == 0) {
            return 1;
        }
        if (PyArray_STRIDES(ap)[i] != sd) {
            return 0;
        }
        sd *= dim;
    }
    return 1;
}

NPY_NO_EXPORT int
_IsAligned(PyArrayObject *ap)
{
    unsigned int i;
    npy_uintp aligned;
    npy_uintp alignment = PyArray_DESCR(ap)->alignment;

    /* alignment 1 types should have a efficient alignment for copy loops */
    if (PyArray_ISFLEXIBLE(ap) || PyArray_ISSTRING(ap)) {
        alignment = 16;
    }

    if (alignment == 1) {
        return 1;
    }
    aligned = (npy_uintp)PyArray_DATA(ap);

    for (i = 0; i < PyArray_NDIM(ap); i++) {
#if NPY_RELAXED_STRIDES_CHECKING
        /* skip dim == 1 as it is not required to have stride 0 */
        if (PyArray_DIM(ap, i) > 1) {
            /* if shape[i] == 1, the stride is never used */
            aligned |= (npy_uintp)PyArray_STRIDES(ap)[i];
        }
        else if (PyArray_DIM(ap, i) == 0) {
            /* an array with zero elements is always aligned */
            return 1;
        }
#else /* not NPY_RELAXED_STRIDES_CHECKING */
        aligned |= (npy_uintp)PyArray_STRIDES(ap)[i];
#endif /* not NPY_RELAXED_STRIDES_CHECKING */
    }
    return npy_is_aligned((void *)aligned, alignment);
}

PyObject* r2k(PyObject *self, PyObject *args)
{
  Py_complex alpha;
  PyArrayObject* a;
  PyArrayObject* b;
  double beta;
  PyArrayObject* c;
  if (!PyArg_ParseTuple(args, "DOOdO", &alpha, &a, &b, &beta, &c))
    return NULL;
  int n = PyArray_DIMS(a)[0];
  int k = PyArray_DIMS(a)[1];
  for (int d = 2; d < PyArray_NDIM(a); d++)
    k *= PyArray_DIMS(a)[d];
  int ldc = PyArray_STRIDES(c)[0] / PyArray_STRIDES(c)[1];
  if (PyArray_DESCR(a)->type_num == NPY_DOUBLE)
    dsyr2k_("u", "t", &n, &k,
            (double*)(&alpha), DOUBLEP(a), &k,
            DOUBLEP(b), &k, &beta,
            DOUBLEP(c), &ldc);
  else
    zher2k_("u", "c", &n, &k,
            (void*)(&alpha), (void*)COMPLEXP(a), &k,
            (void*)COMPLEXP(b), &k, &beta,
            (void*)COMPLEXP(c), &ldc);
  Py_RETURN_NONE;
}

PyObject* gemm(PyObject *self, PyObject *args)
{
  Py_complex alpha;
  PyArrayObject* a;
  PyArrayObject* b;
  Py_complex beta;
  PyArrayObject* c;
  char transa = 'n';
  if (!PyArg_ParseTuple(args, "DOODO|c", &alpha, &a, &b, &beta, &c, &transa))
    return NULL;
  int m, k, lda, ldb, ldc;
  if (transa == 'n')
    {
      m = PyArray_DIMS(a)[1];
      for (int i = 2; i < PyArray_NDIM(a); i++)
        m *= PyArray_DIMS(a)[i];
      k = PyArray_DIMS(a)[0];
      lda = MAX(1, PyArray_STRIDES(a)[0] / PyArray_STRIDES(a)[PyArray_NDIM(a) - 1]);
      ldb = MAX(1, PyArray_STRIDES(b)[0] / PyArray_STRIDES(b)[1]);
      ldc = MAX(1, PyArray_STRIDES(c)[0] / PyArray_STRIDES(c)[PyArray_NDIM(c) - 1]);
    }
  else
    {
      k = PyArray_DIMS(a)[1];
      for (int i = 2; i < PyArray_NDIM(a); i++)
        k *= PyArray_DIMS(a)[i];
      m = PyArray_DIMS(a)[0];
      lda = MAX(1, k);
      ldb = MAX(1, PyArray_STRIDES(b)[0] / PyArray_STRIDES(b)[PyArray_NDIM(b) - 1]);
      ldc = MAX(1, PyArray_STRIDES(c)[0] / PyArray_STRIDES(c)[1]);

    }
  int n = PyArray_DIMS(b)[0];
  if (PyArray_DESCR(a)->type_num == NPY_DOUBLE)
    dgemm_(&transa, "n", &m, &n, &k,
           &(alpha.real),
           DOUBLEP(a), &lda,
           DOUBLEP(b), &ldb,
           &(beta.real),
           DOUBLEP(c), &ldc);
  else
    zgemm_(&transa, "n", &m, &n, &k,
           &alpha,
           (void*)COMPLEXP(a), &lda,
           (void*)COMPLEXP(b), &ldb,
           &beta,
           (void*)COMPLEXP(c), &ldc);
  Py_RETURN_NONE;
}

PyObject *
gdkpixbuf_get_pixels_array(PyObject *pixbuf_pyobject)
{
    GdkPixbuf *pixbuf = GDK_PIXBUF(((PyGObject *)pixbuf_pyobject)->obj);
    PyArrayObject *array;
    npy_intp dims[3] = { 0, 0, 3 };

    dims[0] = gdk_pixbuf_get_height(pixbuf);
    dims[1] = gdk_pixbuf_get_width(pixbuf);

    if (gdk_pixbuf_get_has_alpha(pixbuf))
        dims[2] = 4;
    guchar *pixels = gdk_pixbuf_get_pixels(pixbuf);

    array = (PyArrayObject *)PyArray_SimpleNewFromData(3, dims, NPY_UBYTE,
                                                       pixels);

    if (array == NULL)
        return NULL;

    PyArray_STRIDES(array)[0] = gdk_pixbuf_get_rowstride(pixbuf);
    // the array holds a ref to the pixbuf pixels through this wrapper
    Py_INCREF(pixbuf_pyobject);

#ifdef NPY_1_7_API_VERSION
    PyArray_SetBaseObject(array, (PyObject *)pixbuf_pyobject);
#else
    array->base = (PyObject *)pixbuf_pyobject;
#endif

    return PyArray_Return(array);
}

NPY_NO_EXPORT int
IsAligned(PyArrayObject *ap)
{
    return raw_array_is_aligned(PyArray_NDIM(ap), PyArray_DIMS(ap),
                                PyArray_DATA(ap), PyArray_STRIDES(ap),
                                PyArray_DESCR(ap)->alignment);
}

/* Gets a half-open range [start, end) which contains the array data */
NPY_NO_EXPORT void
get_array_memory_extents(PyArrayObject *arr,
                    npy_uintp *out_start, npy_uintp *out_end)
{
    npy_uintp start, end;
    npy_intp idim, ndim = PyArray_NDIM(arr);
    npy_intp *dimensions = PyArray_DIMS(arr),
            *strides = PyArray_STRIDES(arr);

    /* Calculate with a closed range [start, end] */
    start = end = (npy_uintp)PyArray_DATA(arr);
    for (idim = 0; idim < ndim; ++idim) {
        npy_intp stride = strides[idim], dim = dimensions[idim];
        /* If the array size is zero, return an empty range */
        if (dim == 0) {
            *out_start = *out_end = (npy_uintp)PyArray_DATA(arr);
            return;
        }
        /* Expand either upwards or downwards depending on stride */
        else {
            if (stride > 0) {
                end += stride*(dim-1);
            }
            else if (stride < 0) {
                start += stride*(dim-1);
            }
        }
    }

    /* Return a half-open range */
    *out_start = start;
    *out_end = end + PyArray_DESCR(arr)->elsize;
}

NPY_NO_EXPORT int
IsUintAligned(PyArrayObject *ap)
{
    return raw_array_is_aligned(PyArray_NDIM(ap), PyArray_DIMS(ap),
                                PyArray_DATA(ap), PyArray_STRIDES(ap),
                                npy_uint_alignment(PyArray_DESCR(ap)->elsize));
}

NPY_NO_EXPORT int
_IsAligned(PyArrayObject *ap)
{
    unsigned int i;
    npy_uintp aligned;
    npy_uintp alignment = PyArray_DESCR(ap)->alignment;

    /* alignment 1 types should have a efficient alignment for copy loops */
    if (PyArray_ISFLEXIBLE(ap) || PyArray_ISSTRING(ap)) {
        npy_intp itemsize = PyArray_ITEMSIZE(ap);
        /* power of two sizes may be loaded in larger moves */
        if (((itemsize & (itemsize - 1)) == 0)) {
            alignment = itemsize > NPY_MAX_COPY_ALIGNMENT ?
                NPY_MAX_COPY_ALIGNMENT : itemsize;
        }
        else {
            /* if not power of two it will be accessed bytewise */
            alignment = 1;
        }
    }

    if (alignment == 1) {
        return 1;
    }
    aligned = (npy_uintp)PyArray_DATA(ap);

    for (i = 0; i < PyArray_NDIM(ap); i++) {
#if NPY_RELAXED_STRIDES_CHECKING
        /* skip dim == 1 as it is not required to have stride 0 */
        if (PyArray_DIM(ap, i) > 1) {
            /* if shape[i] == 1, the stride is never used */
            aligned |= (npy_uintp)PyArray_STRIDES(ap)[i];
        }
        else if (PyArray_DIM(ap, i) == 0) {
            /* an array with zero elements is always aligned */
            return 1;
        }
#else /* not NPY_RELAXED_STRIDES_CHECKING */
        aligned |= (npy_uintp)PyArray_STRIDES(ap)[i];
#endif /* not NPY_RELAXED_STRIDES_CHECKING */
    }
    return npy_is_aligned((void *)aligned, alignment);
}

PyObject* gemv(PyObject *self, PyObject *args)
{
  Py_complex alpha;
  PyArrayObject* a;
  PyArrayObject* x;
  Py_complex beta;
  PyArrayObject* y;
  char trans = 't';
  if (!PyArg_ParseTuple(args, "DOODO|c", &alpha, &a, &x, &beta, &y, &trans))
    return NULL;

  int m, n, lda, itemsize, incx, incy;

  if (trans == 'n')
    {
      m = PyArray_DIMS(a)[1];
      for (int i = 2; i < PyArray_NDIM(a); i++)
        m *= PyArray_DIMS(a)[i];
      n = PyArray_DIMS(a)[0];
      lda = MAX(1, m);
    }
  else
    {
      n = PyArray_DIMS(a)[0];
      for (int i = 1; i < PyArray_NDIM(a)-1; i++)
        n *= PyArray_DIMS(a)[i];
      m = PyArray_DIMS(a)[PyArray_NDIM(a)-1];
      lda = MAX(1, m);
    }

  if (PyArray_DESCR(a)->type_num == NPY_DOUBLE)
    itemsize = sizeof(double);
  else
    itemsize = sizeof(double_complex);

  incx = PyArray_STRIDES(x)[0]/itemsize;
  incy = 1;

  if (PyArray_DESCR(a)->type_num == NPY_DOUBLE)
    dgemv_(&trans, &m, &n,
           &(alpha.real),
           DOUBLEP(a), &lda,
           DOUBLEP(x), &incx,
           &(beta.real),
           DOUBLEP(y), &incy);
  else
    zgemv_(&trans, &m, &n,
           &alpha,
           (void*)COMPLEXP(a), &lda,
           (void*)COMPLEXP(x), &incx,
           &beta,
           (void*)COMPLEXP(y), &incy);
  Py_RETURN_NONE;
}

	UMatData* allocate(PyObject* o, int dims, const int* sizes, int type,
			size_t* step) const {
		UMatData* u = new UMatData(this);
		u->data = u->origdata = (uchar*) PyArray_DATA((PyArrayObject*) o);
		npy_intp* _strides = PyArray_STRIDES((PyArrayObject*) o);
		for (int i = 0; i < dims - 1; i++)
			step[i] = (size_t) _strides[i];
		step[dims - 1] = CV_ELEM_SIZE(type);
		u->size = sizes[0] * step[0];
		u->userdata = o;
		return u;
	}

void local_histogram(double* H, 
		     unsigned int clamp, 
		     PyArrayIterObject* iter, 
		     const unsigned int* size)
{
  PyArrayObject *block, *im = iter->ao; 
  PyArrayIterObject* block_iter; 
  unsigned int i, left, right, center, halfsize, dim, offset=0; 
  npy_intp block_dims[3];

  UPDATE_ITERATOR_COORDS(iter); 

  /* Compute block corners */ 
  for (i=0; i<3; i++) {
    center = iter->coordinates[i];
    halfsize = size[i]/2; 
    dim = PyArray_DIM(im, i);
  
    /* Left handside corner */ 
    if (center<halfsize)
      left = 0; 
    else
      left = center-halfsize; 

    /* Right handside corner (plus one)*/ 
    right = center+halfsize+1; 
    if (right>dim) 
      right = dim; 

    /* Block properties */ 
    offset += left*PyArray_STRIDE(im, i); 
    block_dims[i] = right-left;

  }

  /* Create the block as a vew and the block iterator */ 
  block = (PyArrayObject*)PyArray_New(&PyArray_Type, 3, block_dims, 
				      PyArray_TYPE(im), PyArray_STRIDES(im), 
				      (void*)(PyArray_DATA(im)+offset), 
				      PyArray_ITEMSIZE(im),
				      NPY_BEHAVED, NULL);
  block_iter = (PyArrayIterObject*)PyArray_IterNew((PyObject*)block); 

  /* Compute block histogram */ 
  histogram(H, clamp, block_iter); 

  /* Free memory */ 
  Py_XDECREF(block_iter); 
  Py_XDECREF(block); 

  return; 
}		     

//initialization of BallTree object
// argument is a single array of size [D,N]
static int
BallTree_init(BallTreeObject *self, PyObject *args, PyObject *kwds){
  //we use goto statements : all variables should be declared up front
    PyObject *arg=NULL;
    PyObject *arr=NULL;
    long int leaf_size=20;

    static char *kwlist[] = {"x", "leafsize", NULL};

    if (!PyArg_ParseTupleAndKeywords(args, kwds, "O|l", kwlist,
        &arg,&leaf_size) )
        goto fail;

    if(leaf_size <= 0){
        PyErr_SetString(PyExc_ValueError,
            "BallTree : leaf size must be greater than zero");
        goto fail;
    }

  //view this object as an array of doubles,
    arr = PyArray_FROM_OTF(arg,NPY_DOUBLE,0);
    if(arr==NULL)
        goto fail;

  //check that it is 2D
    if( PyArray_NDIM(arr) != 2)
        goto fail;

    if(self != NULL){
    //create the list of points
        self->size = PyArray_DIMS(arr)[0];
        self->dim = PyArray_DIMS(arr)[1];

        int inc = PyArray_STRIDES(arr)[1]/PyArray_DESCR(arr)->elsize;
        self->Points = new std::vector<BallTree_Point*>(self->size);
        for(int i=0;i<self->size;i++)
            self->Points->at(i) = new BallTree_Point(arr,
            (double*)PyArray_GETPTR2(arr,i,0),
            inc, PyArray_DIM(arr,1));
        self->tree = new BallTree<BallTree_Point>(*(self->Points),
            leaf_size);
    }
    self->data = arr;
  //Py_DECREF(arr);
    return 0;

    fail:
    Py_XDECREF(arr);
    return -1;
}

PyObject* rk(PyObject *self, PyObject *args)
{
    double alpha;
    PyArrayObject* a;
    double beta;
    PyArrayObject* c;
    char trans = 'c';
    if (!PyArg_ParseTuple(args, "dOdO|c", &alpha, &a, &beta, &c, &trans))
        return NULL;

    int n = PyArray_DIMS(c)[0];
    
    int k, lda;
    
    if (trans == 'c') {
        k = PyArray_DIMS(a)[1];
        for (int d = 2; d < PyArray_NDIM(a); d++)
            k *= PyArray_DIMS(a)[d];
        lda = k;
    }
    else {
        k = PyArray_DIMS(a)[0];
        lda = n;
    }
    
    int ldc = PyArray_STRIDES(c)[0] / PyArray_STRIDES(c)[1];
    if (PyArray_DESCR(a)->type_num == NPY_DOUBLE)
        dsyrk_("u", &trans, &n, &k,
               &alpha, DOUBLEP(a), &lda, &beta,
               DOUBLEP(c), &ldc);
    else
        zherk_("u", &trans, &n, &k,
               &alpha, (void*)COMPLEXP(a), &lda, &beta,
               (void*)COMPLEXP(c), &ldc);
    Py_RETURN_NONE;
}

/* This also makes sure that the data segment is aligned with
   an itemsize address as well by returning one if not true.
*/
static int
_bad_strides(PyArrayObject *ap)
{
    register int itemsize = PyArray_ITEMSIZE(ap);
    register int i, N=PyArray_NDIM(ap);
    register intp *strides = PyArray_STRIDES(ap);

    if (((intp)(ap->data) % itemsize) != 0)
        return 1;
    for (i=0; i<N; i++) {
        if ((strides[i] < 0) || (strides[i] % itemsize) != 0)
            return 1;
    }

    return 0;
}

/* Fill in the info structure */
static _buffer_info_t*
_buffer_info_new(PyArrayObject *arr)
{
    _buffer_info_t *info;
    _tmp_string_t fmt = {NULL, 0, 0};
    int k;

    info = malloc(sizeof(_buffer_info_t));
    if (info == NULL) {
        goto fail;
    }

    /* Fill in format */
    if (_buffer_format_string(PyArray_DESCR(arr), &fmt, arr, NULL, NULL) != 0) {
        free(fmt.s);
        goto fail;
    }
    _append_char(&fmt, '\0');
    info->format = fmt.s;

    /* Fill in shape and strides */
    info->ndim = PyArray_NDIM(arr);

    if (info->ndim == 0) {
        info->shape = NULL;
        info->strides = NULL;
    }
    else {
        info->shape = malloc(sizeof(Py_ssize_t) * PyArray_NDIM(arr) * 2 + 1);
        if (info->shape == NULL) {
            goto fail;
        }
        info->strides = info->shape + PyArray_NDIM(arr);
        for (k = 0; k < PyArray_NDIM(arr); ++k) {
            info->shape[k] = PyArray_DIMS(arr)[k];
            info->strides[k] = PyArray_STRIDES(arr)[k];
        }
    }

    return info;

fail:
    free(info);
    return NULL;
}

/* Gets a half-open range [start, end) which contains the array data */
static void
get_array_memory_extents(PyArrayObject *arr,
                         npy_uintp *out_start, npy_uintp *out_end,
                         npy_uintp *num_bytes)
{
    npy_intp low, upper;
    int j;
    offset_bounds_from_strides(PyArray_ITEMSIZE(arr), PyArray_NDIM(arr),
                               PyArray_DIMS(arr), PyArray_STRIDES(arr),
                               &low, &upper);
    *out_start = (npy_uintp)PyArray_DATA(arr) + (npy_uintp)low;
    *out_end = (npy_uintp)PyArray_DATA(arr) + (npy_uintp)upper;

    *num_bytes = PyArray_ITEMSIZE(arr);
    for (j = 0; j < PyArray_NDIM(arr); ++j) {
        *num_bytes *= PyArray_DIM(arr, j);
    }
}

static PyObject *
array_repr_builtin(PyArrayObject *self, int repr)
{
    PyObject *ret;
    char *string;
    /* max_n initial value is arbitrary, dump_data will extend it */
    Py_ssize_t n = 0, max_n = PyArray_NBYTES(self) * 4 + 7;

    if ((string = PyArray_malloc(max_n)) == NULL) {
        return PyErr_NoMemory();
    }

    if (dump_data(&string, &n, &max_n, PyArray_DATA(self),
                  PyArray_NDIM(self), PyArray_DIMS(self),
                  PyArray_STRIDES(self), self) < 0) {
        PyArray_free(string);
        return NULL;
    }

    if (repr) {
        if (PyArray_ISEXTENDED(self)) {
            ret = PyUString_FromFormat("array(%s, '%c%d')",
                                       string,
                                       PyArray_DESCR(self)->type,
                                       PyArray_DESCR(self)->elsize);
        }
        else {
            ret = PyUString_FromFormat("array(%s, '%c')",
                                       string,
                                       PyArray_DESCR(self)->type);
        }
    }
    else {
        ret = PyUString_FromStringAndSize(string, n);
    }

    PyArray_free(string);
    return ret;
}

NPY_NO_EXPORT char *
index2ptr(PyArrayObject *mp, npy_intp i)
{
    npy_intp dim0;

    if (PyArray_NDIM(mp) == 0) {
        PyErr_SetString(PyExc_IndexError, "0-d arrays can't be indexed");
        return NULL;
    }
    dim0 = PyArray_DIMS(mp)[0];
    if (i < 0) {
        i += dim0;
    }
    if (i == 0 && dim0 > 0) {
        return PyArray_DATA(mp);
    }
    if (i > 0 && i < dim0) {
        return PyArray_DATA(mp)+i*PyArray_STRIDES(mp)[0];
    }
    PyErr_SetString(PyExc_IndexError,"index out of bounds");
    return NULL;
}