test_interp_method_dnn_parallel.c

Go to the documentation of this file.

// Copyright (c) 2026 The YAC Authors
//
// SPDX-License-Identifier: BSD-3-Clause
 
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
 
#include "tests.h"
#include "interpolation/methods/interp_method.h"
#include "interpolation/methods/interp_method_dnn.h"
#include "interpolation/methods/interp_method_fixed.h"
#include "dist_grid_utils.h"
#include "yac_mpi.h"
#include "geometry.h"
#include "yac_lapack_interface.h"
 
#include <mpi.h>
#include <yaxt.h>
 
static void utest_compute_reference_result(
  struct yac_interp_method_dnn_config const dnn_config,
  double const * src_coordinates_x, double const * src_coordinates_y,
  int const * src_global_mask, size_t const src_size_x, size_t const src_size_y,
  double const * tgt_coordinates_x, double const * tgt_coordinates_y,
  size_t const tgt_size_x, size_t const tgt_size_y,
  double * ref_tgt_results);
 
static double const fixed_value = -1.0;
 
int main (void) {
 
  MPI_Init(NULL, NULL);
 
  xt_initialize(MPI_COMM_WORLD);
 
  int comm_rank, comm_size;
  MPI_Comm_rank(MPI_COMM_WORLD, &comm_rank);
  MPI_Comm_size(MPI_COMM_WORLD, &comm_size);
  MPI_Barrier(MPI_COMM_WORLD);
 
  if (comm_size != 5) {
    PUT_ERR("ERROR: wrong number of processes");
    xt_finalize();
    MPI_Finalize();
    return TEST_EXIT_CODE;
  }
 
  enum yac_interp_weights_reorder_type reorder_type[2] =
    {YAC_MAPPING_ON_SRC, YAC_MAPPING_ON_TGT};
  enum {NUM_REORDER_TYPES = sizeof(reorder_type)/sizeof(reorder_type[0])};
 
  { // Test DNN with fixed search distance, distance weights, and fixed value fallback
 
    // corner and cell ids for a 7 x 7 grid (x == target point position)
    // uses 4 source processes (ranks 0 to 3) and 1 target process (rank 4)
 
    // 56-----57-----58-----59-----60-----61-----62-----63
    //  |     x|     x|     x|     x|     x|     x|     x|
    //  |  42  |  43  |  44  |  45  |  46  |  47  |  48  |
    //  |      |      |      |      |      |      |      |
    // 48-----49-----50-----51-----52-----53-----54-----55
    //  |     x|     x|     x|     x|     x|     x|     x|
    //  |  35  |  36  |  37  |  38  |  39  |  40  |  41  |
    //  |      |      |      |      |      |      |      |
    // 40-----41-----42-----43-----44-----45-----46-----47
    //  |     x|     x|     x|     x|     x|     x|     x|
    //  |  28  |  29  |  30  |  31  |  32  |  33  |  34  |
    //  |      |      |      |      |      |      |      |
    // 32-----33-----34-----35-----36-----37-----38-----39
    //  |     x|     x|     x|     x|     x|     x|     x|
    //  |  21  |  22  |  23  |  24  |  25  |  26  |  27  |
    //  |      |      |      |      |      |      |      |
    // 24-----25-----26-----27-----28-----29-----30-----31
    //  |     x|     x|     x|     x|     x|     x|     x|
    //  |  14  |  15  |  16  |  17  |  18  |  19  |  20  |
    //  |      |      |      |      |      |      |      |
    // 16-----17-----18-----19-----20-----21-----22-----23
    //  |     x|     x|     x|     x|     x|     x|     x|
    //  |  07  |  08  |  09  |  10  |  11  |  12  |  13  |
    //  |      |      |      |      |      |      |      |
    // 08-----09-----10-----11-----12-----13-----14-----15
    //  |     x|     x|     x|     x|     x|     x|     x|
    //  |  00  |  01  |  02  |  03  |  04  |  05  |  06  |
    //  |      |      |      |      |      |      |      |
    // 00-----01-----02-----03-----04-----05-----06-----07
    //
    // the grid is distributed among the processes as follows:
    // (index == process)
    //
    // 3---3---3---3---3---3---3---3
    // | 3 | 3 | 3 | 3 | 3 | 3 | 3 |
    // 3---3---3---3---3---3---3---3
    // | 3 | 3 | 3 | 3 | 3 | 3 | 3 |
    // 3---3---3---1---2---2---3---3
    // | 1 | 1 | 1 | 2 | 2 | 2 | 2 |
    // 1---1---1---2---2---2---2---2
    // | 1 | 1 | 1 | 2 | 2 | 2 | 2 |
    // 1---1---1---1---2---2---2---2
    // | 1 | 1 | 1 | 2 | 2 | 2 | 2 |
    // 1---1---1---0---0---0---2---2
    // | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
    // 0---0---0---0---0---0---0---0
    // | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
    // 0---0---0---0---0---0---0---0
    //
    // the source mask looks as follows (# == masked point)
    //
    // +---+---+---+---+---+---+---#
    // |   |   |   |   |   |   |   |
    // +---+---+---+---+---+---#---+
    // |   |   |   |   |   |   |   |
    // +---+---+---+---+---#---+---+
    // |   |   |   |   |   |   |   |
    // +---+---+---+---#---+---+---+
    // |   |   |   |   |   |   |   |
    // #---#---#---#---+---+---+---+
    // |   |   |   |   |   |   |   |
    // #---#---#---#---+---+---+---+
    // |   |   |   |   |   |   |   |
    // #---#---#---#---+---+---+---+
    // |   |   |   |   |   |   |   |
    // #---#---#---#---+---+---+---+
 
    enum {NUM_CELLS_X = 7, NUM_CELLS_Y = 7, NUM_PROCS = 5, TGT_RANK = 4};
    char const * src_grid_name= "src_grid";
    char const * tgt_grid_name= "tgt_grid";
    double coordinates_x[NUM_CELLS_X+1] = {0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0};
    double coordinates_y[NUM_CELLS_Y+1] = {0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0};
    double cell_coordinates_x[NUM_CELLS_X] = {0.75,1.75,2.75,3.75,4.75,5.75,6.75};
    double cell_coordinates_y[NUM_CELLS_Y] = {0.75,1.75,2.75,3.75,4.75,5.75,6.75};
    size_t const num_cells[2] = {NUM_CELLS_X,NUM_CELLS_Y};
    size_t const local_start[NUM_PROCS][2] = {{0,0},{0,2},{3,2},{0,5}, {0,0}};
    size_t const local_count[NUM_PROCS][2] = {{7,2},{3,3},{4,3},{7,2}, {7,7}};
    int const with_halo = 0;
    int const src_global_corner_mask[NUM_CELLS_X+1][NUM_CELLS_Y+1] = {
      {0,0,0,0,1,1,1,1},
      {0,0,0,0,1,1,1,1},
      {0,0,0,0,1,1,1,1},
      {0,0,0,0,1,1,1,1},
      {1,1,1,1,0,1,1,1},
      {1,1,1,1,1,0,1,1},
      {1,1,1,1,1,1,0,1},
      {1,1,1,1,1,1,1,0}};
    for (size_t i = 0; i <= NUM_CELLS_X; ++i) coordinates_x[i] *= YAC_RAD;
    for (size_t i = 0; i <= NUM_CELLS_Y; ++i) coordinates_y[i] *= YAC_RAD;
    for (size_t i = 0; i < NUM_CELLS_X; ++i) cell_coordinates_x[i] *= YAC_RAD;
    for (size_t i = 0; i < NUM_CELLS_Y; ++i) cell_coordinates_y[i] *= YAC_RAD;
 
    int const is_tgt = comm_rank == TGT_RANK;
 
    struct yac_basic_grid_data grid_data =
      utest_generate_basic_grid_data_reg2d(
        coordinates_x, coordinates_y, num_cells,
        local_start[comm_rank], local_count[comm_rank], with_halo);
 
    struct yac_basic_grid * grid =
      yac_basic_grid_new(is_tgt?tgt_grid_name:src_grid_name, grid_data);
    struct yac_basic_grid * empty_grid =
      yac_basic_grid_empty_new(is_tgt?src_grid_name:tgt_grid_name);
 
    // define custom coordinates and masks in grid data
    if (is_tgt) {
      yac_coordinate_pointer tgt_cell_coordinates =
        xmalloc(NUM_CELLS_X * NUM_CELLS_Y * sizeof(*tgt_cell_coordinates));
      for (size_t i = 0, k = 0; i < NUM_CELLS_Y; ++i)
        for (size_t j = 0; j < NUM_CELLS_X; ++j, ++k)
          LLtoXYZ(
            cell_coordinates_x[j], cell_coordinates_y[i], tgt_cell_coordinates[k]);
      yac_basic_grid_add_coordinates_nocpy(
        grid, YAC_LOC_CELL, tgt_cell_coordinates);
    } else { // is src
      int * src_corner_mask =
        xmalloc(grid_data.num_vertices * sizeof(*src_corner_mask));
      for (size_t i = 0; i < grid_data.num_vertices; ++i)
        src_corner_mask[i] =
          ((int*)(&(src_global_corner_mask[0][0])))[grid_data.vertex_ids[i]];
      yac_basic_grid_add_mask_nocpy(
        grid, YAC_LOC_CORNER, src_corner_mask, NULL);
    }
 
    struct yac_dist_grid_pair * grid_pair =
      yac_dist_grid_pair_new(grid, empty_grid, MPI_COMM_WORLD);
 
    struct yac_interp_field src_fields[] =
      {{.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = 0}};
    size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
    struct yac_interp_field tgt_field =
      {.location = YAC_LOC_CELL, .coordinates_idx = 0, .masks_idx = SIZE_MAX};
 
    struct yac_interp_grid * interp_grid =
      yac_interp_grid_new(
        grid_pair, src_grid_name, tgt_grid_name,
        num_src_fields, src_fields, tgt_field);
 
    // Test different DNN configurations:
    struct yac_interp_method_dnn_config dnn_configs[] = {
      // Fixed search distance (1.3 degrees in radians), AVG weighting
      {.type = YAC_INTERP_DNN_WEIGHT_AVG,
       .search_distance =
         &(struct yac_interp_method_dnn_config_search_distance){
           .type = YAC_INTERP_DNN_SEARCH_DISTANCE_FIXED,
           .search_distance = {.fixed = 1.3 * YAC_RAD}},
       .n_min = YAC_INTERP_DNN_N_MIN_DEFAULT,
       .gauss_scale = YAC_INTERP_DNN_GAUSS_SCALE_DEFAULT,
       .rbf_scale = YAC_INTERP_DNN_RBF_SCALE_DEFAULT},
      // Fixed search distance, DIST weighting
      {.type = YAC_INTERP_DNN_WEIGHT_DIST,
       .search_distance =
         &(struct yac_interp_method_dnn_config_search_distance){
           .type = YAC_INTERP_DNN_SEARCH_DISTANCE_FIXED,
           .search_distance = {.fixed = 1.3 * YAC_RAD}},
       .n_min = YAC_INTERP_DNN_N_MIN_DEFAULT,
       .gauss_scale = YAC_INTERP_DNN_GAUSS_SCALE_DEFAULT,
       .rbf_scale = YAC_INTERP_DNN_RBF_SCALE_DEFAULT},
      // Fixed search distance, GAUSS weighting
      {.type = YAC_INTERP_DNN_WEIGHT_GAUSS,
       .search_distance =
         &(struct yac_interp_method_dnn_config_search_distance){
           .type = YAC_INTERP_DNN_SEARCH_DISTANCE_FIXED,
           .search_distance = {.fixed = 1.3 * YAC_RAD}},
       .n_min = YAC_INTERP_DNN_N_MIN_DEFAULT,
       .gauss_scale = YAC_INTERP_DNN_GAUSS_SCALE_DEFAULT,
       .rbf_scale = YAC_INTERP_DNN_RBF_SCALE_DEFAULT},
      // Smaller fixed search distance (0.4 degrees in radians)
      {.type = YAC_INTERP_DNN_WEIGHT_AVG,
       .search_distance =
         &(struct yac_interp_method_dnn_config_search_distance){
           .type = YAC_INTERP_DNN_SEARCH_DISTANCE_FIXED,
           .search_distance = {.fixed = 0.4 * YAC_RAD}},
       .n_min = YAC_INTERP_DNN_N_MIN_DEFAULT,
       .gauss_scale = YAC_INTERP_DNN_GAUSS_SCALE_DEFAULT,
       .rbf_scale = YAC_INTERP_DNN_RBF_SCALE_DEFAULT},
      // Very smaller fixed search distance (0.1 degrees in radians)
      {.type = YAC_INTERP_DNN_WEIGHT_AVG,
       .search_distance =
         &(struct yac_interp_method_dnn_config_search_distance){
           .type = YAC_INTERP_DNN_SEARCH_DISTANCE_FIXED,
           .search_distance = {.fixed = 0.4 * YAC_RAD}},
       .n_min = YAC_INTERP_DNN_N_MIN_DEFAULT,
       .gauss_scale = YAC_INTERP_DNN_GAUSS_SCALE_DEFAULT,
       .rbf_scale = YAC_INTERP_DNN_RBF_SCALE_DEFAULT},
      // Cell area search distance with scale 1.0
      {.type = YAC_INTERP_DNN_WEIGHT_DIST,
       .search_distance =
         &(struct yac_interp_method_dnn_config_search_distance){
           .type = YAC_INTERP_DNN_SEARCH_DISTANCE_CELL_AREA,
           .search_distance = {.scale = 1.0}},
       .n_min = YAC_INTERP_DNN_N_MIN_DEFAULT,
       .gauss_scale = YAC_INTERP_DNN_GAUSS_SCALE_DEFAULT,
       .rbf_scale = YAC_INTERP_DNN_RBF_SCALE_DEFAULT},
      // Cell area search distance with larger scale
      {.type = YAC_INTERP_DNN_WEIGHT_AVG,
       .search_distance =
         &(struct yac_interp_method_dnn_config_search_distance){
           .type = YAC_INTERP_DNN_SEARCH_DISTANCE_CELL_AREA,
           .search_distance = {.scale = 2.0}},
       .n_min = YAC_INTERP_DNN_N_MIN_DEFAULT,
       .gauss_scale = YAC_INTERP_DNN_GAUSS_SCALE_DEFAULT,
       .rbf_scale = YAC_INTERP_DNN_RBF_SCALE_DEFAULT},
      // RBF weighting with fixed search distance
      {.type = YAC_INTERP_DNN_WEIGHT_RBF,
       .search_distance =
         &(struct yac_interp_method_dnn_config_search_distance){
           .type = YAC_INTERP_DNN_SEARCH_DISTANCE_FIXED,
           .search_distance = {.fixed = 1.3 * YAC_RAD}},
       .n_min = YAC_INTERP_DNN_N_MIN_DEFAULT,
       .gauss_scale = YAC_INTERP_DNN_GAUSS_SCALE_DEFAULT,
       .rbf_scale = YAC_INTERP_DNN_RBF_SCALE_DEFAULT},
      // Fixed search distance with a minimum number of 3 source points
      {.type = YAC_INTERP_DNN_WEIGHT_AVG,
       .search_distance =
         &(struct yac_interp_method_dnn_config_search_distance){
           .type = YAC_INTERP_DNN_SEARCH_DISTANCE_FIXED,
           .search_distance = {.fixed = 1.3 * YAC_RAD}},
       .n_min = 3,
       .gauss_scale = YAC_INTERP_DNN_GAUSS_SCALE_DEFAULT,
       .rbf_scale = YAC_INTERP_DNN_RBF_SCALE_DEFAULT},
      // Fixed search distance with a minimum number of 5 source points
      {.type = YAC_INTERP_DNN_WEIGHT_DIST,
       .search_distance =
         &(struct yac_interp_method_dnn_config_search_distance){
           .type = YAC_INTERP_DNN_SEARCH_DISTANCE_FIXED,
           .search_distance = {.fixed = 1.3 * YAC_RAD}},
       .n_min = 5,
       .gauss_scale = YAC_INTERP_DNN_GAUSS_SCALE_DEFAULT,
       .rbf_scale = YAC_INTERP_DNN_RBF_SCALE_DEFAULT}
    };
    enum {NUM_DNN_CONFIGS = sizeof(dnn_configs)/sizeof(dnn_configs[0])};
 
    for (size_t i = 0; i < NUM_DNN_CONFIGS; ++i) {
 
      double ref_tgt_results[NUM_CELLS_X * NUM_CELLS_Y];
 
      // compute reference data
      if (is_tgt) {
 
        utest_compute_reference_result(
          dnn_configs[i], coordinates_x, coordinates_y,
          (int const *)(&(src_global_corner_mask[0][0])),
          NUM_CELLS_X+1, NUM_CELLS_Y+1,
          cell_coordinates_x, cell_coordinates_y, NUM_CELLS_X, NUM_CELLS_Y,
          ref_tgt_results);
      }
 
      // generate interpolation method stack
      struct interp_method * method_stack[3] =
        {yac_interp_method_dnn_new(dnn_configs[i]),
         yac_interp_method_fixed_new(fixed_value), NULL};
 
      // generate weights
      struct yac_interp_weights * weights =
        yac_interp_method_do_search(method_stack, interp_grid);
      yac_interp_method_delete(method_stack);
 
      for (size_t j = 0; j < NUM_REORDER_TYPES; ++j) {
        enum {COLLECTION_SIZE = 1, COLLECTION_IDX = 0};
 
        struct yac_interpolation * interpolation =
          yac_interp_weights_get_interpolation(
            weights, reorder_type[j], COLLECTION_SIZE,
            YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
 
        // prepare data
        double *** src_data = NULL;
        double ** tgt_data = NULL;
        if (is_tgt) {
 
          tgt_data = xmalloc(COLLECTION_SIZE * sizeof(*tgt_data));
          tgt_data[COLLECTION_IDX] =
            xmalloc(grid_data.num_cells * sizeof(**tgt_data));
          for (size_t k = 0; k < grid_data.num_cells; ++k) {
            tgt_data[COLLECTION_IDX][k] = -999.0;
          }
 
        } else { // is src
 
          src_data = xmalloc(COLLECTION_SIZE * sizeof(*src_data));
 
          // only one field
          src_data[COLLECTION_IDX] = xmalloc(1 * sizeof(**src_data));
          src_data[COLLECTION_IDX][0] =
            xmalloc(grid_data.num_vertices * sizeof(***src_data));
          for (size_t k = 0; k < grid_data.num_vertices; ++k) {
            src_data[COLLECTION_IDX][0][k] =
              (double)(grid_data.vertex_ids[k]);
          }
        }
 
        // do exchange
        yac_interpolation_execute(interpolation, src_data, tgt_data);
 
        // check interpolation result
        if (is_tgt) {
 
          // verify results based on weight type
          for (size_t k = 0; k < grid_data.num_cells; ++k) {
 
            if (fabs(tgt_data[COLLECTION_IDX][k] - ref_tgt_results[k]) > 1e-6) {
              PUT_ERR("wrong result");
            }
          }
        }
 
        // cleanup
        if (is_tgt) {
 
          free(tgt_data[COLLECTION_IDX]);
          free(tgt_data);
 
        } else { // is src
 
          free(src_data[COLLECTION_IDX][0]);
          free(src_data[COLLECTION_IDX]);
          free(src_data);
        }
 
        yac_interpolation_delete(interpolation);
      } // NUM_REORDER_TYPES
      yac_interp_weights_delete(weights);
    } // NUM_STACKS
    yac_interp_grid_delete(interp_grid);
    yac_dist_grid_pair_delete(grid_pair);
    yac_basic_grid_delete(empty_grid);
    yac_basic_grid_delete(grid);
  }
 
  xt_finalize();
  MPI_Finalize();
 
  return TEST_EXIT_CODE;
}
 
static void utest_determine_src_points(
  double const * src_coordinates_x, double const * src_coordinates_y,
  size_t const src_size_x, size_t const src_size_y,
  double const * tgt_coordinates_x, double const * tgt_coordinates_y,
  size_t const tgt_size_x, size_t const tgt_size_y,
  double const inc_angle, int const * src_mask,
  size_t * num_src_points, size_t ** src_points) {
 
  *src_points =
    xmalloc(
      src_size_x * src_size_y * tgt_size_x * tgt_size_y * sizeof(**src_points));
  size_t total_num_src_points = 0;
 
  for (size_t tgt_y = 0, tgt_point_idx = 0; tgt_y < tgt_size_y; ++tgt_y) {
 
    for (size_t tgt_x = 0; tgt_x < tgt_size_x; ++tgt_x, ++tgt_point_idx) {
 
      // compute target point coordinates in 3D
      double tgt_coord[3];
      LLtoXYZ(tgt_coordinates_x[tgt_x], tgt_coordinates_y[tgt_y], tgt_coord);
 
      size_t curr_num_src = 0;
 
      // iterate over all source points
      for (size_t src_y = 0, src_point_idx = 0; src_y < src_size_y; ++src_y) {
        for (size_t src_x = 0; src_x < src_size_x; ++src_x, ++src_point_idx) {
 
          // compute source point coordinates in 3D
          double src_coord[3];
          LLtoXYZ(
            src_coordinates_x[src_x], src_coordinates_y[src_y], src_coord);
 
          // compute angle between target and source point
          double const angle = get_vector_angle(tgt_coord, src_coord);
 
          // check if source point is within the bounding circle and not masked
          if (angle <= inc_angle && src_mask[src_point_idx]) {
            (*src_points)[total_num_src_points] = src_point_idx;
            ++total_num_src_points;
            ++curr_num_src;
          }
        }
      }
 
      num_src_points[tgt_point_idx] = curr_num_src;
    }
  }
 
  *src_points =
    xrealloc(*src_points, total_num_src_points * sizeof(**src_points));
}
 
static double utest_compute_avg_result(
  size_t const n, size_t const * src_indices, size_t const tgt_index,
  double const * src_coordinates_x, double const * src_coordinates_y,
  size_t const src_size_x,
  double const * tgt_coordinates_x, double const * tgt_coordinates_y,
  size_t const tgt_size_x) {
 
  UNUSED(tgt_index);
  UNUSED(src_coordinates_x);
  UNUSED(src_coordinates_y);
  UNUSED(src_size_x);
  UNUSED(tgt_coordinates_x);
  UNUSED(tgt_coordinates_y);
  UNUSED(tgt_size_x);
 
  double sum = 0.0;
  for (size_t i = 0; i < n; ++i) {
    sum += (double)(src_indices[i]);
  }
  return sum / (double)n;
}
 
static double utest_compute_dist_result(
  size_t const n, size_t const * src_indices, size_t const tgt_index,
  double const * src_coordinates_x, double const * src_coordinates_y,
  size_t const src_size_x,
  double const * tgt_coordinates_x, double const * tgt_coordinates_y,
  size_t const tgt_size_x) {
 
  double tgt_coord[3];
  LLtoXYZ(
    tgt_coordinates_x[tgt_index%tgt_size_x],
    tgt_coordinates_y[tgt_index/tgt_size_x], tgt_coord);
 
  double result = 0.0;
  double weights_sum = 0.0;
 
  for (size_t i = 0; i < n; ++i) {
 
    double src_coord[3];
    LLtoXYZ(src_coordinates_x[src_indices[i]%src_size_x],
            src_coordinates_y[src_indices[i]/src_size_x], src_coord);
 
    double weight = 1.0 / get_vector_angle(tgt_coord, src_coord);
    result += weight * (double)src_indices[i];
    weights_sum += weight;
  }
  return result / weights_sum;
}
 
static double utest_compute_gauss_result(
  size_t const n, size_t const * src_indices, size_t const tgt_index,
  double const * src_coordinates_x, double const * src_coordinates_y,
  size_t const src_size_x,
  double const * tgt_coordinates_x, double const * tgt_coordinates_y,
  size_t const tgt_size_x) {
 
  // return early if there is only a single source point
  if (n == 1) {
    return (double)src_indices[0];
  }
 
  double src_coords[n][3];
 
  for (size_t i = 0; i < n; ++i) {
    LLtoXYZ(src_coordinates_x[src_indices[i]%src_size_x],
            src_coordinates_y[src_indices[i]/src_size_x], src_coords[i]);
  }
 
  // "- n" because we do not count the diagonal
  double src_distances_sum = 0.0;
  for (size_t i = 0; i < n; ++i)
    for (size_t j = 0; j < n; ++j)
      src_distances_sum += get_vector_angle(src_coords[i], src_coords[j]);
 
  size_t src_distance_count = n * n - n;
 
  // compute mean distance
  double src_distance_mean = src_distances_sum / (double)src_distance_count;
  double src_distance_mean_sq = src_distance_mean * src_distance_mean;
 
  // compute weights
  double tgt_coord[3];
  LLtoXYZ(tgt_coordinates_x[tgt_index%tgt_size_x],
          tgt_coordinates_y[tgt_index/tgt_size_x], tgt_coord);
  double weights[n];
  for (size_t i = 0; i < n; ++i) {
    double tgt_distance = get_vector_angle(src_coords[i], tgt_coord);
    weights[i] =
      exp(-1.0 * (tgt_distance * tgt_distance) /
                 (YAC_INTERP_DNN_GAUSS_SCALE_DEFAULT * src_distance_mean_sq));
  }
 
  // compute sum of weights
  double weights_sum = 0.0;
  for (size_t i = 0; i < n; ++i) weights_sum += weights[i];
  for (size_t i = 0; i < n; ++i) weights[i] /= weights_sum;
 
  // compute interpolation results
  double result = 0.0;
  for (size_t i = 0; i < n; ++i) {
    result += weights[i] * (double)src_indices[i];
  }
 
  return result;
}
 
static void utest_inverse(double * A, size_t n);
 
static double utest_compute_rbf_result(
  size_t const n, size_t const * src_indices, size_t const tgt_index,
  double const * src_coordinates_x, double const * src_coordinates_y,
  size_t const src_size_x,
  double const * tgt_coordinates_x, double const * tgt_coordinates_y,
  size_t const tgt_size_x) {
 
  double src_coords[n][3];
  for (size_t i = 0; i < n; ++i) {
    LLtoXYZ(src_coordinates_x[src_indices[i]%src_size_x],
            src_coordinates_y[src_indices[i]/src_size_x], src_coords[i]);
  }
 
  double A[n][n];
  double sum_d = 0.0;
 
  // compute distance matrix for all source points
  for (size_t i = 0; i < n; ++i) A[i][i] = 0.0;
  for (size_t i = 0; i < n - 1; ++i) {
    for (size_t j = i + 1; j < n; ++j) {
      double d = get_vector_angle(src_coords[i], src_coords[j]);
      A[i][j] = d;
      A[j][i] = d;
      sum_d += d;
    }
  }
 
  // compute and apply scale factor for distance matrix
  double scale_d = 1.0;
  if (sum_d > 0.0)
    scale_d = ((double)((n - 1) * n)) / (2.0 * sum_d);
  scale_d /= YAC_INTERP_DNN_RBF_SCALE_DEFAULT;
 
  double sq_scale_d = scale_d * scale_d;
 
  // compute matrix A[n][n] with rbf(A)
  for (size_t i = 0; i < n; ++i) {
    for (size_t j = 0; j < n; ++j) {
      double d = A[i][j];
      A[i][j] = exp(-1.0 * d * d * sq_scale_d);
    }
  }
 
  // compute inverse of A
  utest_inverse(&A[0][0], n);
 
  // compute a[n] - distances from target to sources with rbf applied
  double tgt_coord[3];
  LLtoXYZ(tgt_coordinates_x[tgt_index%tgt_size_x],
          tgt_coordinates_y[tgt_index/tgt_size_x], tgt_coord);
 
  double a[n];
  for (size_t i = 0; i < n; ++i) {
    double d = get_vector_angle(tgt_coord, src_coords[i]);
    a[i] = exp(-1.0 * d * d * sq_scale_d);
  }
 
  // compute weights: w_i = SUM(A_inv[i][j]*a[j])
  double weights[n];
  for (size_t i = 0; i < n; ++i) {
    weights[i] = 0.0;
    for (size_t j = 0; j < n; ++j) weights[i] += A[i][j] * a[j];
  }
 
  // compute interpolation result
  double result = 0.0;
  for (size_t i = 0; i < n; ++i) {
    result += weights[i] * (double)src_indices[i];
  }
 
  return result;
}
 
static void utest_inverse(double * A, size_t n) {
 
// LAPACKE_dsytrf_work and LAPACKE_dsytri might not be available
// (see yac_lapack_interface.h).
#ifdef YAC_LAPACK_NO_DSYTR
  lapack_int ipiv[n+1];
  double work[n*n];
 
  for (size_t i = 0; i < n+1; ++i) ipiv[i] = 0;
  for (size_t i = 0; i < n*n; ++i) work[i] = 0;
 
  YAC_ASSERT(
    !LAPACKE_dgetrf(
      LAPACK_COL_MAJOR, (lapack_int) n, (lapack_int) n,
      A, (lapack_int) n, ipiv), "internal ERROR: dgetrf")
 
  YAC_ASSERT(
    !LAPACKE_dgetri_work(
      LAPACK_COL_MAJOR, (lapack_int) n, A, (lapack_int) n,
      ipiv, work, (lapack_int) (n * n)), "internal ERROR: dgetri")
#else
  lapack_int ipiv[n];
  double work[n];
 
  YAC_ASSERT(
    !LAPACKE_dsytrf_work(
      LAPACK_COL_MAJOR, 'L', (lapack_int) n , A,
      (lapack_int) n, ipiv, work, (lapack_int) n), "internal ERROR: dsytrf")
 
  YAC_ASSERT(
    !LAPACKE_dsytri_work(
      LAPACK_COL_MAJOR, 'L', (lapack_int) n , A,
      (lapack_int) n, ipiv, work), "internal ERROR: dsytri")
 
  for (size_t i = 0; i < n; ++i)
    for (size_t j = i + 1; j < n; ++j)
      A[j*n+i] = A[i*n+j];
#endif
}
 
static void utest_compute_reference_result(
  struct yac_interp_method_dnn_config const dnn_config,
  double const * src_coordinates_x, double const * src_coordinates_y,
  int const * src_global_mask, size_t const src_size_x, size_t const src_size_y,
  double const * tgt_coordinates_x, double const * tgt_coordinates_y,
  size_t const tgt_size_x, size_t const tgt_size_y,
  double * ref_tgt_results) {
 
  // determine search distance
  double search_distance;
  YAC_ASSERT(dnn_config.search_distance != NULL, "search_distance is NULL")
  switch (dnn_config.search_distance->type) {
    YAC_UNREACHABLE_DEFAULT("invalid search distance type");
    case (YAC_INTERP_DNN_SEARCH_DISTANCE_FIXED): {
      search_distance = dnn_config.search_distance->search_distance.fixed;
      break;
    }
    case (YAC_INTERP_DNN_SEARCH_DISTANCE_CELL_AREA): {
      // the search distance computation is simplified and assumes a planar cell
      // and bounding circle, but for this test it should be sufficient
      double cell_area = YAC_RAD * YAC_RAD;
      search_distance =
        dnn_config.search_distance->search_distance.scale *
        sqrt(cell_area / M_PI);
      break;
    }
  }
 
  // compute source points overlapping with bounding circles of each
  // target point
  size_t * src_points = NULL;
  size_t num_src_per_tgt[tgt_size_x * tgt_size_y];
  utest_determine_src_points(
    src_coordinates_x, src_coordinates_y, src_size_x, src_size_y,
    tgt_coordinates_x, tgt_coordinates_y, tgt_size_x, tgt_size_y,
    search_distance, src_global_mask, num_src_per_tgt, &src_points);
 
  // get pointer to weight computation function
  // (based on dnn_configs[i].type)
  double (*compute_result)(
    size_t const, size_t const *, size_t const,
    double const *, double const *, size_t const,
    double const *, double const *, size_t const) = NULL;
  switch (dnn_config.type) {
    YAC_UNREACHABLE_DEFAULT("invalid DNN_WEIGHT_TYPE");
    case YAC_INTERP_DNN_WEIGHT_AVG:
      compute_result = utest_compute_avg_result;
      break;
    case YAC_INTERP_DNN_WEIGHT_DIST:
      compute_result = utest_compute_dist_result;
      break;
    case YAC_INTERP_DNN_WEIGHT_GAUSS:
      compute_result = utest_compute_gauss_result;
      break;
    case YAC_INTERP_DNN_WEIGHT_RBF:
      compute_result = utest_compute_rbf_result;
      break;
  }
 
  // determine the minimum number of source points that have to be found
  // within the search distance in order to interpolate a target point
  size_t n_min = (size_t)dnn_config.n_min;
 
  // compute reference tgt results
  size_t src_offset = 0;
  for (size_t k = 0; k < tgt_size_x * tgt_size_y; ++k) {
 
    if (num_src_per_tgt[k] >= n_min) {
      size_t * curr_src_points = &src_points[src_offset];
      ref_tgt_results[k] =
        compute_result(
          num_src_per_tgt[k], curr_src_points, k,
          src_coordinates_x, src_coordinates_y, src_size_x,
          tgt_coordinates_x, tgt_coordinates_y, tgt_size_x);
    } else {
      // No valid source points in search distance, use fallback value
      ref_tgt_results[k] = fixed_value;
    }
    src_offset += num_src_per_tgt[k];
  }
 
  free(src_points);
}