MonotCubicInterpolator.hpp

Go to the documentation of this file.

/* -*-C++-*- */
 
#ifndef _MONOTCUBICINTERPOLATOR_H
#define _MONOTCUBICINTERPOLATOR_H
 
#include <vector>
#include <map>
#include <string>
 
/*
  MonotCubicInterpolator
  Copyright (C) 2006 Statoil ASA
 
  This program is free software; you can redistribute it and/or
  modify it under the terms of the GNU General Public License
  as published by the Free Software Foundation; either version 2
  of the License, or (at your option) any later version.
 
  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
 
  You should have received a copy of the GNU General Public License
  along with this program; if not, write to the Free Software
  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
*/
 
namespace Opm
{
 
class MonotCubicInterpolator {
 public:
 
  explicit MonotCubicInterpolator(const std::string & datafilename)
  {
    if (!read(datafilename)) {
      throw("Unable to constuct MonotCubicInterpolator from file.") ;
    } ;
  }
 
 
 explicit MonotCubicInterpolator(const char* datafilename)
  {
    if (!read(std::string(datafilename))) {
      throw("Unable to constuct MonotCubicInterpolator from file.") ;
    } ;
  }
 
 
   MonotCubicInterpolator(const char* datafilename, int xColumn, int fColumn)
  {
    if (!read(std::string(datafilename),xColumn,fColumn)) {
      throw("Unable to constuct MonotCubicInterpolator from file.") ;
    } ;
  }
 
   MonotCubicInterpolator(const std::string & datafilename, int xColumn, int fColumn)
  {
    if (!read(datafilename,xColumn,fColumn)) {
      throw("Unable to constuct MonotCubicInterpolator from file.") ;
    } ;
  }
 
     MonotCubicInterpolator(const std::vector<double> & x ,
                            const std::vector<double> & f);
 
   MonotCubicInterpolator() { }
 
 
 
  bool read(const std::string & datafilename) {
    return read(datafilename,1,2) ;
  }
 
  bool read(const std::string &  datafilename, int xColumn, int fColumn) ;
 
 
 
  double operator () (double x) const { return evaluate(x) ; }
 
   double evaluate(double x) const;
 
   double evaluate(double x, double & errorestimate_output ) const ;
 
   std::pair<double,double> getMinimumX() const {
       // Easy since the data is sorted on x:
       return *data.begin();
   }
 
   std::pair<double,double> getMaximumX() const {
       // Easy since the data is sorted on x:
       return *data.rbegin();
   }
 
   std::pair<double,double> getMaximumF() const ;
 
   std::pair<double,double> getMinimumF() const  ;
 
 
   std::vector<double> get_xVector() const ;
 
   std::vector<double> get_fVector() const ;
 
   void scaleData(double factor);
 
   bool isStrictlyMonotone() {
 
       /* Use cached value if it can be trusted */
       if (strictlyMonotoneCached) {
           return strictlyMonotone;
       }
       else {
           computeInternalFunctionData();
           return strictlyMonotone;
       }
   }
 
   bool isMonotone() const {
       if (monotoneCached) {
           return monotone;
       }
       else {
           computeInternalFunctionData();
           return monotone;
       }
   }
   bool isStrictlyIncreasing() {
 
       /* Use cached value if it can be trusted */
       if (strictlyMonotoneCached) {
           return (strictlyMonotone && strictlyIncreasing);
       }
       else {
           computeInternalFunctionData();
           return (strictlyMonotone && strictlyIncreasing);
       }
   }
 
   bool isMonotoneIncreasing() const {
       if (monotoneCached) {
           return (monotone && increasing);
       }
       else {
           computeInternalFunctionData();
           return (monotone && increasing);
       }
   }
   bool isStrictlyDecreasing() {
 
       /* Use cached value if it can be trusted */
       if (strictlyMonotoneCached) {
           return (strictlyMonotone && strictlyDecreasing);
       }
       else {
           computeInternalFunctionData();
           return (strictlyMonotone && strictlyDecreasing);
       }
   }
 
   bool isMonotoneDecreasing() const {
       if (monotoneCached) {
           return (monotone && decreasing);
       }
       else {
           computeInternalFunctionData();
           return (monotone && decreasing);
       }
   }
 
 
 
   void addPair(double newx, double newf);
 
   std::pair<double,double> getMissingX() const;
 
   std::string toString() const;
 
   int getSize() const {
       return data.size();
   }
 
    void chopFlatEndpoints(const double);
 
    void chopFlatEndpoints() {
        chopFlatEndpoints(1e-14);
    }
 
    void shrinkFlatAreas(const double);
 
    void shrinkFlatAreas() {
        shrinkFlatAreas(1e-14);
    }
 
 
 
private:
 
   // Data structure to store x- and f-values
   std::map<double, double> data;
 
   // Data structure to store x- and d-values
   mutable std::map<double, double> ddata;
 
 
   // Storage containers for precomputed interpolation data
   //   std::vector<double> dvalues; // derivatives in Hermite interpolation.
 
   // Flag to determine whether the boolean strictlyMonotone can be
   // trusted.
   mutable bool strictlyMonotoneCached;
   mutable bool monotoneCached; /* only monotone, not stricly montone */
 
   mutable bool strictlyMonotone;
   mutable bool monotone;
 
   // if strictlyMonotone is true (and can be trusted), the two next are meaningful
   mutable bool strictlyDecreasing;
   mutable bool strictlyIncreasing;
   mutable bool decreasing;
   mutable bool increasing;
 
 
   /* Hermite basis functions, t \in [0,1] ,
      notation from:
      http://en.wikipedia.org/w/index.php?title=Cubic_Hermite_spline&oldid=84495502
   */
 
   double H00(double t) const {
       return 2*t*t*t - 3*t*t + 1;
   }
   double H10(double t) const {
       return t*t*t - 2*t*t + t;
   }
   double H01(double t) const {
       return -2*t*t*t + 3*t*t;
   }
   double H11(double t) const {
       return t*t*t - t*t;
   }
 
 
   void computeInternalFunctionData() const ;
 
   void computeSimpleDerivatives() const ;
 
 
  void adjustDerivativesForMonotoneness() const ;
 
   bool isMonotoneCoeff(double alpha, double beta) const {
       if ((alpha*alpha + beta*beta) <= 9) {
         return true;
       } else {
         return false;
       }
   }
 
};
 
 
} // namespace Opm
 
#endif