Alexandria 2.31.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"
29#include <cstdlib> // for size_t
30
31namespace Euclid {
32namespace MathUtils {
33
35 : m_x_sampling(x_sampling), m_y_sampling(y_sampling) {
36 if (x_sampling.size() != y_sampling.size()) {
37 throw Elements::Exception("Cumulative: X and Y sampling do not have the same length.");
38 }
39}
40
41Cumulative::Cumulative(const XYDataset::XYDataset& sampling) : m_x_sampling{}, m_y_sampling{} {
42 auto iter = sampling.begin();
43 while (iter != sampling.end()) {
44 m_x_sampling.push_back((*iter).first);
45 m_y_sampling.push_back((*iter).second);
46 ++iter;
47 }
48}
49
53 auto iter = sampling.begin();
54 while (iter != sampling.end()) {
55 xs.push_back((*iter).first);
56 ys.push_back((*iter).second);
57 ++iter;
58 }
59 return Cumulative::fromPdf(xs, ys);
60}
61
75
77 double total = m_y_sampling.back();
79 auto iter = m_y_sampling.begin();
80 while (iter != m_y_sampling.end()) {
81 cumul.push_back(*iter / total);
82 ++iter;
83 }
84
86}
87
88double Cumulative::findValue(double ratio, TrayPosition position) const {
89
90 if (ratio > 1. || ratio < 0.) {
91 throw Elements::Exception("Cumulative::findValue : ratio parameter must be in range [0,1]");
92 }
93
94 double value = m_y_sampling.back() * ratio;
95
96 auto iter_x = m_x_sampling.cbegin();
97 auto iter_y = m_y_sampling.cbegin();
98 while (iter_y != m_y_sampling.cend() && (*iter_y) < value) {
99 ++iter_x;
100 ++iter_y;
101 }
102
103 double begin_value = *iter_x;
104 double tray = *iter_y;
105 while (iter_y != m_y_sampling.cend() && (*iter_y) == tray) {
106 ++iter_x;
107 ++iter_y;
108 }
109
110 double end_value = *(--iter_x);
111
112 double result = 0;
113 switch (position) {
116 break;
118 result = 0.5 * (begin_value + end_value);
119 break;
122 break;
123 }
124
125 return result;
126}
127
129
130 if (rate > 1. || rate <= 0.) {
131 throw Elements::Exception("Cumulative::findMinInterval : rate parameter must be in range ]0,1]");
132 }
133
135
137
139 auto iter_2_x = ++(m_x_sampling.cbegin());
141 auto iter_2_y = ++(m_y_sampling.cbegin());
142 auto min_x = m_x_sampling.cbegin();
143 auto max_x = --(m_x_sampling.cend());
144 while (iter_1_x != m_x_sampling.cend()) {
146 ++iter_2_x;
147 ++iter_2_y;
148 }
149 if (iter_2_x == m_x_sampling.cend()) {
150 break;
151 }
152 if ((*iter_2_x - *iter_1_x) <= (*max_x - *min_x)) {
153 max_x = iter_2_x;
154 min_x = iter_1_x;
155 }
156 ++iter_1_x;
157 ++iter_1_y;
158 first_correction = 0.;
159 }
160
161 return std::make_pair(*min_x, *max_x);
162}
163
165
166 if (rate > 1. || rate <= 0.) {
167 throw Elements::Exception("Cumulative::findCenteredInterval : rate parameter must be in range ]0,1]");
168 }
169
170 double min_value = m_y_sampling.back() * (0.5 - rate / 2.0);
171 double max_value = m_y_sampling.back() * (0.5 + rate / 2.0);
172
173 auto iter_x = m_x_sampling.cbegin();
174 auto iter_y = m_y_sampling.cbegin();
175 while (iter_y != m_y_sampling.cend() && (*iter_y) < min_value) {
176 ++iter_x;
177 ++iter_y;
178 }
179
180 if ((*iter_y) < max_value) {
181 double tray = *iter_y;
182 while (iter_y != m_y_sampling.cend() && (*iter_y) == tray) {
183 ++iter_x;
184 ++iter_y;
185 }
186 double min_x = (iter_x != m_x_sampling.begin()) ? *(iter_x - 1) : *iter_x;
187
188 while (iter_y != m_y_sampling.cend() && (*iter_y) < max_value) {
189 ++iter_x;
190 ++iter_y;
191 }
192 double max_x = *iter_x;
193
194 return std::make_pair(min_x, max_x);
195 } else {
196 double min_x = (iter_x != m_x_sampling.begin()) ? *(iter_x - 1) : *iter_x;
197 double max_x = *iter_x;
198
199 return std::make_pair(min_x, max_x);
200 }
201}
202
203
204double Cumulative::eval(double x_value) const {
206 throw Elements::Exception("Cumulative::eval : provided value is outside of the x range");
207 }
208
210 return (*function)(x_value);
211}
212
213} // namespace MathUtils
214} // namespace Euclid
T back(T... args)
T begin(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:163
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...
std::vector< double > m_y_sampling
Definition Cumulative.h:164
void normalize()
Normalize the Cumulative. After calling this function the last vertical value is 1....
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....
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...
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...
double eval(double x_value) const
return the value of the cumulative at a given value of the horizontal axis. If the value do not match...
This module provides an interface for accessing two dimensional datasets (pairs of (X,...
Definition XYDataset.h:59
T end(T... args)
T front(T... args)
T make_pair(T... args)
T move(T... args)
ELEMENTS_API std::unique_ptr< Function > interpolate(const std::vector< double > &x, const std::vector< double > &y, InterpolationType type, bool extrapolate=false)
std::array< std::vector< double >, N > Coordinates
Used to pass the grid coordinates to interpn. Internally will make a copy of the required values.
T push_back(T... args)
T size(T... args)