Interpolator.cxx

Go to the documentation of this file.

#include "AtlfastEvent/Interpolator.h"
#include <cmath>

namespace Atlfast{
  
  Interpolator::Interpolator(string filename){

    int level = 1;
    int jump = 0;
    vector<int> npoints;
    int nline = 0;
    
    // Make a filestream from the input filename to pass to Setup
    ifstream ifs(filename.c_str());
    Setup(ifs,level,jump,npoints,nline);
    
    // Fill m_limits with the limits in each dimension of all
    // held values
    findLimits(m_limits);

    // Initialise default return values in all dimensions to zero
    setDefault(0);
    
    // Check to see if number of level names (if used)
    // is equal to the number of levels
    int nlevelnames = m_levelnames.size();
    int nlevels = npoints.size();
    if ( nlevelnames && nlevelnames != nlevels ){
      std::cout << nlevelnames << " level identifiers given but "
                << nlevels << " levels detected, exiting" << std::endl;
      exit(1);  
    }
    
    m_contBounds = false;
    
  }
  
  Interpolator::Interpolator(double value, ifstream &ifs, int &level, 
                             int &jump, vector<int> &npoints, int &nline):
    m_value(value){
    ++level;
    Setup(ifs,level,jump,npoints,nline);
  }

  void Interpolator::dump(){
    if (m_interpolators.size()){
      vector<Interpolator>::iterator viitr = m_interpolators.begin();
      for (;viitr!=m_interpolators.end();++viitr){
        std::cout << "Value = " << (*viitr).getValue() << std::endl;
        (*viitr).dump();
      }
    }else{
      vector<pair<double,double> >::iterator vdpitr = m_pairs.begin();
      for (;vdpitr!=m_pairs.end();++vdpitr){
        std::cout << "Value = " << (*vdpitr).first << ", eff = " 
                  << (*vdpitr).second << std::endl;
      }
    }
  }
  
  void Interpolator::Setup(ifstream &ifs,int &level, int &jump,
                           vector<int> &npoints, int &nline){
    
    // This method can certainly be improved - possibly using XML
    // Setup parses the input file and decides to recursively create 
    // more Interpolators at each new dimension
    // At the lowest dimension, a vector of pairs of doubles is held
    // instead, each pair being a point-value pairing
    
    // Add a new entry in for each new level encountered
    if ( npoints.size() < static_cast<std::vector<int>::size_type>(level) ) npoints.push_back(0);
    
    while (!ifs.eof()){
      
      string read_line;
      double read_double;
      
      getline(ifs,read_line); nline++;
      istringstream iss(read_line);
      
      // count number of values in this line
      int numbercount = 0;
      while (iss >> read_double) ++numbercount;

      if (!numbercount){

        // No numbers on this line, two possibilities:
        // 1) It's a blank line, ignore it
        // 2) It's a list of level identifiers at the beginning 
        //    of the file, store these

        if (nline==1){ 

          // These are the level identifiers
          // No new Interpolators have been made yet so this is the main
          // Interpolator and can directly access m_levelnames

          string level_name;
          int nlevel = 0;
          iss.clear(); iss.str(read_line);
          while (iss >> level_name){
            nlevel++;
            m_levelnames[level_name] = nlevel; 
          }
        }
        
        continue;
      }
      
      // only one value - must be a point so need to make an Interpolator
      if (numbercount == 1){
        
        Interpolator w(read_double,ifs,level,jump,npoints,nline);
        m_interpolators.push_back(w);
        
        // Take care of moving down the appropriate number of levels
        if (jump){
          if (!npoints[level-1]) npoints[level-1] = m_interpolators.size();
          else if ( m_interpolators.size() != static_cast<std::vector<Interpolator>::size_type>(npoints[level-1]) ){
            string message("inconsistent number of points");
            exitParseWithMessage(message,nline,read_line);
          }
          --jump; --level;
          return;
        }
        
      }else{

        // More than one value in read line, must the lowest dimension, 
        // read in the points and values from both consecutive lines and
        // store as pairs of doubles
        
        iss.clear(); iss.str(read_line); // Go to beginning of string?
        string read_second_line;
        getline(ifs,read_second_line); nline++;
        istringstream iss2(read_second_line);
        double read_second_double;
        
        // Check to see if the number of points is consistent
        while (iss >> read_double){
          if (!(iss2 >> read_second_double)) 
            exitParseWithMessage("#points > #values, reorganise input file",
                            nline,read_second_line);
          pair<double,double> lowest_dim_value(read_double,read_second_double);
          m_pairs.push_back(lowest_dim_value);
        }
        if (iss2 >> read_second_double)
          exitParseWithMessage("#points < #values, reorganise input file",
                          --nline,read_line);
        if (!npoints[level-1]) npoints[level-1] = m_pairs.size();
        else if ( m_pairs.size() != static_cast<std::vector<pair<double,double> >::size_type>(npoints[level-1]) ){
          exitParseWithMessage("inconsistent number of points in lowest dimension",
                          --nline,read_line);
        }

        // Find number of dimensions to jump
        nDimensionsToJump(ifs,level,jump);

        return;
      }

    }// end of while (!ifs.eof())

   return;
  }
  
  double Interpolator::interpolate(deque<double> values, bool performedchecks){
    
    if ( !performedchecks && values.size() != m_limits.size()){ 
      std::cout << "Interpolator has " << m_limits.size()
                << " dimensions, supplied point exists in "
                << values.size() << " dimensions!" << std::endl;
      exit(1);
    }
    
    // Only check the limits the first time this function is called
    double def;
    if (!m_contBounds && !performedchecks && !inLimits(values,def))
      return def;

    performedchecks = true;
    
    
    double fraction,value1,value2;
    double value = values.front();
    
    if (m_interpolators.size()){
      
      // If m_interpolators is filled, this is not the lowest dimension and interpolate
      // must be called for Interpolators in the next dimension down. Use ClosestInterpolators
      // to find the appropriate two Interpolators
      
      values.pop_front();
      ClosestInterpolators ci(value);
      vector<Interpolator>::iterator viitr = std::adjacent_find(m_interpolators.begin(),
                                                                m_interpolators.end(),
                                                                ci);
      vector<Interpolator>::iterator viitr2;
      
      if (m_contBounds && viitr == m_interpolators.end()){
        if ( value < (*m_interpolators.begin()).getValue() ){
          viitr = m_interpolators.begin();
        }else{
          viitr = (m_interpolators.end()-1);
        }
        fraction = 0;
        value1 = (*viitr).interpolate(values,performedchecks);
        value2 = 0;
        
      }else{
        viitr2 = viitr + 1;
        // Fraction tells us where between the two points we are
        
        fraction = (value - (*viitr).getValue())/((*viitr2).getValue() - (*viitr).getValue());
        value1 = (*viitr).interpolate(values,performedchecks);
        value2 = (*viitr2).interpolate(values,performedchecks);
      }
      
    }else{
      
      // This is the lowest dimension. Use ClosestDoublePairs to find the appropriate two
      // pairs containing the values to interpolate between.
      
      ClosestDoublePairs cdp(value);
      vector<pair<double,double> >::iterator vdpitr = std::adjacent_find(m_pairs.begin(),
                                                                         m_pairs.end(),
                                                                         cdp);
      vector<pair<double,double> >::iterator vdpitr2;
      
      if (m_contBounds && vdpitr == m_pairs.end()){
        if ( value < (*m_pairs.begin()).first ){
          vdpitr = m_pairs.begin();
        }else{
          vdpitr = (m_pairs.end()-1);
        }

        fraction = 0;
        value1 = (*vdpitr).second;
        value2 = 0;
      }else{
        vdpitr2 = vdpitr + 1;
        // p1 and p2 are pointers to the adjacent elements found, replace them later
        fraction = (value - (*vdpitr).first)/((*vdpitr2).first - (*vdpitr).first);
        value1 = (*vdpitr).second;
        value2 = (*vdpitr2).second;
      }
      
    }
    
    // Do interpolation
    return value1 + fraction*(value2-value1);   
  }
  
  void Interpolator::findLimits(deque<pair<double,double> > &limits_deque){
    
    // Get limits in this quantity first
    // The points are assumed to be in ascending numerical order
    pair<double,double> limits;
    
     limits.first = (m_interpolators.size()) ? 
       m_interpolators.front().getValue() :
       m_pairs.front().first;      

     //std::cout << "Lower limit = " << limits.first << std::endl;

     limits.second = (m_interpolators.size()) ?
       m_interpolators.back().getValue():
       m_pairs.back().first;
    
     //std::cout << "Higher limit = " << limits.second << std::endl;

     limits_deque.push_back(limits);
    
     // Now call for next dimension
     if (m_interpolators.size())
       m_interpolators.front().findLimits(limits_deque);
     
     return;
     
   }
  
  void Interpolator::nDimensionsToJump(ifstream &ifs,int &level, int &jump){

    // How many dimensions to jump up?
    
    // Current position in file buffer, we're just looking ahead
    // so will reset to this position afterwards
    int position = ifs.tellg();
    
    int nsingleentries = 0;
    bool singleentry = true;
    
    // Count number of lines with just a single value
    while (singleentry && !ifs.eof()){
      
      string read_line;
      getline(ifs,read_line);
      istringstream iss(read_line);
      int numbercount = 0;
      double read_double;
      while (iss >> read_double) ++numbercount;
      if (!numbercount) continue;
      if (numbercount == 1) ++nsingleentries;
      else singleentry = false;
      
    }
    
    // If we've hit the end of the input file, jump up 
    // the total number of dimensions to finish off the
    // setup
    if (ifs.eof()) nsingleentries = level;
    
    --level;
    
    // Reset the buffer position
    ifs.seekg(position);

    // Record number of dimensions to jump up
    jump = nsingleentries - 1;

    return;
    
  }
  
  bool Interpolator::inLimits(deque<double> values, double& def){

    // Make local copy of held point limits so we can safely
    // delete entries
    deque<double> values_copy = values;
    deque<pair<double,double> > limits = m_limits;
    deque<pair<double,double> > defaults = m_defaults;
    
    // Check each value is within the limits of a given dimension
    while (values.size()){
      if ( values.front() < limits.front().first ){
        def = m_contBounds ? getNearestValue(values_copy) : defaults.front().first;
        return false;
      }
      if ( values.front() > limits.front().second ){
        def = m_contBounds ? getNearestValue(values_copy) : defaults.front().second;
        return false;
      }
      values.pop_front();
      limits.pop_front();
      defaults.pop_front();
    }
    return true;
  }

  void Interpolator::setDefault(double def){

    // Find number of levels to which to apply this default
    int nlevels = m_limits.size();

    // Get rid of previous defaults, if any
    m_defaults.clear();

    // For each level, add def as the default return value
    // for being outside limits in either direction
    while (nlevels--){
      pair<double,double> default_pair(def,def);
      m_defaults.push_back(default_pair);
    }
    
    return;
  }

  void Interpolator::setDefault(int level, double def){

    setDefault(level,def,def);
    return;
    
  }
  
  void Interpolator::setDefault(int level, double def_below, double def_above){
    
    if (level < 1 || static_cast<std::deque<pair<double,double> >::size_type>(level) > m_limits.size()){
      std::cout << "Cannot set default for level " << level << std::endl;
      exit(1);
    }

    pair<double,double> default_pair(def_below,def_above);
    m_defaults[level-1] = default_pair;
    return;
    
  }
  
  void Interpolator::setDefault(string levelname, double def){
    setDefault(levelname,def,def);
    return;
  }
  
  void Interpolator::setDefault(string levelname, double def_below, double def_above){

    if ( !m_levelnames.size()){
      std::cout << "No level names specified in input file, cannot set default with name" << std::endl;
      exit(1);
    }

    if (!m_levelnames[levelname]){
      std::cout << "Level name " << levelname << " not found (typo?)" << std::endl;
      exit(1);
    }

    setDefault(m_levelnames[levelname],def_below,def_above);

    return;
  }
  
  double Interpolator::getNearestValue(deque<double> values){
    
    double value = values.front();

    if (m_interpolators.size()){
      
      values.pop_front();
      vector<Interpolator>::iterator closestPoint;
      
      for (vector<Interpolator>::iterator intItr=m_interpolators.begin();
           intItr!=m_interpolators.end();++intItr){
        if (intItr==m_interpolators.begin() && 
            fabs(value-(*intItr).getValue()) < fabs(value-(*(intItr+1)).getValue()))
          closestPoint = intItr;
        else if (intItr==m_interpolators.end()-1 && 
                 fabs(value-(*intItr).getValue()) < fabs(value-(*(intItr-1)).getValue()))
          closestPoint = intItr;
        else if (fabs(value-(*intItr).getValue()) < fabs(value-(*(intItr-1)).getValue()) &&
                 fabs(value-(*intItr).getValue()) < fabs(value-(*(intItr+1)).getValue()))
          closestPoint = intItr;
      }
      return (*closestPoint).getNearestValue(values);      
      
    }else{
      
      vector<pair<double,double> >::iterator closestPoint;
      
      for (vector<pair<double,double> >::iterator vpdItr=m_pairs.begin();
           vpdItr!=m_pairs.end();++vpdItr){
        if (vpdItr==m_pairs.begin() && 
            fabs(value-(*vpdItr).first) < fabs(value-(*(vpdItr+1)).first))
          closestPoint = vpdItr;
        else if (vpdItr==m_pairs.end()-1 && 
                 fabs(value-(*vpdItr).first) < fabs(value-(*(vpdItr-1)).first))
          closestPoint = vpdItr;
        else if (fabs(value-(*vpdItr).first) < fabs(value-(*(vpdItr-1)).first) &&
                 fabs(value-(*vpdItr).first) < fabs(value-(*(vpdItr+1)).first))
          closestPoint = vpdItr;
      }
      return (*closestPoint).second;      
      
    }
  }
  
  void Interpolator::exitParseWithMessage(string message, int nline, string line){
    std::cout << "Exiting Interpolator, " << message << std::endl;
    std::cout << "Problem occurred at about line " << nline << " of input file:" << std::endl;
    std::cout << line << std::endl;
    exit(1);
  }

  //=======================================================
  // ClosestInterpolator methods
  //=======================================================

  ClosestInterpolators::ClosestInterpolators(double value): m_value(value){}
  
  bool ClosestInterpolators::operator()(Interpolator &i1, Interpolator &i2){
    return (i1.getValue() <= m_value && i2.getValue() >= m_value) ? true : false;
  }
  
  //=======================================================
  // ClosestDoublePairs methods
  //=======================================================
  
  ClosestDoublePairs::ClosestDoublePairs(double value): m_value(value){}
  
  bool ClosestDoublePairs::operator()(pair<double,double> &p1, 
                                      pair<double,double> &p2){
    return (p1.first <= m_value && p2.first >= m_value) ? true : false;
  }
  
}

---