tinterpolation.cc - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
HTML git clone git://src.adamsgaard.dk/pism
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
tinterpolation.cc (10139B)
---
1 /* Copyright (C) 2015, 2016, 2017, 2018, 2019 PISM Authors
2 *
3 * This file is part of PISM.
4 *
5 * PISM is free software; you can redistribute it and/or modify it under the
6 * terms of the GNU General Public License as published by the Free Software
7 * Foundation; either version 3 of the License, or (at your option) any later
8 * version.
9 *
10 * PISM is distributed in the hope that it will be useful, but WITHOUT ANY
11 * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
12 * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
13 * details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with PISM; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
18 */
19
20 #include <gsl/gsl_interp.h>
21 #include <cassert>
22
23 #include "interpolation.hh"
24 #include "error_handling.hh"
25
26 namespace pism {
27
28 Interpolation::Interpolation(InterpolationType type,
29 const std::vector<double> &input_x,
30 const std::vector<double> &output_x,
31 double period)
32 : Interpolation(type, input_x.data(), input_x.size(),
33 output_x.data(), output_x.size(), period) {
34 // empty
35 }
36
37 Interpolation::Interpolation(InterpolationType type,
38 const double *input_x, unsigned int input_x_size,
39 const double *output_x, unsigned int output_x_size,
40 double period) {
41 switch (type) {
42 case LINEAR:
43 init_linear(input_x, input_x_size, output_x, output_x_size);
44 break;
45 case NEAREST:
46 init_nearest(input_x, input_x_size, output_x, output_x_size);
47 break;
48 case PIECEWISE_CONSTANT:
49 init_piecewise_constant(input_x, input_x_size, output_x, output_x_size);
50 break;
51 case LINEAR_PERIODIC:
52 init_linear_periodic(input_x, input_x_size, output_x, output_x_size, period);
53 break;
54 default:
55 throw RuntimeError(PISM_ERROR_LOCATION, "invalid interpolation type");
56 }
57 }
58
59 /**
60 * Compute linear interpolation indexes and weights.
61 *
62 * @param[in] input_x coordinates of the input grid
63 * @param[in] input_x_size number of points in the input grid
64 * @param[in] output_x coordinates of the output grid
65 * @param[in] output_x_size number of points in the output grid
66 */
67 void Interpolation::init_linear(const double *input_x, unsigned int input_x_size,
68 const double *output_x, unsigned int output_x_size) {
69
70 m_left.resize(output_x_size);
71 m_right.resize(output_x_size);
72 m_alpha.resize(output_x_size);
73
74 // the trivial case (the code below requires input_x_size >= 2)
75 if (input_x_size < 2) {
76 for (unsigned int k = 0; k < output_x_size; ++k) {
77 m_left[k] = 0.0;
78 m_right[k] = 0.0;
79 m_alpha[k] = 0.0;
80 }
81 return;
82 }
83
84 // input grid points have to be stored in the increasing order
85 for (unsigned int i = 0; i < input_x_size - 1; ++i) {
86 if (input_x[i] >= input_x[i + 1]) {
87 throw RuntimeError(PISM_ERROR_LOCATION, "an input grid for linear interpolation has to be "
88 "strictly increasing");
89 }
90 }
91
92 // compute indexes and weights
93 for (unsigned int i = 0; i < output_x_size; ++i) {
94 double x = output_x[i];
95
96 unsigned int
97 L = gsl_interp_bsearch(input_x, x, 0, input_x_size),
98 R = L + 1;
99
100 double alpha = 0.0;
101 if (x >= input_x[L] and R < input_x_size) {
102 // regular case
103 alpha = (x - input_x[L]) / (input_x[R] - input_x[L]);
104 } else {
105 // extrapolation
106 alpha = 0.0;
107 R = L;
108 }
109
110 assert(L < input_x_size);
111 assert(R < input_x_size);
112 assert(alpha >= 0.0 and alpha <= 1.0);
113
114 m_left[i] = L;
115 m_right[i] = R;
116 m_alpha[i] = alpha;
117 }
118
119 init_weights_linear(input_x, input_x_size, output_x, output_x_size);
120 }
121
122 const std::vector<int>& Interpolation::left() const {
123 return m_left;
124 }
125
126 const std::vector<int>& Interpolation::right() const {
127 return m_right;
128 }
129
130 const std::vector<double>& Interpolation::alpha() const {
131 return m_alpha;
132 }
133
134 int Interpolation::left(size_t j) const {
135 return m_left[j];
136 }
137
138 int Interpolation::right(size_t j) const {
139 return m_right[j];
140 }
141
142 double Interpolation::alpha(size_t j) const {
143 return m_alpha[j];
144 }
145
146 std::vector<double> Interpolation::interpolate(const std::vector<double> &input_values) const {
147 std::vector<double> result(m_alpha.size());
148
149 interpolate(input_values.data(), result.data());
150
151 return result;
152 }
153
154 void Interpolation::interpolate(const double *input, double *output) const {
155 size_t n = m_alpha.size();
156 for (size_t k = 0; k < n; ++k) {
157 const int
158 L = m_left[k],
159 R = m_right[k];
160 output[k] = input[L] + m_alpha[k] * (input[R] - input[L]);
161 }
162 }
163
164 void Interpolation::init_nearest(const double *input_x, unsigned int input_x_size,
165 const double *output_x, unsigned int output_x_size) {
166
167 init_linear(input_x, input_x_size, output_x, output_x_size);
168
169 for (unsigned int j = 0; j < m_alpha.size(); ++j) {
170 m_alpha[j] = m_alpha[j] > 0.5 ? 1.0 : 0.0;
171 }
172 }
173
174 void Interpolation::init_piecewise_constant(const double *input_x, unsigned int input_x_size,
175 const double *output_x, unsigned int output_x_size) {
176
177 m_left.resize(output_x_size);
178 m_right.resize(output_x_size);
179 m_alpha.resize(output_x_size);
180
181 // the trivial case
182 if (input_x_size < 2) {
183 for (unsigned int i = 0; i < output_x_size; ++i) {
184 m_left[i] = 0;
185 m_right[i] = 0;
186 m_alpha[i] = 0.0;
187 }
188 return;
189 }
190
191 // input grid points have to be stored in the increasing order
192 for (unsigned int i = 0; i < input_x_size - 1; ++i) {
193 if (input_x[i] >= input_x[i + 1]) {
194 throw RuntimeError(PISM_ERROR_LOCATION, "an input grid for interpolation has to be "
195 "strictly increasing");
196 }
197 }
198
199 // compute indexes and weights
200 for (unsigned int i = 0; i < output_x_size; ++i) {
201
202 size_t L = gsl_interp_bsearch(input_x, output_x[i], 0, input_x_size);
203
204 m_left[i] = L;
205 m_right[i] = L;
206 m_alpha[i] = 0.0;
207
208 assert(m_left[i] >= 0 and m_left[i] < (int)input_x_size);
209 assert(m_right[i] >= 0 and m_right[i] < (int)input_x_size);
210 assert(m_alpha[i] >= 0.0 and m_alpha[i] <= 1.0);
211 }
212 }
213
214 void Interpolation::init_linear_periodic(const double *input_x, unsigned int input_x_size,
215 const double *output_x, unsigned int output_x_size,
216 double period) {
217
218 assert(period > 0);
219
220 m_left.resize(output_x_size);
221 m_right.resize(output_x_size);
222 m_alpha.resize(output_x_size);
223
224 // the trivial case
225 if (input_x_size < 2) {
226 for (unsigned int i = 0; i < output_x_size; ++i) {
227 m_left[i] = 0;
228 m_right[i] = 0;
229 m_alpha[i] = 0.0;
230 }
231 return;
232 }
233
234 // input grid points have to be stored in the increasing order
235 for (unsigned int i = 0; i < input_x_size - 1; ++i) {
236 if (input_x[i] >= input_x[i + 1]) {
237 throw RuntimeError(PISM_ERROR_LOCATION, "an input grid for interpolation has to be "
238 "strictly increasing");
239 }
240 }
241
242 // compute indexes and weights
243 for (unsigned int i = 0; i < output_x_size; ++i) {
244 double x = output_x[i];
245
246 unsigned int L = 0, R = 0;
247 if (x < input_x[0]) {
248 L = input_x_size - 1;
249 R = 0.0;
250 } else {
251 L = gsl_interp_bsearch(input_x, x, 0, input_x_size);
252 R = L + 1 < input_x_size ? L + 1 : 0;
253 }
254
255 double
256 x_l = input_x[L],
257 x_r = input_x[R],
258 alpha = 0.0;
259 if (L < R) {
260 // regular case
261 alpha = (x - x_l) / (x_r - x_l);
262 } else {
263 double
264 x0 = input_x[0],
265 dx = (period - x_l) + x0;
266 assert(dx > 0);
267 if (x > x0) {
268 // interval from the last point of the input grid to the period
269 alpha = (x - x_l) / dx;
270 } else {
271 // interval from 0 to the first point of the input grid
272 alpha = 1.0 - (x_r - x) / dx;
273 }
274 }
275
276 assert(L < input_x_size);
277 assert(R < input_x_size);
278 assert(alpha >= 0.0 and alpha <= 1.0);
279
280 m_left[i] = L;
281 m_right[i] = R;
282 m_alpha[i] = alpha;
283 }
284 }
285
286 /*!
287 * Initialize integration weights (trapezoid rule).
288 */
289 void Interpolation::init_weights_linear(const double *x,
290 unsigned int x_size,
291 const double *output_x,
292 unsigned int output_x_size) {
293
294 if (output_x_size == 1) {
295 m_w = {0.0};
296 m_interval_length = 0.0;
297 return;
298 }
299
300 int
301 N = output_x_size - 1,
302 al = m_left[0],
303 ar = m_right[0],
304 bl = m_left[N],
305 br = m_right[N];
306
307 double
308 a = output_x[0],
309 b = output_x[N],
310 alpha_a = m_alpha[0],
311 alpha_b = m_alpha[N];
312
313 if (a >= b) {
314 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "invalid interval: (%f, %f)", a, b);
315 }
316
317 m_interval_length = b - a;
318
319 m_w.resize(x_size);
320 for (unsigned int k = 0; k < x_size; ++k) {
321 m_w[k] = 0.0;
322 }
323
324 if (al == bl and ar == br) {
325 // both end points are in the same interval
326
327 m_w[al] += 0.5 * (2.0 - alpha_a - alpha_b) * (b - a);
328 m_w[ar] += 0.5 * (alpha_a + alpha_b) * (b - a);
329 } else {
330 // first interval
331 m_w[al] += 0.5 * (1.0 - alpha_a) * (x[ar] - a);
332 m_w[ar] += 0.5 * (1.0 + alpha_a) * (x[ar] - a);
333
334 // intermediate intervals
335 for (int k = ar; k < bl; ++k) {
336 int
337 L = k,
338 R = k + 1;
339 m_w[L] += 0.5 * (x[R] - x[L]);
340 m_w[R] += 0.5 * (x[R] - x[L]);
341 }
342
343 // last interval
344 m_w[bl] += 0.5 * (2.0 - alpha_b) * (b - x[bl]);
345 m_w[br] += 0.5 * alpha_b * (b - x[bl]);
346 }
347 }
348
349 double Interpolation::integral(const double *input) const {
350 double result = 0.0;
351 for (size_t k = 0; k < m_w.size(); ++k) {
352 result += input[k] * m_w[k];
353 }
354 return result;
355 }
356
357 double Interpolation::integral(const std::vector<double> &input) const {
358 return integral(input.data());
359 }
360
361 double Interpolation::interval_length() const {
362 return m_interval_length;
363 }
364
365
366 } // end of namespace pism