Alexandria  2.27.0
SDC-CH common library for the Euclid project
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
InverseCumulative.icpp
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2012-2021 Euclid Science Ground Segment
3  *
4  * This library is free software; you can redistribute it and/or modify it under
5  * the terms of the GNU Lesser General Public License as published by the Free
6  * Software Foundation; either version 3.0 of the License, or (at your option)
7  * any later version.
8  *
9  * This library is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11  * FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
12  * details.
13  *
14  * You should have received a copy of the GNU Lesser General Public License
15  * along with this library; if not, write to the Free Software Foundation, Inc.,
16  * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
17  */
18 
19 #ifdef INVERSE_CUMULATIVE_IMPL
20 
22 #include <algorithm>
23 #include <cassert>
24 
25 namespace Euclid {
26 namespace MathUtils {
27 
31 template <typename TKnot>
32 class InverseCumulative<TKnot, typename std::enable_if<std::is_floating_point<TKnot>::value>::type> {
33 public:
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();
38  }
39  if (!std::is_sorted(m_knots.begin(), m_knots.end())) {
40  throw Elements::Exception() << "Knots must be sorted";
41  }
42 
43  m_cdf[0] = m_pdf[0];
44  for (std::size_t i = 1; i < m_cdf.size(); ++i) {
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];
47  }
48 
49  // Remove trailing knots with no probability
50  while (m_cdf.size() > 1 && m_cdf.back() == *(m_cdf.end() - 2)) {
51  m_cdf.pop_back();
52  m_knots.pop_back();
53  }
54 
55  m_min = m_cdf.front();
56  m_range = m_cdf.back() - m_min;
57  }
58 
59  double operator()(double p) const {
60  if (p < 0. || p > 1.) {
61  throw Elements::Exception() << "Cumulative::findInterpolatedValue : p parameter must be in the range [0,1]";
62  }
63 
64  const double unnormed_p = p * m_range + m_min;
65 
66  // Find segment
67  if (unnormed_p <= m_cdf.front())
68  return m_knots.front();
69  if (unnormed_p >= m_cdf.back())
70  return m_knots.back();
71 
72  std::size_t i = std::upper_bound(m_cdf.begin(), m_cdf.end(), unnormed_p) - m_cdf.begin() - 1;
73 
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];
77 
78  // If p0 == p1 we are on a uniform area
79  if (p0 == p1) {
80  double interval_p = (unnormed_p - cdf0) / (cdf1 - cdf0);
81  return (x1 - x0) * interval_p + x0;
82  }
83 
84  // If both cdf are the same, then the interval has 0 probability.
85  // If both x are the same, this is probably a discontinuity.
86  // In those cases, return the lower bound, which might be at least defined (or has a probability of exactly p).
87  if (x1 == x0 || cdf0 == cdf1 || cdf0 >= unnormed_p) {
88  return x0;
89  }
90 
91  // Since we assume a linear interpolation, we know that
92  // p0 = a * x0 + b
93  // p1 = a * x1 + b
94  // So we can solve for a and b
95  const double a = (p1 - p0) / (x1 - x0);
96  const double b = p0 - a * x0;
97 
98  assert(std::abs(a * x0 + b - p0) < 1e-8);
99  assert(std::abs(a * x1 + b - p1) < 1e-8);
100 
101  // The CDF is the integral, so we also know that
102  // cdf0 = a/2 * x0^2 + b * x0 + c
103  // cdf1 = a/2 * x1^2 + b * x1 + c
104  // We already know a and b, so it is easy to solve for c
105  const double c = cdf0 - (a / 2.) * (x0 * x0) - b * x0;
106 
107  // Double check that the equation passes through the CDF at the limits
108  // assert(std::abs((a / 2.) * (x0 * x0) + b * x0 + c - cdf0) < 1e-4);
109  // assert(std::abs((a / 2.) * (x1 * x1) + b * x1 + c - cdf1) < 1e-4);
110 
111  // We have the equation, so we now need to solve x for p
112  double s0, s1;
113  std::tie(s0, s1) = solveSquare(a / 2, b, c, unnormed_p);
114 
115  // Pick the possible result that lies within [x0, x1]
116  if (s0 >= x0 - 1e-8 && s0 <= x1 + 1e-8) {
117  return s0;
118  }
119  assert(s1 >= x0 - 1e-8 && s1 <= x1 + 1e-8);
120  return s1;
121  }
122 
123 private:
124  std::vector<TKnot> m_knots;
125  std::vector<double> m_pdf, m_cdf;
126  double m_min, m_range;
127 };
128 
132 template <typename TKnot>
133 class InverseCumulative<TKnot, typename std::enable_if<!std::is_floating_point<TKnot>::value>::type> {
134 public:
135  InverseCumulative(std::vector<TKnot> knots, std::vector<double> pdf) : m_knots(std::move(knots)), m_cdf(std::move(pdf)) {
136  // Compute cumulative and normalize
137  for (std::size_t i = 1; i < m_cdf.size(); ++i) {
138  m_cdf[i] += m_cdf[i - 1];
139  }
140  for (auto& v : m_cdf) {
141  v /= m_cdf.back();
142  }
143  }
144 
145  TKnot operator()(double p) const {
146  if (p < 0. || p > 1.) {
147  throw Elements::Exception() << "Cumulative::findInterpolatedValue : p parameter must be in the range [0,1]";
148  }
149  std::size_t i = std::lower_bound(m_cdf.begin(), m_cdf.end(), p) - m_cdf.begin();
150  return m_knots[i];
151  }
152 
153 private:
154  std::vector<TKnot> m_knots;
155  std::vector<double> m_cdf;
156 };
157 
158 } // namespace MathUtils
159 } // namespace Euclid
160 
161 #endif
T tie(T...args)
T upper_bound(T...args)
T lower_bound(T...args)
std::pair< double, double > solveSquare(double a, double b, double c, double y)
Definition: Solvers.cpp:25
TKnot operator()(double p) const
T move(T...args)
T size(T...args)
STL class.
T is_sorted(T...args)
InverseCumulative(std::vector< TKnot > knots, std::vector< double > pdf)