Alexandria  2.18
Please provide a description of the project.
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
interpolation.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 
26 #include "ElementsKernel/Logging.h"
28 #include "implementations.h"
30 
31 namespace Euclid {
32 namespace MathUtils {
33 
34 namespace {
35 
37 
38 } // Anonymous namespace
39 
41  bool extrapolate) {
42 
43  if (x.size() != y.size()) {
44  throw InterpolationException() << "Interpolation using vectors of incompatible "
45  << "size: X=" << x.size() << ", Y=" << y.size();
46  }
47 
48  if (x.size() == 1 && extrapolate) {
49  auto c = y.front();
50  return make_unique<FunctionAdapter>([c](double) { return c; });
51  }
52 
53  // We remove any duplicate lines and we check that we have only increasing
54  // X values and no step functions
55  std::vector<double> final_x;
56  std::vector<double> final_y;
57 
58  final_x.reserve(x.size());
59  final_y.reserve(x.size());
60 
61  final_x.emplace_back(x[0]);
62  final_y.emplace_back(y[0]);
63  for (std::size_t i = 1; i < x.size(); ++i) {
64  if (x[i] == x[i - 1]) {
65  if (y[i] == y[i - 1]) {
66  logger.warn() << "Ignoring duplicate pair (" << x[i] << ", " << y[i] << ") during interpolation";
67  continue;
68  } else {
69  throw InterpolationException() << "Interpolation of step functions is not "
70  << "supported. Entries: (" << x[i] << ", " << y[i] << ") and (" << x[i - 1] << ", "
71  << y[i - 1] << ")";
72  }
73  }
74  if (x[i] < x[i - 1]) {
75  throw InterpolationException() << "Only increasing X values are supported "
76  << "but found " << x[i] << " after " << x[i - 1];
77  }
78  final_x.emplace_back(x[i]);
79  final_y.emplace_back(y[i]);
80  }
81 
82  switch (type) {
84  return linearInterpolation(final_x, final_y, extrapolate);
86  return splineInterpolation(final_x, final_y, extrapolate);
87  }
88  return nullptr;
89 }
90 
94  x.reserve(dataset.size());
95  y.reserve(dataset.size());
96  for (auto& pair : dataset) {
97  x.emplace_back(pair.first);
98  y.emplace_back(pair.second);
99  }
100  return interpolate(x, y, type, extrapolate);
101 }
102 
103 } // namespace MathUtils
104 } // end of namespace Euclid
size_t size() const
Get the size of the vector container.
Definition: XYDataset.h:149
T front(T...args)
void warn(const std::string &logMessage)
static Elements::Logging logger
std::unique_ptr< Function > splineInterpolation(const std::vector< double > &x, const std::vector< double > &y, bool extrapolate)
Performs cubic spline interpolation for the given set of data points.
Definition: spline.cpp:34
T size(T...args)
STL class.
std::unique_ptr< Function > linearInterpolation(const std::vector< double > &x, const std::vector< double > &y, bool extrapolate)
Performs linear interpolation for the given set of data points.
Definition: linear.cpp:34
This module provides an interface for accessing two dimensional datasets (pairs of (X...
Definition: XYDataset.h:59
InterpolationType
Enumeration of the different supported interpolation types.
Definition: interpolation.h:41
ELEMENTS_API std::unique_ptr< Function > interpolate(const std::vector< double > &x, const std::vector< double > &y, InterpolationType type, bool extrapolate=false)
static Logging getLogger(const std::string &name="")
T reserve(T...args)
T emplace_back(T...args)