forked from npshub/mantid
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Interpolation.cpp
194 lines (159 loc) · 5.39 KB
/
Interpolation.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
// NScD Oak Ridge National Laboratory, European Spallation Source,
// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
// SPDX - License - Identifier: GPL - 3.0 +
#include "MantidKernel/Interpolation.h"
#include "MantidKernel/Logger.h"
#include "MantidKernel/StringTokenizer.h"
#include "MantidKernel/UnitFactory.h"
#include <algorithm>
#include <iterator>
#include <stdexcept>
namespace Mantid {
namespace Kernel {
namespace {
/// static logger
Logger g_log("Interpolation");
} // namespace
/* Functor used in std::lower_bound to replicate the original behavior.
*/
struct LessOrEqualFunctor {
bool operator()(const DataXY &lhs, const DataXY &rhs) { return lhs.first <= rhs.first; }
};
/** Constructor default to linear interpolation and x-unit set to TOF
*/
Interpolation::Interpolation() : m_method("linear") {
m_xUnit = UnitFactory::Instance().create("TOF");
m_yUnit = UnitFactory::Instance().create("TOF");
}
/** Get iterator of item that is next larger than the supplied x value
* @param key :: the x value to base the search on
* @return iterator of the next largest x value
*/
std::vector<DataXY>::const_iterator Interpolation::findIndexOfNextLargerValue(double key) const {
return std::lower_bound(m_data.begin(), m_data.end(), DataXY(key, 0), LessOrEqualFunctor());
}
std::vector<DataXY>::const_iterator Interpolation::cbegin() const { return m_data.cbegin(); }
std::vector<DataXY>::const_iterator Interpolation::cend() const { return m_data.cend(); }
void Interpolation::setXUnit(const std::string &unit) { m_xUnit = UnitFactory::Instance().create(unit); }
void Interpolation::setYUnit(const std::string &unit) { m_yUnit = UnitFactory::Instance().create(unit); }
/** Get interpolated value at location at
* @param at :: Location where to get interpolated value
* @return the value
*/
double Interpolation::value(const double &at) const {
size_t N = m_data.size();
if (N == 0) {
g_log.error() << "Need at least one value for interpolation. Return "
"interpolation value zero.";
return 0.0;
}
if (N == 1) {
return m_data[0].second;
}
// check first if at is within the limits of interpolation interval
if (at < m_data[0].first) {
return m_data[0].second -
(m_data[0].first - at) * (m_data[1].second - m_data[0].second) / (m_data[1].first - m_data[0].first);
}
if (at > m_data[N - 1].first) {
return m_data[N - 1].second + (at - m_data[N - 1].first) * (m_data[N - 1].second - m_data[N - 2].second) /
(m_data[N - 1].first - m_data[N - 2].first);
}
try {
// otherwise
// General case. Find index of next largest value by std::lower_bound.
auto pos = findIndexOfNextLargerValue(at);
auto posBefore = std::prev(pos);
if (posBefore->first == at) {
return posBefore->second;
} else {
double interpolatedY = posBefore->second + (at - posBefore->first) * (pos->second - posBefore->second) /
(pos->first - posBefore->first);
return interpolatedY;
}
} catch (const std::range_error &) {
return 0.0;
}
}
/** Add point in the interpolation.
*
* @param xx :: x-value
* @param yy :: y-value
*/
void Interpolation::addPoint(const double &xx, const double &yy) {
size_t N = m_data.size();
std::vector<DataXY>::const_iterator it;
if (N == 0) {
DataXY newpair(xx, yy);
m_data.emplace_back(newpair);
return;
}
// check first if xx is within the limits of interpolation interval
if (xx < m_data[0].first) {
it = m_data.begin();
it = m_data.insert(it, DataXY(xx, yy));
return;
}
if (xx > m_data[N - 1].first) {
m_data.emplace_back(xx, yy);
return;
}
// otherwise
it = findIndexOfNextLargerValue(xx);
auto posBefore = std::prev(it);
if (posBefore->first != xx) {
m_data.insert(it, DataXY(xx, yy));
}
}
/**
Prints object to stream
@param os :: the Stream to output to
*/
void Interpolation::printSelf(std::ostream &os) const {
os << m_method << " ; " << m_xUnit->unitID() << " ; " << m_yUnit->unitID();
for (unsigned int i = 0; i < m_data.size(); i++) {
os << " ; " << m_data[i].first << " " << m_data[i].second;
}
}
/**
Resets interpolation data by clearing the internal storage for x- and y-values
*/
void Interpolation::resetData() { m_data.clear(); }
/**
Prints the value of parameter
@param os :: the Stream to output to
@param f :: the FitParameter to output
@return the output stream
*/
std::ostream &operator<<(std::ostream &os, const Interpolation &f) {
f.printSelf(os);
return os;
}
/**
Reads in parameter value
@param in :: Input Stream
@param f :: FitParameter to write to
@return Current state of stream
*/
std::istream &operator>>(std::istream &in, Interpolation &f) {
using tokenizer = Mantid::Kernel::StringTokenizer;
std::string str;
getline(in, str);
tokenizer values(str, ";", tokenizer::TOK_TRIM);
f.setMethod(values[0]);
f.setXUnit(values[1]);
f.setYUnit(values[2]);
f.resetData(); // Reset data, in case the interpolation table is not empty
for (unsigned int i = 3; i < values.count(); i++) {
std::stringstream strstream(values[i]);
double x, y;
strstream >> x >> y;
f.addPoint(x, y);
}
return in;
}
} // namespace Kernel
} // namespace Mantid