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:

  MonotCubicInterpolator(const std::string & datafilename) throw (const char*)
  {
    if (!read(datafilename)) {
      throw("Unable to constuct MonotCubicInterpolator from file.") ;
    } ;
  }


 MonotCubicInterpolator(const char* datafilename) throw (const char*)
  {
    if (!read(std::string(datafilename))) {
      throw("Unable to constuct MonotCubicInterpolator from file.") ;
    } ;
  }


   MonotCubicInterpolator(const char* datafilename, int xColumn, int fColumn) throw (const char*)
  {
    if (!read(std::string(datafilename),xColumn,fColumn)) {
      throw("Unable to constuct MonotCubicInterpolator from file.") ;
    } ;
  }

   MonotCubicInterpolator(const std::string & datafilename, int xColumn, int fColumn) throw (const char*)
  {
    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 throw(const char*);

   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 throw(const char*) ;

   std::pair<double,double> getMinimumF() const throw(const char*)  ;


   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) throw(const char*);

   std::pair<double,double> getMissingX() const throw(const char*) ;

   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