test_point_sphere_part.c

Go to the documentation of this file.

// Copyright (c) 2024 The YAC Authors
//
// SPDX-License-Identifier: BSD-3-Clause
 
#include <stdlib.h>
#include <stdio.h>
#include <float.h>
#include <string.h>
 
#include "tests.h"
#include "test_common.h"
#include "grids/basic_grid.h"
#include "sphere_part.h"
#include "geometry.h"
 
static void utest_random_test(size_t const num_points);
static void utest_zero_point_test();
static void utest_single_point_test();
static void utest_empty_branch_test(size_t num_points);
static void utest_zero_balance_point();
static void utest_duplicated_point_test();
static void utest_bnd_circle_search();
static void utest_bnd_circle_random_test(size_t num_points, size_t num_circles);
 
int main(void) {
 
  srand(1365441);
 
  for (size_t i = 0; i < 10; ++i) {
    utest_random_test(1);
    utest_random_test(10);
    utest_random_test(100);
    utest_random_test(650);
  }
  utest_zero_point_test();
  utest_single_point_test();
  utest_empty_branch_test(4);
  utest_empty_branch_test(128);
  utest_zero_balance_point();
  utest_duplicated_point_test();
  utest_bnd_circle_search();
  for (size_t i = 0; i < 10; ++i) {
    utest_bnd_circle_random_test(10, 5);
    utest_bnd_circle_random_test(100, 20);
    utest_bnd_circle_random_test(500, 50);
  }
 
  return TEST_EXIT_CODE;
}
 
struct distance_index {
  double distance;
  size_t index;
};
 
static int utest_compare_distance_index(const void * a, const void * b) {
 
  int ret = (((struct distance_index*)a)->distance >
             ((struct distance_index*)b)->distance) -
            (((struct distance_index*)a)->distance <
             ((struct distance_index*)b)->distance);
  if (ret == 0)
    ret = (((struct distance_index*)a)->index >
           ((struct distance_index*)b)->index) -
          (((struct distance_index*)a)->index <
           ((struct distance_index*)b)->index);
 
  return ret;
}
 
static void utest_point_sphere_test_NNN(
  struct point_sphere_part_search * search, size_t n,
  double (*sphere_part_coordinates_xyz)[3],
  double (*search_coordinates_xyz)[3], size_t num_points,
  size_t ** local_point_ids, size_t * local_point_ids_array_size,
  size_t * num_local_point_ids, struct distance_index * distances) {
 
  for (int with_result_coords = 0; with_result_coords < 2;
    ++with_result_coords) {
 
    double (*result_coordinates_xyz)[3] = NULL;
    size_t result_coordinates_xyz_array_size = 0;
 
    yac_point_sphere_part_search_NNN(
      search, num_points, search_coordinates_xyz, n, NULL, NULL,
      with_result_coords?&result_coordinates_xyz:NULL,
      with_result_coords?&result_coordinates_xyz_array_size:NULL,
      local_point_ids, local_point_ids_array_size, num_local_point_ids);
 
    // check results
    for (size_t j = 0, offset = 0; j < num_points; ++j) {
 
      struct distance_index * curr_distances = distances + j * num_points;
 
      size_t ref_num_local_points = 1, m = 0;
      for (; (ref_num_local_points < num_points); ++ref_num_local_points)
        if (utest_double_are_unequal(
              curr_distances[ref_num_local_points-1].distance,
              curr_distances[ref_num_local_points].distance))
          if (++m == n) break;
 
      if (num_local_point_ids[j] != ref_num_local_points) {
        PUT_ERR("wrong number of local points found\n");
        continue;
      }
 
      size_t num_matching_points = 0;
      for (size_t k = 0; k < ref_num_local_points; ++k) {
        for (size_t l = 0; l < ref_num_local_points; ++l) {
          if ((curr_distances[k].index == (*local_point_ids)[offset+l]) &&
              (!with_result_coords ||
               (points_are_identically(
                  sphere_part_coordinates_xyz[curr_distances[k].index],
                  result_coordinates_xyz[offset+l])))) {
            ++num_matching_points;
          }
        }
      }
 
      if (num_matching_points != ref_num_local_points)
        PUT_ERR("wrong results\n");
 
      offset += ref_num_local_points;
    }
 
    free(result_coordinates_xyz);
  }
}
 
static void utest_point_sphere_test_NNN_ubound(
  struct point_sphere_part_search * search, size_t num_points,
  yac_coordinate_pointer xyz_coordinates, size_t n,
  struct distance_index * distances) {
 
  struct sin_cos_angle angles[num_points];
 
  yac_point_sphere_part_search_NNN_ubound(
    search, num_points, xyz_coordinates, n, angles);
 
  for (size_t i = 0; i < num_points; ++i) {
    double ref_ubound =
      distances[i * num_points + n - 1].distance;
 
    if (compare_angles(
          sin_cos_angle_new(sin(ref_ubound), cos(ref_ubound)),
          angles[i]) > 0)
      PUT_ERR("wrong result");
  }
}
 
static void utest_point_sphere_test_NNN_bnd_circle(
  struct point_sphere_part_search * search, size_t num_points,
  yac_coordinate_pointer xyz_coordinates, size_t n,
  struct distance_index * distances) {
 
  if (n == 0) return;
 
  struct bounding_circle bnd_circles[num_points];
 
  size_t cut_off = n / 2;
 
  for (size_t i = 0; i < num_points; ++i) {
 
    memcpy(
      bnd_circles[i].base_vector, xyz_coordinates[i],
      sizeof(*xyz_coordinates));
    double distance =
      distances[
        i * num_points + ((i < num_points / 2)?cut_off:(n-1))].distance +
        yac_angle_tol;
    bnd_circles[i].inc_angle =
      sin_cos_angle_new(sin(distance), cos(distance));
  }
 
  size_t * local_point_ids = NULL;
  size_t local_point_ids_array_size = 0;
  size_t * num_local_point_ids =
    malloc(num_points * sizeof(*num_local_point_ids));
 
  yac_point_sphere_part_search_NNN_bnd_circle(
    search, num_points, bnd_circles, n, &local_point_ids,
    &local_point_ids_array_size, num_local_point_ids);
 
  size_t offset = 0;
  for (size_t i = 0; i < num_points; ++i) {
 
    size_t ref_num_local_point_ids = (i < num_points / 2)?(n / 2 + 1):n;
    double ref_distance =
      distances[i * num_points + ref_num_local_point_ids - 1].distance;
    for (; (ref_num_local_point_ids < MIN(n, num_points)) &&
         (distances[
           i * num_points + ref_num_local_point_ids].distance == ref_distance);
         ++ref_num_local_point_ids);
 
    if (ref_num_local_point_ids != num_local_point_ids[i])
      PUT_ERR("wrong result");
 
    for (size_t j = 0; j < num_local_point_ids[i]; ++j) {
      int match_flag = 0;
      for (size_t k = 0; (k < num_local_point_ids[i]) && !match_flag; ++k)
        match_flag =
          distances[i * num_points + k].index == local_point_ids[offset + j];
      if (!match_flag) PUT_ERR("wrong result");
    }
    offset += num_local_point_ids[i];
  }
 
  free(local_point_ids);
  free(num_local_point_ids);
}
 
static void utest_gen_random_point(double * point) {
 
  double x_coordinate = 2.0 * M_PI * (((double)rand()) / ((double)RAND_MAX));
  double y_coordinate = M_PI * (((double)rand()) / ((double)RAND_MAX)) - M_PI_2;
  LLtoXYZ(x_coordinate, y_coordinate, point);
}
 
static void utest_gen_random_bnd_circle(struct bounding_circle * circle) {
 
  utest_gen_random_point(circle->base_vector);
  double angle = M_PI_2 * (((double)rand()) / ((double)RAND_MAX));
  circle->inc_angle = sin_cos_angle_new(sin(angle), cos(angle));
  circle->sq_crd = DBL_MAX;
}
 
static void utest_random_test(size_t const num_points) {
 
  double (*coordinates_xyz)[3];
  double (*coordinates_xyz_a)[3];
  double (*coordinates_xyz_b)[3];
  struct bounding_circle * bnd_circles =
    xmalloc(num_points * sizeof(*bnd_circles));
  yac_int * global_ids = xmalloc(num_points * sizeof(*global_ids));
 
  coordinates_xyz = xmalloc(2 * num_points * sizeof(*coordinates_xyz));
  coordinates_xyz_a = coordinates_xyz;
  coordinates_xyz_b = coordinates_xyz + num_points;
 
  size_t * local_point_ids = NULL;
  size_t local_point_ids_array_size = 0;
  size_t * num_local_point_ids =
    xmalloc(num_points * sizeof(*num_local_point_ids));
  struct distance_index * distances =
    xmalloc(num_points * num_points * sizeof(*distances));
 
  for (size_t i = 0; i < num_points; ++i) {
    utest_gen_random_bnd_circle(bnd_circles + i);
    global_ids[i] = (yac_int)i;
  }
  for (size_t k = 0; k < 2 * num_points; ++k)
    utest_gen_random_point(coordinates_xyz[k]);
 
  for (size_t j = 0; j < num_points; ++j) {
    for (size_t k = 0; k < num_points; ++k) {
 
      distances[j * num_points + k].distance =
        get_vector_angle(coordinates_xyz_a[k], coordinates_xyz_b[j]);
      distances[j * num_points + k].index = k;
    }
    qsort(
      distances + j * num_points, num_points, sizeof(*distances),
      utest_compare_distance_index);
  }
 
  struct point_sphere_part_search * search =
    yac_point_sphere_part_search_new(
      num_points, (yac_const_coordinate_pointer)coordinates_xyz_a,
      global_ids);
 
  for (size_t n = 1; n < 16; ++n)
    utest_point_sphere_test_NNN(
      search, n, coordinates_xyz_a, coordinates_xyz_b, num_points,
      &local_point_ids, &local_point_ids_array_size, num_local_point_ids,
      distances);
 
  for (size_t n = 1; n < MIN(num_points,16); ++n)
    utest_point_sphere_test_NNN_ubound(
      search, num_points, coordinates_xyz_b, n, distances);
 
  for (size_t n = 1; n < MIN(num_points,16); ++n)
    utest_point_sphere_test_NNN_bnd_circle(
      search, num_points, coordinates_xyz_b, n, distances);
 
  yac_delete_point_sphere_part_search(search);
 
  free(distances);
  free(local_point_ids);
  free(num_local_point_ids);
  free(global_ids);
  free(bnd_circles);
  free(coordinates_xyz);
}
 
static void utest_single_point_test() {
 
  struct point_sphere_part_search * search;
  double xyz_coordinate[1][3];
  {
    double x_coordinate = 0, y_coordinate = 0;
    yac_int global_id[1] = {0};
    LLtoXYZ(x_coordinate, y_coordinate, xyz_coordinate[0]);
    search =
      yac_point_sphere_part_search_new(
        1, (yac_const_coordinate_pointer)xyz_coordinate, global_id);
  }
 
  double x_coordinates[] = {0.0,90.0,180.0,-90.0, 0.0,  0.0};
  double y_coordinates[] = {0.0, 0.0,  0.0,  0.0,90.0,-90.0};
  size_t num_points = sizeof(x_coordinates) / sizeof(x_coordinates[0]);
  double xyz_coordinates[num_points][3];
  for (size_t i = 0; i < num_points; ++i)
    LLtoXYZ_deg(x_coordinates[i], y_coordinates[i], xyz_coordinates[i]);
 
  size_t * local_point_ids = NULL;
  size_t local_point_ids_array_size = 0;
  size_t * num_local_point_ids =
    xmalloc(num_points * sizeof(*num_local_point_ids));
 
  yac_point_sphere_part_search_NN(
    search, num_points, xyz_coordinates, NULL, NULL, NULL, &local_point_ids,
    &local_point_ids_array_size, num_local_point_ids);
 
  for (size_t i = 0; i < num_points; ++i)
    if ((num_local_point_ids[i] != 1) || (local_point_ids[i] != 0))
      PUT_ERR("utest_single_point_test: wrong result\n");
 
  struct sin_cos_angle angles[num_points];
  yac_point_sphere_part_search_NNN_ubound(
    search, num_points, xyz_coordinates, 1, angles);
 
  for (size_t i = 0; i < num_points; ++i)
    if (compare_angles(
          angles[i],
          get_vector_angle_2(
            xyz_coordinate[0], xyz_coordinates[i])) != 0)
      PUT_ERR("wrong ubound");
 
  free(local_point_ids);
  free(num_local_point_ids);
  yac_delete_point_sphere_part_search(search);
}
 
static void utest_zero_point_test() {
 
  struct point_sphere_part_search * search =
    yac_point_sphere_part_search_new(0, NULL, NULL);
 
  double x_coordinates[] = {0.0,90.0,180.0,-90.0, 0.0,  0.0};
  double y_coordinates[] = {0.0, 0.0,  0.0,  0.0,90.0,-90.0};
  enum {NUM_POINTS = sizeof(x_coordinates) / sizeof(x_coordinates[0])};
  double xyz_coordinates[NUM_POINTS][3];
  for (size_t i = 0; i < NUM_POINTS; ++i)
    LLtoXYZ_deg(x_coordinates[i], y_coordinates[i], xyz_coordinates[i]);
 
  struct bounding_circle bnd_circles[NUM_POINTS];
  for (size_t i = 0; i < NUM_POINTS; ++i) {
    LLtoXYZ_deg(
      x_coordinates[i], y_coordinates[i], bnd_circles[i].base_vector);
    bnd_circles[i].inc_angle = sin_cos_angle_new(sin(0.1), cos(0.1));
  }
 
  size_t * local_point_ids = NULL;
  size_t local_point_ids_array_size = 0;
  size_t num_local_point_ids[NUM_POINTS];
 
  for (size_t i = 0; i < NUM_POINTS; ++i)
    num_local_point_ids[i] = (size_t)-1;
  yac_point_sphere_part_search_NN(
    search, NUM_POINTS, xyz_coordinates, NULL, NULL, NULL, &local_point_ids,
    &local_point_ids_array_size, num_local_point_ids);
 
  for (size_t i = 0; i < NUM_POINTS; ++i)
    if (num_local_point_ids[i] != 0) PUT_ERR("utest_zero_point_test: wrong result\n");
 
  for (size_t i = 0; i < NUM_POINTS; ++i)
    num_local_point_ids[i] = (size_t)-1;
  yac_point_sphere_part_search_NNN(
    search, NUM_POINTS, xyz_coordinates, 16, NULL, NULL, NULL, NULL,
    &local_point_ids, &local_point_ids_array_size, num_local_point_ids);
 
  for (size_t i = 0; i < NUM_POINTS; ++i)
    if (num_local_point_ids[i] != 0) PUT_ERR("utest_zero_point_test: wrong result\n");
 
  for (size_t i = 0; i < NUM_POINTS; ++i)
    num_local_point_ids[i] = (size_t)-1;
  yac_point_sphere_part_search_NNN_bnd_circle(
    search, NUM_POINTS, bnd_circles, 1, &local_point_ids,
    &local_point_ids_array_size, num_local_point_ids);
 
  for (size_t i = 0; i < NUM_POINTS; ++i)
    if (num_local_point_ids[i] != 0) PUT_ERR("utest_zero_point_test: wrong result\n");
 
  free(local_point_ids);
  yac_delete_point_sphere_part_search(search);
}
 
static void utest_empty_branch_test(size_t num_points_a) {
 
  yac_coordinate_pointer xyz_coordinates_a =
    malloc(num_points_a * sizeof(*xyz_coordinates_a));
  yac_int * global_ids = malloc(num_points_a * sizeof(*global_ids));
  for (size_t i = 0; i < num_points_a; ++i) {
    LLtoXYZ_deg(0.0, -5.0 + (10.0 / (double)num_points_a) * (double)i,
                xyz_coordinates_a[i]);
    global_ids[i] = (yac_int)i;
  }
 
  struct point_sphere_part_search * search =
    yac_point_sphere_part_search_new(
      num_points_a, (yac_const_coordinate_pointer)xyz_coordinates_a,
      global_ids);
  free(global_ids);
 
  double x_coordinates_b[] = {-1.0,-1.0, 1.0,1.0};
  double y_coordinates_b[] = {-1.0, 1.0,-1.0,1.0};
  size_t num_points_b = sizeof(x_coordinates_b) / sizeof(x_coordinates_b[0]);
  double xyz_coordinates_b[num_points_b][3];
  for (size_t i = 0; i < num_points_b; ++i)
    LLtoXYZ_deg(x_coordinates_b[i], y_coordinates_b[i], xyz_coordinates_b[i]);
 
  size_t * local_point_ids = NULL;
  size_t local_point_ids_array_size = 0;
  size_t * num_local_point_ids =
    xmalloc(num_points_b * sizeof(*num_local_point_ids));
 
  struct distance_index distances[num_points_a * num_points_b];
 
  for (size_t i = 0; i < num_points_b; ++i) {
 
    for (size_t j = 0; j < num_points_a; ++j) {
 
      distances[i * num_points_a + j].distance =
        get_vector_angle(xyz_coordinates_a[j], xyz_coordinates_b[i]);
      distances[i * num_points_a + j].index = j;
    }
    qsort(distances + i * num_points_a, num_points_a, sizeof(*distances),
          utest_compare_distance_index);
  }
 
  for (size_t i = 0; i < num_points_b; ++i)
    num_local_point_ids[i] = (size_t)-1;
  yac_point_sphere_part_search_NN(
    search, num_points_b, xyz_coordinates_b, NULL, NULL, NULL,
    &local_point_ids, &local_point_ids_array_size, num_local_point_ids);
 
    // check results
  for (size_t j = 0, offset = 0; j < num_points_b; ++j) {
 
    struct distance_index * curr_distances = distances + j * num_points_a;
 
    size_t ref_num_local_points;
    for (ref_num_local_points = 1;
         (ref_num_local_points < num_points_a); ++ref_num_local_points)
      if (utest_double_are_unequal(
            curr_distances[ref_num_local_points-1].distance,
            curr_distances[ref_num_local_points].distance)) break;
 
    if (num_local_point_ids[j] != ref_num_local_points) {
      PUT_ERR("wrong number of local points found\n");
      continue;
    }
 
    size_t num_matching_points = 0;
    for (size_t k = 0; k < ref_num_local_points; ++k)
      for (size_t l = 0; l < ref_num_local_points; ++l)
        if (curr_distances[k].index == local_point_ids[offset+l])
          ++num_matching_points;
 
    if (num_matching_points != ref_num_local_points)
      PUT_ERR("wrong results\n");
 
    offset += ref_num_local_points;
  }
  free(local_point_ids);
  free(num_local_point_ids);
 
  struct sin_cos_angle angles[num_points_b];
 
  yac_point_sphere_part_search_NNN_ubound(
    search, num_points_b, xyz_coordinates_b, 1, angles);
 
  // check results
  for (size_t j = 0; j < num_points_b; ++j) {
 
    struct sin_cos_angle min_angle = SIN_COS_M_PI;
 
    for (size_t i = 0; i < num_points_a; ++i) {
 
      struct sin_cos_angle curr_angle =
        get_vector_angle_2(
          xyz_coordinates_a[i], xyz_coordinates_b[j]);
 
      if (compare_angles(curr_angle, min_angle) < 0)
        min_angle = curr_angle;
    }
 
    if (compare_angles(min_angle, angles[j]))
      PUT_ERR("wrong results");
  }
 
  free(xyz_coordinates_a);
 
  yac_delete_point_sphere_part_search(search);
}
 
static void utest_zero_balance_point() {
 
  struct point_sphere_part_search * search;
  {
    double x_coordinate = 0, y_coordinate = 0;
    double xyz_coordinate[6][3] =
      {{1,0,0}, {-1,0,0},
       {0,1,0}, {0,-1,0},
       {0,0,1}, {0,0,-1}};
    yac_int global_ids[6] = {0,1,2,3,4,5};
    LLtoXYZ(x_coordinate, y_coordinate, xyz_coordinate[0]);
    search =
      yac_point_sphere_part_search_new(
        6, (yac_const_coordinate_pointer)xyz_coordinate, global_ids);
  }
  yac_delete_point_sphere_part_search(search);
}
 
static void utest_duplicated_point_test() {
 
  struct point_sphere_part_search * search;
  {
    double x_coordinates[] = {0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
                              1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,
                              2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,
                              3.0,3.0,3.0,3.0,3.0,3.0,3.0,3.0,3.0,3.0,
                              4.0,4.0,4.0,4.0,4.0,4.0,4.0,4.0,4.0,4.0,
                              2.0,2.0,2.0};
    double y_coordinates[] = {0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,
                              0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,
                              0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,
                              0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,
                              0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,
                              4.0,4.0,4.0};
    enum {NUM_POINTS = sizeof(x_coordinates) / sizeof(x_coordinates[0])};
    double xyz_coordinates[NUM_POINTS][3];
    yac_int global_ids[NUM_POINTS] = {
      0,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,
      100,100,101};
    for (size_t i = 0; i < NUM_POINTS; ++i) {
      LLtoXYZ_deg(x_coordinates[i], y_coordinates[i], xyz_coordinates[i]);
    }
    search =
      yac_point_sphere_part_search_new(
        NUM_POINTS, (yac_const_coordinate_pointer)xyz_coordinates,
        global_ids);
  }
 
  double x_coordinates[] = {1.0,1.0,2.0,2.0,2.0,3.0,3.0,2.1,0.1,2.5};
  double y_coordinates[] = {1.0,8.0,3.0,4.0,5.0,1.0,8.0,4.0,0.0,3.0};
  enum {
    NUM_POINTS = sizeof(x_coordinates) / sizeof(x_coordinates[0]),
    MAX_NUM_RESULTS = 2};
  double xyz_coordinates[NUM_POINTS][3];
  for (size_t i = 0; i < NUM_POINTS; ++i)
    LLtoXYZ_deg(x_coordinates[i], y_coordinates[i], xyz_coordinates[i]);
 
  size_t * local_point_ids = NULL;
  size_t local_point_ids_array_size = 0;
  size_t num_local_point_ids[NUM_POINTS];
 
  double (*result_coordinates_xyz)[3] = NULL;
  size_t result_coordinates_xyz_array_size = 0;
 
  yac_point_sphere_part_search_NN(
    search, NUM_POINTS, xyz_coordinates, NULL,
    &result_coordinates_xyz, &result_coordinates_xyz_array_size,
    &local_point_ids, &local_point_ids_array_size, num_local_point_ids);
 
  size_t ref_local_point_ids[NUM_POINTS][MAX_NUM_RESULTS] =
    {{11}, {18}, {23}, {52}, {25}, {31}, {38}, {52}, {0}, {23,33}};
  size_t ref_num_local_point_ids[NUM_POINTS] = {1,1,1,1,1,1,1,1,1,2};
  double ref_x_coordinates[NUM_POINTS][MAX_NUM_RESULTS] =
    {{1.0},{1.0},{2.0},{2.0},{2.0},{3.0},{3.0},{2.0},{0.0},{2.0,3.0}};
  double ref_y_coordinates[NUM_POINTS][MAX_NUM_RESULTS] =
    {{1.0},{8.0},{3.0},{4.0},{5.0},{1.0},{8.0},{4.0},{0.0},{3.0,3.0}};
  double ref_result_coordinates_xyz[NUM_POINTS][MAX_NUM_RESULTS][3];
 
  for (size_t i = 0; i < NUM_POINTS; ++i)
    for (size_t j = 0; j < ref_num_local_point_ids[i]; ++j)
      LLtoXYZ_deg(
        ref_x_coordinates[i][j], ref_y_coordinates[i][j],
        ref_result_coordinates_xyz[i][j]);
 
  for (size_t i = 0, offset = 0; i < NUM_POINTS;
       offset += num_local_point_ids[i++]) {
 
    if (num_local_point_ids[i] != ref_num_local_point_ids[i]) {
      PUT_ERR("ERROR in yac_point_sphere_part_search_NN (num_local_point_ids)");
    } else {
 
      int match_flag[ref_num_local_point_ids[i]];
      size_t match_count = 0;
      memset(match_flag, 0, ref_num_local_point_ids[i] * sizeof(match_flag[0]));
 
      for (size_t j = 0; j < ref_num_local_point_ids[i]; ++j) {
 
        for (size_t k = 0; k < ref_num_local_point_ids[i]; ++k) {
 
          if (points_are_identically(
                result_coordinates_xyz[offset+j],
                ref_result_coordinates_xyz[i][k]) &&
              (local_point_ids[offset+j] == ref_local_point_ids[i][k])) {
 
            if (match_flag[j]) {
              PUT_ERR("ERROR in yac_point_sphere_part_search_NN "
                      "duplicated result");
            } else {
              match_flag[j] = 1;
              match_count++;
            }
          }
        }
      }
      if (match_count != ref_num_local_point_ids[i]) {
        PUT_ERR("ERROR in yac_point_sphere_part_search_NN result mismatch");
      }
    }
  }
 
  free(local_point_ids);
  free(result_coordinates_xyz);
  yac_delete_point_sphere_part_search(search);
}
 
static void utest_bnd_circle_search() {
 
  // Test 1: basic bounding circle search with known points
  {
    double x_coordinates[] = {0.0, 10.0, 20.0, 30.0, 40.0};
    double y_coordinates[] = {0.0,  0.0,  0.0,  0.0,  0.0};
    enum {NUM_POINTS = sizeof(x_coordinates) / sizeof(x_coordinates[0])};
    double xyz_coordinates[NUM_POINTS][3];
    yac_int global_ids[NUM_POINTS];
    for (size_t i = 0; i < NUM_POINTS; ++i) {
      LLtoXYZ_deg(x_coordinates[i], y_coordinates[i], xyz_coordinates[i]);
      global_ids[i] = (yac_int)i;
    }
 
    struct point_sphere_part_search * search =
      yac_point_sphere_part_search_new(
        NUM_POINTS, (yac_const_coordinate_pointer)xyz_coordinates, global_ids);
 
    // Create a bounding circle centered at (15, 0)/(35, 0) with angle
    // covering points at (10, 0)/(30, 0) and (20, 0)/(40, 0)
    struct bounding_circle bnd_circles[2];
    double inc_angle = 6.0 * YAC_RAD; // ~6 degrees
    LLtoXYZ_deg(15.0, 0.0, bnd_circles[0].base_vector);
    bnd_circles[0].inc_angle =
      sin_cos_angle_new(sin(inc_angle), cos(inc_angle));
    bnd_circles[0].sq_crd = DBL_MAX;
    LLtoXYZ_deg(35.0, 0.0, bnd_circles[1].base_vector);
    bnd_circles[1].inc_angle =
      sin_cos_angle_new(sin(inc_angle), cos(inc_angle));
    bnd_circles[1].sq_crd = DBL_MAX;
 
    size_t * local_point_ids = NULL;
    size_t local_point_ids_array_size = 0;
    size_t num_local_point_ids[2];
 
    // Perform search
    yac_point_sphere_part_search_bnd_circle(
      search, 2, bnd_circles,
      &local_point_ids, &local_point_ids_array_size, num_local_point_ids);
 
    struct {
      size_t count;
      size_t * ids;
    } ref_results[2] =
      {{.count = 2, .ids = (size_t[]){1,2}},
       {.count = 2, .ids = (size_t[]){3,4}}};
 
    // Check against reference data
    for (size_t i = 0, offset = 0; i < 2; ++i) {
 
      if (num_local_point_ids[i] != ref_results[i].count){
        PUT_ERR("utest_bnd_circle_search: wrong number of points");
      }
 
      int match_flag[ref_results[i].count];
      size_t match_count = 0;
      memset(match_flag, 0, ref_results[i].count * sizeof(match_flag[0]));
 
      size_t * curr_local_point_ids = local_point_ids + offset;
      for (size_t j = 0; j < num_local_point_ids[i]; ++j) {
        for (size_t k = 0; k < ref_results[i].count; ++k) {
          if (curr_local_point_ids[j] == ref_results[i].ids[k]) {
 
            if (match_flag[k]) {
              PUT_ERR("utest_bnd_circle_search: multiple matches");
            } else {
              match_count++;
            }
          }
        }
      }
 
      if (match_count != ref_results[i].count) {
        PUT_ERR("utest_bnd_circle_search: wrong points");
      }
 
      offset += num_local_point_ids[i];
    }
 
    free(local_point_ids);
    yac_delete_point_sphere_part_search(search);
  }
 
  // Test 2: empty search (null search)
  {
    struct point_sphere_part_search * search =
      yac_point_sphere_part_search_new(0, NULL, NULL);
 
    struct bounding_circle bnd_circle;
    LLtoXYZ_deg(0.0, 0.0, bnd_circle.base_vector);
    double inc_angle = 10.0 * YAC_RAD;
    bnd_circle.inc_angle = sin_cos_angle_new(sin(inc_angle), cos(inc_angle));
    bnd_circle.sq_crd = DBL_MAX;
 
    size_t * local_point_ids = NULL;
    size_t local_point_ids_array_size = 0;
    size_t num_local_point_ids[1] = {(size_t)-1};
 
    yac_point_sphere_part_search_bnd_circle(
      search, 1, &bnd_circle,
      &local_point_ids, &local_point_ids_array_size, num_local_point_ids);
 
    if (num_local_point_ids[0] != 0)
      PUT_ERR("utest_bnd_circle_search: empty search should find 0 points");
 
    free(local_point_ids);
    yac_delete_point_sphere_part_search(search);
  }
 
  // Test 3: zero bounding circles
  {
    double xyz_coordinates[3][3] = {{1,0,0}, {0,1,0}, {0,0,1}};
    yac_int global_ids[3] = {0, 1, 2};
 
    struct point_sphere_part_search * search =
      yac_point_sphere_part_search_new(
        3, (yac_const_coordinate_pointer)xyz_coordinates, global_ids);
 
    size_t * local_point_ids = NULL;
    size_t local_point_ids_array_size = 0;
 
    // Should handle zero circles gracefully
    yac_point_sphere_part_search_bnd_circle(
      search, 0, NULL,
      &local_point_ids, &local_point_ids_array_size, NULL);
 
    free(local_point_ids);
    yac_delete_point_sphere_part_search(search);
  }
 
  // Test 4: single point within bounding circle
  {
    double xyz_coordinate[1][3];
    LLtoXYZ_deg(45.0, 30.0, xyz_coordinate[0]);
    yac_int global_id[1] = {0};
 
    struct point_sphere_part_search * search =
      yac_point_sphere_part_search_new(
        1, (yac_const_coordinate_pointer)xyz_coordinate, global_id);
 
    // Bounding circle centered at same location
    struct bounding_circle bnd_circle;
    memcpy(bnd_circle.base_vector, xyz_coordinate[0], sizeof(xyz_coordinate[0]));
    double inc_angle = 1.0 * YAC_RAD;
    bnd_circle.inc_angle = sin_cos_angle_new(sin(inc_angle), cos(inc_angle));
    bnd_circle.sq_crd = DBL_MAX;
 
    size_t * local_point_ids = NULL;
    size_t local_point_ids_array_size = 0;
    size_t num_local_point_ids[1];
 
    yac_point_sphere_part_search_bnd_circle(
      search, 1, &bnd_circle,
      &local_point_ids, &local_point_ids_array_size, num_local_point_ids);
 
    if (num_local_point_ids[0] != 1)
      PUT_ERR("utest_bnd_circle_search: single point search failed (count)");
    if (local_point_ids[0] != 0)
      PUT_ERR("utest_bnd_circle_search: single point search failed (id)");
 
    free(local_point_ids);
    yac_delete_point_sphere_part_search(search);
  }
 
  // Test 5: no points within bounding circle
  {
    double x_coordinates[] = {0.0, 90.0, 180.0, -90.0};
    double y_coordinates[] = {0.0, 0.0, 0.0, 0.0};
    enum {NUM_POINTS = sizeof(x_coordinates) / sizeof(x_coordinates[0])};
    double xyz_coordinates[NUM_POINTS][3];
    yac_int global_ids[NUM_POINTS];
    for (size_t i = 0; i < NUM_POINTS; ++i) {
      LLtoXYZ_deg(x_coordinates[i], y_coordinates[i], xyz_coordinates[i]);
      global_ids[i] = (yac_int)i;
    }
 
    struct point_sphere_part_search * search =
      yac_point_sphere_part_search_new(
        NUM_POINTS, (yac_const_coordinate_pointer)xyz_coordinates, global_ids);
 
    // Bounding circle at north pole with small angle
    struct bounding_circle bnd_circle;
    LLtoXYZ_deg(0.0, 90.0, bnd_circle.base_vector);
    double inc_angle = 0.5 * YAC_RAD; // Very small, won't contain any point
    bnd_circle.inc_angle = sin_cos_angle_new(sin(inc_angle), cos(inc_angle));
    bnd_circle.sq_crd = DBL_MAX;
 
    size_t * local_point_ids = NULL;
    size_t local_point_ids_array_size = 0;
    size_t num_local_point_ids[1];
 
    yac_point_sphere_part_search_bnd_circle(
      search, 1, &bnd_circle,
      &local_point_ids, &local_point_ids_array_size, num_local_point_ids);
 
    if (num_local_point_ids[0] != 0)
      PUT_ERR("utest_bnd_circle_search: should find 0 points");
 
    free(local_point_ids);
    yac_delete_point_sphere_part_search(search);
  }
 
  // Test 6: large bounding circle containing all points
  {
    double x_coordinates[] = {0.0, 1.0, 2.0, 3.0, 4.0};
    double y_coordinates[] = {0.0, 0.0, 0.0, 0.0, 0.0};
    enum {NUM_POINTS = sizeof(x_coordinates) / sizeof(x_coordinates[0])};
    double xyz_coordinates[NUM_POINTS][3];
    yac_int global_ids[NUM_POINTS];
    for (size_t i = 0; i < NUM_POINTS; ++i) {
      LLtoXYZ_deg(x_coordinates[i], y_coordinates[i], xyz_coordinates[i]);
      global_ids[i] = (yac_int)i;
    }
 
    struct point_sphere_part_search * search =
      yac_point_sphere_part_search_new(
        NUM_POINTS, (yac_const_coordinate_pointer)xyz_coordinates, global_ids);
 
    // Large bounding circle centered at (2, 0) containing all points
    struct bounding_circle bnd_circle;
    LLtoXYZ_deg(2.0, 0.0, bnd_circle.base_vector);
    double inc_angle = 10.0 * YAC_RAD; // Large enough to contain all
    bnd_circle.inc_angle = sin_cos_angle_new(sin(inc_angle), cos(inc_angle));
    bnd_circle.sq_crd = DBL_MAX;
 
    size_t * local_point_ids = NULL;
    size_t local_point_ids_array_size = 0;
    size_t num_local_point_ids[1];
 
    yac_point_sphere_part_search_bnd_circle(
      search, 1, &bnd_circle,
      &local_point_ids, &local_point_ids_array_size, num_local_point_ids);
 
    if (num_local_point_ids[0] != NUM_POINTS)
      PUT_ERR("utest_bnd_circle_search: should find all points");
 
    free(local_point_ids);
    yac_delete_point_sphere_part_search(search);
  }
}
 
static void utest_bnd_circle_random_test(
  size_t num_points, size_t num_circles) {
 
  double (*xyz_coordinates)[3] = xmalloc(num_points * sizeof(*xyz_coordinates));
  yac_int * global_ids = xmalloc(num_points * sizeof(*global_ids));
  struct bounding_circle * bnd_circles =
    xmalloc(num_circles * sizeof(*bnd_circles));
 
  // Generate random points
  for (size_t i = 0; i < num_points; ++i) {
    utest_gen_random_point(xyz_coordinates[i]);
    global_ids[i] = (yac_int)i;
  }
 
  // Generate random bounding circles
  for (size_t i = 0; i < num_circles; ++i)
    utest_gen_random_bnd_circle(&bnd_circles[i]);
 
  struct point_sphere_part_search * search =
    yac_point_sphere_part_search_new(
      num_points, (yac_const_coordinate_pointer)xyz_coordinates, global_ids);
 
  size_t * local_point_ids = NULL;
  size_t local_point_ids_array_size = 0;
  size_t * num_local_point_ids = xmalloc(num_circles * sizeof(*num_local_point_ids));
 
  yac_point_sphere_part_search_bnd_circle(
    search, num_circles, bnd_circles,
    &local_point_ids, &local_point_ids_array_size, num_local_point_ids);
 
  // Verify results by brute-force checking
  size_t offset = 0;
  for (size_t i = 0; i < num_circles; ++i) {
 
    // Compute reference: find all points within bounding circle
    size_t ref_count = 0;
    size_t * ref_ids = xmalloc(num_points * sizeof(*ref_ids));
 
    for (size_t j = 0; j < num_points; ++j) {
      struct sin_cos_angle angle_to_point =
        get_vector_angle_2(bnd_circles[i].base_vector, xyz_coordinates[j]);
      if (compare_angles(angle_to_point, bnd_circles[i].inc_angle) <= 0) {
        ref_ids[ref_count++] = j;
      }
    }
 
    // Check count matches
    if (num_local_point_ids[i] != ref_count) {
      PUT_ERR("utest_bnd_circle_random_test: wrong number of points found");
    } else {
      // Check all returned points are in reference set
      for (size_t j = 0; j < num_local_point_ids[i]; ++j) {
        int found = 0;
        for (size_t k = 0; k < ref_count; ++k) {
          if (local_point_ids[offset + j] == ref_ids[k]) {
            found = 1;
            break;
          }
        }
        if (!found)
          PUT_ERR("utest_bnd_circle_random_test: unexpected point in result");
      }
 
      // Check all reference points are in returned set
      for (size_t k = 0; k < ref_count; ++k) {
        int found = 0;
        for (size_t j = 0; j < num_local_point_ids[i]; ++j) {
          if (local_point_ids[offset + j] == ref_ids[k]) {
            found = 1;
            break;
          }
        }
        if (!found)
          PUT_ERR("utest_bnd_circle_random_test: missing point in result");
      }
    }
 
    offset += num_local_point_ids[i];
    free(ref_ids);
  }
 
  free(local_point_ids);
  free(num_local_point_ids);
  free(bnd_circles);
  free(global_ids);
  free(xyz_coordinates);
  yac_delete_point_sphere_part_search(search);
}