Herwig/Utilities/Interpolator.tcc

0001 // -*- C++ -*-
0002 //
0003 // Interpolator.tcc 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 //
0010 // This is the implementation of the non-inlined, non-templated member
0011 // functions of the Interpolator class.
0012 //
0013
0014 #include "Interpolator.h"
0015 #include "ThePEG/Interface/ClassDocumentation.h"
0016 #include "ThePEG/Interface/ParVector.h"
0017 #include "ThePEG/Interface/Parameter.h"
0018 #include "ThePEG/Persistency/PersistentOStream.h"
0019 #include "ThePEG/Persistency/PersistentIStream.h"
0020
0021 using namespace Herwig;
0022
0023 template <typename ValT, typename ArgT>
0024 void Interpolator<ValT,ArgT>::persistentOutput(PersistentOStream & os) const {
0025   os << _xval << _fun << _order
0026      << ounit(_funit,TypeTraits<ValT>::baseunit())
0027      << ounit(_xunit,TypeTraits<ArgT>::baseunit());
0028 }
0029
0030 template <typename ValT, typename ArgT>
0031 void Interpolator<ValT,ArgT>::persistentInput(PersistentIStream & is, int) {
0032   is >> _xval >> _fun >> _order
0033      >> iunit(_funit,TypeTraits<ValT>::baseunit())
0034      >> iunit(_xunit,TypeTraits<ArgT>::baseunit());
0035 }
0036
0037 /**
0038  * Macro for Interpolator user classes to use. Only they know what the
0039  * template arguments are going to be.
0040  */
0041 #define HERWIG_INTERPOLATOR_CLASSDESC(Name,ValT,ArgT) \
0042 /**                                                   \
0043 * Register the Interpolator with ThePEG               \
0044 */                                                    \
0045 DescribeClass<Interpolator<ValT,ArgT>,Interfaced>     \
0046 describeHerwigInterpolatorFor##Name("Herwig::Interpolator<"#ValT","#ArgT">","");\
0047
0048
0049
0050 template <typename ValT, typename ArgT>
0051 void Interpolator<ValT,ArgT>::Init() {
0052
0053   static ClassDocumentation<Interpolator<ValT,ArgT> > documentation
0054     ("The Interpolator class is design to interpolate a table of values");
0055
0056   static Parameter<Interpolator<ValT,ArgT>,unsigned int> interfaceOrder
0057     ("Order",
0058      "Order of the interpolation",
0059      &Interpolator::_order, 3, 1, 10,
0060      false, false, Interface::limited);
0061
0062   static ParVector<Interpolator<ValT,ArgT>,double> interfaceXValues
0063     ("XValues",
0064      "The x values for the interpolation",
0065      &Interpolator::_xval, -1, 0., 0, 0,
0066      false, false, Interface::nolimits);
0067
0068   static ParVector<Interpolator<ValT,ArgT>,double> interfaceFunctionValues
0069     ("FunctionValues",
0070      "The function values for the interpolation",
0071      &Interpolator::_fun, -1, 0., 0, 0,
0072      false, false, Interface::nolimits);
0073
0074   static Parameter<Interpolator<ValT,ArgT>,ValT> interfaceValueType
0075     ("ValueType",
0076      "The unit of the function values",
0077      &Interpolator<ValT,ArgT>::_funit,
0078      TypeTraits<ValT>::baseunit(),
0079      1.0*TypeTraits<ValT>::baseunit(),
0080      0*TypeTraits<ValT>::baseunit(),
0081      0*TypeTraits<ValT>::baseunit(),
0082      false, true, Interface::nolimits);
0083
0084   static Parameter<Interpolator<ValT,ArgT>,ArgT> interfaceArgType
0085     ("ArgType",
0086      "The unit of the function arguments",
0087      &Interpolator<ValT,ArgT>::_xunit,
0088      TypeTraits<ArgT>::baseunit(),
0089      1.0*TypeTraits<ArgT>::baseunit(),
0090      0*TypeTraits<ArgT>::baseunit(),
0091      0*TypeTraits<ArgT>::baseunit(),
0092      false, true, Interface::nolimits);
0093
0094 }
0095
0096 template <typename ValT, typename ArgT>
0097 ValT Interpolator<ValT,ArgT>::operator ()(ArgT xpt) const {
0098   const double xpoint = xpt / _xunit;
0099   // size of the vectors
0100   unsigned int isize(_xval.size());
0101   // workout the numer of points we need
0102   unsigned int m(std::min(_order,isize)),mp(m+1),ix,iy;
0103   // search for the point if the function increases
0104   int mid,iupp=isize,ilow=0;
0105   if(_xval[0]>_xval[_xval.size()-1]) {
0106     do {
0107       mid=(iupp+ilow)/2;
0108       if(xpoint>_xval[mid]){iupp=mid;}
0109       else{ilow=mid;}
0110     }
0111     while(iupp-ilow>1);
0112   }
0113   // search for the point if the function decreases
0114   else {
0115     do {
0116       mid=(iupp+ilow)/2;
0117       if(xpoint<_xval[mid]){iupp=mid;}
0118       else{ilow=mid;}
0119     }
0120     while(iupp-ilow>1);
0121   }
0122   // ilow is now the midpoint
0123   mid=ilow;
0124   // copy the re-ordered interpolation points
0125   // number of points
0126   unsigned int npoints(_order+2-_order%2),icopy,j(0);
0127   int iloc(0),i(0);
0128   do {
0129     icopy=mid+iloc;
0130     if(icopy>isize-1) npoints=mp;
0131     else {
0132       _copyx[j]   = _xval[icopy];
0133       _copyfun[j] = _fun [icopy];
0134       ++j;
0135     }
0136     iloc=-iloc;
0137     if(iloc>=0){++iloc;}
0138   }
0139   while(j<npoints);
0140   // do this interpolation
0141   bool extra(npoints!=mp);
0142   for(ix=0;ix<m;++ix) {
0143     if(extra) {
0144       icopy=m-ix-1;
0145       _copyfun[m+1]=(_copyfun[m+1]-_copyfun[m-1])/(_copyx[m+1]-_copyx[icopy]);
0146     }
0147     i=m;
0148     for(iy=ix;iy<m;++iy) {
0149       icopy=i-ix-1;
0150       _copyfun[i]=(_copyfun[i]-_copyfun[i-1])/(_copyx[i]-_copyx[icopy]);
0151       --i;
0152     }
0153   }
0154   double sum(_copyfun[m]);
0155   if(extra) sum=0.5*(sum+_copyfun[m+1]);
0156   i=m-1;
0157   for(ix=0;ix<m;++ix) {
0158     sum=_copyfun[i]+(xpoint-_copyx[i])*sum;
0159     --i;
0160   }
0161   return sum * _funit;
0162 }