Alexandria 2.31.0
SDC-CH common library for the Euclid project
Loading...
Searching...
No Matches
Piecewise.cpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2012-2022 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
29#include <algorithm>
30
31namespace Euclid {
32namespace MathUtils {
33
35 : PiecewiseBase(std::move(knots)) {
36 if (m_knots.size() - functions.size() != 1) {
37 throw Elements::Exception() << "Invalid number of knots(" << m_knots.size() << ")-functions(" << m_functions.size()
38 << ")";
39 }
40
41 m_functions.reserve(functions.size());
43 [](const std::shared_ptr<Function>& f) { return f->clone(); });
44
45 auto knotsIter = m_knots.begin();
46 while (++knotsIter != m_knots.end()) {
47 if (*knotsIter < *(knotsIter - 1)) {
48 throw Elements::Exception("knots must be increasing");
49 }
50 }
51}
52
54 : PiecewiseBase(std::move(knots)), m_functions{std::move(functions)} {
55 if (m_knots.size() - m_functions.size() != 1) {
56 throw Elements::Exception() << "Invalid number of knots(" << m_knots.size() << ")-functions(" << m_functions.size()
57 << ")";
58 }
59 auto knotsIter = m_knots.begin();
60 while (++knotsIter != m_knots.end()) {
61 if (*knotsIter < *(knotsIter - 1)) {
62 throw Elements::Exception("knots must be increasing");
63 }
64 }
65}
66
70
71double Piecewise::operator()(const double x) const {
72 auto i = findKnot(x);
73 if (i < 0 || i >= static_cast<ssize_t>(m_knots.size())) {
74 return 0.;
75 }
76 if (i == 0) {
77 return (*m_functions[0])(x);
78 }
79 return (*m_functions[i - 1])(x);
80}
81
83 out.resize(xs.size());
84 std::transform(xs.begin(), xs.end(), out.begin(), std::cref(*this));
85}
86
94
95double Piecewise::integrate(const double x1, const double x2) const {
96 if (x1 == x2) {
97 return 0;
98 }
99 int direction = 1;
100 double min = x1;
101 double max = x2;
102 if (min > max) {
103 direction = -1;
104 min = x2;
105 max = x1;
106 }
107 double result = 0;
109 if (knotIter != m_knots.begin()) {
110 --knotIter;
111 }
113 while (++knotIter != m_knots.end()) {
114 auto prevKnotIter = knotIter - 1;
115 if (max <= *prevKnotIter) {
116 break;
117 }
118 if (min < *knotIter) {
119 double down = (min > *prevKnotIter) ? min : *prevKnotIter;
120 double up = (max < *knotIter) ? max : *knotIter;
122 }
123 ++functionIter;
124 }
125 return direction * result;
126}
127
128} // namespace MathUtils
129} // end of namespace Euclid
T back_inserter(T... args)
T begin(T... args)
Represents a piecewise function.
Definition Piecewise.h:48
ssize_t findKnot(double x) const
Definition Piecewise.h:60
std::vector< double > m_knots
A vector where the knots are kept.
Definition Piecewise.h:74
Piecewise(std::vector< double > knots, std::vector< std::shared_ptr< Function > > functions)
Definition Piecewise.cpp:34
double integrate(const double x1, const double x2) const override
Definition Piecewise.cpp:95
double operator()(const double) const override
Definition Piecewise.cpp:71
std::vector< std::unique_ptr< Function > > m_functions
A vector where the sub-functions are kept.
Definition Piecewise.h:136
const std::vector< std::unique_ptr< Function > > & getFunctions() const
Returns the functions in the ranges between the knots.
Definition Piecewise.cpp:67
std::unique_ptr< Function > clone() const override
Definition Piecewise.cpp:87
T end(T... args)
T move(T... args)
ELEMENTS_API double integrate(const Function &function, const double min, const double max, std::unique_ptr< NumericalIntegrationScheme > numericalIntegrationScheme=nullptr)
std::array< std::vector< double >, N > Coordinates
Used to pass the grid coordinates to interpn. Internally will make a copy of the required values.
STL namespace.
T cref(T... args)
T resize(T... args)
T size(T... args)
T transform(T... args)
T upper_bound(T... args)