19#ifdef INVERSE_CUMULATIVE_IMPL
22#include "ElementsKernel/Logging.h"
34template <
typename TKnot>
35class InverseCumulative<TKnot, typename
std::enable_if<std::is_floating_point<TKnot>::value>::type> {
38 : m_knots(
std::
move(knots)), m_pdf(
std::
move(pdf)), m_cdf(m_pdf.size()) {
39 if (m_knots.size() != m_knots.size()) {
40 throw Elements::Exception() <<
"PDF and knots dimensionality do not match: " << m_knots.size() <<
" != " << knots.size();
43 throw Elements::Exception() <<
"Knots must be sorted";
48 m_cdf[i] = (m_knots[i] - m_knots[i - 1]) * (m_pdf[i] + m_pdf[i - 1]) / 2.;
49 m_cdf[i] += m_cdf[i - 1];
53 while (m_cdf.size() > 1 && m_cdf.back() == *(m_cdf.end() - 2)) {
58 m_min = m_cdf.front();
59 m_range = m_cdf.back() - m_min;
63 if (p < 0. || p > 1.) {
64 throw Elements::Exception() <<
"Cumulative::findInterpolatedValue : p parameter must be in the range [0,1]";
67 const double unnormed_p = p * m_range + m_min;
70 if (unnormed_p <= m_cdf.front()) {
71 return m_knots.front();
73 if (unnormed_p >= m_cdf.back()) {
74 return m_knots.back();
79 const double x0 = m_knots[i], x1 = m_knots[i + 1];
80 const double cdf0 = m_cdf[i], cdf1 = m_cdf[i + 1];
81 const double p0 = m_pdf[i], p1 = m_pdf[i + 1];
85 double interval_p = (unnormed_p - cdf0) / (cdf1 - cdf0);
86 return (x1 - x0) * interval_p + x0;
92 if (x1 == x0 || cdf0 == cdf1 || cdf0 >= unnormed_p) {
100 const double a = (p1 - p0) / (x1 - x0);
101 const double b = p0 - a * x0;
103 assert(std::abs(a * x0 + b - p0) < 1e-8);
104 assert(std::abs(a * x1 + b - p1) < 1e-8);
110 const double c = cdf0 - (a / 2.) * (x0 * x0) - b * x0;
121 if (s0 >= x0 - 1e-8 && s0 <= x1 + 1e-8) {
124 assert(s1 >= x0 - 1e-8 && s1 <= x1 + 1e-8);
128 loggerinverseCumul.
warn()<<
"Computation of the inverse Cumulative is not finite: use the median value.";
139 double m_min, m_range;
145template <
typename TKnot>
146class InverseCumulative<TKnot, typename
std::enable_if<!std::is_floating_point<TKnot>::value>::type> {
151 m_cdf[i] += m_cdf[i - 1];
153 for (
auto& v : m_cdf) {
159 if (p < 0. || p > 1.) {
160 throw Elements::Exception() <<
"Cumulative::findInterpolatedValue : p parameter must be in the range [0,1]";
static Logging getLogger(const std::string &name="")
void warn(const std::string &logMessage)
InverseCumulative(std::vector< TKnot > knots, std::vector< double > pdf)
TKnot operator()(double p) const
std::pair< double, double > solveSquare(double a, double b, double c, double y)