Sampling/exsample/linear_interpolator.icc

0001 // -*- C++ -*-
0002 //
0003 // linear_interpolator.icc 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 namespace exsample {
0012
0013   inline linear_interpolator::linear_interpolator()
0014     : interpolation_(), range_() {}
0015
0016   inline linear_interpolator::linear_interpolator(const std::map<double,double>& points)
0017     : interpolation_(points) {
0018     reset();
0019   }
0020
0021
0022   inline void linear_interpolator::set_interpolation(double point, double value) {
0023     interpolation_[point] = value;
0024     if (value > range_.second)
0025       range_.second = value;
0026     if(value < range_.first)
0027       range_.first = value;
0028   }
0029
0030
0031   inline void linear_interpolator::reset() {
0032     range_.first = interpolation_.begin()->second;
0033     range_.second = interpolation_.begin()->second;
0034     for (std::map<double,double>::const_iterator c = interpolation_.begin();
0035      c != interpolation_.end(); ++c) {
0036       if (c->second < range_.first)
0037     range_.first = c->second;
0038       if (c->second > range_.second)
0039     range_.second = c->second;
0040     }
0041   }
0042
0043
0044   inline double linear_interpolator::operator()(double x) const {
0045     std::map<double, double>::const_iterator upper =
0046       interpolation_.upper_bound(x);
0047     if (upper == interpolation_.end()) {
0048       upper = interpolation_.upper_bound(x-1e-10);
0049     }
0050     if (upper == interpolation_.end()) {
0051       upper = interpolation_.upper_bound(x+1e-10);
0052     }
0053     assert(upper != interpolation_.begin() &&
0054        upper != interpolation_.end());
0055     return ((upper->second-std::prev(upper)->second)*x +
0056         std::prev(upper)->second*upper->first - upper->second*std::prev(upper)->first)/
0057       (upper->first - std::prev(upper)->first);
0058   }
0059
0060
0061   inline double linear_interpolator::unique_inverse(double f) const {
0062     if(!invertible(f)) throw inversion_has_no_solution();
0063     std::map<double, double>::const_iterator lower = interpolation_.begin();
0064     bool gotone = false;
0065     for (; lower != --interpolation_.end(); ++lower)
0066       if ((lower->second >= f && std::next(lower)->second <= f) ||
0067       (lower->second <= f && std::next(lower)->second >= f)) {
0068     gotone = true;
0069     break;
0070       }
0071     if(!gotone) throw inversion_has_no_solution();
0072     if (lower->second == std::next(lower)->second) {
0073       throw constant_interpolation(lower->first,std::next(lower)->first,lower->second);
0074     }
0075     double xdiff = std::next(lower)->first - lower->first;
0076     double wdiff = std::next(lower)->second - lower->second;
0077     double woffset = lower->second * std::next(lower)->first - std::next(lower)->second * lower->first;
0078     return (xdiff/wdiff)*(f - woffset/xdiff);
0079   }
0080
0081
0082   template<class OStream>
0083   void linear_interpolator::put(OStream& os) const {
0084     os << interpolation_.size();
0085     ostream_traits<OStream>::separator(os);
0086     for (std::map<double, double>::const_iterator p
0087        = interpolation_.begin(); p != interpolation_.end(); ++p) {
0088       os << p->first;
0089       ostream_traits<OStream>::separator(os);
0090       os << p->second;
0091       ostream_traits<OStream>::separator(os);
0092     }
0093     os << range_.first;
0094     ostream_traits<OStream>::separator(os);
0095     os << range_.second;
0096     ostream_traits<OStream>::separator(os);
0097   }
0098
0099   template<class IStream>
0100   void linear_interpolator::get(IStream& is) {
0101     std::size_t size;
0102     is >> size;
0103     std::pair<double,double> point;
0104     for (std::size_t k = 0; k < size; ++k) {
0105       is >> point.first >> point.second;
0106       interpolation_.insert(point);
0107     }
0108     is >> range_.first >> range_.second;
0109   }
0110
0111 }