YetAnotherCoupler 3.2.0
Loading...
Searching...
No Matches
test_point_sphere_part.c

This contains a test of the proc_sphere_part grid search algorithm.

// Copyright (c) 2024 The YAC Authors
//
// SPDX-License-Identifier: BSD-3-Clause
#include <stdlib.h>
#include <stdio.h>
#include <float.h>
#include "tests.h"
#include "test_common.h"
#include "basic_grid.h"
#include "sphere_part.h"
#include "geometry.h"
static void random_test(size_t const num_points);
static void zero_point_test();
static void single_point_test();
static void empty_branch_test(size_t num_points);
static void zero_balance_point();
static void duplicated_point_test();
int main(void) {
srand(1365441);
for (size_t i = 0; i < 10; ++i) {
random_test(1);
random_test(10);
random_test(100);
random_test(650);
}
zero_point_test();
single_point_test();
empty_branch_test(4);
empty_branch_test(128);
zero_balance_point();
duplicated_point_test();
return TEST_EXIT_CODE;
}
struct distance_index {
double distance;
size_t index;
};
static int 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 point_sphere_test_NNN(
struct point_sphere_part_search * search, size_t n,
double (*xyz_coordinates)[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) {
search, num_points, xyz_coordinates, n, NULL, 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; ++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 (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])
++num_matching_points;
if (num_matching_points != ref_num_local_points)
PUT_ERR("wrong results\n");
offset += ref_num_local_points;
}
}
static void 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 gen_random_bnd_circle(struct bounding_circle * circle) {
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 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) {
gen_random_bnd_circle(bnd_circles + i);
global_ids[i] = (yac_int)i;
}
for (size_t k = 0; k < 2 * num_points; ++k)
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),
compare_distance_index);
}
struct point_sphere_part_search * search =
global_ids);
for (size_t n = 1; n < 16; ++n)
point_sphere_test_NNN(
search, n, coordinates_xyz_b, num_points, &local_point_ids,
&local_point_ids_array_size, num_local_point_ids, distances);
{ // yac_point_sphere_part_search_bnd_circle_contains_points
for (size_t i = 0; i < num_points; ++i) {
int ref_flag = 0;
for (size_t j = 0; (j < num_points) && !ref_flag; ++j)
ref_flag =
coordinates_xyz_a[j], bnd_circles + i);
if (ref_flag !=
search, bnd_circles[i]))
PUT_ERR(
"ERROR in yac_point_sphere_part_search_bnd_circle_contains_points");
}
}
free(distances);
free(local_point_ids);
free(num_local_point_ids);
free(global_ids);
free(bnd_circles);
free(coordinates_xyz);
}
static void single_point_test() {
struct point_sphere_part_search * search;
{
double x_coordinate = 0, y_coordinate = 0;
double xyz_coordinate[1][3];
yac_int global_id[1] = {0};
LLtoXYZ(x_coordinate, y_coordinate, xyz_coordinate[0]);
search =
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));
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("single_point_test: wrong result\n");
free(local_point_ids);
free(num_local_point_ids);
}
static void zero_point_test() {
struct point_sphere_part_search * search =
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));
for (size_t i = 0; i < num_points; ++i)
num_local_point_ids[i] = (size_t)-1;
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("zero_point_test: wrong result\n");
for (size_t i = 0; i < num_points; ++i)
num_local_point_ids[i] = (size_t)-1;
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("zero_point_test: wrong result\n");
free(local_point_ids);
free(num_local_point_ids);
}
static void 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 =
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),
compare_distance_index);
}
for (size_t i = 0; i < num_points_b; ++i)
num_local_point_ids[i] = (size_t)-1;
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 (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(xyz_coordinates_a);
free(local_point_ids);
free(num_local_point_ids);
}
static void 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 =
6, (yac_const_coordinate_pointer)xyz_coordinate, global_ids);
}
}
static void 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]);
global_ids[i] = (yac_int)i;
}
search =
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};
double y_coordinates[] = {1.0,8.0,3.0,4.0,5.0,1.0,8.0,4.0,0.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]);
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;
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][3] =
{{11}, {18}, {23}, {52}, {25}, {31}, {38}, {52}, {0}};
size_t ref_num_local_point_ids[NUM_POINTS] = {1,1,1,1,1,1,1,1,1};
double ref_x_coordinates[NUM_POINTS] = {1.0,1.0,2.0,2.0,2.0,3.0,3.0,2.0,0.0};
double ref_y_coordinates[NUM_POINTS] = {1.0,8.0,3.0,4.0,5.0,1.0,8.0,4.0,0.0};
double ref_result_coordinates_xyz[NUM_POINTS][3][3];
for (size_t i = 0; i < NUM_POINTS; ++i)
for (size_t j = 0; j < ref_num_local_point_ids[i]; ++j)
ref_x_coordinates[i], ref_y_coordinates[i],
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 {
for (size_t j = 0; j < ref_num_local_point_ids[i]; ++j) {
result_coordinates_xyz[offset+j],
ref_result_coordinates_xyz[i][j]))
PUT_ERR("ERROR in yac_point_sphere_part_search_NN "
"(result_coordinates_xyz)");
if (local_point_ids[offset+j] != ref_local_point_ids[i][j])
PUT_ERR("ERROR in yac_point_sphere_part_search_NN (local_point_ids)");
}
}
}
free(local_point_ids);
free(result_coordinates_xyz);
}
static int points_are_identically(double const *a, double const *b)
Definition geometry.h:648
static int yac_point_in_bounding_circle_vec(double point_vector[3], struct bounding_circle *bnd_circle)
Definition geometry.h:266
static void LLtoXYZ_deg(double lon, double lat, double p_out[])
Definition geometry.h:304
static void LLtoXYZ(double lon, double lat, double p_out[])
Definition geometry.h:287
static double get_vector_angle(double const a[3], double const b[3])
Definition geometry.h:364
static struct sin_cos_angle sin_cos_angle_new(double sin, double cos)
Definition geometry.h:405
#define xmalloc(size)
Definition ppm_xfuncs.h:66
void yac_delete_point_sphere_part_search(struct point_sphere_part_search *search)
struct point_sphere_part_search * yac_point_sphere_part_search_new(size_t num_points, yac_const_coordinate_pointer coordinates_xyz, yac_int const *ids)
int yac_point_sphere_part_search_bnd_circle_contains_points(struct point_sphere_part_search *search, struct bounding_circle circle)
void yac_point_sphere_part_search_NN(struct point_sphere_part_search *search, size_t num_points, double(*coordinates_xyz)[3], double *cos_angles, double(**result_coordinates_xyz)[3], size_t *result_coordinates_xyz_array_size, size_t **local_point_ids, size_t *local_point_ids_array_size, size_t *num_local_point_ids)
void yac_point_sphere_part_search_NNN(struct point_sphere_part_search *search, size_t num_points, double(*coordinates_xyz)[3], size_t n, double **cos_angles, size_t *cos_angles_array_size, double(**result_coordinates_xyz)[3], size_t *result_coordinates_xyz_array_size, size_t **local_point_ids, size_t *local_point_ids_array_size, size_t *num_local_point_ids)
algorithm for searching cells and points on a grid
struct sin_cos_angle inc_angle
angle between the middle point and the boundary of the spherical cap
Definition geometry.h:61
double base_vector[3]
Definition geometry.h:59
double sq_crd
Definition geometry.h:64
int main(int argc, char **argv)
static size_t num_points
Definition yac.c:121
Xt_int yac_int
Definition yac_types.h:15
double const (* yac_const_coordinate_pointer)[3]
Definition yac_types.h:20
double(* yac_coordinate_pointer)[3]
Definition yac_types.h:19