Alexandria  2.27.0
SDC-CH common library for the Euclid project
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 
31 namespace Euclid {
32 namespace 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());
42  std::transform(functions.begin(), functions.end(), std::back_inserter(m_functions),
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 
68  return m_functions;
69 }
70 
71 double 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 
88  std::vector<std::unique_ptr<Function>> cloned_functions;
89  cloned_functions.reserve(m_functions.size());
90  std::transform(m_functions.begin(), m_functions.end(), std::back_inserter(cloned_functions),
91  [](const std::unique_ptr<Function>& f) { return f->clone(); });
92  return make_unique<Piecewise>(m_knots, std::move(cloned_functions));
93 }
94 
95 double 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;
108  auto knotIter = std::upper_bound(m_knots.begin(), m_knots.end(), min);
109  if (knotIter != m_knots.begin()) {
110  --knotIter;
111  }
112  auto functionIter = m_functions.begin() + (knotIter - m_knots.begin());
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;
121  result += Euclid::MathUtils::integrate(**functionIter, down, up);
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
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
Piecewise(std::vector< double > knots, std::vector< std::shared_ptr< Function >> functions)
Definition: Piecewise.cpp:34
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)
STL namespace.
T cref(T... args)
T reserve(T... args)
T resize(T... args)
T size(T... args)
T transform(T... args)
T upper_bound(T... args)