TabulatedFunctionComplexLgAmpPhase.cc

Go to the documentation of this file.

#include <iomanip>
#include <sstream>
#include <cmath>
#if !defined(__clang_major__) || __clang_major__ < 5
#  include <ext/algorithm>
#endif
#include <utl/TabulatedFunctionComplexLgAmpPhase.h>
#include <utl/Test.h>
#include <utl/AugerException.h>
#include <utl/ErrorLogger.h>
#include <limits>

using namespace std;
using namespace utl;

typedef TabulatedFunctionComplexLgAmpPhase::Array Array;
typedef TabulatedFunctionComplexLgAmpPhase::ArrayComplex ArrayComplex;

void
TabulatedFunctionComplexLgAmpPhase::FillTable(const Array& xValues,
                                    const ArrayComplex& yValues)
{
  const unsigned int xSize = xValues.size();

  if (xSize != yValues.size()) {
    ostringstream err;
    err << "Size of x and y arrays are different.";
    ERROR(err);
    throw OutOfBoundException(err.str());
  }

  fX = xValues;
  fY = yValues;

  if (is_sorted(fX.begin(), fX.end()))
    return;

  // insertion sort
  for (unsigned int i = 0; i < xSize; ++i) {
    unsigned int mini = i;
    double minx = fX[i];
    for (unsigned int j = i+1; j < xSize; ++j) {
      const double x = fX[j];
      if (x < minx) {
        minx = x;
        mini = j;
      }
    }
    if (mini != i) {
      swap(fX[i], fX[mini]);
      swap(fY[i], fY[mini]);
    }
  }
}


void
TabulatedFunctionComplexLgAmpPhase::PushBack(const double x, const ComplexLgAmpPhase& y)
{
  // Check that same X value is not accidentally pushed twice in a row
  // (a common error)
  if (fX.size() && abs(fX.back() - x) < 1./numeric_limits<float>::max()) {
    WARNING("Refusing to push the same x value twice.");
    return;
  }

  fX.push_back(x);
  fY.push_back(y);
}


ComplexLgAmpPhase
TabulatedFunctionComplexLgAmpPhase::Y(double x)
  const
{
  const unsigned int size = fX.size();

  if (!size) {
    ostringstream msg;
    msg << "No points in table.";
    ERROR(msg);
    throw OutOfBoundException(msg.str());
  }

  if (size == 1)
    return fY[0];

  if (x < fX[0] || x > fX[size-1]) {

    if (x < fX[0] && (fX[0] - x) / (fX[1] - fX[2]) < fBoundaryTolerance)
      x = fX[0];
    else if (x > fX[size-1] &&
            (x - fX[size-1]) / (fX[size-1] - fX[size-2]) < fBoundaryTolerance)
      x = fX[size-1];
    else {
      boost::format err("Value x=%.10f is out of table bounds [%.10f, %.12f].");
      err % x % fX[0] % fX[size-1];
      ERROR(err);
      throw OutOfBoundException(err.str());
    }

  }

  const Array::const_iterator xBegin = fX.begin();
  const Array::const_iterator xIt = upper_bound(xBegin, fX.end(), x);
  int index = xIt - xBegin - 1;
  const double x1 = fX[index];
  const ComplexLgAmpPhase y1 = fY[index];
  if (x1 == x)
    return y1;
  ++index;
  const double x2 = fX[index];
  const ComplexLgAmpPhase y2 = fY[index];

  // linear interpolation done separately in log10(amplitude) and phase
  return Interpolate(x1,y1,x2,y2,x);
}

ConstTabulatedFunctionComplexLgAmpPhaseIterator
TabulatedFunctionComplexLgAmpPhase::FindX(const double x)
  const
{
  const Array::const_iterator xBegin = fX.begin();
  const Array::const_iterator xEnd = fX.end();
  const Array::const_iterator xIt = lower_bound(xBegin, xEnd, x);
  if (xIt != xEnd && x == *xIt) {
    const ArrayComplex::const_iterator yIt = fY.begin() + (xIt - xBegin);
    return ConstTabulatedFunctionComplexLgAmpPhaseIterator(xIt, yIt);
  } else
    return ConstTabulatedFunctionComplexLgAmpPhaseIterator(xEnd, fY.end());
}


ConstTabulatedFunctionComplexLgAmpPhaseIterator
TabulatedFunctionComplexLgAmpPhase::FindY(const ComplexLgAmpPhase& y)
  const
{
  const ArrayComplex::const_iterator yIt = find(fY.begin(), fY.end(), y);
  const Array::const_iterator xIt = (yIt != fY.end()) ?
    fX.begin() + (yIt - fY.begin()) : fX.end();
  return ConstTabulatedFunctionComplexLgAmpPhaseIterator(xIt, yIt);
}


TabulatedFunctionComplexLgAmpPhaseIterator
TabulatedFunctionComplexLgAmpPhase::Insert(const double x, const ComplexLgAmpPhase& y)
{
  const Array::iterator xIt = lower_bound(fX.begin(), fX.end(), x);
  const ArrayComplex::iterator yIt = fY.begin() + (xIt - fX.begin());
  const Array::iterator ixIt = fX.insert(xIt, x);
  const ArrayComplex::iterator iyIt = fY.insert(yIt, y);
  return TabulatedFunctionComplexLgAmpPhaseIterator(ixIt, iyIt);
}


bool
TabulatedFunctionComplexLgAmpPhase::IsInValidRange(const double x)
  const
{
  return fX.front() <= x && x <= fX.back();
}


/*
TabulatedFunctionComplexLgAmpPhase
TabulatedFunctionComplexLgAmpPhase::operator*(const TabulatedFunctionComplexLgAmpPhase& rhs)
  const
{
  // to be written
}
*/

void
TabulatedFunctionComplexLgAmpPhase::Pow(const double power)
{
  for (ArrayComplex::iterator it = YBegin(); it != YEnd(); ++it)
    it->Pow(power);
}


TabulatedFunctionComplexLgAmpPhase&
TabulatedFunctionComplexLgAmpPhase::operator*=(const TabulatedFunctionComplexLgAmpPhase& rhs)
{

  // check if both tabulated functions have at least two elements
  if (this->GetNPoints()<2) {
    ostringstream err;
    err << "LHS TabulatedFunctionComplexLgAmpPhase has less than 2 elements, cannot be multiplied.";
    ERROR(err);
    throw OutOfBoundException(err.str());
  }

  if (rhs.GetNPoints()<2) {
    ostringstream err;
    err << "RHS TabulatedFunctionComplexLgAmpPhase has less than 2 elements, cannot be multiplied.";
    ERROR(err);
    throw OutOfBoundException(err.str());
  }

  // determine the common range of the two TabulatedFunctionComplexLgAmpPhase
  const double commonstart = max(*XBegin(),*rhs.XBegin());
  const double commonend = min(this->GetX(this->GetNPoints()-1),rhs.GetX(rhs.GetNPoints()-1));

  /* comment out warning, it is not really useful and produces a lot of verbosity
  // warn if the valid range of the combined TabulatedFunctionComplexLgAmpPhase is smaller than those of the two individual TabluatedFunctionComplex
  if (  (Test<Not<CloseTo> >(*XBegin(),*rhs.XBegin()) )
     || (Test<Not<CloseTo> >(this->GetX(this->GetNPoints()-1),rhs.GetX(rhs.GetNPoints()-1))) ) {
    ostringstream warn;
    warn << "Range of combined TabulatedFunctionComplexLgAmpPhase has been decreased to ["
         << commonstart
         << ":"
         << commonend
         << "]";
    WARNING(warn);
  }
  */

  // initialize iterator for lhs TabulatedFunctionComplexLgAmpPhase and propagate it to the first element in the common range
  Array::const_iterator lxIt = XBegin();
  while (*lxIt < commonstart)
    ++lxIt;

  // initialize iterator for lhs TabulatedFunctionComplexLgAmpPhase and propagate it to the first element in the common range
  Array::const_iterator rxIt = rhs.XBegin();
  while (*rxIt < commonstart)
    ++rxIt;

  Array tempX;
  ArrayComplex tempY;

  // create a union of the two x-arrays
  do
  {
    // test if point is present in both tabulated functions
    if (Test<CloseTo>(*lxIt,*rxIt)) {
      tempX.push_back(*lxIt);
      ++lxIt; // point has been used, go to next one
      ++rxIt; // point has been used, go to next one
      continue;
    }

    // points were not the same, so test if lhs point is the next to be used
    if (*lxIt < *rxIt) {
      tempX.push_back(*lxIt);
      ++lxIt; // point has been used, go to next one
      continue;
    }

    // if we arrive here, then the rhs point is the next to be used
    tempX.push_back(*rxIt);
    ++rxIt; // point has been used, go to next one

  } while ((*lxIt < commonend) || (*rxIt < commonend));
  // can this comparison fail because of rounding errors? then we could read past the end of an xArray!

  // check if last point has to be added manually
  if (Test<Not<CloseTo> >(tempX.back(),commonend))
    tempX.push_back(commonend);

  // now create the corresponding array of Y values
  tempY.reserve(tempX.size());
  for (Array::const_iterator xIt = tempX.begin(); xIt != tempX.end(); ++xIt)
    tempY.push_back(this->Y(*xIt)*rhs.Y(*xIt));

  // set the internal arrays to the temporary ones (swap is most efficient)
  fX.swap(tempX);
  fY.swap(tempY);

  return *this;
}


// Configure (x)emacs for this file ...
// Local Variables:
// mode: c++
// compile-command: "make -C .. -k"
// End: