Alexandria  2.16
Please provide a description of the project.
interpolation.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 
26 #include "ElementsKernel/Logging.h"
29 #include "implementations.h"
30 
31 namespace Euclid {
32 namespace MathUtils {
33 
34 namespace {
35 
37 
38 } // Anonymous namespace
39 
41  InterpolationType type, 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]
67  << ") during interpolation";
68  continue;
69  } else {
70  throw InterpolationException() << "Interpolation of step functions is not "
71  << "supported. Entries: (" << x[i] << ", " << y[i] << ") and ("
72  << x[i-1] << ", " << y[i-1] << ")";
73  }
74  }
75  if (x[i] < x[i-1]) {
76  throw InterpolationException() << "Only increasing X values are supported "
77  << "but found " << x[i] << " after " << x[i-1];
78  }
79  final_x.emplace_back(x[i]);
80  final_y.emplace_back(y[i]);
81  }
82 
83  switch (type) {
85  return linearInterpolation(final_x, final_y, extrapolate);
87  return splineInterpolation(final_x, final_y, extrapolate);
88  }
89  return nullptr;
90 }
91 
93  bool extrapolate) {
96  x.reserve(dataset.size());
97  y.reserve(dataset.size());
98  for (auto& pair : dataset) {
99  x.emplace_back(pair.first);
100  y.emplace_back(pair.second);
101  }
102  return interpolate(x, y, type, extrapolate);
103 }
104 
105 } // End of MathUtils
106 } // end of namespace Euclid
T front(T... args)
size_t size() const
Get the size of the vector container.
Definition: XYDataset.h:152
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)