YetAnotherCoupler 3.2.0
Loading...
Searching...
No Matches
interp_method_nnn.c
Go to the documentation of this file.
1// Copyright (c) 2024 The YAC Authors
2//
3// SPDX-License-Identifier: BSD-3-Clause
4
5#ifdef HAVE_CONFIG_H
6// Get the definition of the 'restrict' keyword.
7#include "config.h"
8#endif
9
10#include <string.h>
11
13#include "interp_method_nnn.h"
15
16static size_t do_search_nnn(struct interp_method * method,
17 struct yac_interp_grid * interp_grid,
18 size_t * tgt_points, size_t count,
19 struct yac_interp_weights * weights);
20static size_t do_search_1nn(struct interp_method * method,
21 struct yac_interp_grid * interp_grid,
22 size_t * tgt_points, size_t count,
23 struct yac_interp_weights * weights);
24static void delete_nnn(struct interp_method * method);
25
26static struct interp_method_vtable
33
35
38 double[3], yac_coordinate_pointer, size_t const, double *, double const);
39 size_t n;
40 double scale;
41};
42
44
47 struct yac_interp_weights * weights, struct remote_points * tgts,
48 struct remote_point * srcs);
49};
50
52 double tgt_coord[3], yac_coordinate_pointer src_coords, size_t const n,
53 double * weights, double const dummy) {
54
55 UNUSED(tgt_coord);
56 UNUSED(src_coords);
57 UNUSED(dummy);
58
59 double weight = 1.0 / (double)n;
60 for (size_t i = 0; i < n; ++i) weights[i] = weight;
61}
62
64 double tgt_coord[3], yac_coordinate_pointer src_coords, size_t const n,
65 double * weights, double const dummy) {
66
67 UNUSED(dummy);
68
69 for (size_t i = 0; i < n; ++i) {
70
71 double distance = get_vector_angle(tgt_coord, (double*)(src_coords[i]));
72
73 // if the target and source point are nearly identical
74 if (distance < yac_angle_tol) {
75 for (size_t j = 0; j < n; ++j) weights[j] = 0.0;
76 weights[i] = 1.0;
77 return;
78 }
79
80 weights[i] = 1.0 / distance;
81 }
82
83 // compute scaling factor for the weights
84 double inv_distance_sum = 0.0;
85 for (size_t i = 0; i < n; ++i) inv_distance_sum += weights[i];
86 double scale = 1.0 / inv_distance_sum;
87
88 for (size_t i = 0; i < n; ++i) weights[i] *= scale;
89}
90
92 double tgt_coord[3], yac_coordinate_pointer src_coords, size_t const n,
93 double * weights, double const gauss_scale) {
94
95 for (size_t i = 0; i < n; ++i) {
96
97 double distance = get_vector_angle(tgt_coord, src_coords[i]);
98
99 // if the target and source point are nearly identical
100 if (distance < yac_angle_tol) {
101 for (size_t j = 0; j < n; ++j) weights[j] = 0.0;
102 weights[i] = 1.0;
103 return;
104 }
105 weights[i] = distance * distance;
106 }
107
108 // a) compute sum of source point distances
109 double src_distances_sum = 0.0;
110 double src_distances_count = 0.5 * (double)(n * n - n);
111 for (size_t i = 0; i < n-1; ++i)
112 for (size_t j = i + 1; j < n; ++j)
113 src_distances_sum += get_vector_angle(src_coords[i], src_coords[j]);
114
115 // b) C = -1 / (c * d_mean^2)
116 double src_distances_mean = src_distances_sum / src_distances_count;
117 double scale = -1.0 / (gauss_scale * src_distances_mean * src_distances_mean);
118
119 // c) calculate weights
120 // w_i = e^(-d_i^2/(c*s^2))
121 // w_i = e^(C * d_i^2)
122 double weights_sum = 0.0;
123 for (size_t i = 0; i < n; ++i) {
124 weights[i] = exp(scale * weights[i]);
125 weights_sum += weights[i];
126 }
127
128 // If the sum of the weights is very low, which can happen in case
129 // the target point is very far away from the group source points.
130 if (fabs(weights_sum) < 1e-9) {
131
132 // Due to limited accuracy the exact contribution of each source
133 // point cannot be computed. Therefore, the normalisation would
134 // generate NaN's. Hence we fall back to inverse distance weighted
135 // averge for this target point.
137 tgt_coord, src_coords, n, weights, gauss_scale);
138 return;
139 }
140
141 // d) scale weights such that SUM(w_i) == 1
142 for (size_t i = 0; i < n; ++i) weights[i] /= weights_sum;
143}
144
145static void inverse(double * A, size_t n) {
146
147// LAPACKE_dsytrf_work and LAPACKE_dsytri might not be available (see yac_lapack_interface.h).
148#ifdef YAC_LAPACK_NO_DSYTR
149 lapack_int ipiv[n+1];
150 double work[n*n];
151
152 for (size_t i = 0; i < n+1; ++i) ipiv[i] = 0;
153 for (size_t i = 0; i < n*n; ++i) work[i] = 0;
154
156 !LAPACKE_dgetrf(
157 LAPACK_COL_MAJOR, (lapack_int) n, (lapack_int) n,
158 A, (lapack_int) n, ipiv), "internal ERROR: dgetrf")
159
161 !LAPACKE_dgetri_work(
162 LAPACK_COL_MAJOR, (lapack_int) n, A, (lapack_int) n,
163 ipiv, work, (lapack_int) (n * n)), "internal ERROR: dgetri")
164#else
165 lapack_int ipiv[n];
166 double work[n];
167
169 !LAPACKE_dsytrf_work(
170 LAPACK_COL_MAJOR, 'L', (lapack_int) n , A,
171 (lapack_int) n, ipiv, work, (lapack_int) n), "internal ERROR: dsytrf")
172
174 !LAPACKE_dsytri_work(
175 LAPACK_COL_MAJOR, 'L', (lapack_int) n , A,
176 (lapack_int) n, ipiv, work), "internal ERROR: dsytri")
177
178 for (size_t i = 0; i < n; ++i)
179 for (size_t j = i + 1; j < n; ++j)
180 A[j*n+i] = A[i*n+j];
181#endif
182}
183
185 double tgt_coord[3], yac_coordinate_pointer src_coords, size_t const n,
186 double * weights, double const rbf_scale) {
187
188 double A[n][n];
189 double a[n];
190
191 double sum_d = 0.0, scale_d = 1.0;
192
193 // compute distance matrix for all found source points
194 for (size_t i = 0; i < n; ++i) A[i][i] = 0.0;
195 for (size_t i = 0; i < n - 1; ++i) {
196 for (size_t j = i + 1; j < n; ++j) {
197 double d = get_vector_angle(src_coords[i], src_coords[j]);
198 A[i][j] = d;
199 A[j][i] = d;
200 sum_d += d;
201 }
202 }
203
204 // compute and apply scale factor for distance matrix
205 if (sum_d > 0.0)
206 scale_d = ((double)((n - 1) * n)) / (2.0 * sum_d);
207 scale_d /= rbf_scale;
208
209 double sq_scale_d = scale_d * scale_d;
210
211 // compute matrix A[n][n]
212 // with A = rbf(A)
213 // rbf(a) => Radial basis function
214 for (size_t i = 0; i < n; ++i) {
215 for (size_t j = 0; j < n; ++j) {
216 double d = A[i][j];
217 // gauß rbf kernel
218 A[i][j] = exp(-1.0 * d * d * sq_scale_d);
219 }
220 }
221
222 // compute inverse of A
223 inverse(&A[0][0], n);
224
225 // compute a[NUM_NEIGH]
226 // with a_i = rbf(d(x_i, y))
227 // x => vector containing the coordinates of the
228 // n nearest neighbours of the current
229 // target point
230 // y => coordinates of target point
231 // d(a, b) => great circle distance between point a and b
232 // rbf(a) => Radial basis function
233 for (size_t i = 0; i < n; ++i) {
234 double d = get_vector_angle(tgt_coord, src_coords[i]);
235 // gauß rbf kernel
236 a[i] = exp(-1.0 * d * d * sq_scale_d);
237 }
238
239 // compute weights
240 // with w_i = SUM(A_inv[i][j]*a[j]) // j = 0..n-1
241 for (size_t i = 0; i < n; ++i) {
242 weights[i] = 0.0;
243 for (size_t j = 0; j < n; ++j) weights[i] += A[i][j] * a[j];
244 }
245}
246
247static size_t do_search_nnn (struct interp_method * method,
248 struct yac_interp_grid * interp_grid,
249 size_t * tgt_points, size_t count,
250 struct yac_interp_weights * weights) {
251
252 struct interp_method_nnn * method_nnn = (struct interp_method_nnn *)method;
253
255 yac_interp_grid_get_num_src_fields(interp_grid) == 1,
256 "ERROR(do_search_nnn): invalid number of source fields")
257
258 // get coordinates of target points
259 yac_coordinate_pointer tgt_coords = xmalloc(count * sizeof(*tgt_coords));
261 interp_grid, tgt_points, count, tgt_coords);
262
263 size_t n = method_nnn->n;
264
265 size_t * size_t_buffer = xmalloc((2 * n + 1) * count * sizeof(*size_t_buffer));
266 size_t * src_points = size_t_buffer;
267 size_t * src_points_reorder = size_t_buffer + n * count;
268 size_t * reorder_idx = size_t_buffer + 2 * n * count;
269 int * successful_flag = xmalloc(count * sizeof(*successful_flag));
270
271 // search for the n nearest source points
273 interp_grid, tgt_coords, count, n, src_points);
274
275 size_t result_count = 0;
276
277 // check for which cells the search was successful
278 for (size_t i = 0, k = 0; i < count; ++i) {
279 int flag = 1;
280 for (size_t j = 0; j < n; ++j, ++k) flag &= src_points[k] != SIZE_MAX;
281 result_count += (size_t)(successful_flag[i] = flag);
282 reorder_idx[i] = i;
283 }
284
285 // sort target points, for which the search was successful to the beginning of
286 // the array
287 for (size_t i = 0; i < count; ++i) reorder_idx[i] = i;
288 yac_quicksort_index_int_size_t(successful_flag, count, reorder_idx);
289 free(successful_flag);
290
291 void (*compute_weights)(
292 double[3], yac_coordinate_pointer, size_t, double *, double const) =
293 method_nnn->compute_weights;
294 double scale = method_nnn->scale;
295 double * w = xmalloc(n * result_count * sizeof(*w));
296 size_t * num_weights_per_tgt =
297 xmalloc(result_count * sizeof(*num_weights_per_tgt));
298 yac_coordinate_pointer src_coord_buffer = xmalloc(n * sizeof(*src_coord_buffer));
299 yac_const_coordinate_pointer src_field_coordinates =
301
302 // compute the weights
303 for (size_t i = 0; i < result_count; ++i) {
304 size_t curr_reorder_idx = reorder_idx[i];
305 for (size_t j = 0; j < n; ++j) {
306 size_t curr_src_point = src_points[curr_reorder_idx * n + j];
307 memcpy(src_coord_buffer[j], src_field_coordinates[curr_src_point],
308 3 * sizeof(double));
309 src_points_reorder[i * n + j] = curr_src_point;
310 }
311
313 tgt_coords[curr_reorder_idx], src_coord_buffer, n, w + i * n, scale);
314 num_weights_per_tgt[i] = n;
315 }
316
317 // move the non-interpolated target points to the end
318 for (size_t i = 0; i < count; ++i) src_points[reorder_idx[i]] = i;
319 yac_quicksort_index_size_t_size_t(src_points, count, tgt_points);
320
321 free(src_coord_buffer);
322 free(tgt_coords);
323
324 struct remote_point * srcs =
326 interp_grid, 0, src_points_reorder, result_count * n);
327
328 struct remote_points tgts = {
329 .data =
331 interp_grid, tgt_points, result_count),
332 .count = result_count};
333
334 // store weights
336 weights, &tgts, num_weights_per_tgt, srcs, w);
337
338 free(size_t_buffer);
339 free(tgts.data);
340 free(srcs);
341 free(w);
342 free(num_weights_per_tgt);
343
344 return result_count;
345}
346
347static size_t do_search_1nn(struct interp_method * method,
348 struct yac_interp_grid * interp_grid,
349 size_t * tgt_points, size_t count,
350 struct yac_interp_weights * weights) {
351
353 yac_interp_grid_get_num_src_fields(interp_grid) == 1,
354 "ERROR(do_search_1nn): invalid number of source fields")
355
356 // get coordinates of target points
357 yac_coordinate_pointer tgt_coords = xmalloc(count * sizeof(*tgt_coords));
359 interp_grid, tgt_points, count, tgt_coords);
360
361 size_t * src_points = xmalloc(count * sizeof(*src_points));
362
363 // search for the nearest source points
365 interp_grid, tgt_coords, count, 1, src_points);
366
367 // move the non-interpolated target points to the end
368 // (source points == SIZE_MAX if no source point was found)
369 yac_quicksort_index_size_t_size_t(src_points, count, tgt_points);
370 free(tgt_coords);
371
372 // count the number of target points for which we found a source point
373 size_t result_count = 0;
374 for (; (result_count < count) && (src_points[result_count] != SIZE_MAX);
375 ++result_count);
376
377 struct remote_point * srcs =
379 interp_grid, 0, src_points, result_count);
380
381 struct remote_points tgts = {
382 .data =
384 interp_grid, tgt_points, result_count),
385 .count = result_count};
386
387 // store weights
388 ((struct interp_method_1nn *)method)->interp_weights_add(
389 weights, &tgts, srcs);
390
391 free(src_points);
392 free(tgts.data);
393 free(srcs);
394
395 return result_count;
396}
397
399 struct yac_interp_weights * weights, struct remote_points * tgts,
400 struct remote_point * srcs) {
401
402 size_t * num_src_per_tgt = xmalloc(tgts->count * sizeof(*num_src_per_tgt));
403 double * w = xmalloc(tgts->count * sizeof(*w));
404
405 for (size_t i = 0; i < tgts->count; ++i) {
406 num_src_per_tgt[i] = 1;
407 w[i] = 0.0;
408 }
409
410 yac_interp_weights_add_wsum(weights, tgts, num_src_per_tgt, srcs, w);
411
412 free(w);
413 free(num_src_per_tgt);
414}
415
417
418 if (config.type == YAC_INTERP_NNN_ZERO) {
419
421 config.n == 1,
422 "ERROR(yac_interp_method_nnn_new): n != 1 does not make sense for zero")
423
424 struct interp_method_1nn * method = xmalloc(1 * sizeof(*method));
427 return (struct interp_method*)method;
428 }
429
430 if (config.n == 1) {
431
432 struct interp_method_1nn * method = xmalloc(1 * sizeof(*method));
435 return (struct interp_method*)method;
436 }
437
438 struct interp_method_nnn * method = xmalloc(1 * sizeof(*method));
439
441
442 method->n = config.n;
443
445 (config.type == YAC_INTERP_NNN_AVG) ||
446 (config.type == YAC_INTERP_NNN_DIST) ||
447 (config.type == YAC_INTERP_NNN_GAUSS) ||
448 (config.type == YAC_INTERP_NNN_RBF),
449 "ERROR(yac_interp_method_nnn_new): invalid NNN type")
450
451 switch (config.type) {
452 case(YAC_INTERP_NNN_AVG):
454 config.n >= 1,
455 "ERROR(yac_interp_method_nnn_new): n needs to be at least 1 for "
456 "distance weighted average")
457 method->scale = 0.0;
459 break;
462 config.n >= 2,
463 "ERROR(yac_interp_method_nnn_new): n needs to be at least 2 for "
464 "distance weighted average")
465 method->scale = 0.0;
467 break;
470 config.n >= 2,
471 "ERROR(yac_interp_method_nnn_new): n needs to be at least 2 for gauss")
472 method->scale = config.data.gauss_scale;
474 break;
475 default:
476 case(YAC_INTERP_NNN_RBF):
478 config.n >= 2,
479 "ERROR(yac_interp_method_nnn_new): n needs to be at least 2 for RBF")
480 method->scale = config.data.rbf_scale;
482 break;
483 };
484
485 return (struct interp_method*)method;
486}
487
488static void delete_nnn(struct interp_method * method) {
489 free(method);
490}
#define UNUSED(x)
Definition core.h:71
#define yac_angle_tol
Definition geometry.h:34
static double get_vector_angle(double const a[3], double const b[3])
Definition geometry.h:364
size_t yac_interp_grid_get_num_src_fields(struct yac_interp_grid *interp_grid)
struct remote_point * yac_interp_grid_get_tgt_remote_points(struct yac_interp_grid *interp_grid, size_t *tgt_points, size_t count)
struct remote_point * yac_interp_grid_get_src_remote_points(struct yac_interp_grid *interp_grid, size_t src_field_idx, size_t *src_points, size_t count)
yac_const_coordinate_pointer yac_interp_grid_get_src_field_coords(struct yac_interp_grid *interp_grid, size_t src_field_idx)
void yac_interp_grid_do_nnn_search_src(struct yac_interp_grid *interp_grid, yac_coordinate_pointer search_coords, size_t count, size_t n, size_t *src_points)
void yac_interp_grid_get_tgt_coordinates(struct yac_interp_grid *interp_grid, size_t *tgt_points, size_t count, yac_coordinate_pointer tgt_coordinates)
static void compute_weights(struct tgt_point_search_data *tgt_point_data, size_t num_tgt_points, struct edge_interp_data *edge_data, size_t num_edges, struct triangle_interp_data *triangle_data, size_t num_triangles, struct weight_vector_data **weights, size_t **num_weights_per_tgt, size_t *total_num_weights)
static void compute_weights_gauss(double tgt_coord[3], yac_coordinate_pointer src_coords, size_t const n, double *weights, double const gauss_scale)
static size_t do_search_nnn(struct interp_method *method, struct yac_interp_grid *interp_grid, size_t *tgt_points, size_t count, struct yac_interp_weights *weights)
static void compute_weights_rbf(double tgt_coord[3], yac_coordinate_pointer src_coords, size_t const n, double *weights, double const rbf_scale)
static void delete_nnn(struct interp_method *method)
static struct interp_method_vtable interp_method_1nn_vtable
static void compute_weights_dist(double tgt_coord[3], yac_coordinate_pointer src_coords, size_t const n, double *weights, double const dummy)
static void interp_weights_add_zero(struct yac_interp_weights *weights, struct remote_points *tgts, struct remote_point *srcs)
static void inverse(double *A, size_t n)
static size_t do_search_1nn(struct interp_method *method, struct yac_interp_grid *interp_grid, size_t *tgt_points, size_t count, struct yac_interp_weights *weights)
static struct interp_method_vtable interp_method_nnn_vtable
struct interp_method * yac_interp_method_nnn_new(struct yac_nnn_config config)
static void compute_weights_avg(double tgt_coord[3], yac_coordinate_pointer src_coords, size_t const n, double *weights, double const dummy)
@ YAC_INTERP_NNN_GAUSS
distance with Gauss weights of n source points
@ YAC_INTERP_NNN_RBF
radial basis functions
@ YAC_INTERP_NNN_AVG
average of n source points
@ YAC_INTERP_NNN_DIST
distance weighted average of n source points
@ YAC_INTERP_NNN_ZERO
all weights are set to zero
void yac_interp_weights_add_wsum(struct yac_interp_weights *weights, struct remote_points *tgts, size_t *num_src_per_tgt, struct remote_point *srcs, double *w)
void yac_interp_weights_add_direct(struct yac_interp_weights *weights, struct remote_points *tgts, struct remote_point *srcs)
#define xmalloc(size)
Definition ppm_xfuncs.h:66
void(* interp_weights_add)(struct yac_interp_weights *weights, struct remote_points *tgts, struct remote_point *srcs)
struct interp_method_vtable * vtable
void(* compute_weights)(double[3], yac_coordinate_pointer, size_t const, double *, double const)
struct interp_method_vtable * vtable
size_t(* do_search)(struct interp_method *method, struct yac_interp_grid *grid, size_t *tgt_points, size_t count, struct yac_interp_weights *weights)
struct remote_point * data
union yac_nnn_config::@22 data
enum yac_interp_nnn_weight_type type
void yac_quicksort_index_int_size_t(int *a, size_t n, size_t *idx)
void yac_quicksort_index_size_t_size_t(size_t *a, size_t n, size_t *idx)
#define YAC_ASSERT(exp, msg)
Definition yac_assert.h:16
double const (* yac_const_coordinate_pointer)[3]
Definition yac_types.h:20
double(* yac_coordinate_pointer)[3]
Definition yac_types.h:19