include/Professor/Ipol.h

0001 #ifndef PROF_IPOL_H
0002 #define PROF_IPOL_H
0003
0004 #include "Professor/ParamPoints.h"
0005 #include <string>
0006 #include <vector>
0007 #include <sstream>
0008 #include <iostream>
0009 #include <stdexcept>
0010
0011 namespace Professor {
0012
0013
0014   /// Throwable error
0015   struct IpolError : public std::runtime_error {
0016     IpolError(const std::string& reason) : std::runtime_error(reason) { }
0017   };
0018
0019
0020   /// @name Calculator functions for parameterisation elements
0021   //@{
0022
0023   /// Calculate the number of coefficients for the given parameter space dimension and polynomial order
0024   /// @todo Deal in uints
0025   int calcnumCoeffs(int dim, int order);
0026
0027   /// Calculate parametrisation coefficients
0028   /// @note structure is the pre-calculated algebraic structure of the polynomial
0029   /// @todo Provide an (unscaled) arrays-only version with a vector<vector<double>> in place of ParamPoints
0030   std::vector<double> calcCoeffs(const ParamPoints& pts, const std::vector<double>& vals, int order,
0031                                  double threshold, const std::vector<std::vector<int> >& structure);
0032
0033   /// Calculate an interpolated value
0034   double calcValue(const std::vector<double>& params,
0035                    const std::vector<double>& coeffs, int order,
0036                    const std::vector<std::vector<int> >& structure);
0037
0038   /// Calculate an interpolated value
0039   double calcValue(const std::vector<double>& paramslongvector,
0040                    const std::vector<double>& coeffs);
0041
0042   /// @brief Make the algebraic coefficient structure
0043   ///
0044   /// The structure is a nested vector of ints, with each entry in the outer vector
0045   /// representing a term in the polynomial and the inner vectors being the powers
0046   /// to which each param should be raised in that term.
0047   ///
0048   /// @note In a given param space, the structure can be reused between any values
0049   ///   or errors parameterised at the same polynomial order. So pre-computing
0050   ///   structures at 3rd, 4th and 5th order would cover everything needed...
0051   std::vector< std::vector<int> > mkStructure(int dim, int order);
0052
0053   /// Make the vector of polynomial terms to which the coeffs are to be applied, at the given order
0054   ///
0055   /// @note The same long vector can be used in any parameterised value calculation
0056   ///    in the "active" parameter space, between values & errors and multiple bins.
0057   ///    Just need to "dot" it with the appropriate fixed coefficient vector.
0058   std::vector<double> mkLongVector(const std::vector<double>& params, int order,
0059                                    const std::vector< std::vector<int> >& structure);
0060
0061   // vector<double> mkLongVectorDerivative(const vector<double>& params, int order,
0062   //                                       const vector<double>& minPV, const vector<double>& maxPV,
0063   //                                       const vector<vector<int> >& structure);
0064
0065   // vector<double> mkLongVectorGradient(const vector<double>& params, int coord, int order,
0066   //                                     const vector<double>& minPV, const vector<double>& maxPV,
0067   //                                     const vector<vector<int> >& structure);
0068
0069   //@}
0070
0071
0072
0073   /// The heart of Professor: the interpolation of a single numerical value through the parameter space
0074   class Ipol {
0075   public:
0076
0077     /// @brief Constructor for calculation of coefficients
0078     ///
0079     /// The @a pts list of N-dimensional parameter points must correspond to the @a ptvals
0080     /// list of values at those points, to be interpolated at the given polynomial @a order.
0081     /// A name may optionally be given.
0082     ///
0083     /// @note Expert settings: The stability of the SVD operation is controlled
0084     /// by the @a svdthreshold parameter, which should not normally be
0085     /// touched. The stability is normally ensured by internally scaling
0086     /// parameter points into unit ranges within the sampled hypercube defined
0087     /// by @a pts; changing @doscaling to false will disable this scaling, which
0088     /// simplifies Ipol I/O (no PMin/Max metadata is needed) but risks SVD
0089     /// instability.
0090     ///
0091     Ipol(const ParamPoints& pts, const std::vector<double>& ptvals, int order,
0092          const std::string& name="", double svdthreshold=1e-20, bool doscaling=true) {
0093       _dim = pts.dim();
0094       _order = order;
0095       _name = name;
0096       _structure = mkStructure(_dim, _order);
0097       if (doscaling) {
0098         _minPV = pts.ptmins();
0099         _maxPV = pts.ptmaxs();
0100       }
0101       _coeffs = calcCoeffs(pts, ptvals, _order, svdthreshold, _structure);
0102     };
0103
0104     /// Constructor to read ipol from file (one string for each object)
0105     /// @todo Also allow optional passing of pmins, pmaxs vectors for the case where the string includes scaling?
0106     Ipol(const std::string& s) {
0107       fromString(s);
0108     };
0109
0110     ~Ipol() {
0111       _structure.clear();
0112       _coeffs.clear();
0113       if (_minPV.size()) _minPV.clear();
0114       if (_maxPV.size()) _maxPV.clear();
0115     };
0116
0117
0118     /// @name String representations
0119     //@{
0120
0121     /// Get a persistent string representation of this Ipol
0122     std::string toString(const std::string& name="") const;
0123
0124     /// Read and set coefficients (name), order from string
0125     /// @todo Also allow optional passing of pmins, pmaxs vectors for the case where the string includes scaling?
0126     void fromString(const std::string& s);
0127
0128     /// Get a string representation of the mathematical expression
0129     std::string exprString() const;
0130
0131     //@}
0132
0133
0134     /// @name Calculations
0135     //@{
0136
0137     /// Get the value of the parametrisation at point p
0138     double value(const std::vector<double>& p) const;
0139
0140     /// Get the value of the derivative of the parametrisation at point p
0141     /// @todo Expose as a standalone calcDerivative function, cf. calcValue
0142     double derivative(const std::vector<double>& p) const;
0143
0144     /// Get the gradient of the parametrisation at point p
0145     /// @todo Expose as a standalone calcGradient function, cf. calcValue
0146     std::vector<double> gradient(const std::vector<double>& p) const;
0147
0148     //@}
0149
0150
0151     /// @name Coefficient access
0152     //@{
0153
0154     /// Get the number of coefficients
0155     const size_t numCoeffs() const {
0156       // ::numCoeffs(dim(),order())
0157       return _coeffs.size();
0158     }
0159
0160     /// Get the vector of coefficients by const reference
0161     const std::vector<double>& coeffs() const { return _coeffs; }
0162
0163     /// Get a single coefficient
0164     const double& coeff(size_t i) const { return coeffs()[i]; }
0165
0166     //@}
0167
0168
0169     /// @name Polynomial term structure
0170     //@{
0171
0172     /// Get the polynomial term exponent structure
0173     const std::vector< std::vector<int> >& structure() const { return _structure; }
0174
0175     /// Convert params to scaled params
0176     std::vector<double> sparams(const std::vector<double>& params) const;
0177
0178     /// Get a long vector of polynomial terms (sans coefficients) for the given param point
0179     std::vector<double> longVector(const std::vector<double>& params) const {
0180       return mkLongVector(sparams(params), order(), structure());
0181     }
0182
0183     //@}
0184
0185
0186     /// @name Basic ipol properties
0187     //@{
0188
0189     /// Accessors to the dimension of the param points
0190     int dim() const { return _dim; }
0191     int numParams() const { return dim(); }
0192
0193     /// Get the order of the parametrisation
0194     int order() const { return _order; }
0195
0196     /// Get the name of the parametrised object
0197     std::string name() const { return _name; }
0198
0199     //@}
0200
0201
0202     /// @name Limit-setting
0203     //@{
0204
0205     void setParamLimits(const std::vector<double>& minpvs, const std::vector<double>& maxpvs) {
0206       setMinParamVals(minpvs);
0207       setMaxParamVals(maxpvs);
0208     }
0209
0210     const std::vector<double>& minParamVals() { return _minPV; }
0211     const std::vector<double>& maxParamVals() { return _maxPV; }
0212
0213     // The if statements here guarantee backwards compatibility with the format 'binned 2' where the
0214     // ranges were read from the ipol meta. In binned 3 they are read for each ipol from the coefficient
0215     // output
0216     void setMinParamVals(const std::vector<double>& minpvs) { if (_minPV.empty()) _minPV = minpvs; }
0217     void setMaxParamVals(const std::vector<double>& maxpvs) { if (_maxPV.empty()) _maxPV = maxpvs; }
0218
0219     //@}
0220
0221
0222   private:
0223
0224     int _dim, _order;
0225     std::vector<std::vector<int> > _structure;
0226     std::string _name;
0227     std::vector<double> _coeffs, _minPV, _maxPV;
0228
0229   };
0230
0231
0232 }
0233
0234 #endif