Herwig/Utilities/Interpolator.h

0001 // -*- C++ -*-
0002 //
0003 // Interpolator.h is a part of Herwig - A multi-purpose Monte Carlo event generator
0004 // Copyright (C) 2002-2019 The Herwig Collaboration
0005 //
0006 // Herwig is licenced under version 3 of the GPL, see COPYING for details.
0007 // Please respect the MCnet academic guidelines, see GUIDELINES for details.
0008 //
0009 #ifndef HERWIG_Interpolator_H
0010 #define HERWIG_Interpolator_H
0011 //
0012 // This is the declaration of the Interpolator class.
0013 //
0014
0015 #include "ThePEG/Interface/Interfaced.h"
0016 #include <cassert>
0017
0018 namespace Herwig {
0019
0020 using namespace ThePEG;
0021
0022 /** \ingroup Utilities
0023  *  \author Peter Richardson
0024  *
0025  *  This class implments a polynominal interpolation of a table of values, it is
0026  *  based on the interpolation code in FORTRAN HERWIG.
0027  *
0028  */
0029
0030 template <typename ValT, typename ArgT>
0031 class Interpolator: public Interfaced {
0032
0033 public:
0034
0035   /**
0036    *  Pointer to an Interpolator
0037    */
0038   typedef typename Ptr<Interpolator<ValT,ArgT> >::pointer Ptr;
0039
0040   /** @name Standard constructors and destructors. */
0041   //@{
0042   /**
0043    * The default constructor.
0044    */
0045   Interpolator() : _order(3), _copyx(5),_copyfun(5) {}
0046
0047   /**
0048    * Constructor with data as vectors.
0049    */
0050   Interpolator(const vector<ValT> & f,
0051          const vector<ArgT> & x,
0052          unsigned int order)
0053     : _fun(f.size(),0.0),_xval(x.size(),0.0),_order(order),
0054       _funit(TypeTraits<ValT>::baseunit()),
0055       _xunit(TypeTraits<ArgT>::baseunit()),
0056       _copyx(order+2),_copyfun(order+2) {
0057     assert(_order>0);
0058     assert(x.size() == f.size());
0059     for (size_t i = 0; i < f.size(); ++i) {
0060       _fun [i] = f[i] / _funit;
0061       _xval[i] = x[i] / _xunit;
0062     }
0063   }
0064
0065   /**
0066    * Constructor with data as arrays.
0067    */
0068   template <size_t N>
0069   Interpolator(const array<ValT,N> & f,
0070          const array<ArgT,N> & x,
0071          unsigned int order)
0072     : _fun(N,0.0),_xval(N,0.0),_order(order),
0073       _funit(TypeTraits<ValT>::baseunit()),
0074       _xunit(TypeTraits<ArgT>::baseunit()),
0075       _copyx(order+2),_copyfun(order+2) {
0076     assert(_order>0);
0077     for (size_t i = 0; i < N; ++i) {
0078       _fun [i] = f[i] / _funit;
0079       _xval[i] = x[i] / _xunit;
0080     }
0081   }
0082   //@}
0083
0084   /**
0085    * Constructor from bare arrays
0086    */
0087   Interpolator(size_t size,
0088            const double f[], ValT funit,
0089            const double x[], ArgT xunit,
0090            unsigned int order)
0091     : _fun(size,0.0),_xval(size,0.0),_order(order),
0092       _funit(funit),_xunit(xunit), _copyx(order+2),_copyfun(order+2) {
0093     assert(_order>0);
0094     for (size_t i = 0; i < size; ++i) {
0095       _fun [i] = f[i];
0096       _xval[i] = x[i];
0097     }
0098   }
0099   //@}
0100
0101   /**
0102    *  Return the interpolated value
0103    */
0104   ValT operator () (ArgT) const;
0105   /** Return type for GaussianIntegrator */
0106   typedef ValT ValType;
0107   /** Argument type for GaussianIntegrator */
0108   typedef ArgT ArgType;
0109
0110 public:
0111
0112   /** @name Functions used by the persistent I/O system. */
0113   //@{
0114   /**
0115    * Function used to write out object persistently.
0116    * @param os the persistent output stream written to.
0117    */
0118   void persistentOutput(PersistentOStream & os) const;
0119
0120   /**
0121    * Function used to read in object persistently.
0122    * @param is the persistent input stream read from.
0123    * @param version the version number of the object when written.
0124    */
0125   void persistentInput(PersistentIStream & is, int version);
0126   //@}
0127
0128   /**
0129    * The standard Init function used to initialize the interfaces.
0130    * Called exactly once for each class by the class description system
0131    * before the main function starts or
0132    * when this class is dynamically loaded.
0133    */
0134   static void Init();
0135
0136 protected:
0137
0138   /** @name Clone Methods. */
0139   //@{
0140   /**
0141    * Make a simple clone of this object.
0142    * @return a pointer to the new object.
0143    */
0144   virtual IBPtr clone() const { return new_ptr(*this); }
0145
0146   /** Make a clone of this object, possibly modifying the cloned object
0147    * to make it sane.
0148    * @return a pointer to the new object.
0149    */
0150   virtual IBPtr fullclone() const { return new_ptr(*this); }
0151   //@}
0152
0153 private:
0154
0155   /**
0156    * The assignment operator is private and must never be called.
0157    * In fact, it should not even be implemented.
0158    */
0159   Interpolator & operator=(const Interpolator &) = delete;
0160
0161 private:
0162
0163   /**
0164    * the function values.
0165    */
0166   vector<double> _fun;
0167
0168   /**
0169    * The x values.
0170    */
0171   vector<double> _xval;
0172
0173   /**
0174    * the order of interpolation.
0175    */
0176   unsigned int _order;
0177
0178   /**
0179    * The Unit of the function values
0180    */
0181   ValT _funit;
0182
0183   /**
0184    * The Unit of the argument values
0185    */
0186   ArgT _xunit;
0187
0188   /**
0189    *  Temporary storage vector
0190    */
0191   mutable vector<double> _copyx;
0192
0193   /**
0194    *  Temporary storage vector
0195    */
0196   mutable vector<double> _copyfun;
0197
0198 };
0199
0200 /**
0201  * helper function to create InterpolatorPtr easily
0202  * from bare arrays (analogous to make_pair() )
0203  */
0204 template <typename ValT, typename ArgT>
0205 inline typename Interpolator<ValT,ArgT>::Ptr
0206 make_InterpolatorPtr(size_t size,
0207          const double f[], ValT funit,
0208          const double x[], ArgT xunit,
0209          unsigned int order)
0210 {
0211   return new_ptr(Interpolator<ValT,ArgT>(size,
0212            f,funit,
0213            x,xunit,
0214            order));
0215 }
0216
0217 /**
0218  * helper function to create InterpolatorPtr easily
0219  * from vectors (analogous to make_pair() )
0220  */
0221 template <typename ValT, typename ArgT>
0222 inline typename Interpolator<ValT,ArgT>::Ptr
0223 make_InterpolatorPtr(const typename std::vector<ValT> & f,
0224          const typename std::vector<ArgT> & x,
0225          unsigned int order)
0226 {
0227   return new_ptr(Interpolator<ValT,ArgT>(f,x,order));
0228 }
0229
0230 /**
0231  * helper function to create InterpolatorPtr easily
0232  * from arrays (analogous to make_pair() )
0233  */
0234 template <typename ValT, typename ArgT, size_t N>
0235 inline typename Interpolator<ValT,ArgT>::Ptr
0236 make_InterpolatorPtr(const typename std::array<ValT,N> & f,
0237          const typename std::array<ArgT,N> & x,
0238          unsigned int order)
0239 {
0240   return new_ptr(Interpolator<ValT,ArgT>(f,x,order));
0241 }
0242
0243 }
0244
0245 #ifndef ThePEG_TEMPLATES_IN_CC_FILE
0246 #include "Interpolator.tcc"
0247 #endif
0248
0249 #endif /* HERWIG_Interpolator_H */