Branch data Line data Source code
1 : : /**
2 : : *
3 : : *
4 : : * @author: William Kingsford (namespace created by Abed Yassine)
5 : : * @data: October 18th, 2019
6 : : *
7 : : * @file interpolation_engine.cxx
8 : : * @brief: contains implementations for the interpolation_engine methods
9 : : */
10 : :
11 : :
12 : : #include "interpolation_engine.h"
13 : : #include "profiler.h"
14 : :
15 : : #include <cmath>
16 : :
17 : : namespace SP_IE
18 : : {
19 : : RaisedCosLUT *RaisedCosLUT::_s_instance = 0;
20 : :
21 : 0 : double sinc(double x) {
22 [ # # ]: 0 : if (x == 0) {
23 : 0 : return 1;
24 : : }
25 : : else {
26 : 0 : return sin(M_PI * x) / (M_PI * x);
27 : : }
28 : : }
29 : :
30 : 0 : double raised_cosine(double x, double beta) {
31 [ # # ]: 0 : if (abs(x) == 1.0 / (2.0 * beta)) {
32 : 0 : return (M_PI / (4.0)) * sinc(1/(2.0*beta));
33 : : }
34 : : else {
35 : 0 : return sinc(x) * cos(M_PI*beta*x) / (1.0 - 4.0*beta*beta*x*x);
36 : : }
37 : : }
38 : :
39 : 0 : double integrate_midpoint(std::function<double(double)> func, double l, double r, unsigned num_pts)
40 : : {
41 : : /*
42 : : * This integration assumes that the delta_x is fixed. i.e, all points along the line are equidistant
43 : : *
44 : : */
45 : :
46 : 0 : double integration = 0.0;
47 : :
48 : 0 : double temp_r = r;
49 : 0 : r = std::max(r, l);
50 : 0 : l = std::min(temp_r, l);
51 : :
52 [ # # ]: 0 : if (num_pts == 0)
53 : : {
54 [ # # ]: 0 : fprintf(stderr, "\033[1;%dmintegrate_midpoint::num_points is zero ==> cannot integrate.\033[0;m\n", 31);
55 : 0 : std::exit(-1);
56 : : }
57 : :
58 : 0 : double delta_x = (r - l)/num_pts;
59 : :
60 : 0 : double x = l + 0.5*delta_x;
61 [ # # ]: 0 : for (unsigned i = 0; i < num_pts; i++)
62 : : {
63 [ # # ]: 0 : integration += func(x);
64 : 0 : x += delta_x;
65 : : }
66 : :
67 : 0 : return integration*delta_x;
68 : : }
69 : :
70 : 0 : Eigen::RowVectorXd integrate_midpoint(std::function<Eigen::RowVectorXd(Eigen::RowVectorXd)> func,
71 : : Eigen::RowVectorXd l, Eigen::RowVectorXd r, unsigned num_pts)
72 : : {
73 [ # # ][ # # ]: 0 : Eigen::RowVectorXd vec = r - l;
74 [ # # ]: 0 : double delta_x = vec.norm()/num_pts;
75 : :
76 [ # # ][ # # ]: 0 : std::cout << "l: " << l << std::endl;
[ # # ]
77 [ # # ]: 0 : vec.normalize();
78 : :
79 [ # # ][ # # ]: 0 : Eigen::RowVectorXd x = l + 0.5*delta_x*vec;
[ # # ]
80 : :
81 [ # # ][ # # ]: 0 : Eigen::RowVectorXd integration = func(x);
82 : :
83 [ # # ]: 0 : for (unsigned i = 1; i < num_pts; i++)
84 : : {
85 [ # # ][ # # ]: 0 : std::cout << "x: " << x << std::endl;
[ # # ]
86 [ # # ][ # # ]: 0 : x += delta_x*vec;
87 [ # # ][ # # ]: 0 : integration += func(x);
[ # # ]
88 : : }
89 : :
90 [ # # ][ # # ]: 0 : return integration*delta_x;
91 : : }
92 : :
93 : 0 : double evaluate_interpolation(double x, double x_min, double x_max, double delta_x,
94 : : vector<double> &f, int N, bool use_raised_cosine, double beta) {
95 : : // determine nearest sample point index to the left of x
96 : 0 : int idx = (int) floor((x - x_min)/delta_x);
97 : : // convolution will include N sample points on each side of x
98 : 0 : int begin = idx - N + 1;
99 : 0 : int end = idx + N;
100 : :
101 : : // ensure sum is within the range of the sample function
102 [ # # ]: 0 : if (begin < 0) {
103 : 0 : std::cout << "WARNING: Sum in evaluate_interpolation is being pruned " <<
104 : 0 : " due to going outside of range (x = " << x <<
105 : 0 : ", x_min = " << x_min << ", x_max = " << x_max << ")" << std::endl;
106 : 0 : begin = 0;
107 : : }
108 [ # # ]: 0 : if (end >= static_cast<int>(f.size())) {
109 : 0 : std::cout << "WARNING: Sum in evaluate_interpolation is being " <<
110 : 0 : " pruned due to going outside of range (x = " << x <<
111 : 0 : ", x_min = " << x_min << ", x_max = " << x_max << ")" << std::endl;
112 : 0 : end = static_cast<int>(f.size()) - 1;
113 : : }
114 : :
115 : 0 : double ret_val = 0;
116 : 0 : double a = (x - x_min)/delta_x; // pre-computing for efficiency
117 [ # # ]: 0 : if (!use_raised_cosine) { // use sinc
118 [ # # ]: 0 : for (int i=begin; i<=end; i++) {
119 : 0 : ret_val += f[i] * sinc(a - (double)i);
120 : : }
121 : : }
122 : : else { // use raised cosine
123 [ # # ]: 0 : for (int i=begin; i<=end; i++) {
124 : 0 : ret_val += f[i] * raised_cosine(a - (double)i, beta);
125 : : }
126 : : }
127 : :
128 : 0 : return ret_val;
129 : : }
130 : :
131 : 0 : double evaluate_interpolation(double x[3], double x_min[3], double x_max[3], double delta_x,
132 : : vector<vector<vector<double>>> &f, int N,
133 : : bool use_raised_cosine, double beta) {
134 : : // determine nearest sample point index to the left of x
135 : : int idx[3], begin[3], end[3];
136 : :
137 : : // number of entries in f array across each dimension
138 : : int length[3];
139 : 0 : length[0] = static_cast<int>(f.size());
140 : 0 : length[1] = static_cast<int>(f[0].size());
141 : 0 : length[2] = static_cast<int>(f[0][0].size());
142 [ # # ]: 0 : for (int i=0; i<3; i++) {
143 : 0 : idx[i] = (int) floor((x[i] - x_min[i])/delta_x);
144 : :
145 : : // convolution will include N sample points on each side of x
146 : 0 : begin[i] = idx[i] - N + 1;
147 : 0 : end[i] = idx[i] + N;
148 : :
149 : : // ensure sum is within the range of the sample function
150 [ # # ]: 0 : if (begin[i] < 0) {
151 [ # # ]: 0 : std::cout << "WARNING: Sum in evaluate_interpolation is being " <<
152 [ # # ][ # # ]: 0 : "pruned due to going outside of range (x[" << i << "] = " << x[i] <<
[ # # ][ # # ]
153 [ # # ][ # # ]: 0 : ", x_min[" << i << "] = " << x_min[i] << ", x_max[" << i << "] = " << x_max[i] << ")" << std::endl;
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
154 : 0 : begin[i] = 0;
155 : : }
156 [ # # ]: 0 : if (end[i] >= static_cast<int>(f.size())) {
157 [ # # ]: 0 : std::cout << "WARNING: Sum in evaluate_interpolation is being " <<
158 [ # # ][ # # ]: 0 : "pruned due to going outside of range (x[" << i << "] = " << x[i] <<
[ # # ][ # # ]
159 [ # # ][ # # ]: 0 : ", x_min[" << i << "] = " << x_min[i] << ", x_max[" << i << "] = " << x_max[i] << ")" << std::endl;
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
160 : 0 : end[i] = length[i] - 1;
161 : : }
162 : : }
163 : :
164 : 0 : double ret_val = 0;
165 : : // pre-computing for efficiency
166 : 0 : double a = (x[0] - x_min[0])/delta_x;
167 : 0 : double b = (x[1] - x_min[1])/delta_x;
168 : 0 : double c = (x[2] - x_min[2])/delta_x;
169 : :
170 [ # # ]: 0 : RaisedCosLUT *lut = RaisedCosLUT::instance();
171 : : // check if initialized already; if not, initialize with default settings
172 [ # # ]: 0 : if (!lut->isInitialized()) {
173 [ # # ]: 0 : lut->initialize((double)(N+1), 100.0, beta); // range of [-(N+1),N+1], resolution of 100
174 : : }
175 : :
176 [ # # ]: 0 : if (!use_raised_cosine) { // use sinc
177 [ # # ]: 0 : for (int i=begin[0]; i<=end[0]; i++) {
178 [ # # ]: 0 : for (int j=begin[1]; j<=end[1]; j++) {
179 [ # # ]: 0 : for (int k=begin[2]; k<=end[2]; k++) {
180 [ # # ]: 0 : ret_val += f[i][j][k] * sinc(a - (double)i) *
181 [ # # ]: 0 : sinc(b - (double)j) *
182 [ # # ]: 0 : sinc(c - (double)k);
183 : : }
184 : : }
185 : : }
186 : : }
187 : : else { // use raised cosine
188 [ # # ]: 0 : for (int i=begin[0]; i<=end[0]; i++) {
189 [ # # ]: 0 : for (int j=begin[1]; j<=end[1]; j++) {
190 [ # # ]: 0 : for (int k=begin[2]; k<=end[2]; k++) {
191 : : // we know that abs(a-(double)i) < N+1, so we don't need to worry about
192 : : // going out of the array's bounds
193 : 0 : ret_val += f[i][j][k] * lut->get(a-(double)i) *
194 : 0 : lut->get(b-(double)j) *
195 : 0 : lut->get(c-(double)k);
196 : : //ret_val += f[i][j][k] * raised_cosine(a-(double)i, beta) *
197 : : // raised_cosine(b-(double)j, beta) *
198 : : // raised_cosine(c-(double)k, beta);
199 : : }
200 : : }
201 : : }
202 : : }
203 : :
204 : 0 : return ret_val;
205 : : }
206 : :
207 : : // overloaded version for vector-valued functions of 2 variables
208 : 0 : Eigen::RowVectorXd evaluate_interpolation(double x[2], double x_min[2],
209 : : double delta_x,
210 : : source_cache &grid_sources, int N,
211 : : bool use_raised_cosine, double beta) {
212 : : // determine nearest sample point index to the left of x
213 : : int idx[2], begin[2], end[2];
214 : :
215 : :
216 : : // number of entries in f array across each dimension
217 : : //int length[2];
218 : : //length[0] = f.size();
219 : : //length[1] = f[0].size();
220 [ # # ]: 0 : for (int i=0; i<2; i++) {
221 : 0 : idx[i] = (int) floor((x[i] - x_min[i])/delta_x);
222 : :
223 : : // convolution will include N sample points on each side of x
224 : 0 : begin[i] = idx[i] - N + 1;
225 : 0 : end[i] = idx[i] + N;
226 : :
227 : : // TODO: no longer need to check bounds now that we're using source_cache
228 : : // ensure sum is within the range of the sample function
229 : : /*if (begin[i] < 0) {
230 : : std::cout << "WARNING: Sum in evaluate_interpolation is being "
231 : : "pruned due to going outside of range (x[" << i << "] = " <<
232 : : x[i] << ", x_min[" << i << "] = " << x_min[i] << ", x_max[" << i << "] = " << x_max[i] << ")" << std::endl;
233 : : begin[i] = 0;
234 : : }
235 : : if (end[i] >= f.size()) {
236 : : std::cout << "WARNING: Sum in evaluate_interpolation is being "
237 : : "pruned due to going outside of range (x[" << i << "] = " <<
238 : : x[i] << ", x_min[" << i << "] = " << x_min[i] << ", x_max[" << i << "] = " << x_max[i] << ")" << std::endl;
239 : : end[i] = length[i] - 1;
240 : : }*/
241 : : }
242 : :
243 : : // initialize return value
244 [ # # ][ # # ]: 0 : Eigen::RowVectorXd ret_val = Eigen::RowVectorXd::Zero(grid_sources.get_num_tetras());
245 : : // pre-computing for efficiency
246 : 0 : double a = (x[0] - x_min[0])/delta_x;
247 : 0 : double b = (x[1] - x_min[1])/delta_x;
248 : : double weight;
249 : :
250 [ # # ]: 0 : RaisedCosLUT *lut = RaisedCosLUT::instance();
251 : : // check if initialized already; if not, initialize with default settings
252 [ # # ]: 0 : if (!lut->isInitialized()) {
253 : : //TODO: currently making larger than we'll ever need, should adjust this so it can
254 : : // check if the size is too small (for when N changes)
255 : : //lut->initialize((double)(N+1), 100.0, beta); // range of [-(N+1),N+1], resolution of 100
256 [ # # ]: 0 : lut->initialize((double)(11), 100.0, beta); // range of [-(11),11], resolution of 100
257 : : }
258 : :
259 : : // TODO: using only linear interpolation for now
260 : :
261 : : // parameters for linear interpolation
262 : : double t_x, t_y;
263 : : int n_x, n_y;
264 : 0 : n_x = static_cast<int>(floor(a));
265 : 0 : n_y = static_cast<int>(floor(b));
266 : 0 : t_x = a - (double)n_x;
267 : 0 : t_y = b - (double)n_y;
268 : :
269 : : //cout << "Evaluating linear interpolation: t_x = " << t_x << ", t_y = " << t_y << endl;
270 : : // TODO: no longer need to cast n_x, n_y to int
271 : :
272 [ # # ]: 0 : if (!use_raised_cosine) { // use linear interpolation
273 [ # # ]: 0 : ret_val = (1.0 - t_y) * (
274 [ # # ][ # # ]: 0 : (1.0 - t_x) * grid_sources.get_source((int)n_x, (int)n_y) +
[ # # ]
275 [ # # ][ # # ]: 0 : t_x * grid_sources.get_source((int)n_x + 1, (int)n_y) ) +
[ # # ]
276 [ # # ]: 0 : t_y * (
277 [ # # ][ # # ]: 0 : (1.0 - t_x) * grid_sources.get_source((int)n_x, (int)n_y+1) +
[ # # ]
278 [ # # ][ # # ]: 0 : t_x * grid_sources.get_source((int)n_x + 1, (int)n_y+1) );
[ # # ]
279 : : // TODO: old code for sinc
280 : : /*
281 : : for (int i=begin[0]; i<=end[0]; i++) {
282 : : for (int j=begin[1]; j<=end[1]; j++) {
283 : : for (int k=begin[2]; k<=end[2]; k++) {
284 : : weight = sinc(a - (double)i) * sinc(b - (double)j) * sinc(c - (double)k);
285 : : ret_val += weight * f[i][j][k];
286 : : }
287 : : }
288 : : }*/
289 : : }
290 : : else { // use raised cosine
291 [ # # ]: 0 : for (int i=begin[0]; i<=end[0]; i++) {
292 [ # # ]: 0 : for (int j=begin[1]; j<=end[1]; j++) {
293 : : // we know that abs(a-(double)i) < N+1, so we don't need to worry about
294 : : // going out of the array's bounds
295 : 0 : weight = lut->get(a-(double)i) *
296 : 0 : lut->get(b-(double)j);
297 : : //ret_val += weight * f[i][j][k];
298 [ # # ][ # # ]: 0 : ret_val += weight * grid_sources.get_source(i, j);
[ # # ]
299 : : }
300 : : }
301 : : }
302 : :
303 : 0 : return ret_val;
304 : : }
305 : :
306 : : // TODO: replace double[2] with std::vector
307 : 0 : vector<Eigen::RowVectorXd> evaluate_interpolation_gradient(double x[2], double x_min[2],
308 : : double delta_x,
309 : : source_cache &grid_sources, int N,
310 : : bool /*use_raised_cosine*/, double beta) {
311 : : // TODO: don't need the below now that we're using source_cache
312 : : /*
313 : : // determine nearest sample point index to the left of x
314 : : int idx[2], begin[2], end[2];
315 : :
316 : : // number of entries in f array across each dimension
317 : : int length[2];
318 : : length[0] = f.size();
319 : : length[1] = f[0].size();
320 : : for (int i=0; i<2; i++) {
321 : : idx[i] = (int) floor((x[i] - x_min[i])/delta_x);
322 : :
323 : : // convolution will include N sample points on each side of x
324 : : begin[i] = idx[i] - N + 1;
325 : : end[i] = idx[i] + N;
326 : :
327 : : // ensure sum is within the range of the sample function
328 : : if (begin[i] < 0) {
329 : : std::cout << "WARNING: Sum in evaluate_interpolation is being " <<
330 : : "pruned due to going outside of range (x[" << i << "] = " <<
331 : : x[i] << ", x_min[" << i << "] = " << x_min[i] << ", x_max[" << i << "] = " << x_max[i] << ")" << std::endl;
332 : : begin[i] = 0;
333 : : }
334 : : if (end[i] >= f.size()) {
335 : : std::cout << "WARNING: Sum in evaluate_interpolation is being " <<
336 : : "pruned due to going outside of range (x[" << i << "] = " <<
337 : : x[i] << ", x_min[" << i << "] = " << x_min[i] << ", x_max[" << i << "] = " << x_max[i] << ")" << std::endl;
338 : : end[i] = length[i] - 1;
339 : : }
340 : : }*/
341 : :
342 : : // initialize return value
343 : : // TODO: replace 2 by the number of components in f
344 [ # # ][ # # ]: 0 : vector<Eigen::RowVectorXd> ret_val(2, Eigen::RowVectorXd::Zero(grid_sources.get_num_tetras()));
[ # # ]
345 : : // pre-computing for efficiency
346 : 0 : double a = (x[0] - x_min[0])/delta_x;
347 : 0 : double b = (x[1] - x_min[1])/delta_x;
348 : :
349 [ # # ]: 0 : RaisedCosLUT *lut = RaisedCosLUT::instance();
350 : : // check if initialized already; if not, initialize with default settings
351 [ # # ]: 0 : if (!lut->isInitialized()) {
352 [ # # ]: 0 : lut->initialize((double)(N+1), 100.0, beta); // range of [-(N+1),N+1], resolution of 100
353 : : }
354 : :
355 : : // TODO: using only linear interpolation for now
356 : :
357 : : // parameters for linear interpolation
358 : : double t_x, t_y;
359 : : int n_x, n_y;
360 : 0 : n_x = static_cast<int>(floor(a));
361 : 0 : n_y = static_cast<int>(floor(b));
362 : 0 : t_x = a - (double)n_x;
363 : 0 : t_y = b - (double)n_y;
364 : :
365 : : //cout << "Evaluating linear interpolation: t_x = " << t_x << ", t_y = " << t_y << endl;
366 : : // TODO: use intermediate variables to make this less ugly (e.g. variables so we don't call
367 : : // grid_sources.get_source so often)
368 [ # # ][ # # ]: 0 : ret_val[0] = (1.0 / delta_x) * ( (1.0 - t_y) * (
369 [ # # ][ # # ]: 0 : grid_sources.get_source(n_x + 1, n_y) - grid_sources.get_source(n_x, n_y) ) +
[ # # ][ # # ]
370 [ # # ]: 0 : t_y * (
371 [ # # ][ # # ]: 0 : grid_sources.get_source(n_x + 1, n_y+1) - grid_sources.get_source(n_x, n_y+1) )
[ # # ]
372 [ # # ]: 0 : );
373 [ # # ]: 0 : ret_val[1] = (1.0 / delta_x) * (
374 [ # # ][ # # ]: 0 : ( (1.0 - t_x) * grid_sources.get_source(n_x, n_y+1) +
[ # # ]
375 [ # # ][ # # ]: 0 : t_x * grid_sources.get_source(n_x + 1, n_y+1) )
376 [ # # ]: 0 : -
377 [ # # ][ # # ]: 0 : ( (1.0 - t_x) * grid_sources.get_source(n_x, n_y) +
[ # # ]
378 [ # # ][ # # ]: 0 : t_x * grid_sources.get_source(n_x + 1, n_y) )
379 [ # # ]: 0 : );
380 : : /*ret_val[0] = (1.0 / delta_x) * ( (1.0 - t_y) * (
381 : : f[n_x + 1][n_y] - f[n_x][n_y] ) +
382 : : t_y * (
383 : : f[n_x + 1][n_y+1] - f[n_x][n_y+1] )
384 : : );
385 : : ret_val[1] = (1.0 / delta_x) * (
386 : : ( (1.0 - t_x) * f[n_x][n_y+1] +
387 : : t_x * f[n_x + 1][n_y+1] )
388 : : -
389 : : ( (1.0 - t_x) * f[n_x][n_y] +
390 : : t_x * f[n_x + 1][n_y] )
391 : : );*/
392 : :
393 : : // TODO: uncomment, add an option for linear interpolation
394 : : /*
395 : : double weight;
396 : : if (!use_raised_cosine) { // use sinc
397 : : for (int i=begin[0]; i<=end[0]; i++) {
398 : : for (int j=begin[1]; j<=end[1]; j++) {
399 : : for (int k=begin[2]; k<=end[2]; k++) {
400 : : weight = sinc(a - (double)i) * sinc(b - (double)j) * sinc(c - (double)k);
401 : : ret_val += weight * f[i][j][k];
402 : : }
403 : : }
404 : : }
405 : : }
406 : : else { // use raised cosine
407 : : for (int i=begin[0]; i<=end[0]; i++) {
408 : : for (int j=begin[1]; j<=end[1]; j++) {
409 : : for (int k=begin[2]; k<=end[2]; k++) {
410 : : // we know that abs(a-(double)i) < N+1, so we don't need to worry about
411 : : // going out of the array's bounds
412 : : weight = lut->get(a-(double)i) *
413 : : lut->get(b-(double)j) *
414 : : lut->get(c-(double)k);
415 : : ret_val += weight * f[i][j][k];
416 : : }
417 : : }
418 : : }
419 : : }*/
420 : :
421 : 0 : return ret_val;
422 : : }
423 : :
424 : : // overloaded version for vector-valued functions of 3 variables
425 : 0 : Eigen::RowVectorXd evaluate_interpolation(double x[3], double x_min[3], double x_max[3],
426 : : double delta_x,
427 : : vector<vector<vector<Eigen::RowVectorXd>>> &f, int N,
428 : : bool /*use_raised_cosine*/, double beta) {
429 : : // determine nearest sample point index to the left of x
430 : : int idx[3], begin[3], end[3];
431 : :
432 : : // number of entries in f array across each dimension
433 : : int length[3];
434 : 0 : length[0] = static_cast<int>(f.size());
435 : 0 : length[1] = static_cast<int>(f[0].size());
436 : 0 : length[2] = static_cast<int>(f[0][0].size());
437 [ # # ]: 0 : for (int i=0; i<3; i++) {
438 : 0 : idx[i] = (int) floor((x[i] - x_min[i])/delta_x);
439 : :
440 : : // convolution will include N sample points on each side of x
441 : 0 : begin[i] = idx[i] - N + 1;
442 : 0 : end[i] = idx[i] + N;
443 : :
444 : : // ensure sum is within the range of the sample function
445 [ # # ]: 0 : if (begin[i] < 0) {
446 [ # # ]: 0 : std::cout << "WARNING: Sum in evaluate_interpolation is being " <<
447 [ # # ][ # # ]: 0 : "pruned due to going outside of range (x[" << i << "] = " <<
[ # # ]
448 [ # # ][ # # ]: 0 : x[i] << ", x_min[" << i << "] = " << x_min[i] << ", x_max[" << i << "] = " << x_max[i] << ")" << std::endl;
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
449 : 0 : begin[i] = 0;
450 : : }
451 [ # # ]: 0 : if (end[i] >= static_cast<int>(f.size())) {
452 [ # # ]: 0 : std::cout << "WARNING: Sum in evaluate_interpolation is being " <<
453 [ # # ][ # # ]: 0 : "pruned due to going outside of range (x[" << i << "] = " <<
[ # # ]
454 [ # # ][ # # ]: 0 : x[i] << ", x_min[" << i << "] = " << x_min[i] << ", x_max[" << i << "] = " << x_max[i] << ")" << std::endl;
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
455 : 0 : end[i] = length[i] - 1;
456 : : }
457 : : }
458 : :
459 : : // initialize return value
460 [ # # ][ # # ]: 0 : Eigen::RowVectorXd ret_val = Eigen::RowVectorXd::Zero(f[0][0][0].size());
[ # # ]
461 : : // pre-computing for efficiency
462 : 0 : double a = (x[0] - x_min[0])/delta_x;
463 : 0 : double b = (x[1] - x_min[1])/delta_x;
464 : 0 : double c = (x[2] - x_min[2])/delta_x;
465 : :
466 [ # # ]: 0 : RaisedCosLUT *lut = RaisedCosLUT::instance();
467 : : // check if initialized already; if not, initialize with default settings
468 [ # # ]: 0 : if (!lut->isInitialized()) {
469 [ # # ]: 0 : lut->initialize((double)(N+1), 100.0, beta); // range of [-(N+1),N+1], resolution of 100
470 : : }
471 : :
472 : : // TODO: using only linear interpolation for now
473 : :
474 : : // parameters for linear interpolation
475 : : double t_x, t_y, t_z;
476 : : double n_x, n_y, n_z;
477 : 0 : n_x = floor(a);
478 : 0 : n_y = floor(b);
479 : 0 : n_z = floor(c);
480 : 0 : t_x = a - n_x;
481 : 0 : t_y = b - n_y;
482 : 0 : t_z = c - n_z;
483 : :
484 [ # # ][ # # ]: 0 : cout << "Evaluating linear interpolation..." << endl;
485 [ # # ]: 0 : ret_val = (1.0 - t_z) * (
486 [ # # ]: 0 : (1.0 - t_y) * (
487 [ # # ][ # # ]: 0 : (1.0 - t_x) * f[(int)n_x][(int)n_y][(int)n_z] +
488 [ # # ][ # # ]: 0 : t_x * f[(int)n_x + 1][(int)n_y][(int)n_z] ) +
489 [ # # ]: 0 : t_y * (
490 [ # # ][ # # ]: 0 : (1.0 - t_x) * f[(int)n_x][(int)n_y+1][(int)n_z] +
491 [ # # ][ # # ]: 0 : t_x * f[(int)n_x + 1][(int)n_y+1][(int)n_z] ) ) +
492 [ # # ]: 0 : t_z * (
493 [ # # ]: 0 : (1.0 - t_y) * (
494 [ # # ][ # # ]: 0 : (1.0 - t_x) * f[(int)n_x][(int)n_y][(int)n_z+1] +
495 [ # # ][ # # ]: 0 : t_x * f[(int)n_x + 1][(int)n_y][(int)n_z+1] ) +
496 [ # # ]: 0 : t_y * (
497 [ # # ][ # # ]: 0 : (1.0 - t_x) * f[(int)n_x][(int)n_y+1][(int)n_z+1] +
498 [ # # ][ # # ]: 0 : t_x * f[(int)n_x + 1][(int)n_y+1][(int)n_z+1] ) );
499 : :
500 : : // TODO: uncomment, add an option for linear interpolation
501 : : /*
502 : : double weight;
503 : : if (!use_raised_cosine) { // use sinc
504 : : for (int i=begin[0]; i<=end[0]; i++) {
505 : : for (int j=begin[1]; j<=end[1]; j++) {
506 : : for (int k=begin[2]; k<=end[2]; k++) {
507 : : weight = sinc(a - (double)i) * sinc(b - (double)j) * sinc(c - (double)k);
508 : : ret_val += weight * f[i][j][k];
509 : : }
510 : : }
511 : : }
512 : : }
513 : : else { // use raised cosine
514 : : for (int i=begin[0]; i<=end[0]; i++) {
515 : : for (int j=begin[1]; j<=end[1]; j++) {
516 : : for (int k=begin[2]; k<=end[2]; k++) {
517 : : // we know that abs(a-(double)i) < N+1, so we don't need to worry about
518 : : // going out of the array's bounds
519 : : weight = lut->get(a-(double)i) *
520 : : lut->get(b-(double)j) *
521 : : lut->get(c-(double)k);
522 : : ret_val += weight * f[i][j][k];
523 : : }
524 : : }
525 : : }
526 : : }*/
527 : :
528 : 0 : return ret_val;
529 : : }
530 : :
531 : 0 : Eigen::RowVectorXd evaluate_interpolation(const Eigen::RowVectorXd & pt,
532 : : const Eigen::RowVectorXd & x_min,
533 : : double delta_x,
534 : : source_cache& grid_sources)
535 : : {
536 : : // initialize return value
537 [ # # ][ # # ]: 0 : Eigen::RowVectorXd ret_val = Eigen::RowVectorXd::Zero(grid_sources.get_num_tetras());
538 : : // pre-computing for efficiency
539 [ # # ][ # # ]: 0 : double a = (pt[0] - x_min[0])/delta_x;
540 [ # # ][ # # ]: 0 : double b = (pt[1] - x_min[1])/delta_x;
541 [ # # ][ # # ]: 0 : double c = (pt[2] - x_min[2])/delta_x;
542 : :
543 : : // parameters for linear interpolation
544 : : double t_x, t_y, t_z;
545 : : int n_x, n_y, n_z;
546 : 0 : n_x = static_cast<int>(floor(a));
547 : 0 : n_y = static_cast<int>(floor(b));
548 : 0 : n_z = static_cast<int>(floor(c));
549 : 0 : t_x = a - (double)n_x;
550 : 0 : t_y = b - (double)n_y;
551 : 0 : t_z = c - (double)n_z;
552 : :
553 [ # # ][ # # ]: 0 : cout << "Evaluating linear interpolation..." << endl;
554 [ # # ]: 0 : ret_val = (1.0 - t_z) * (
555 [ # # ]: 0 : (1.0 - t_y) * (
556 [ # # ][ # # ]: 0 : (1.0 - t_x) * grid_sources.get_source(n_x, n_y, n_z) +
[ # # ]
557 [ # # ][ # # ]: 0 : t_x * grid_sources.get_source(n_x + 1, n_y, n_z) ) +
[ # # ]
558 [ # # ]: 0 : t_y * (
559 [ # # ][ # # ]: 0 : (1.0 - t_x) * grid_sources.get_source(n_x, n_y+1, n_z) +
[ # # ]
560 [ # # ][ # # ]: 0 : t_x * grid_sources.get_source(n_x + 1, n_y+1, n_z) ) ) +
[ # # ]
561 [ # # ]: 0 : t_z * (
562 [ # # ]: 0 : (1.0 - t_y) * (
563 [ # # ][ # # ]: 0 : (1.0 - t_x) * grid_sources.get_source(n_x, n_y, n_z+1) +
[ # # ]
564 [ # # ][ # # ]: 0 : t_x * grid_sources.get_source(n_x + 1, n_y, n_z+1) ) +
[ # # ]
565 [ # # ]: 0 : t_y * (
566 [ # # ][ # # ]: 0 : (1.0 - t_x) * grid_sources.get_source(n_x, n_y+1, n_z+1) +
[ # # ]
567 [ # # ][ # # ]: 0 : t_x * grid_sources.get_source(n_x + 1, n_y+1, n_z+1) ) );
[ # # ]
568 : :
569 : 0 : return ret_val;
570 : : }
571 : :
572 : 0 : void sample_function(std::function<double(double)> &f,
573 : : double x_min,
574 : : double x_max,
575 : : double step,
576 : : vector<double> & func_values) {
577 : : // allocate memory
578 : 0 : int n_samples = (int)floor((x_max - x_min)/step);
579 : 0 : func_values.resize(n_samples);
580 : :
581 : : // sample function
582 [ # # ]: 0 : for (int i = 0; i<n_samples; i++) {
583 : 0 : func_values[i] = f(x_min + step*i);
584 : : }
585 : 0 : }
586 : :
587 : : // overloaded 3d version
588 : 0 : void sample_function(std::function<double(double, double, double)> &f,
589 : : double x_min[3],
590 : : double x_max[3],
591 : : double step,
592 : : vector<vector<vector<double>>> & func_values) {
593 : : // allocate memory
594 : : int n_samples[3];
595 [ # # ]: 0 : for (int i=0; i<3; i++) {
596 : 0 : n_samples[i] = (int)floor((x_max[i] - x_min[i])/step);
597 : : }
598 [ # # ][ # # ]: 0 : func_values.resize(n_samples[0], vector<vector<double>>(n_samples[1],
599 [ # # ]: 0 : vector<double>(n_samples[2])));
600 : : // sample function
601 [ # # ]: 0 : for (int i = 0; i<n_samples[0]; i++) {
602 [ # # ]: 0 : for (int j = 0; j<n_samples[1]; j++) {
603 [ # # ]: 0 : for (int k = 0; k<n_samples[2]; k++) {
604 [ # # ]: 0 : func_values[i][j][k] = f(x_min[0] + step*i, x_min[1] + step*j, x_min[2] + step*k);
605 : : }
606 : : }
607 : : }
608 : 0 : }
609 : :
610 : 0 : vector<float> interpolate_line_source(const SP_Line& line, const vector<double> & x_min,
611 : : source_cache & grid_sources, unsigned num_pts, double delta_x)
612 : : {
613 : : // first define the function that would compute the fluence from a point source
614 : 0 : std::function<Eigen::RowVectorXd(Eigen::RowVectorXd)> fluence_pt = [&](Eigen::RowVectorXd t)->Eigen::RowVectorXd
615 : : {
616 : : // 1- define delta_x TODO
617 [ # # ]: 0 : Eigen::RowVectorXd xmin(3);
618 [ # # ]: 0 : for (unsigned i = 0; i < 3; i++) {
619 [ # # ]: 0 : xmin(i) = x_min[i];
620 : : }
621 : : // 2- Call evaluate interpolation for point sources that make use of grid_sources
622 [ # # ]: 0 : return SP_IE::evaluate_interpolation(t, xmin, delta_x, grid_sources);
623 [ # # ]: 0 : };
624 : :
625 [ # # ][ # # ]: 0 : Eigen::RowVectorXd l(3), r(3);
626 [ # # ][ # # ]: 0 : l(0) = static_cast<double>(line.end0.get_xcoord());
627 [ # # ][ # # ]: 0 : l(1) = static_cast<double>(line.end0.get_ycoord());
628 [ # # ][ # # ]: 0 : l(2) = static_cast<double>(line.end0.get_zcoord());
629 [ # # ][ # # ]: 0 : r(0) = static_cast<double>(line.end1.get_xcoord());
630 [ # # ][ # # ]: 0 : r(1) = static_cast<double>(line.end1.get_ycoord());
631 [ # # ][ # # ]: 0 : r(2) = static_cast<double>(line.end1.get_zcoord());
632 [ # # ][ # # ]: 0 : Eigen::RowVectorXd fluence_vals = integrate_midpoint(fluence_pt, l, r, num_pts);
[ # # ][ # # ]
633 : :
634 : : // TODO: might be expensive
635 : 0 : vector<float> fluence;
636 [ # # ][ # # ]: 0 : for (unsigned i = 0; i < fluence_vals.size(); i++) {
637 [ # # ][ # # ]: 0 : fluence.push_back(static_cast<float>(fluence_vals(i)));
638 : : }
639 : :
640 : 0 : return fluence;
641 : : }
642 : :
643 [ + - ][ + - ]: 44 : };
|