Alexandria  2.16
Please provide a description of the project.
Piecewise.cpp
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2012-2020 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 
25 #include <algorithm>
29 
30 namespace Euclid {
31 namespace MathUtils {
32 
34  : m_knots{std::move(knots)} {
35  if (m_knots.size() - functions.size() != 1) {
36  throw Elements::Exception() << "Invalid number of knots(" << m_knots.size()
37  << ")-functions(" << m_functions.size() << ")";
38  }
39 
40  m_functions.reserve(functions.size());
41  std::transform(functions.begin(), functions.end(), std::back_inserter(m_functions),
42  [](std::shared_ptr<Function>& f) { return f->clone(); });
43 
44  auto knotsIter = m_knots.begin();
45  while (++knotsIter != m_knots.end()) {
46  if (*knotsIter <= *(knotsIter-1)) {
47  throw Elements::Exception("knots must be strictly increasing");
48  }
49  }
50 }
51 
52 Piecewise::Piecewise(std::vector<double> knots, std::vector<std::unique_ptr<Function>>&& functions)
53  : m_knots{std::move(knots)}, m_functions{std::move(functions)} {
54  if (m_knots.size() - m_functions.size() != 1) {
55  throw Elements::Exception() << "Invalid number of knots(" << m_knots.size()
56  << ")-functions(" << m_functions.size() << ")";
57  }
58  auto knotsIter = m_knots.begin();
59  while (++knotsIter != m_knots.end()) {
60  if (*knotsIter <= *(knotsIter-1)) {
61  throw Elements::Exception("knots must be strictly increasing");
62  }
63  }
64 }
65 
67  return m_knots;
68 }
69 
71  return m_functions;
72 }
73 
74 double Piecewise::operator()(const double x) const {
75  auto knotsBegin = m_knots.begin();
76  if (x < *knotsBegin) {
77  return 0;
78  }
79  if (x == *knotsBegin) {
80  return (*m_functions[0])(x);
81  }
82  auto knotsEnd = m_knots.end();
83  auto findX = std::lower_bound(knotsBegin, knotsEnd, x);
84  if (findX == knotsEnd) {
85  return 0;
86  }
87  return (*m_functions[findX-knotsBegin-1])(x);
88 }
89 
91  std::vector<std::unique_ptr<Function>> cloned_functions;
92  cloned_functions.reserve(m_functions.size());
93  for (auto& f: m_functions) {
94  cloned_functions.emplace_back(f->clone());
95  }
96  return std::unique_ptr<Function> {new Piecewise(m_knots, std::move(cloned_functions))};
97 }
98 
99 double Piecewise::integrate(const double x1, const double x2) const {
100  if (x1 == x2) {
101  return 0;
102  }
103  int direction = 1;
104  double min = x1;
105  double max = x2;
106  if (min > max) {
107  direction = -1;
108  min = x2;
109  max = x1;
110  }
111  double result = 0;
112  auto knotIter = std::upper_bound(m_knots.begin(), m_knots.end(), min);
113  if (knotIter != m_knots.begin()) {
114  --knotIter;
115  }
116  auto functionIter = m_functions.begin() + (knotIter - m_knots.begin());
117  while (++knotIter != m_knots.end()) {
118  auto prevKnotIter = knotIter - 1;
119  if (max <= *prevKnotIter) {
120  break;
121  }
122  if (min < *knotIter) {
123  double down = (min > *prevKnotIter) ? min : *prevKnotIter;
124  double up = (max < *knotIter) ? max : *knotIter;
125  result += Euclid::MathUtils::integrate(**functionIter, down, up);
126  }
127  ++functionIter;
128  }
129  return direction * result;
130 }
131 
132 } // End of MathUtils
133 } // end of namespace Euclid
std::shared_ptr
STL class.
Euclid::MathUtils::Piecewise::m_functions
std::vector< std::unique_ptr< Function > > m_functions
A vector where the sub-functions are kept.
Definition: Piecewise.h:100
std::move
T move(T... args)
std::vector::reserve
T reserve(T... args)
std::vector< double >
std::back_inserter
T back_inserter(T... args)
Piecewise.h
Euclid::MathUtils::integrate
ELEMENTS_API double integrate(const Function &function, const double min, const double max, std::unique_ptr< NumericalIntegrationScheme > numericalIntegrationScheme=nullptr)
Definition: function_tools.cpp:33
Exception.h
Elements::Exception
Euclid::MathUtils::Piecewise::Piecewise
Piecewise(std::vector< double > knots, std::vector< std::shared_ptr< Function >> functions)
Euclid::MathUtils::Piecewise::operator()
double operator()(const double) const override
Definition: Piecewise.cpp:74
std::transform
T transform(T... args)
std::upper_bound
T upper_bound(T... args)
Euclid::MathUtils::Piecewise::m_knots
std::vector< double > m_knots
A vector where the knots are kept.
Definition: Piecewise.h:98
std::vector::emplace_back
T emplace_back(T... args)
Euclid::MathUtils::Piecewise::integrate
double integrate(const double x1, const double x2) const override
Definition: Piecewise.cpp:99
std::lower_bound
T lower_bound(T... args)
Euclid::MathUtils::Piecewise::getFunctions
const std::vector< std::unique_ptr< Function > > & getFunctions() const
Returns the functions in the ranges between the knots.
Definition: Piecewise.cpp:70
std::vector::begin
T begin(T... args)
function_tools.h
Euclid::MathUtils::Piecewise::getKnots
const std::vector< double > & getKnots() const
Returns the knots of the piecewise function.
Definition: Piecewise.cpp:66
std::vector::end
T end(T... args)
std::unique_ptr
STL class.
Euclid
Definition: InstOrRefHolder.h:29
Euclid::MathUtils::Piecewise::clone
std::unique_ptr< Function > clone() const override
Definition: Piecewise.cpp:90