Optimize continuum

Author: kormang

README.md

@@ -4,13 +4,14 @@
 ./startlab.sh
 ```
 
-## Python depenencies
+## Python dependencies
 ```
 spectral
 scipy
 numpy
 matplotlib
 sklearn
+numba
 ```
 
 ## To build C module use

ccontinuum.c

@@ -1,13 +1,13 @@
 #include "ccontinuum.h"
 
 #include <stdio.h>
+#include <math.h>
 
 // Used to represent wavelengths in case we switch
 // to algorithm that includes wavelengths.
 #define WL(x) (wavelengths[x])
 
-void continuum(double input[], double output[], double wavelengths[], size_t n) {
-
+void continuum_old(double input[], double output[], double wavelengths[], size_t n) {
 	output[0] = input[0];
 	// i points to the last point that belongs to the curve.
 	// j points to the current potential point.
@@ -44,9 +44,74 @@ void continuum(double input[], double output[], double wavelengths[], size_t n)
 	}
 }
 
-void continuum_removed(double input[], double output[], double wavelengths[], size_t n) {
-	continuum(input, output, wavelengths, n);
+void continuum_removed(double input[], double output[], double wavelengths[], size_t indices[], size_t n) {
+	continuum(input, output, wavelengths, indices, n);
 	for (size_t i = 0; i < n; ++i) {
 		output[i] = input[i] / output[i];
 	}
 }
+
+struct data_t {
+	double* spectrum;
+	double* wavelengths;
+	size_t* indices;
+	size_t n;
+	size_t ind_fill;
+};
+
+static void find_indices(struct data_t* data, size_t ibegin, size_t iend) {
+	double* spectrum = data->spectrum;
+	double* wavelengths = data->wavelengths;
+	size_t iendi = iend - 1;
+	double naxis_y = wavelengths[iendi] - wavelengths[ibegin];
+	double naxis_x = spectrum[ibegin] - spectrum[iendi];
+	double maxval = -INFINITY;
+	double imax = ibegin;
+
+	for(size_t i = ibegin; i < iendi; ++i) {
+		double newval = wavelengths[i] * naxis_x + spectrum[i] * naxis_y;
+		if (newval > maxval) {
+			maxval = newval;
+			imax = i;
+		}
+	}
+
+	if (imax == ibegin) {
+		return;
+	}
+
+	if (imax > ibegin + 1) {
+		find_indices(data, ibegin, imax + 1);
+	}
+
+	data->indices[data->ind_fill++] = imax;
+
+	if (imax < iend - 2) {
+		find_indices(data, imax, iend);
+	}
+}
+
+void continuum(double spectrum[], double output[], double wavelengths[], size_t indices[], size_t n) {
+	struct data_t data = { .spectrum = spectrum, .wavelengths = wavelengths, .indices = indices, .n = n, .ind_fill = 1 };
+
+	// Find indices of points that belong to convex hull.
+	indices[0] = 0;
+	find_indices(&data, 0, n);
+	indices[data.ind_fill] = n - 1; // Didn't increase ind_fill on purpose.
+
+	// Linear interpolation of points.
+	for (size_t i = 0; i < data.ind_fill; ++i) {
+		size_t ibegin = indices[i];
+		size_t iend = indices[i + 1];
+		// Put exact values where possible.
+		output[ibegin] = spectrum[ibegin];
+		// Calculate line parameters.
+		double a = (spectrum[iend] - spectrum[ibegin]) / (wavelengths[iend] - wavelengths[ibegin]);
+		double b = spectrum[ibegin] - a * wavelengths[ibegin];
+		// Fill.
+		for (size_t j = ibegin + 1; j < iend; ++j) {
+			output[j] = a * wavelengths[j] + b;
+		}
+	}
+	output[n - 1] = spectrum[n - 1];
+}

ccontinuum.h

@@ -1,5 +1,6 @@
 
 #include <stddef.h>
 
-void continuum(double input[], double output[], double wavelength[], size_t n);
-void continuum_removed(double input[], double output[], double wavelength[], size_t n);
+void continuum_old(double input[], double output[], double wavelength[], size_t n);
+void continuum_removed(double input[], double output[], double wavelength[], size_t indices[], size_t n);
+void continuum(double spectrum[], double output[], double wavelengths[], size_t indices[], size_t n);

ccontinuummodule.c

@@ -23,7 +23,7 @@
 static PyObject *ErrorObject;
 
 /*
-  Check that PyArrayObject is a double (Float) type and a 1d, or 2d array.
+  Check that PyArrayObject is a double (Float) type and a 1d, or 2d, or 3d array.
   Return 1 if an error and raise exception.
 */
 static int check_type(PyArrayObject* a, int min_dim, int max_dim)  {
@@ -41,11 +41,15 @@ static int check_type(PyArrayObject* a, int min_dim, int max_dim)  {
 */
 static int check_same_shapes(PyArrayObject* a, PyArrayObject* b) {
   if (PyArray_NDIM(a) != PyArray_NDIM(b)) {
+    PyErr_SetString(PyExc_ValueError,
+	"In check_same_shapes: arrays must have same number of dimensions.");
     return 1;
   }
   Py_ssize_t ndims = PyArray_NDIM(a);
   for (int i = 0; i < ndims; ++i) {
     if (PyArray_SHAPE(a)[i] != PyArray_SHAPE(b)[i]) {
+      PyErr_SetString(PyExc_ValueError,
+        "In check_same_shapes: arrays must have same number of dimensions.");
       return 1;
     }
   }
@@ -55,7 +59,7 @@ static int check_same_shapes(PyArrayObject* a, PyArrayObject* b) {
 
 static PyObject *
 continuum_generic_impl(PyObject *self, PyObject *args,
-  void (*continuum_processing_f)(double*, double*, double*, size_t))
+  void (*continuum_processing_f)(double*, double*, double*, size_t*, size_t))
 {
     PyArrayObject* ain;
     PyArrayObject* aout;
@@ -69,8 +73,8 @@ continuum_generic_impl(PyObject *self, PyObject *args,
     if (NULL == aout)  return Py_None;
     if (NULL == awl)  return Py_None;
 
-    if (check_type(ain, 1, 2)) return Py_None;
-    if (check_type(aout, 1, 2)) return Py_None;
+    if (check_type(ain, 1, 3)) return Py_None;
+    if (check_type(aout, 1, 3)) return Py_None;
     if (check_type(awl, 1, 1)) return Py_None;
 
     if (check_same_shapes(ain, aout)) return Py_None;
@@ -86,8 +90,10 @@ continuum_generic_impl(PyObject *self, PyObject *args,
               "In continuum_generic_impl: wavelengths array has incorrect length.");
           return Py_None;
       }
-      continuum_processing_f(datain, dataout, dataawl, spectrum_length);
-    } else {
+      size_t* indices = malloc(sizeof(size_t) * spectrum_length);
+      continuum_processing_f(datain, dataout, dataawl, indices, spectrum_length);
+      free(indices);
+    } else if (PyArray_NDIM(ain) == 2) {
       Py_ssize_t num_spectra = PyArray_SHAPE(ain)[0];
       Py_ssize_t spectrum_length = PyArray_SHAPE(ain)[1];
       if (spectrum_length != PyArray_SHAPE(awl)[0]) {
@@ -96,14 +102,49 @@ continuum_generic_impl(PyObject *self, PyObject *args,
           return Py_None;
       }
 
-      //#pragma omp parallel for
+      #pragma omp parallel
+      {
+      size_t* indices = malloc(sizeof(size_t) * spectrum_length);
+
+      #pragma omp for
       for (Py_ssize_t i = 0; i < num_spectra; ++i) {
         double* datain = (double*)(PyArray_DATA(ain)
           + i * PyArray_STRIDES(ain)[0]);
         double* dataout = (double*)(PyArray_DATA(aout)
           + i * PyArray_STRIDES(aout)[0]);
-        continuum_processing_f(datain, dataout, dataawl, spectrum_length);
+        continuum_processing_f(datain, dataout, dataawl, indices, spectrum_length);
+      }
+
+      free(indices);
+      } // pragma omp parallel
+    } else {
+      Py_ssize_t num_rows = PyArray_SHAPE(ain)[0];
+      Py_ssize_t num_cols = PyArray_SHAPE(ain)[1];
+      Py_ssize_t spectrum_length = PyArray_SHAPE(ain)[2];
+      if (spectrum_length != PyArray_SHAPE(awl)[0]) {
+          PyErr_SetString(PyExc_ValueError,
+              "In continuum_generic_impl: wavelengths array has incorrect length.");
+          return Py_None;
       }
+
+      //TODO: This could be optimized for better work sharing.
+      #pragma omp parallel
+      {
+      size_t* indices = malloc(sizeof(size_t) * spectrum_length);
+
+      #pragma omp for
+      for (Py_ssize_t i = 0; i < num_rows; ++i) {
+        for (Py_ssize_t j = 0; j < num_cols; ++j) {
+	  double* datain = (double*)(PyArray_DATA(ain)
+            + i * PyArray_STRIDES(ain)[0] + j * PyArray_STRIDES(ain)[1]);
+          double* dataout = (double*)(PyArray_DATA(aout)
+            + i * PyArray_STRIDES(aout)[0] + j * PyArray_STRIDES(aout)[1]);
+          continuum_processing_f(datain, dataout, dataawl, indices, spectrum_length);
+	}
+      }
+
+      free(indices);
+      } // pragma omp parallel
     }
 
     return Py_None;
@@ -114,7 +155,7 @@ continuum_generic_impl(PyObject *self, PyObject *args,
 PyDoc_STRVAR(ccontinuum_continuum_doc,
 "continuum(spectrum)\n\
 \n\
-Return continuum of the spectrum.");
+Return continuum of the spectrum or each spectrum in image.");
 
 static PyObject *
 ccontinuum_continuum(PyObject *self, PyObject *args)
@@ -127,7 +168,7 @@ ccontinuum_continuum(PyObject *self, PyObject *args)
 PyDoc_STRVAR(ccontinuum_continuum_removed_doc,
 "continuum_removed(spectrum)\n\
 \n\
-Return continuum removed spectrum.");
+Return continuum removed spectrum or image.");
 
 static PyObject *
 ccontinuum_continuum_removed(PyObject *self, PyObject *args)