Index: /trunk/anuga_core/source/anuga/fit_interpolate/fitsmooth.c
===================================================================
--- /trunk/anuga_core/source/anuga/fit_interpolate/fitsmooth.c	(revision 8692)
+++ /trunk/anuga_core/source/anuga/fit_interpolate/fitsmooth.c	(revision 8692)
@@ -0,0 +1,1022 @@
+
+#include "Python.h"
+
+#include <stdio.h>   /* gets */
+#include <stdlib.h>  /* atoi, malloc */
+#include <string.h>  /* strcpy */
+#include "numpy/arrayobject.h"
+#include "math.h"
+
+#include "util_ext.h" /* in utilities */
+#include "sparse_dok.h" /*in utilities */
+#include "quad_tree.h"
+#include <netcdf.h>
+#include "omp.h"
+
+// Errors defined for netcdf reading
+#define ERRCODE 2
+#define ERR(e) {printf("Error: %s\n", nc_strerror(e)); exit(ERRCODE);}
+
+//------------------------ NET CDF READING ----------------------------
+
+// Note Padarn 06/12/12: I am removing the use of these functions from
+// the workflow of the new code. I think they should be reintroduced, but
+// to make sure all the unit tests work properly more checking of filetypes
+// etc needs to be done.
+// Note Padarn 12/12/12: These have been moved back into the workflow, 
+// but still there is no way of handling .csv or .txt files efficiently.
+
+
+// Reads a block of a netcdf file into the arrays 'point_coordiantes' and 'point_values'. The
+// block is specified by a start point (np_start) and and an end point (np_end)
+// NOTE PADARN: FIX - only reads elevation at the moment
+int _read_net_cdf_points_block(char * filename, char * attribute_name, int np_start,int np_end,double *point_coordinates,double *point_values){
+
+    //point_coordinates and _values are intitialised outside
+    int np=np_end-np_start+1;
+    int status;                       /* error status */
+    int ncid;                          /* netCDF ID */
+    int rh_id;                         /* variable ID */
+    size_t start_points[] = {np_start, 0}; /* start at first value */
+    size_t count_points[] = {np, 2};
+    size_t start[] = {np_start};
+    size_t count[] = {np};
+    
+    if ((status = nc_open(filename, NC_NOWRITE, &ncid))) ERR(status);
+
+    if ((status = nc_inq_varid(ncid, "points", &rh_id))) ERR(status);
+
+    /* read coordinates from netCDF variable */
+    if ((status = nc_get_vara_double(ncid, rh_id, start_points, count_points, point_coordinates))) ERR(status);
+
+    if ((status = nc_inq_varid(ncid, attribute_name, &rh_id))) ERR(status);
+
+    /* read values from netCDF variable */
+    if ((status = nc_get_vara_double(ncid, rh_id, start, count, point_values))) ERR(status);
+
+    if ((status = nc_close(ncid))) ERR(status);
+
+    return 0;
+
+}
+
+// Takes a netcdf file and gets the number of points stored in it (also the number of values for
+// attributes specified)
+// NOTE PADARN: FIX - uses 'points' to get points, this may be a problem.
+int _read_net_cdf_entries(char * filename,int * x){
+
+    int status;                       /* error status */
+    int ncid;
+    int p_id;
+    int ndims;
+    size_t length = {0};
+
+    if ((status = nc_open(filename, NC_NOWRITE, &ncid))) ERR(status);
+
+    if ((status = nc_inq_varid(ncid, "points", &p_id))) ERR(status);
+    
+    if((status = nc_inq_varndims(ncid, p_id, &ndims))) ERR(status);
+    
+    int * ndimsize = malloc(sizeof(int)*ndims);
+
+    if((status = nc_inq_vardimid(ncid, p_id, ndimsize))) ERR(status);
+
+    if((status = nc_inq_dimlen(ncid, ndimsize[0], &length))) ERR(status);
+
+    *x = (int)length;
+
+    if ((status = nc_close(ncid))) ERR(status);
+
+    return 0;
+
+}
+
+//--------------------------------------------------------------------------
+
+//-------------------------- QUANTITY FITTING ------------------------------
+
+
+// Builds the matrix D used to smooth the interpolation 
+// of a variables from scattered data points to a mesh. See fit.py for more details.s
+int _build_smoothing_matrix(int n,
+                      long* triangles,
+        		      double* areas,
+                      double* vertex_coordinates,
+                      int* strides,
+                      sparse_dok * smoothing_mat)
+		      {
+
+
+    int k;
+    int k3,k6;
+    int err = 0;
+    edge_key_t key;
+
+    double det,area,x0,x1,x2,y0,y1,y2;
+    double a0,b0,a1,b1,a2,b2,e01,e12,e20;
+    int v0,v1,v2;
+    double smoothing_val;
+
+    
+    for(k=0; k<n; k++) {
+        // multiple k by 3 to give index for triangles which store in 3-blocks
+        k3=k*3;
+        k6=k*6;
+        // store the area for the current triangle
+        area = areas[k];
+        // store current triangles global vertex indicies
+        v0 = triangles[k3];
+        v1 = triangles[k3+1];
+        v2 = triangles[k3+2];
+        // store the locations of the three verticies
+        x0 = vertex_coordinates[k6];
+        y0 = vertex_coordinates[k6+1];
+        x1 = vertex_coordinates[k6+2];
+        y1 = vertex_coordinates[k6+3];
+        x2 = vertex_coordinates[k6+4];
+        y2 = vertex_coordinates[k6+5];
+
+        // calculate gradients (move to external function?)
+
+        det = (y2-y0)*(x1-x0) - (y1-y0)*(x2-x0);
+        a0 = (y2-y0)*(0-1) - (y1-y0)*(0-1);
+        a0 /= det;
+
+        b0 = (x1-x0)*(0-1) - (x2-x0)*(0-1);
+        b0 /= det;
+
+        a1 = (y2-y0)*(1-0) - (y1-y0)*(0-0);
+        a1 /= det;
+
+        b1 = (x1-x0)*(0-0) - (x2-x0)*(1-0);
+        b1 /= det;
+
+        a2 = (y2-y0)*(0-0) - (y1-y0)*(1-0);
+        a2 /= det;
+
+        b2 = (x1-x0)*(1-0) - (x2-x0)*(0-0);
+        b2 /= det;
+        
+        // insert diagonal contributions
+
+        // v0,v0
+        key.i = v0;
+        key.j = v0;
+        smoothing_val = (a0*a0 + b0*b0)*area;
+        add_dok_entry(smoothing_mat,key,smoothing_val);
+
+        // v1,v1
+        key.i = v1;
+        key.j = v1;
+        smoothing_val = (a1*a1 + b1*b1)*area;
+        add_dok_entry(smoothing_mat,key,smoothing_val);
+
+        // v2,v2
+        key.i = v2;
+        key.j = v2;
+        smoothing_val = (a2*a2 + b2*b2)*area;
+        add_dok_entry(smoothing_mat,key,smoothing_val);
+
+
+        // insert off diagonal contributions
+        e01 = (a0*a1 + b0*b1)*area;
+        // v0,v1 (v1,v0)
+        key.i = v0;
+        key.j = v1;
+        add_dok_entry(smoothing_mat,key,e01);
+        key.i = v1;
+        key.j = v0;
+        add_dok_entry(smoothing_mat,key,e01);
+
+        e12 = (a1*a2 + b1*b2)*area;
+        // v1,v2 (v2,v1)
+        key.i = v1;
+        key.j = v2;
+        add_dok_entry(smoothing_mat,key,e12);
+        key.i = v2;
+        key.j = v1;
+        add_dok_entry(smoothing_mat,key,e12);
+
+        e20 = (a2*a0 + b2*b0)*area;
+        // v2,v0 (v0,v2)
+        key.i = v2;
+        key.j = v0;
+        add_dok_entry(smoothing_mat,key,e20);
+        key.i = v0;
+        key.j = v2;
+        add_dok_entry(smoothing_mat,key,e20);
+    }
+
+    return err;
+
+}
+
+// Builds a quad tree out of a list of triangles for quick 
+// searching. 
+quad_tree * _build_quad_tree(int n,
+                      long* triangles,
+                      double* vertex_coordinates,
+                      double* extents)               
+{   
+    
+    int k,k6;
+    double x0,y0,x1,y1,x2,y2;
+
+    // set up quad tree and allocate memory
+    quad_tree * tree = new_quad_tree(extents[0],extents[1],extents[2],extents[3]);
+    
+    // iterate through triangles
+    for(k=0; k<n; k++) {
+        // multiple k by 3 to give index for triangles which store in 3-blocks
+        k6=k*6;
+        // store the locations of the three verticies
+        x0 = vertex_coordinates[k6];
+        y0 = vertex_coordinates[k6 + 1];
+        x1 = vertex_coordinates[k6 + 2];
+        y1 = vertex_coordinates[k6 + 3];
+        x2 = vertex_coordinates[k6 + 4];
+        y2 = vertex_coordinates[k6 + 5];
+        triangle * T = new_triangle(k,x0,y0,x1,y1,x2,y2);
+        quad_tree_insert_triangle(tree,T);
+    }
+  
+    // return pointer to new tree struct
+    return tree;
+    
+}
+
+// Builds the AtA and Atz interpolation matrix
+// and residual. Uses a quad_tree for fast access to the triangles of the mesh.
+// This function takes a filename and an attribute name, in netcdf format .pts, 
+// and reads the attribute value and point location from this file.
+int _build_matrix_AtA_Atz_fileread(int N, long * triangles,
+                      char * filename,
+                      char * attribute_name,
+                      int blocksize,
+                      sparse_dok * AtA,
+                      double * Atz,quad_tree * quadtree)
+              {
+
+
+    int np;
+    _read_net_cdf_entries(filename,&np);
+    //n=10;
+
+    double * point_coordinates = malloc(2*sizeof(double)*blocksize);
+    double * point_values = malloc(sizeof(double)*blocksize);
+    int k;
+    int left;
+
+    int i;
+    for(i=0;i<N;i++){
+        Atz[i]=0;
+    } 
+
+    edge_key_t key;
+    int w;
+    int start = 0;
+    left = np;
+
+    
+    while(left > 0){
+        if(left-blocksize<=0){
+            blocksize = left;
+        }
+
+        // Read data into arrays
+        _read_net_cdf_points_block(filename,attribute_name,start,
+        start+blocksize-1,point_coordinates,point_values);
+        
+        start = start + blocksize;
+
+            #pragma omp parallel for private(k,i,key,w)
+            for(k=0;k<blocksize;k++){
+
+
+                double x = point_coordinates[2*k];
+                double y = point_coordinates[2*k+1];
+                triangle * T = search(quadtree,x,y);
+
+                if(T!=NULL){
+                    double * sigma = calculate_sigma(T,x,y);
+                    int js[3];
+                    for(i=0;i<3;i++){
+                        js[i]=triangles[3*(T->index)+i];
+                    }
+                    
+                    #pragma omp critical
+                    { 
+                    for(i=0;i<3;i++){
+
+                        Atz[js[i]] += sigma[i]*point_values[k];
+
+                        for(w=0;w<3;w++){
+                            
+                            key.i=js[i];
+                            key.j=js[w];
+
+                            add_dok_entry(AtA,key,sigma[i]*sigma[w]);
+                        }                        
+                    }
+                    }
+                    free(sigma);
+                    sigma=NULL;
+                }
+
+            } 
+        
+        left=left-blocksize;
+    }
+    
+    return 0;
+}
+
+// Builds the AtA and Atz interpolation matrix
+// and residual. Uses a quad_tree for fast access to the triangles of the mesh.
+// This function takes a list of point coordinates, and associated point values
+// (for any number of attributes).
+int _build_matrix_AtA_Atz_points(int N, long * triangles,
+                      double * point_coordinates, double * point_values,
+                      int zdims, int npts,
+                      sparse_dok * AtA,
+                      double ** Atz,quad_tree * quadtree)
+              {
+
+
+
+    int k;
+
+
+
+    int i,w;
+    for(w=0;w<zdims;w++){
+        for(i=0;i<N;i++){
+            Atz[w][i]=0;
+        } 
+    }
+
+    edge_key_t key;
+
+
+
+
+    //#pragma omp parallel for private(k,i,key,w)
+    for(k=0;k<npts;k++){
+
+
+        double x = point_coordinates[2*k];
+        double y = point_coordinates[2*k+1];
+        triangle * T = search(quadtree,x,y);
+
+        if(T!=NULL){
+            double * sigma = calculate_sigma(T,x,y);
+            int js[3];
+            for(i=0;i<3;i++){
+                js[i]=triangles[3*(T->index)+i];
+            }
+            
+            #pragma omp critical
+            { 
+            for(i=0;i<3;i++){
+
+               for(w=0;w<zdims;w++){
+                    Atz[w][js[i]] += sigma[i]*point_values[zdims*k+w];
+               }
+               
+               for(w=0;w<3;w++){
+                    
+                    key.i=js[i];
+                    key.j=js[w];
+
+                   add_dok_entry(AtA,key,sigma[i]*sigma[w]);
+                }                        
+            }
+            }
+            free(sigma);
+            sigma=NULL;
+
+       } 
+    }
+
+    return 0;
+}
+
+// Combines two sparse_dok matricies and two vectors of doubles. 
+void _combine_partial_AtA_Atz(sparse_dok * dok_AtA1,sparse_dok * dok_AtA2,
+                             double* Atz1,
+                             double* Atz2,
+                             int n, int zdim){
+
+    add_sparse_dok(dok_AtA1,1,dok_AtA2,1);
+
+    int i;
+    for(i=0;i<n*zdim;i++){
+        Atz1[i]+=Atz2[i];
+    }
+}
+
+// --------------------------- Utilities -----------------------------------
+
+// Converts a double array into a PyList object for return to python.
+static PyObject *c_double_array_to_list(double * mat,int cols){
+
+    int j;
+
+    PyObject *lst;
+
+    lst = PyList_New(cols);
+    if (!lst) return NULL;
+    for (j = 0; j < cols; j++) {
+        PyObject *num = PyFloat_FromDouble(mat[j]);
+        if (!num) {
+            Py_DECREF(lst);
+            return NULL;
+        }
+        PyList_SET_ITEM(lst, j, num);   // reference to num stolen
+    }
+ 
+    return lst;
+}
+
+// Converts a int array into a PyList object for return to python.
+static PyObject *c_int_array_to_list(int * mat,int cols){
+
+    int j;
+
+    PyObject *lst;
+
+    lst = PyList_New(cols);
+    if (!lst) return NULL;
+    for (j = 0; j < cols; j++) {
+        PyObject *num = PyInt_FromLong(mat[j]);
+        if (!num) {
+            Py_DECREF(lst);
+            return NULL;
+        }
+        PyList_SET_ITEM(lst, j, num);   // reference to num stolen
+    }
+ 
+    return lst;
+}
+
+// ----------------------------------------------------------------------------
+
+// --------------------- Capsule destructor methods ---------------------------
+
+// PADARN NOTE 11/12/12: The following two 'destructor' functions will free the 
+// memory allocated to a quad_tree or sparse_dok. An error should be thrown 
+// if the pointer is already Null, as this shouldn't happen.
+
+// Delete capsule containing a quad tree - name of capsule must be exactly
+// "quad tree". 
+void delete_quad_tree_cap(PyObject * cap){
+    quad_tree * kill = (quad_tree*) PyCapsule_GetPointer(cap,"quad tree");
+    if(kill!=NULL){
+        delete_quad_tree(kill);
+    } else{
+    }
+}
+
+// Delete capsule containing a sparse_dok - name of capsule must be exactly
+// "sparse dok".
+void delete_dok_cap(PyObject *cap){
+
+    sparse_dok * kill = (sparse_dok*) PyCapsule_GetPointer(cap,"sparse dok");
+    if(kill!=NULL){
+        delete_dok_matrix(kill);
+    }
+
+}
+
+//----------------------- PYTHON WRAPPER FUNCTION -----------------------------
+
+// Parses python information about triangles
+// and vertex coordiantes to build a quad tree for efficient location of points in
+// the triangle mesh.
+PyObject *build_quad_tree(PyObject *self, PyObject *args) {
+
+
+
+    PyArrayObject *triangles;
+    PyArrayObject *vertex_coordinates;
+    PyArrayObject *extents;
+
+
+    // Convert Python arguments to C
+    if (!PyArg_ParseTuple(args, "OOO",
+                                            &triangles,
+                                            &vertex_coordinates,
+                                            &extents
+                                            )) {
+      PyErr_SetString(PyExc_RuntimeError,
+              "fitsmooth.c: could not parse input");
+      return NULL;
+    }
+
+    CHECK_C_CONTIG(triangles);
+    CHECK_C_CONTIG(vertex_coordinates);
+    CHECK_C_CONTIG(extents);
+
+    int n = triangles->dimensions[0];
+
+    return  PyCapsule_New((void*) _build_quad_tree(n,
+                      (long*) triangles -> data,
+                      (double*) vertex_coordinates -> data,
+                      (double*) extents -> data),
+                      "quad tree",
+                      &delete_quad_tree_cap); 
+
+}
+
+// Builds the smoothing matrix D. Parses 
+// input mesh information from python and then builds a sparse_dok, returning
+// a capsule object wrapping a pointer to this struct.
+PyObject *build_smoothing_matrix(PyObject *self, PyObject *args) {
+
+
+    PyArrayObject *triangles;
+    PyArrayObject *areas;
+    PyArrayObject *vertex_coordinates;
+    int err;
+
+    // Convert Python arguments to C
+    if (!PyArg_ParseTuple(args, "OOO" ,
+                                            &triangles,
+                                            &areas,
+                                            &vertex_coordinates
+                                            )) {
+      PyErr_SetString(PyExc_RuntimeError,
+              "fitsmooth.c: could not parse input");
+      return NULL;
+    }
+
+    CHECK_C_CONTIG(triangles);
+    CHECK_C_CONTIG(areas);
+    CHECK_C_CONTIG(vertex_coordinates);
+
+    int n = triangles->dimensions[0];
+
+
+    sparse_dok * smoothing_mat; // Should be an input argument?
+    smoothing_mat = make_dok();
+
+    err = _build_smoothing_matrix(n,
+                      (long*) triangles  -> data,
+                      (double*) areas -> data,
+                      (double*) vertex_coordinates -> data,
+                      (int *) vertex_coordinates -> strides,
+                      smoothing_mat);
+
+   
+    if (err != 0) {
+      PyErr_SetString(PyExc_RuntimeError,
+              "Unknown Error");
+      return NULL;
+    }
+
+    return  PyCapsule_New((void*) smoothing_mat,
+                  "sparse dok",
+                  &delete_dok_cap); 
+    
+
+}
+
+// Read a data file (netcdf format) to build
+// AtA and Atz. Returns a pointer to the sparse dok matrix AtA wrapped in a capsule
+// object and a python list for the array Atz.
+PyObject *build_matrix_AtA_Atz_fileread(PyObject *self, PyObject *args) {
+
+
+    PyArrayObject *triangles;
+    char * file_name;
+    char * attribute_name;
+    int block_size;
+    int N; // Number of triangles
+    int err;
+    PyObject *tree;
+
+    // Convert Python arguments to C
+    if (!PyArg_ParseTuple(args, "OiOssi",&tree, &N,
+                                            &triangles,
+                                            &file_name,
+                                            &attribute_name,
+                                            &block_size
+                                            )) {
+      PyErr_SetString(PyExc_RuntimeError,
+              "fitsmooth.c: could not parse input");
+      return NULL;
+    }
+
+    
+
+    CHECK_C_CONTIG(triangles);
+
+    quad_tree * quadtree = (quad_tree*) PyCapsule_GetPointer(tree,"quad tree");
+
+
+    sparse_dok * dok_AtA; // Should be an input argument?
+    dok_AtA = make_dok();
+    double * Atz = malloc(sizeof(double)*N);
+
+    err = _build_matrix_AtA_Atz_fileread(N,(long*) triangles->data,
+                      file_name,
+                      attribute_name,
+                      block_size,
+                      dok_AtA,
+                      Atz,
+                      quadtree);
+
+    if (err != 0) {
+      PyErr_SetString(PyExc_RuntimeError,
+              "Unknown Error");
+      return NULL;
+    }
+
+    PyObject * AtA_cap =  PyCapsule_New((void*) dok_AtA,
+                  "sparse dok",
+                  &delete_dok_cap);
+    PyObject *Atz_ret = c_double_array_to_list(Atz,N); 
+
+    int npts;
+    err = _read_net_cdf_entries(file_name, &npts);
+
+    PyObject *lst = PyList_New(3);
+    PyList_SET_ITEM(lst, 0, AtA_cap);
+    PyList_SET_ITEM(lst, 1, Atz_ret);
+    PyList_SET_ITEM(lst, 2, PyInt_FromLong((long)npts));
+    return lst;
+
+}
+
+// Builds AtA and Atz directly from a list of 
+// points and values Returns a pointer to the sparse dok matrix AtA wrapped in 
+// a capsule object and a python list for the array Atz.
+PyObject *build_matrix_AtA_Atz_points(PyObject *self, PyObject *args) {
+
+
+    PyArrayObject *triangles;
+    PyArrayObject *point_coordinates;
+    PyArrayObject *z;
+    int N; // Number of triangles
+    int err;
+    int npts;
+    int zdims;
+    PyObject *tree;
+
+    // Convert Python arguments to C
+    if (!PyArg_ParseTuple(args, "OiOOOii",&tree, &N,
+                                            &triangles,
+                                            &point_coordinates,
+                                            &z,
+                                            &zdims,
+                                            &npts
+                                            )) {
+      PyErr_SetString(PyExc_RuntimeError,
+              "fitsmooth.c: could not parse input");
+      return NULL;
+    }
+
+    
+    CHECK_C_CONTIG(triangles);
+    CHECK_C_CONTIG(point_coordinates);
+    CHECK_C_CONTIG(z);
+
+    quad_tree * quadtree = (quad_tree*) PyCapsule_GetPointer(tree,"quad tree");
+
+
+    sparse_dok * dok_AtA; // Should be an input argument?
+    dok_AtA = make_dok();
+    double ** Atz = malloc(sizeof(double*)*zdims);
+    int i;
+    for(i=0;i<zdims;i++){
+        Atz[i] = malloc(sizeof(double)*N);
+    }
+
+    err = _build_matrix_AtA_Atz_points(N,(long*) triangles->data,
+                      (double*) point_coordinates->data,
+                      (double*) z->data,
+                      zdims,
+                      npts,
+                      dok_AtA,
+                      Atz,
+                      quadtree);
+
+
+    if (err != 0) {
+      PyErr_SetString(PyExc_RuntimeError,
+              "Unknown Error");
+      return NULL;
+    }
+
+    PyObject * AtA_cap =  PyCapsule_New((void*) dok_AtA,
+                  "sparse dok",
+                  &delete_dok_cap); 
+    PyObject * Atz_ret = PyList_New(zdims);
+    for(i=0;i<zdims;i++){
+        PyList_SET_ITEM(Atz_ret,i,c_double_array_to_list(Atz[i],N));
+        free(Atz[i]);
+    }
+    free(Atz);
+
+    PyObject *lst = PyList_New(2);
+    PyList_SET_ITEM(lst, 0, AtA_cap);
+    PyList_SET_ITEM(lst, 1, Atz_ret);
+    return lst;
+
+}
+
+// Searches a quad tree struct for the triangle containing a given point,
+// returns the sigma values produced by this point and the triangle found,
+// and the triangle index. Found is returned as 0 if no triangle is found
+// and 1 otherwise.
+PyObject *individual_tree_search(PyObject *self, PyObject *args) {
+
+    // Setting up variables to parse input
+    PyObject *tree;
+    PyArrayObject *point;
+
+    // Convert Python arguments to C
+    if (!PyArg_ParseTuple(args, "OO",&tree, &point
+                                            )) {
+      PyErr_SetString(PyExc_RuntimeError,
+              "fitsmooth.c: could not parse input");
+      return NULL;
+    }
+
+    CHECK_C_CONTIG(point);
+
+    quad_tree * quadtree = (quad_tree*) PyCapsule_GetPointer(tree,"quad tree");
+    double xp,yp;
+    if(PyArray_TYPE(point)==7){
+        int *pointd = (int*)point->data;
+        xp = (double)(pointd[0]);
+        yp = (double)(pointd[1]);
+    }
+    else{
+        double *pointd = (double*)point->data;
+        xp = pointd[0];
+        yp = pointd[1];
+    }
+
+    triangle * T = search(quadtree,xp,yp);
+    PyObject * sigmalist = PyList_New(3);
+    long found;
+    long index;
+    
+    if(T!=NULL){
+            double * sigma = calculate_sigma(T,xp,yp);
+            sigmalist = c_double_array_to_list(sigma,3);
+            free(sigma);
+            found = 1;
+            index = (long)T->index;
+    }else{
+            double sigma[3];
+            sigma[0]=sigma[1]=sigma[2]=-1;
+            sigmalist = c_double_array_to_list(sigma,3);
+            index = -10;
+            found = 0;
+    }
+    PyObject * retlist = PyList_New(3);
+    PyList_SET_ITEM(retlist,0,PyInt_FromLong(found));
+    PyList_SET_ITEM(retlist,1,sigmalist);
+    PyList_SET_ITEM(retlist,2,PyInt_FromLong(index));
+    return retlist;
+}
+
+// Returns the total number of triangles stored in quad_tree.
+// Takes a capsule object holding a pointer to the quad_tree as input.
+//
+// PADARN NOTE: This function was only added for the purpose of passing unit
+// tests. Perhaps this, and other tree functions, should be moved to another 
+// file to provide a more comprehensive python interface to the tree structure. 
+PyObject *items_in_tree(PyObject *self, PyObject *args) {
+
+    // Setting up variables to parse input
+    PyObject *tree; // capsule to hold quad_tree pointer
+
+    // Convert Python arguments to C
+    if (!PyArg_ParseTuple(args, "O",&tree
+                                            )) {
+      PyErr_SetString(PyExc_RuntimeError,
+              "fitsmooth.c: could not parse input");
+      return NULL;
+    }
+
+    // Extract quad_tree pointer from capsule
+    quad_tree * quadtree = (quad_tree*) PyCapsule_GetPointer(tree,"quad tree");
+    
+    // Return the number of elements in the tree (stored in struct)
+    return PyInt_FromLong((long)quadtree->count);
+    
+}
+
+// Returns the total number of nodes stored in quad_tree.
+// Takes a capsule object holding a pointer to the quad_tree as input. 
+PyObject *nodes_in_tree(PyObject *self, PyObject *args) {
+
+    // Setting up variables to parse input
+    PyObject *tree; // capsule to hold quad_tree pointer
+
+    // Convert Python arguments to C
+    if (!PyArg_ParseTuple(args, "O",&tree
+                                            )) {
+      PyErr_SetString(PyExc_RuntimeError,
+              "fitsmooth.c: could not parse input");
+      return NULL;
+    }
+
+    // Extract quad_tree pointer from capsule
+    quad_tree * quadtree = (quad_tree*) PyCapsule_GetPointer(tree,"quad tree");
+    
+    // Return the number of elements in the tree (stored in struct)
+    return PyInt_FromLong((long)quad_tree_node_count(quadtree));
+    
+}
+
+// Combines two sparse_dok and two double arrays together
+// which represent partial completion of the AtA and Atz matrices, when they are build
+// in parts due to points being read in blocks in python. Result is stored in the first
+// sparse_dok and double array.
+// Function takes as arguments two capsule objects holding pointers to the sparse_dok
+// structs, and two numpy double arrays. Aslo the size of the double arrays, and the number
+// of columns (variables) for the array Atz are passed as arguments.
+//
+// PADARN NOTE: Blocking the points in python is far slower than reading them directly
+// in c and bypassing the need for this.
+PyObject *combine_partial_AtA_Atz(PyObject *self, PyObject *args) {
+
+    // Setting up variables to parse input
+    PyObject *AtA_cap1, *AtA_cap2;
+    PyArrayObject *Atz1, *Atz2;
+    int n,zdim; // number of nodes
+
+    // Convert Python arguments to C
+    if (!PyArg_ParseTuple(args, "OOOOii",&AtA_cap1, &AtA_cap2,
+                                        &Atz1, &Atz2, &zdim, &n
+                                            )) {
+      PyErr_SetString(PyExc_RuntimeError,
+              "fitsmooth.c: could not parse input");
+      return NULL;
+    }
+
+    // Get pointers to sparse_dok objects from capsules
+    sparse_dok * dok_AtA1 = (sparse_dok*) PyCapsule_GetPointer(AtA_cap1,"sparse dok");
+    sparse_dok * dok_AtA2 = (sparse_dok*) PyCapsule_GetPointer(AtA_cap2,"sparse dok");
+
+    // Combine the partial AtA and Atz
+    _combine_partial_AtA_Atz(dok_AtA1,dok_AtA2,
+                             (double*) Atz1->data,
+                             (double*) Atz2->data,
+                             n, zdim);
+
+    // Return nothing interesting
+    return Py_BuildValue("");
+}
+
+// Converts a sparse_dok matrix to a full non-compressed matrix expressed
+// as a list of lists (python). Takes as input a capsule object containing a pointer to the
+// sparse_dok object. Also takes an integer n as input, specifying the (n x n) size of the 
+// matrix. Matrix is assumed square.
+//
+// PADARN NOTE: This function does not seem particularly sensible, but was required for the 
+// purposes of passing unit tests. 
+PyObject *return_full_D(PyObject *self, PyObject *args) {
+
+    // Setting up variables to parse input
+    PyObject *D_cap; // capsule object holding sparse_dok pointer
+    int n; // number of matrix columns/rows
+
+    // Convert Python arguments to C
+    if (!PyArg_ParseTuple(args, "Oi",&D_cap, &n
+                                            )) {
+      PyErr_SetString(PyExc_RuntimeError,
+              "fitsmooth.return_full_D: could not parse input");
+      return NULL;
+    }
+
+    // Get pointer to spars_dok struct
+    sparse_dok * D_mat = (sparse_dok*) PyCapsule_GetPointer(D_cap,"sparse dok");
+    
+    // Check to make sure that the specified size of the matrix is at least as big
+    // as the entries stored in the sparse_dok.
+    if (D_mat->num_rows>n || get_dok_rows(D_mat)>n){
+        PyErr_SetString(PyExc_RuntimeError,
+              "fitsmooth.return_full_D: sparse_dok is bigger than size specified for return.");
+        return NULL;
+    }
+
+    // Build new python list containing the full matrix to return
+    PyObject *ret_D = PyList_New(n);
+    int i,j;
+    edge_key_t key;
+    edge_t *s;
+    for(i=0;i<n;i++)
+    {
+        PyObject *temp = PyList_New(n);
+        for(j=0;j<n;j++)
+        {
+            key.i=i;
+            key.j=j;
+            s = find_dok_entry(D_mat,key);
+            if (s){
+                PyList_SET_ITEM(temp,j,PyFloat_FromDouble(s->entry));
+            }else{
+                PyList_SET_ITEM(temp,j,PyFloat_FromDouble(0));
+            }
+        PyList_SET_ITEM(ret_D,i,temp);
+        }
+    }
+
+    return ret_D;
+}
+
+// Takes as input two capsule objects corresponding to the two matricies
+// D and AtA, along with a regularization coefficient alpha.
+// Returns three matricies corresponding to the CSR format storage of the matrix resulting
+// from the computation B = AtA + alpha * D.
+//
+// Capsule objects are not freed and are left alone to be freed as they naturally go out of 
+// scope in Python. The temporary CSR matrix created is cleaned up.
+PyObject *build_matrix_B(PyObject *self, PyObject *args) {
+
+
+    // Setting up variables to parse input
+    double alpha; // Regularization coefficient
+    PyObject *smoothing_mat_cap; // Capsule storing D pointer (sparse_dok struct)
+    PyObject *AtA_cap; // Capsule storing AtA pointer (sparse_dok struct)
+
+    // Convert Python arguments to C
+    if (!PyArg_ParseTuple(args, "OOd",&smoothing_mat_cap, &AtA_cap, &alpha
+                                            )) {
+      PyErr_SetString(PyExc_RuntimeError,
+              "fitsmooth.build_matrix_B: could not parse input");
+      return NULL;
+    }
+
+    // Extract pointers to c structs from capsules
+    sparse_dok * smoothing_mat = (sparse_dok*) PyCapsule_GetPointer(smoothing_mat_cap,"sparse dok");
+    sparse_dok * dok_AtA = (sparse_dok*) PyCapsule_GetPointer(AtA_cap,"sparse dok");
+
+    // Add two sparse_dok matrices
+    add_sparse_dok(smoothing_mat,alpha,dok_AtA,1);
+    
+    // Create sparse_csr matrix and convert result to this format
+    sparse_csr * B;
+    B = make_csr();
+    convert_to_csr_ptr(B,smoothing_mat);
+    
+    // Extract the sparse_csr data to be returned as python lists
+    PyObject *data = c_double_array_to_list(B->data,
+                                B->num_entries);
+    PyObject *colind = c_int_array_to_list(B->colind,
+                                B->num_entries);
+    PyObject *row_ptr = c_int_array_to_list(B->row_ptr,
+                                B->num_rows);
+
+    // Build python list for return
+    PyObject *lst = PyList_New(3);
+    PyList_SET_ITEM(lst, 0, data);
+    PyList_SET_ITEM(lst, 1, colind);
+    PyList_SET_ITEM(lst, 2, row_ptr);
+
+    // Clean up
+    delete_csr_matrix(B);
+    
+    return lst;
+   
+
+
+}
+//------------------------------------------------------------------------------
+
+
+// ------------------------------ PYTHON GLUE ----------------------------------
+
+//==============================================================================
+// Structures to allow calling from python
+//==============================================================================
+
+// Method table for python module
+static struct PyMethodDef MethodTable[] = {
+    {"items_in_tree",items_in_tree, METH_VARARGS, "Print out"},
+    {"nodes_in_tree",nodes_in_tree, METH_VARARGS, "Print out"},
+    {"return_full_D",return_full_D, METH_VARARGS, "Print out"},
+	{"build_matrix_B",build_matrix_B, METH_VARARGS, "Print out"},
+    {"build_quad_tree",build_quad_tree, METH_VARARGS, "Print out"},
+    {"build_smoothing_matrix",build_smoothing_matrix, METH_VARARGS, "Print out"},
+    {"build_matrix_AtA_Atz_fileread",build_matrix_AtA_Atz_fileread, METH_VARARGS, "Print out"},
+    {"build_matrix_AtA_Atz_points",build_matrix_AtA_Atz_points, METH_VARARGS, "Print out"},
+    {"combine_partial_AtA_Atz",combine_partial_AtA_Atz, METH_VARARGS, "Print out"},
+    {"individual_tree_search",individual_tree_search, METH_VARARGS, "Print out"},
+	{NULL, NULL, 0, NULL}   // sentinel
+};
+
+
+// Module initialisation
+void initfitsmooth(void){
+  Py_InitModule("fitsmooth", MethodTable);
+  import_array(); // Necessary for handling of NumPY structures
+
+}
+
+// --------------------------------------------------------------------------------
