Alexandria 2.25.0
SDC-CH common library for the Euclid project
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
30namespace Euclid {
31namespace 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
42Cumulative::Cumulative(const Cumulative& other) : m_x_sampling(other.m_x_sampling), m_y_sampling(other.m_y_sampling) {}
43
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
57Cumulative::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
104double 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
T back(T... args)
T cbegin(T... args)
Class for build cumulative from PDF and extract feature out of it.
Definition: Cumulative.h:41
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
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
std::vector< double > m_y_sampling
Definition: Cumulative.h:155
Cumulative(Cumulative &&other)
move constructor
Definition: Cumulative.cpp:33
Cumulative & operator=(Cumulative &&other)
move assignation operator
Definition: Cumulative.cpp:36
void normalize()
Normalize the Cumulative. After calling this function the last vertical value is 1....
Definition: Cumulative.cpp:92
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
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
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
This module provides an interface for accessing two dimensional datasets (pairs of (X,...
Definition: XYDataset.h:59
const_iterator begin() const
Returns a const iterator to the first pair of the dataset.
Definition: XYDataset.cpp:36
const_iterator end() const
Returns a const iterator to the one after last pair dataset.
Definition: XYDataset.cpp:40
T cend(T... args)
T front(T... args)
T make_pair(T... args)
T move(T... args)
STL namespace.
T push_back(T... args)
T size(T... args)