bp::list	wrapped_sparse_hess_no_repeat(short tape_tag, bpn::ndarray &bpn_x, bpn::ndarray &bpn_options){
	size_t tape_stats[STAT_SIZE];
	tapestats(tape_tag, tape_stats);
	int N = tape_stats[NUM_INDEPENDENTS];

	double* x = (double*) nu::data(bpn_x);
	int* options  = (int*) nu::data(bpn_options);
// 	int options[2] = {0,0};

	int nnz=-1;
	unsigned int *rind = NULL;
	unsigned int *cind = NULL;
	double   *values   = NULL;
	sparse_hess(tape_tag,  N, 0, x, &nnz, &rind, &cind, &values, options);

	npy_intp ret_nnz = static_cast<npy_intp>(nnz);
	bp::object bp_rind   ( bp::handle<>(PyArray_SimpleNewFromData(1, &ret_nnz, NPY_UINT, (char*) rind )));
	bp::object bp_cind   ( bp::handle<>(PyArray_SimpleNewFromData(1, &ret_nnz, NPY_UINT, (char*) cind )));
	bp::object bp_values ( bp::handle<>(PyArray_SimpleNewFromData(1, &ret_nnz, NPY_DOUBLE, (char*) values )));

	bpn::ndarray ret_rind   = boost::python::extract<boost::python::numpy::ndarray>(bp_rind);
	bpn::ndarray ret_cind   = boost::python::extract<boost::python::numpy::ndarray>(bp_cind);
	bpn::ndarray ret_values = boost::python::extract<boost::python::numpy::ndarray>(bp_values);

	bp::list retvals;
	retvals.append(ret_nnz);
	retvals.append(ret_rind);
	retvals.append(ret_cind);
	retvals.append(ret_values);

	return retvals;

}

void HIDDEN gadget_initkernel(PyObject * m) {
    fill_koverlap();
    npy_intp kline_dims[] = {KLINE_BINS};
    PyArrayObject * kline_a = (PyArrayObject *)PyArray_SimpleNewFromData(1, kline_dims, NPY_FLOAT, kline);
    PyArrayObject * klinesq_a = (PyArrayObject *)PyArray_SimpleNewFromData(1, kline_dims, NPY_FLOAT, klinesq);
    npy_intp koverlap_dims[] = {KOVERLAP_BINS, KOVERLAP_BINS, KOVERLAP_BINS, KOVERLAP_BINS};
    PyArrayObject * koverlap_a = (PyArrayObject *)PyArray_SimpleNewFromData(4, koverlap_dims, NPY_FLOAT, koverlap);

    generic_functions[0] =  PyUFunc_f_f;
    generic_functions[1] =  PyUFunc_d_d;
    koverlap_functions[0] =  PyUFunc_ffff_f;
    koverlap_functions[1] =  PyUFunc_dddd_d;
    PyObject * k0_u = PyUFunc_FromFuncAndData(generic_functions,
                      k0_data, generic_signatures,
                      2, 1, 1,
                      PyUFunc_None,
                      "k0", k0_doc_string, 0);
    PyObject * kline_u = PyUFunc_FromFuncAndData(generic_functions,
                         kline_data, generic_signatures,
                         2, 1, 1,
                         PyUFunc_None,
                         "kline", kline_doc_string, 0);
    PyObject * koverlap_u = PyUFunc_FromFuncAndData(koverlap_functions,
                            koverlap_data, koverlap_signatures,
                            2, 4, 1, PyUFunc_None,
                            "koverlap", koverlap_doc_string, 0);

    PyModule_AddObject(m, "k0", k0_u);
    PyModule_AddObject(m, "kline", kline_u);
    PyModule_AddObject(m, "koverlap", koverlap_u);
    PyModule_AddObject(m, "akline", (PyObject*)kline_a);
    PyModule_AddObject(m, "aklinesq", (PyObject*)klinesq_a);
    PyModule_AddObject(m, "akoverlap", (PyObject*)koverlap_a);
}

static PyObject* 
SparseMatrix_getCSRrepresentation(SparseMatrix *self, 
				  PyObject *args)
{
  PyObject *rowout,*colout,*aout;
  int nnz,nr;
  npy_intp dims[1];
  int_t *rowptr,*colind;
  double *a;
  a = (double*)self->A.nzval;
  rowptr = self->A.rowptr;
  colind = self->A.colind;
  nr  = self->dim[0];
  nnz = self->A.nnz;

  dims[0] = nr+1;
  rowout = PyArray_SimpleNewFromData(1,dims,NPY_INT,
				   (void*)rowptr);
  dims[0] = nnz;
  colout = PyArray_SimpleNewFromData(1,dims,NPY_INT,
				     (void*)colind);

  dims[0] = nnz;
  aout = PyArray_SimpleNewFromData(1,dims,NPY_DOUBLE,
				   (void*)a);

  return Py_BuildValue("OOO",
		       rowout,
		       colout,
		       aout);
}

PyObject* VideoCapture_get_array(VideoCapture *self)
{
  switch(self->camptr->format)
  {
  case FG_BINARY:
  case FG_GRAY:
    {
      npy_intp dims[2] = {self->camptr->height,self->camptr->width};
      return PyArray_SimpleNewFromData(2,dims,NPY_UINT8,self->camptr->ImgPtr);
    }
  case FG_GRAY16:
    {
      npy_intp dims[2] = {self->camptr->height,self->camptr->width};
      return PyArray_SimpleNewFromData(2,dims,NPY_UINT16,self->camptr->ImgPtr);
    }
  case FG_COL24:
    {
      cout << "24 bits!" << endl;
      npy_intp dims[3] = {self->camptr->height,self->camptr->width,3};
      return PyArray_SimpleNewFromData(3,dims,NPY_UINT8,self->camptr->ImgPtr);
    }
  default:
    {
      cout << "VideoCapture.get_array(): Unsupported data type!" << endl;
      Py_RETURN_NONE;
    }
  }
}

void _pylm_callback(PyObject *func, double *p, double *hx, int m, int n, int  jacobian) {
    int i;
	PyObject *args		= NULL;
    PyObject *result	= NULL;
	
	// marshall parameters from c -> python
	// construct numpy arrays from c 
	
	npy_intp dims_m[1] = {m};
	npy_intp dims_n[1] = {n};
	
    PyObject *estimate = PyArray_SimpleNewFromData(1, dims_m, PyArray_DOUBLE, p);
    PyObject *measurement = PyArray_SimpleNewFromData(1, dims_n , PyArray_DOUBLE, hx);
	
    args = Py_BuildValue("(OO)", estimate, measurement);
    if (!args) {
        goto cleanup;
    }

    // call func
    result = PyObject_CallObject(func, args);
    if (result == NULL) {
        PyErr_Print();
        goto cleanup;
    }
		
	if(!PyArray_Check(result)){
		PyErr_SetString(PyExc_TypeError, 
						"Return value from callback "
						"should be of numpy array type");
		goto cleanup;		
	}	

    // marshall results from python -> c
    
	npy_intp result_size = PyArray_DIM(result, 0);
	
	if ((!jacobian && (result_size == n)) ||
        (jacobian &&  (result_size == m*n))) {

        for (i = 0; i < result_size; i++) {
            double *j = PyArray_GETPTR1(result, i);
			hx[i] = *j;
		}
		
    } else {
        PyErr_SetString(PyExc_TypeError, 
						"Return value from callback "
                        "should be same size as measurement");
    }

 cleanup:
    
    Py_XDECREF(args);
    Py_XDECREF(estimate);
    Py_XDECREF(measurement);
    Py_XDECREF(result);
    return;
}

/* rle_decode(byte_array, lbrow, lbnpt, mdi) */
static PyObject *rle_decode_py(PyObject *self, PyObject *args)
{
    char *bytes_in=NULL;
    PyArrayObject *npy_array_out=NULL;
    int bytes_in_len;
    npy_intp dims[2];
    int lbrow, lbnpt, npts;
    float mdi;

    if (!PyArg_ParseTuple(args, "s#iif", &bytes_in, &bytes_in_len, &lbrow, &lbnpt, &mdi)) return NULL;

    // Unpacking algorithm accepts an int - so assert that lbrow*lbnpt does not overflow
    if (lbrow > 0 && lbnpt >= INT_MAX / (lbrow+1)) {
        PyErr_SetString(PyExc_ValueError, "Resulting unpacked PP field is larger than PP supports.");
        return NULL;
    } else {
        npts = lbnpt*lbrow;
    }

    // We can't use the macros Py_BEGIN_ALLOW_THREADS / Py_END_ALLOW_THREADS
    // because they declare a new scope block, but we want multiple exits.
    PyThreadState *_save;
    _save = PyEval_SaveThread();

    float *dataout = (float*)calloc(npts, sizeof(float));

    if (dataout == NULL) {
        PyEval_RestoreThread(_save);
        PyErr_SetString(PyExc_ValueError, "Unable to allocate memory for wgdos_unpacking.");
        return NULL;
    }

    function func;  // function is defined by wgdosstuff.
    set_function_name(__func__, &func, 0);
    int status = unpack_ppfield(mdi, (bytes_in_len/BYTES_PER_INT_UNPACK_PPFIELD), bytes_in, LBPACK_RLE_PACKED, npts, dataout, &func);

    /* Raise an exception if there was a problem with the REL algorithm */
    if (status != 0) {
        free(dataout);
        PyEval_RestoreThread(_save);
        PyErr_SetString(PyExc_ValueError, "RLE decode encountered an error.");
        return NULL;
    }
    else {
        /* The data came back fine, so make a Numpy array and return it */
        dims[0]=lbrow;
        dims[1]=lbnpt;
        PyEval_RestoreThread(_save);
        npy_array_out=(PyArrayObject *) PyArray_SimpleNewFromData(2, dims, NPY_FLOAT, dataout);

        if (npy_array_out == NULL) {
            PyErr_SetString(PyExc_ValueError, "Failed to make the numpy array for the packed data.");
            return NULL;
        }

        // give ownership of dataout to the Numpy array - Numpy will then deal with memory cleanup.
        npy_array_out->flags = npy_array_out->flags | NPY_OWNDATA;
        return (PyObject *)npy_array_out;
    }
}

static PyObject *__getattro__(PyObject *self, PyObject *attr_name)
{
	const char *name = PyString_AsString(attr_name);
	pylal_COMPLEX16FrequencySeries *obj = (pylal_COMPLEX16FrequencySeries *) self;

	if(!strcmp(name, "name"))
		return PyString_FromString(obj->series->name);
	if(!strcmp(name, "epoch"))
		return pylal_LIGOTimeGPS_new(obj->series->epoch);
	if(!strcmp(name, "f0"))
		return PyFloat_FromDouble(obj->series->f0);
	if(!strcmp(name, "deltaF"))
		return PyFloat_FromDouble(obj->series->deltaF);
	if(!strcmp(name, "sampleUnits"))
		return pylal_LALUnit_new(0, obj->series->sampleUnits);
	if(!strcmp(name, "data")) {
		npy_intp dims[] = {obj->series->data->length};
		PyObject *array = PyArray_SimpleNewFromData(1, dims, NPY_CDOUBLE, obj->series->data->data);
		if(!array)
			return NULL;
		/* incref self to prevent data from disappearing while
		 * array is still in use, and tell numpy to decref self
		 * when the array is deallocated */
		Py_INCREF(self);
		PyArray_SetBaseObject((PyArrayObject *) array, self);
		return array;
	}
	PyErr_SetString(PyExc_AttributeError, name);
	return NULL;
}

// Function to extract the node information
PyObject *
px_GetNodes(PyObject *self, PyObject *args)
{
    xf_Mesh *Mesh = NULL;
    PyObject *np_Coord;
    npy_intp dims[2];

    // Get the pointer to the xf_Mesh.
    if (!PyArg_ParseTuple(args, "n", &Mesh))
        return NULL;

    // Get dimensions
    dims[0] = Mesh->nNode;
    dims[1] = Mesh->Dim;

    // Error checking
    if (Mesh->nNode <= 0) return NULL;

    // Make the mesh.
    np_Coord = PyArray_SimpleNewFromData( \
                                          2, dims, NPY_DOUBLE, *Mesh->Coord);

    // Output (Dim, nNode, Coord).
    return Py_BuildValue("iiO", Mesh->Dim, Mesh->nNode, np_Coord);
}

// Function to read the BFaceGroup
PyObject *
px_ElemGroup(PyObject *self, PyObject *args)
{
    xf_Mesh *Mesh = NULL;
    xf_ElemGroup *EG = NULL;
    int i, nElem, nNode, QOrder;
    const char *QBasis;
    PyObject *Node;
    npy_intp dims[2];

    // Get the pointer to the xf_BFaceGroup
    if (!PyArg_ParseTuple(args, "ni", &Mesh, &i))
        return NULL;

    // Assign the element group
    EG = Mesh->ElemGroup;

    // Read the data
    nElem  = EG[i].nElem;
    nNode  = EG[i].nNode;
    QOrder = EG[i].QOrder;
    // Convert the enumeration to a string.
    QBasis = xfe_BasisName[EG[i].QBasis];

    // Dimensions of the nodes array
    dims[0] = nElem;
    dims[1] = nNode;
    // Read the nodes to a numpy array.
    Node = PyArray_SimpleNewFromData( \
                                      2, dims, NPY_INT, *EG[i].Node);

    // Output: (nElem, nNode, QOrder, QBasis, Node)
    return Py_BuildValue("iiisO", nElem, nNode, QOrder, QBasis, Node);
}

static PyObject*
get_rgb_screen(PyObject *self, PyObject *args) {

    int x, y, dx, dy;

    int good_args = PyArg_ParseTuple(args, "IIII", &x, &y, &dx, &dy);
    if (!good_args) {
        return NULL;
    }

    Display *display = XOpenDisplay(NULL);
    Window window = RootWindow(display, DefaultScreen(display));
	XImage *img = XGetImage(display, window, x, y, dx, dy, AllPlanes, ZPixmap);

    int w = img->width;
    int h = img->height;

    unsigned char *rgb_arr = (unsigned char*) malloc(4*w*h);
    if (rgb_arr == NULL) {
        PyErr_NoMemory();
        return NULL;
    }

    npy_intp dims[] = {h, w, 4};
    PyObject *np_arr = PyArray_SimpleNewFromData(3, dims, NPY_UINT8, img->data);
    return np_arr;
}

/*
  Export a fff_matrix to a PyArray, and delete it. This function is a
  fff_matrix destructor compatible with any of the following
  constructors: fff_matrix_new and fff_matrix_fromPyArray.
*/ 
PyArrayObject* fff_matrix_toPyArray(fff_matrix* y) 
{
  PyArrayObject* x; 
  size_t size1;
  size_t size2;
  size_t tda;
  npy_intp dims[2]; 
  if (y == NULL) 
    return NULL;
  size1 = y->size1;
  size2 = y->size2;
  tda = y->tda; 

  dims[0] = (npy_intp) size1;  
  dims[1] = (npy_intp) size2;  
  
  /* If the fff_matrix is contiguous and owner, just pass the
     buffer to Python and transfer ownership */ 
  if ((tda == size2) && (y->owner)) {
    x = (PyArrayObject*) PyArray_SimpleNewFromData(2, dims, NPY_DOUBLE, (void*)y->data);
    x->flags = (x->flags) | NPY_OWNDATA; 
  }
  /* Otherwise, create PyArray from scratch. Note, the input
     fff_matrix is necessarily in row-major order. */ 
  else 
    x = fff_matrix_const_toPyArray(y); 
  
  /* Ciao bella */ 
  free(y);
   
  return x;
}

/*
  Export a fff_vector to a PyArray, and delete it. This function is a
  fff_vector destructor compatible with any either fff_vector_new or
  _fff_vector_new_from_buffer.
*/ 
PyArrayObject* fff_vector_toPyArray(fff_vector* y) 
{
  PyArrayObject* x; 
  size_t size;
  npy_intp dims[1]; 
   if (y == NULL) 
    return NULL;
   size = y->size;

  dims[0] = (npy_intp) size;  
 
  /* If the fff_vector is owner (hence contiguous), just pass the
     buffer to Python and transfer ownership */ 
  if (y->owner) {
    x = (PyArrayObject*) PyArray_SimpleNewFromData(1, dims, NPY_DOUBLE, (void*)y->data);
    x->flags = (x->flags) | NPY_OWNDATA; 
  }
  /* Otherwise, create Python array from scratch */ 
  else 
    x = fff_vector_const_toPyArray(y); 
 
  /* Ciao bella */ 
  free(y);

  return x;
}

static PyObject* auto_correlations(PyObject* self, PyObject* args)
{
    oskar_VisBlock* h = 0;
    oskar_Mem* m = 0;
    PyObject *capsule = 0;
    PyArrayObject *array = 0;
    npy_intp dims[4];
    if (!PyArg_ParseTuple(args, "O", &capsule)) return 0;
    if (!(h = (oskar_VisBlock*) get_handle(capsule, name))) return 0;

    /* Check that auto-correlations exist. */
    if (!oskar_vis_block_has_auto_correlations(h))
    {
        PyErr_SetString(PyExc_RuntimeError, "No auto-correlations in block.");
        return 0;
    }

    /* Return an array reference to Python. */
    m = oskar_vis_block_auto_correlations(h);
    dims[0] = oskar_vis_block_num_times(h);
    dims[1] = oskar_vis_block_num_channels(h);
    dims[2] = oskar_vis_block_num_stations(h);
    dims[3] = oskar_vis_block_num_pols(h);
    array = (PyArrayObject*)PyArray_SimpleNewFromData(4, dims,
            (oskar_mem_is_double(m) ? NPY_CDOUBLE : NPY_CFLOAT),
            oskar_mem_void(m));
    return Py_BuildValue("N", array); /* Don't increment refcount. */
}

static PyObject * 
creategrid(PyObject *self, PyObject *args)
{
npy_intp dimensions[1];

int i,j, n;
double** input_params;
double* data;
PyObject *input;
PyArrayObject *input_arr, *result;

if (!PyArg_ParseTuple(args, "O", &input))
	return NULL;
input_arr = (PyArrayObject *)PyArray_ContiguousFromObject(input, PyArray_DOUBLE, 2, 2);

n=input_arr->dimensions[0];

input_params = (double**)malloc(sizeof(double*)*n);
for (i=0; i<n; i++)
	{
	input_params[i] = (double *)malloc(sizeof(double)*3);
	for (j=0; j<3; j++)
		input_params[i][j] = *(double *)(input_arr->data+ i*input_arr->strides[0] + j*input_arr->strides[1]);
	}

int size;
data = creategrid_c(n, input_params, &size);
dimensions[0] = size*n;
result = (PyArrayObject *)PyArray_SimpleNewFromData(1, dimensions, PyArray_DOUBLE, (char*)data);
return PyArray_Return(result);
}

    PyObject* createNumpyArray(const cv::Mat& mat) const
    {
      assert(mat.type() == CV_64F && mat.isContinuous() && mat.data != nullptr);

      std::array<npy_intp, 2> dims {{ mat.rows-1, mat.cols }};
      return PyArray_SimpleNewFromData(dims.size(), std::begin(dims), NPY_FLOAT, reinterpret_cast<void*>(mat.data));
    }

static PyObject *premalloced_npy_double_array(double *data, npy_intp len)
{
    PyObject *premalloced = premalloced_new(data);
    if (!premalloced)
        return NULL;

    npy_intp dims[1] = {len};

    PyArrayObject *out = (PyArrayObject *)
        PyArray_SimpleNewFromData(1, dims, NPY_DOUBLE, data);
    if (!out)
    {
        Py_DECREF(premalloced);
        return NULL;
    }

#ifdef PyArray_BASE
    /* FIXME: PyArray_BASE changed from a macro to a getter function in
     * Numpy 1.7. When we drop Numpy 1.6 support, remove this #ifdef block. */
    PyArray_BASE(out) = premalloced;
#else
    if (PyArray_SetBaseObject(out, premalloced))
    {
        Py_DECREF(out);
        return NULL;
    }
#endif

    return (PyObject *)out;
}

static PyObject *
pixbuf_get_pixels_array(PyObject *self, PyObject *args)
{
    /* 1) read in Python pixbuf, get the underlying gdk_pixbuf */
    PyGObject *py_pixbuf;
    GdkPixbuf *gdk_pixbuf;
    PyArrayObject *array;
    npy_intp dims[3] = { 0, 0, 3 };

    if (!PyArg_ParseTuple(args, "O!:pixbuf_get_pixels_array",
			  &PyGdkPixbuf_Type, &py_pixbuf))
	return NULL;

    gdk_pixbuf = GDK_PIXBUF(py_pixbuf->obj);

    /* 2) same as pygtk/gtk/gdk.c _wrap_gdk_pixbuf_get_pixels_array()
     * with 'self' changed to py_pixbuf
     */

    dims[0] = gdk_pixbuf_get_height(gdk_pixbuf);
    dims[1] = gdk_pixbuf_get_width(gdk_pixbuf);
    if (gdk_pixbuf_get_has_alpha(gdk_pixbuf))
        dims[2] = 4;

    array = (PyArrayObject *)PyArray_SimpleNewFromData(3, dims, NPY_UBYTE,
			     (char *)gdk_pixbuf_get_pixels(gdk_pixbuf));
    if (array == NULL)
        return NULL;

    array->strides[0] = gdk_pixbuf_get_rowstride(gdk_pixbuf);
    /* the array holds a ref to the pixbuf pixels through this wrapper*/
    Py_INCREF(py_pixbuf);
    array->base = (PyObject *)py_pixbuf;
    return PyArray_Return(array);
}

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);
}

PyObject* ICoinPackedMatrix::np_getIndices(){

	npy_intp dims = this->getNumElements();
	_import_array();
	PyObject *Arr = PyArray_SimpleNewFromData( 1, &dims, PyArray_INT32, this->IgetIndices() );
	return Arr;
}

void py_coulomb_set_hubbard_u(PyObject *self, void *p, double *U, int *error)
{
  particles_t *py_p;
  PyObject *py_U, *r;
  npy_intp dims[1];

  INIT_ERROR(error);

#ifdef DEBUG
  printf("[py_coulomb_set_Hubbard_U] %s %p %p\n",
	 PyString_AsString(PyObject_Str(self)), U, error);
#endif

  f_particles_get_tag(p, (void**) &py_p);
  assert(py_p->f90obj == p);

  dims[0] = f_particles_get_nel(py_p->f90obj);
  py_U = PyArray_SimpleNewFromData(1, dims, NPY_DOUBLE, U);

  r = PyObject_CallMethod(self, "set_Hubbard_U", "(OO)", py_p, py_U);

  Py_DECREF(py_U);
  PASS_PYTHON_ERROR(error, r);
  Py_DECREF(r);
}