From de35eeabd803f586dcedc0ab79cff4402cb4859c Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Thu, 8 Jan 2009 09:18:03 -0600 Subject: [PATCH] Generic interpolation --- src/interpolate.cpp | 212 ++++++++++++++++++++++++++++++++++++++++++++ src/interpolate.hpp | 44 +++++++++ 2 files changed, 256 insertions(+) create mode 100644 src/interpolate.cpp create mode 100644 src/interpolate.hpp diff --git a/src/interpolate.cpp b/src/interpolate.cpp new file mode 100644 index 0000000..932e3be --- /dev/null +++ b/src/interpolate.cpp @@ -0,0 +1,212 @@ +#include +#include "interpolate.hpp" +#include +#include +#include +#include + +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 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 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; +} + diff --git a/src/interpolate.hpp b/src/interpolate.hpp new file mode 100644 index 0000000..748c6c8 --- /dev/null +++ b/src/interpolate.hpp @@ -0,0 +1,44 @@ +#ifndef __MAK_INTERPOLATE_HPP +#define __MAK_INTERPOLATE_HPP + +#include "config.hpp" +#include +#include +#include +#include +#include + +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 > 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