Generic interpolation

src/interpolate.cpp

@@ -0,0 +1,212 @@
+#include <algorithm>
+#include "interpolate.hpp"
+#include <vector>
+#include <fstream>
+#include <inttypes.h>
+#include <sstream>
+
+using namespace std;
+using namespace CosmoTool;
+  
+Interpolate::Interpolate(double *xs, double *values, uint32_t N, bool autofree)
+{
+  spline = gsl_spline_alloc (gsl_interp_linear, N);
+  gsl_spline_init (spline, xs, values, N);
+  accel_interp = gsl_interp_accel_alloc();
+
+  this->autoFree = autofree;
+}
+
+Interpolate::~Interpolate()
+{
+  if (spline != 0)
+    gsl_spline_free (spline);
+  if (accel_interp != 0)
+    gsl_interp_accel_free (accel_interp);
+}
+
+double Interpolate::compute(double x)
+  throw (InvalidRangeException)
+{
+  double y;
+  int err = gsl_spline_eval_e(spline, x, accel_interp, &y);
+
+  if (err)
+    throw InvalidRangeException("Interpolate argument outside range");
+}
+
+uint32_t Interpolate::getNumPoints() const
+{
+  return spline->size;
+}
+
+void Interpolate::fillWithXY(double *x, double *y) const
+{
+  if (x != 0)
+    memcpy(x, spline->x, sizeof(double)*spline->size);
+  if (y != 0)
+    memcpy(y, spline->y, sizeof(double)*spline->size);
+}
+
+const Interpolate& Interpolate::operator=(const Interpolate& a)
+{
+  double *x, *y;
+
+  if (spline != NULL)
+    {
+      gsl_spline_free (spline);
+      gsl_interp_accel_free (accel_interp);
+    }
+
+  autoFree = true;
+  spline = gsl_spline_alloc (gsl_interp_linear, a.spline->size);
+  accel_interp = gsl_interp_accel_alloc ();
+  gsl_spline_init(spline, a.spline->x, a.spline->y, a.spline->size);
+}
+
+double Interpolate::getMaxX() const
+{
+  return spline->x[spline->size-1];
+}
+
+typedef struct {
+  double a, b;
+} MyPair;
+
+bool operator<(const MyPair& a, const MyPair& b)
+{
+  return a.a < b.a;
+}
+
+Interpolate CosmoTool::buildFromVector(const InterpolatePairs& v)
+{
+  double *x = new double[v.size()];
+  double *y = new double[v.size()];
+  
+  for (uint32_t i = 0; i < v.size(); i++)
+    {
+      x[i] = v[i].first;
+      y[i] = v[i].second;
+    }
+
+  Interpolate inter = Interpolate(x, y, v.size(), true);
+
+  delete[] x;
+  delete[] y;
+
+  return inter;  
+}
+
+Interpolate CosmoTool::buildInterpolateFromFile(const char *fname)
+  throw (NoSuchFileException)
+{
+  vector<MyPair> allData;
+  ifstream f(fname);
+
+  if (!f)
+    throw NoSuchFileException(fname);
+
+  do 
+    {
+      MyPair m;
+
+      if ((f >> m.a >> m.b).eof())
+	break;
+
+      allData.push_back(m);
+    }
+  while (1);
+
+  sort(allData.begin(), allData.end());
+
+  double *x = new double[allData.size()];
+  double *y = new double[allData.size()];
+  
+  for (uint32_t i = 0; i < allData.size(); i++)
+    {
+      x[i] = allData[i].a;
+      y[i] = allData[i].b;
+    }
+
+  Interpolate inter = Interpolate(x, y, allData.size(), true);
+
+  delete[] x;
+  delete[] y;
+
+  return inter;
+}
+
+Interpolate CosmoTool::buildInterpolateFromColumns(const char *fname, uint32_t col1, uint32_t col2)
+  throw (NoSuchFileException,InvalidRangeException)
+{
+  vector<MyPair> allData;
+  ifstream f(fname);
+
+  if (!f)
+    throw NoSuchFileException(fname);
+  
+  bool swapped = (col1 > col2);
+  uint32_t colMin = min(col1, col2);
+  uint32_t colMax = max(col1, col2);
+
+  do 
+    {
+      MyPair m;
+      string line;
+
+      if (getline(f, line).eof())
+	break;
+
+      istringstream iss(line);
+      double dummy;
+      double val1, val2;
+
+      for (uint32_t i = 0; i < colMin; i++)
+	iss >> dummy;
+      if (!(iss >> val1))
+	throw InvalidRangeException("Invalid first column");
+
+      if (col2 != col1)
+	{
+	  for (uint32_t i = 0; i < (colMax-colMin-1); i++)
+	    iss >> dummy;
+	  if (!(iss >> val2))
+	    throw InvalidRangeException("Invalid second column");
+	}
+      else
+	val2 = val2;
+
+      if (!swapped)
+	{
+	  m.a = val1;
+	  m.b = val2;
+	}
+      else
+	{
+	  m.a = val2;
+	  m.b = val1;
+	}
+
+      allData.push_back(m);
+    }
+  while (1);
+
+  sort(allData.begin(), allData.end());
+
+  double *x = new double[allData.size()];
+  double *y = new double[allData.size()];
+  
+  for (uint32_t i = 0; i < allData.size(); i++)
+    {
+      x[i] = allData[i].a;
+      y[i] = allData[i].b;
+    }
+
+  Interpolate inter = Interpolate(x, y, allData.size(), true);
+
+  delete[] x;
+  delete[] y;
+  
+  return inter;
+}
+

src/interpolate.hpp

@@ -0,0 +1,44 @@
+#ifndef __MAK_INTERPOLATE_HPP
+#define __MAK_INTERPOLATE_HPP
+
+#include "config.hpp"
+#include <inttypes.h>
+#include <gsl/gsl_interp.h>
+#include <gsl/gsl_spline.h>
+#include <vector>
+#include <utility>
+
+namespace CosmoTool
+{
+  
+  class Interpolate
+  {
+  public:
+    Interpolate() : spline(0), accel_interp(0) { }
+    Interpolate(double *xs, double *values, uint32_t N, bool autofree = false);
+    ~Interpolate();
+
+    double compute(double x)
+      throw (InvalidRangeException);
+    const Interpolate& operator=(const Interpolate& a);
+
+    uint32_t getNumPoints() const;
+    void fillWithXY(double *x, double *y) const;
+    double getMaxX() const;
+  protected:
+    gsl_interp_accel *accel_interp;
+    gsl_spline *spline;
+    bool autoFree;
+  };
+
+  typedef std::vector< std::pair<double,double> > InterpolatePairs;
+
+  Interpolate buildInterpolateFromFile(const char *fname)
+    throw (NoSuchFileException);
+  Interpolate buildInterpolateFromColumns(const char *fname, uint32_t col1, uint32_t col2)
+    throw (NoSuchFileException,InvalidRangeException);
+  Interpolate buildFromVector(const InterpolatePairs& v);
+  
+};
+
+#endif