HcalCubicInterpolator.cc

Go to the documentation of this file.

#include <algorithm>
#include "FWCore/Utilities/interface/Exception.h"

#include "CondFormats/HcalObjects/interface/HcalCubicInterpolator.h"

HcalCubicInterpolator::HcalCubicInterpolator()
{
}

HcalCubicInterpolator::HcalCubicInterpolator(
    const std::vector<Triple>& points)
{
    const std::size_t sz = points.size();
    if (sz)
    {
        std::vector<Triple> tmp(points);
        std::sort(tmp.begin(), tmp.end());
        abscissae_.reserve(sz);
        values_.reserve(sz);
        derivatives_.reserve(sz);
        for (std::size_t i=0; i<sz; ++i)
        {
            const Triple& t(tmp[i]);
            abscissae_.push_back(std::get<0>(t));
            values_.push_back(std::get<1>(t));
            derivatives_.push_back(std::get<2>(t));
        }
        const std::size_t szm1 = sz - 1;
        for (std::size_t i=0; i<szm1; ++i)
            if (abscissae_[i] == abscissae_[i+1])
                throw cms::Exception("In HcalCubicInterpolator constructor:"
                                            " abscissae must not coincide");
    }
}

double HcalCubicInterpolator::operator()(const double x) const
{
    double result = 0.0;
    const std::size_t sz = abscissae_.size();
    if (sz)
    {
        if (sz > 1)
        {
            const std::size_t szm1 = sz - 1;
            if (x >= abscissae_[szm1])
                result = values_[szm1] + derivatives_[szm1]*(x - abscissae_[szm1]);
            else if (x <= abscissae_[0])
                result = values_[0] + derivatives_[0]*(x - abscissae_[0]);
            else
            {
                const std::size_t cell = std::upper_bound(abscissae_.begin(), abscissae_.end(), x) -
                                         abscissae_.begin() - 1;
                const std::size_t cellp1 = cell + 1;
                const double dx = abscissae_[cellp1] - abscissae_[cell];
                const double t = (x - abscissae_[cell])/dx;
                const double onemt = 1.0 - t;
                const double h00 = onemt*onemt*(1.0 + 2.0*t);
                const double h10 = onemt*onemt*t;
                const double h01 = t*t*(3.0 - 2.0*t);
                const double h11 = t*t*onemt;
                result = h00*values_[cell] + h10*dx*derivatives_[cell] +
                         h01*values_[cellp1] - h11*dx*derivatives_[cellp1];
            }
        }
        else
            result = values_[0] + derivatives_[0]*(x - abscissae_[0]);
    }
    return result;
}

double HcalCubicInterpolator::xmin() const
{
    double result = 0.0;
    if (!abscissae_.empty())
        result = abscissae_[0];
    return result;
}

double HcalCubicInterpolator::xmax() const
{
    double result = 0.0;
    if (!abscissae_.empty())
        result = abscissae_.back();
    return result;
}

HcalCubicInterpolator HcalCubicInterpolator::approximateInverse() const
{
    const bool monotonous = isStrictlyIncreasing(values_.begin(), values_.end()) ||
                            isStrictlyDecreasing(values_.begin(), values_.end());
    if (!monotonous)
        throw cms::Exception("In HcalCubicInterpolator::inverse:"
                                    " can't invert non-monotonous functor");
    const std::size_t sz = abscissae_.size();
    std::vector<Triple> points;
    points.reserve(sz);
    for (std::size_t i=0; i<sz; ++i)
    {
        const double dydx = derivatives_[i];
        if (dydx == 0.0)
            throw cms::Exception("In HcalCubicInterpolator::inverse:"
                                        " can't invert functor with derivatives of 0");
        points.push_back(Triple(values_[i], abscissae_[i], 1.0/dydx));
    }
    return HcalCubicInterpolator(points);
}

BOOST_CLASS_EXPORT_IMPLEMENT(HcalCubicInterpolator)