Alexandria  2.25.0
SDC-CH common library for the Euclid project
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Cumulative.cpp
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 
27 #include "XYDataset/XYDataset.h"
28 #include <cstdlib> // for size_t
29 
30 namespace Euclid {
31 namespace MathUtils {
32 
34  : m_x_sampling(std::move(other.m_x_sampling)), m_y_sampling(std::move(other.m_y_sampling)) {}
35 
37  m_x_sampling = std::move(other.m_x_sampling);
38  m_y_sampling = std::move(other.m_y_sampling);
39  return *this;
40 }
41 
42 Cumulative::Cumulative(const Cumulative& other) : m_x_sampling(other.m_x_sampling), m_y_sampling(other.m_y_sampling) {}
43 
45  m_x_sampling = other.m_x_sampling;
46  m_y_sampling = other.m_y_sampling;
47  return *this;
48 }
49 
51  : m_x_sampling(x_sampling), m_y_sampling(y_sampling) {
52  if (x_sampling.size() != y_sampling.size()) {
53  throw Elements::Exception("Cumulative: X and Y sampling do not have the same length.");
54  }
55 }
56 
57 Cumulative::Cumulative(const XYDataset::XYDataset& sampling) : m_x_sampling{}, m_y_sampling{} {
58  auto iter = sampling.begin();
59  while (iter != sampling.end()) {
60  m_x_sampling.push_back((*iter).first);
61  m_y_sampling.push_back((*iter).second);
62  ++iter;
63  }
64 }
65 
69  auto iter = sampling.begin();
70  while (iter != sampling.end()) {
71  xs.push_back((*iter).first);
72  ys.push_back((*iter).second);
73  ++iter;
74  }
75  return Cumulative::fromPdf(xs, ys);
76 }
77 
79  double total = 0.;
80  std::vector<double> cumul{};
81  auto iter_pdf = pdf_sampling.cbegin();
82 
83  while (iter_pdf != pdf_sampling.cend()) {
84  total += *(iter_pdf);
85  cumul.push_back(total);
86  ++iter_pdf;
87  }
88 
89  return Cumulative(x_sampling, cumul);
90 }
91 
93  double total = m_y_sampling.back();
94  std::vector<double> cumul{};
95  auto iter = m_y_sampling.begin();
96  while (iter != m_y_sampling.end()) {
97  cumul.push_back(*iter / total);
98  ++iter;
99  }
100 
101  m_y_sampling = std::move(cumul);
102 }
103 
104 double Cumulative::findValue(double ratio, TrayPosition position) const {
105 
106  if (ratio > 1. || ratio < 0.) {
107  throw Elements::Exception("Cumulative::findValue : ratio parameter must be in range [0,1]");
108  }
109 
110  double value = m_y_sampling.back() * ratio;
111 
112  auto iter_x = m_x_sampling.cbegin();
113  auto iter_y = m_y_sampling.cbegin();
114  while (iter_y != m_y_sampling.cend() && (*iter_y) < value) {
115  ++iter_x;
116  ++iter_y;
117  }
118 
119  double begin_value = *iter_x;
120  double tray = *iter_y;
121  while (iter_y != m_y_sampling.cend() && (*iter_y) == tray) {
122  ++iter_x;
123  ++iter_y;
124  }
125 
126  double end_value = *(--iter_x);
127 
128  double result = 0;
129  switch (position) {
130  case TrayPosition::begin:
131  result = begin_value;
132  break;
133  case TrayPosition::middle:
134  result = 0.5 * (begin_value + end_value);
135  break;
136  case TrayPosition::end:
137  result = end_value;
138  break;
139  }
140 
141  return result;
142 }
143 
145 
146  if (rate > 1. || rate <= 0.) {
147  throw Elements::Exception("Cumulative::findMinInterval : rate parameter must be in range ]0,1]");
148  }
149 
150  rate *= m_y_sampling.back();
151 
152  double first_correction = m_y_sampling.front();
153 
154  auto iter_1_x = m_x_sampling.cbegin();
155  auto iter_2_x = ++(m_x_sampling.cbegin());
156  auto iter_1_y = m_y_sampling.cbegin();
157  auto iter_2_y = ++(m_y_sampling.cbegin());
158  auto min_x = m_x_sampling.cbegin();
159  auto max_x = --(m_x_sampling.cend());
160  while (iter_1_x != m_x_sampling.cend()) {
161  while (iter_2_x != m_x_sampling.cend() && (*iter_2_y - *iter_1_y + first_correction) < rate) {
162  ++iter_2_x;
163  ++iter_2_y;
164  }
165  if (iter_2_x == m_x_sampling.cend()) {
166  break;
167  }
168  if ((*iter_2_x - *iter_1_x) <= (*max_x - *min_x)) {
169  max_x = iter_2_x;
170  min_x = iter_1_x;
171  }
172  ++iter_1_x;
173  ++iter_1_y;
174  first_correction = 0.;
175  }
176 
177  return std::make_pair(*min_x, *max_x);
178 }
179 
181 
182  if (rate > 1. || rate <= 0.) {
183  throw Elements::Exception("Cumulative::findCenteredInterval : rate parameter must be in range ]0,1]");
184  }
185 
186  double min_value = m_y_sampling.back() * (0.5 - rate / 2.0);
187  double max_value = m_y_sampling.back() * (0.5 + rate / 2.0);
188 
189  auto iter_x = m_x_sampling.cbegin();
190  auto iter_y = m_y_sampling.cbegin();
191  while (iter_y != m_y_sampling.cend() && (*iter_y) < min_value) {
192  ++iter_x;
193  ++iter_y;
194  }
195 
196  if ((*iter_y) < max_value) {
197  double tray = *iter_y;
198  while (iter_y != m_y_sampling.cend() && (*iter_y) == tray) {
199  ++iter_x;
200  ++iter_y;
201  }
202  double min_x = (iter_x != m_x_sampling.begin()) ? *(iter_x - 1) : *iter_x;
203 
204  while (iter_y != m_y_sampling.cend() && (*iter_y) < max_value) {
205  ++iter_x;
206  ++iter_y;
207  }
208  double max_x = *iter_x;
209 
210  return std::make_pair(min_x, max_x);
211  } else {
212  double min_x = (iter_x != m_x_sampling.begin()) ? *(iter_x - 1) : *iter_x;
213  double max_x = *iter_x;
214 
215  return std::make_pair(min_x, max_x);
216  }
217 }
218 
219 } // namespace MathUtils
220 } // namespace Euclid
TrayPosition
when looking for the position having a given value, one may encounter tray where the value is constan...
Definition: Cumulative.h:51
std::vector< double > m_x_sampling
Definition: Cumulative.h:154
Class for build cumulative from PDF and extract feature out of it.
Definition: Cumulative.h:41
double findValue(double ratio, TrayPosition position=TrayPosition::middle) const
Find the first horizontal sample which vertical value is bigger or equal to the ratio value...
Definition: Cumulative.cpp:104
T front(T...args)
static Cumulative fromPdf(std::vector< double > &x_sampling, std::vector< double > &pdf_sampling)
Factory from the sampling of a PDF. The Cumulative vertical samples are build as the sum of the the p...
Definition: Cumulative.cpp:78
std::vector< double > m_y_sampling
Definition: Cumulative.h:155
T cend(T...args)
Cumulative & operator=(Cumulative &&other)
move assignation operator
Definition: Cumulative.cpp:36
std::pair< double, double > findCenteredInterval(double rate) const
return the horizontal interval starting where the Cumulative has value (1-ratio)/2 and ending where t...
Definition: Cumulative.cpp:180
void normalize()
Normalize the Cumulative. After calling this function the last vertical value is 1.0.
Definition: Cumulative.cpp:92
const_iterator begin() const
Returns a const iterator to the first pair of the dataset.
Definition: XYDataset.cpp:36
T make_pair(T...args)
T move(T...args)
T size(T...args)
T cbegin(T...args)
std::pair< double, double > findMinInterval(double rate) const
Scan the horizontal axis looking for the smallest x-interval for which the vertical interval is at le...
Definition: Cumulative.cpp:144
Cumulative(Cumulative &&other)
move constructor
Definition: Cumulative.cpp:33
This module provides an interface for accessing two dimensional datasets (pairs of (X...
Definition: XYDataset.h:59
T back(T...args)
const_iterator end() const
Returns a const iterator to the one after last pair dataset.
Definition: XYDataset.cpp:40