Alexandria 2.27.0
SDC-CH common library for the Euclid project
Loading...
Searching...
No Matches
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(x_sampling), m_y_sampling(y_sampling) {
35 if (x_sampling.size() != y_sampling.size()) {
36 throw Elements::Exception("Cumulative: X and Y sampling do not have the same length.");
37 }
38}
39
40Cumulative::Cumulative(const XYDataset::XYDataset& sampling) : m_x_sampling{}, m_y_sampling{} {
41 auto iter = sampling.begin();
42 while (iter != sampling.end()) {
43 m_x_sampling.push_back((*iter).first);
44 m_y_sampling.push_back((*iter).second);
45 ++iter;
46 }
47}
48
52 auto iter = sampling.begin();
53 while (iter != sampling.end()) {
54 xs.push_back((*iter).first);
55 ys.push_back((*iter).second);
56 ++iter;
57 }
58 return Cumulative::fromPdf(xs, ys);
59}
60
62 double total = 0.;
63 std::vector<double> cumul{};
64 auto iter_pdf = pdf_sampling.cbegin();
65
66 while (iter_pdf != pdf_sampling.cend()) {
67 total += *(iter_pdf);
68 cumul.push_back(total);
69 ++iter_pdf;
70 }
71
72 return Cumulative(x_sampling, cumul);
73}
74
76 double total = m_y_sampling.back();
77 std::vector<double> cumul{};
78 auto iter = m_y_sampling.begin();
79 while (iter != m_y_sampling.end()) {
80 cumul.push_back(*iter / total);
81 ++iter;
82 }
83
84 m_y_sampling = std::move(cumul);
85}
86
87double Cumulative::findValue(double ratio, TrayPosition position) const {
88
89 if (ratio > 1. || ratio < 0.) {
90 throw Elements::Exception("Cumulative::findValue : ratio parameter must be in range [0,1]");
91 }
92
93 double value = m_y_sampling.back() * ratio;
94
95 auto iter_x = m_x_sampling.cbegin();
96 auto iter_y = m_y_sampling.cbegin();
97 while (iter_y != m_y_sampling.cend() && (*iter_y) < value) {
98 ++iter_x;
99 ++iter_y;
100 }
101
102 double begin_value = *iter_x;
103 double tray = *iter_y;
104 while (iter_y != m_y_sampling.cend() && (*iter_y) == tray) {
105 ++iter_x;
106 ++iter_y;
107 }
108
109 double end_value = *(--iter_x);
110
111 double result = 0;
112 switch (position) {
114 result = begin_value;
115 break;
117 result = 0.5 * (begin_value + end_value);
118 break;
120 result = end_value;
121 break;
122 }
123
124 return result;
125}
126
128
129 if (rate > 1. || rate <= 0.) {
130 throw Elements::Exception("Cumulative::findMinInterval : rate parameter must be in range ]0,1]");
131 }
132
133 rate *= m_y_sampling.back();
134
135 double first_correction = m_y_sampling.front();
136
137 auto iter_1_x = m_x_sampling.cbegin();
138 auto iter_2_x = ++(m_x_sampling.cbegin());
139 auto iter_1_y = m_y_sampling.cbegin();
140 auto iter_2_y = ++(m_y_sampling.cbegin());
141 auto min_x = m_x_sampling.cbegin();
142 auto max_x = --(m_x_sampling.cend());
143 while (iter_1_x != m_x_sampling.cend()) {
144 while (iter_2_x != m_x_sampling.cend() && (*iter_2_y - *iter_1_y + first_correction) < rate) {
145 ++iter_2_x;
146 ++iter_2_y;
147 }
148 if (iter_2_x == m_x_sampling.cend()) {
149 break;
150 }
151 if ((*iter_2_x - *iter_1_x) <= (*max_x - *min_x)) {
152 max_x = iter_2_x;
153 min_x = iter_1_x;
154 }
155 ++iter_1_x;
156 ++iter_1_y;
157 first_correction = 0.;
158 }
159
160 return std::make_pair(*min_x, *max_x);
161}
162
164
165 if (rate > 1. || rate <= 0.) {
166 throw Elements::Exception("Cumulative::findCenteredInterval : rate parameter must be in range ]0,1]");
167 }
168
169 double min_value = m_y_sampling.back() * (0.5 - rate / 2.0);
170 double max_value = m_y_sampling.back() * (0.5 + rate / 2.0);
171
172 auto iter_x = m_x_sampling.cbegin();
173 auto iter_y = m_y_sampling.cbegin();
174 while (iter_y != m_y_sampling.cend() && (*iter_y) < min_value) {
175 ++iter_x;
176 ++iter_y;
177 }
178
179 if ((*iter_y) < max_value) {
180 double tray = *iter_y;
181 while (iter_y != m_y_sampling.cend() && (*iter_y) == tray) {
182 ++iter_x;
183 ++iter_y;
184 }
185 double min_x = (iter_x != m_x_sampling.begin()) ? *(iter_x - 1) : *iter_x;
186
187 while (iter_y != m_y_sampling.cend() && (*iter_y) < max_value) {
188 ++iter_x;
189 ++iter_y;
190 }
191 double max_x = *iter_x;
192
193 return std::make_pair(min_x, max_x);
194 } else {
195 double min_x = (iter_x != m_x_sampling.begin()) ? *(iter_x - 1) : *iter_x;
196 double max_x = *iter_x;
197
198 return std::make_pair(min_x, max_x);
199 }
200}
201
202} // namespace MathUtils
203} // 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
Cumulative(Cumulative &&other)=default
move constructor
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:163
std::vector< double > m_y_sampling
Definition: Cumulative.h:155
void normalize()
Normalize the Cumulative. After calling this function the last vertical value is 1....
Definition: Cumulative.cpp:75
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:87
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:127
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:61
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)
T push_back(T... args)
T size(T... args)