19 #ifdef INVERSE_CUMULATIVE_IMPL
31 template <
typename TKnot>
32 class InverseCumulative<TKnot, typename std::enable_if<std::is_floating_point<TKnot>::value>::type> {
35 : m_knots(std::
move(knots)), m_pdf(std::
move(pdf)), m_cdf(m_pdf.size()) {
36 if (m_knots.size() != m_knots.size()) {
37 throw Elements::Exception() <<
"PDF and knots dimensionality do not match: " << m_knots.size() <<
" != " << knots.
size();
45 m_cdf[i] = (m_knots[i] - m_knots[i - 1]) * (m_pdf[i] + m_pdf[i - 1]) / 2.;
46 m_cdf[i] += m_cdf[i - 1];
50 while (m_cdf.size() > 1 && m_cdf.back() == *(m_cdf.end() - 2)) {
55 m_min = m_cdf.front();
56 m_range = m_cdf.back() - m_min;
60 if (p < 0. || p > 1.) {
61 throw Elements::Exception() <<
"Cumulative::findInterpolatedValue : p parameter must be in the range [0,1]";
64 const double unnormed_p = p * m_range + m_min;
67 if (unnormed_p <= m_cdf.front())
68 return m_knots.front();
69 if (unnormed_p >= m_cdf.back())
70 return m_knots.back();
74 const double x0 = m_knots[i], x1 = m_knots[i + 1];
75 const double cdf0 = m_cdf[i], cdf1 = m_cdf[i + 1];
76 const double p0 = m_pdf[i], p1 = m_pdf[i + 1];
80 double interval_p = (unnormed_p - cdf0) / (cdf1 - cdf0);
81 return (x1 - x0) * interval_p + x0;
87 if (x1 == x0 || cdf0 == cdf1 || cdf0 >= unnormed_p) {
95 const double a = (p1 - p0) / (x1 - x0);
96 const double b = p0 - a * x0;
98 assert(std::abs(a * x0 + b - p0) < 1e-8);
99 assert(std::abs(a * x1 + b - p1) < 1e-8);
105 const double c = cdf0 - (a / 2.) * (x0 * x0) - b * x0;
116 if (s0 >= x0 - 1e-8 && s0 <= x1 + 1e-8) {
119 assert(s1 >= x0 - 1e-8 && s1 <= x1 + 1e-8);
126 double m_min, m_range;
132 template <
typename TKnot>
133 class InverseCumulative<TKnot, typename std::enable_if<!std::is_floating_point<TKnot>::value>::type> {
138 m_cdf[i] += m_cdf[i - 1];
140 for (
auto& v : m_cdf) {
146 if (p < 0. || p > 1.) {
147 throw Elements::Exception() <<
"Cumulative::findInterpolatedValue : p parameter must be in the range [0,1]";
std::pair< double, double > solveSquare(double a, double b, double c, double y)
TKnot operator()(double p) const
InverseCumulative(std::vector< TKnot > knots, std::vector< double > pdf)