Sampling/exsample/linear_interpolator.h

0001 // -*- C++ -*-
0002 //
0003 // linear_interpolator.h is part of ExSample -- A Library for Sampling Sudakov-Type Distributions
0004 //
0005 // Copyright (C) 2008-2019 Simon Platzer -- simon.plaetzer@desy.de, The Herwig Collaboration
0006 //
0007 // ExSample is licenced under version 3 of the GPL, see COPYING for details.
0008 // Please respect the MCnet academic guidelines, see GUIDELINES for details.
0009 //
0010 //
0011 #ifndef EXSAMPLE_linear_interpolator_h_included
0012 #define EXSAMPLE_linear_interpolator_h_included
0013
0014 #include "utility.h"
0015
0016 namespace exsample {
0017
0018   /// \brief Exception thrown, if inversion of the interpolation has
0019   /// no solution
0020   struct inversion_has_no_solution { };
0021
0022   /// \brief Exception thrown, if a constant piece of the
0023   /// interpolation has been hit
0024   struct constant_interpolation {
0025
0026     /// standard constructor
0027     constant_interpolation(double xlow,
0028                double xhigh,
0029                double h)
0030       : range(xlow,xhigh), value(h) {}
0031
0032     /// The range in which the interpolation is constant
0033     std::pair<double,double> range;
0034
0035     /// The value of the interpolation in the range
0036     double value;
0037
0038   };
0039
0040   /// \brief A linear interpolation allowing for inversion of the
0041   /// linear interpolation
0042   class linear_interpolator {
0043
0044   public:
0045
0046     /// default constructor
0047     linear_interpolator();
0048
0049     /// construct a linear interpolator given the
0050     /// map of interpolation points and values
0051     explicit linear_interpolator(const std::map<double,double>& points);
0052
0053   public:
0054
0055     /// return the interpolations value at the given point
0056     double operator()(double x) const;
0057
0058     /// return the range of the interpolation
0059     std::pair<double,double> range() const { return range_; }
0060
0061   public:
0062
0063     /// return true, if an inverse exists for
0064     /// the given value
0065     bool invertible(double f) const {
0066       return (range_.first <= f &&
0067           range_.second >= f);
0068     }
0069
0070     /// return the inverse, assuming the inverse
0071     /// is unique
0072     double unique_inverse(double f) const;
0073
0074     /// access the interpolation map
0075     std::map<double,double>& interpolation() {
0076       return interpolation_;
0077     }
0078
0079     /// set the value at the given point
0080     void set_interpolation(double point, double value);
0081
0082     /// reset after interpolation points have been changed
0083     void reset();
0084
0085   public:
0086
0087     /// put to ostream
0088     template<class OStream>
0089     void put(OStream& os) const;
0090
0091     /// get from istream
0092     template<class IStream>
0093     void get(IStream& is);
0094
0095   private:
0096
0097     /// map points to values
0098     std::map<double, double> interpolation_;
0099
0100     /// the range of the interpolation
0101     std::pair<double,double> range_;
0102
0103   };
0104
0105   /// ostream operator for linear_interpolator objects
0106   template<class OStream>
0107   inline OStream& operator<<(OStream& os, const linear_interpolator& ip) {
0108     ip.put(os);
0109     return os;
0110   }
0111
0112   /// istream operator for linear_interpolator objects
0113   template<class IStream>
0114   inline IStream& operator>>(IStream& is, linear_interpolator& ip) {
0115     ip.get(is);
0116     return is;
0117   }
0118
0119 }
0120
0121 #include "linear_interpolator.icc"
0122
0123 #endif // EXSAMPLE_linear_interpolator_h_included