From 35c3ff9a0590dadc27216dd04e06cd3d028d5b1d Mon Sep 17 00:00:00 2001 From: Cleveland Date: Thu, 26 Aug 2021 11:32:02 -0600 Subject: [PATCH 1/4] rely on quick_index for spherical reconstruction --- src/kde/kde.cc | 135 +--- src/kde/kde.hh | 33 +- src/kde/kde.i.hh | 67 +- src/kde/quick_index.cc | 814 ++++---------------- src/kde/quick_index.hh | 73 +- src/kde/test/tstkde.cc | 234 ++++-- src/kde/test/tstquick_index.cc | 1267 ++++++++++++++++++++++++-------- 7 files changed, 1370 insertions(+), 1253 deletions(-) diff --git a/src/kde/kde.cc b/src/kde/kde.cc index 1cd3b80ac..df5963d5d 100644 --- a/src/kde/kde.cc +++ b/src/kde/kde.cc @@ -15,7 +15,7 @@ * of the spatial domain. Other approaches that could be considered are quadrature based * approaches that fully sample the Kernel space reducing the need for the normalization. * - * \note Copyright (C) 2020-2021 Triad National Security, LLC., All rights reserved. + * \note Copyright (C) 2021-2021 Triad National Security, LLC., All rights reserved. */ //------------------------------------------------------------------------------------------------// @@ -28,9 +28,9 @@ namespace rtt_kde { //------------------------------------------------------------------------------------------------// /*! - * \brief Calculate Cartesian Weight + * \brief Calculate Weight * - * \pre Calculate the effective weight in Cartesian geometry from a given location to the current + * \pre Calculate the effective weight in Cartesian and Spherical reconstructions from a given location to the current * kernel * * \param[in] r0 current kernel center location @@ -45,21 +45,33 @@ namespace rtt_kde { * * \post the local reconstruction of the original data is returned. */ -double kde::calc_cartesian_weight(const std::array &r0, - const std::array &one_over_h0, - const std::array &r, - const std::array &one_over_h, - const quick_index &qindex, - const double &discontinuity_cutoff) const { +double kde::calc_weight(const std::array &r0, const std::array &one_over_h0, + const std::array &r, const std::array &one_over_h, + const quick_index &qindex, const double &discontinuity_cutoff) const { Require(one_over_h0[0] > 0.0); Require(qindex.dim > 1 ? one_over_h0[1] > 0.0 : true); Require(qindex.dim > 2 ? one_over_h0[2] > 0.0 : true); Require(one_over_h[0] > 0.0); Require(qindex.dim > 1 ? one_over_h[1] > 0.0 : true); Require(qindex.dim > 2 ? one_over_h[2] > 0.0 : true); + // do not allow spherical reflection of the radial direction + Require(qindex.spherical ? !reflect_boundary[0] : true); + Require(qindex.spherical ? !reflect_boundary[1] : true); double weight = 1.0; + // set the radius at which all arch_lengths should be computed + // this value is ignored by calc_orthogonal_distance in non-spherical geometry + double arch_radius = r0[0]; + std::array distance = qindex.calc_orthogonal_distance(r0, r, arch_radius); + std::array low_reflect_r0_distance = + qindex.calc_orthogonal_distance(qindex.bounding_box_min, r0, arch_radius); + std::array low_reflect_r_distance = + qindex.calc_orthogonal_distance(qindex.bounding_box_min, r, arch_radius); + std::array high_reflect_r0_distance = + qindex.calc_orthogonal_distance(r0, qindex.bounding_box_max, arch_radius); + std::array high_reflect_r_distance = + qindex.calc_orthogonal_distance(r, qindex.bounding_box_max, arch_radius); for (size_t d = 0; d < qindex.dim; d++) { - const double u = (r0[d] - r[d]) * one_over_h0[d]; + const double u = distance[d] * one_over_h0[d]; const double scale = fabs(one_over_h0[d] - one_over_h[d]) / std::max(one_over_h0[d], one_over_h[d]) > discontinuity_cutoff @@ -71,14 +83,12 @@ double kde::calc_cartesian_weight(const std::array &r0, const bool high_reflect = reflect_boundary[d * 2 + 1]; if (low_reflect) { const double low_u = - ((r0[d] - qindex.bounding_box_min[d]) + (r[d] - qindex.bounding_box_min[d])) * - one_over_h0[d]; + (low_reflect_r0_distance[d] + low_reflect_r_distance[d]) * one_over_h0[d]; bc_weight += epan_kernel(low_u); } if (high_reflect) { const double high_u = - ((qindex.bounding_box_max[d] - r0[d]) + (qindex.bounding_box_max[d] - r[d])) * - one_over_h0[d]; + (high_reflect_r0_distance[d] + high_reflect_r_distance[d]) * one_over_h0[d]; bc_weight += epan_kernel(high_u); } weight *= scale * bc_weight * epan_kernel(u) * one_over_h0[d]; @@ -87,78 +97,6 @@ double kde::calc_cartesian_weight(const std::array &r0, return weight; } -//------------------------------------------------------------------------------------------------// -/*! - * \brief calculate spherical weight - * - * \pre Calculate the effective weight from a given location to the current kernel - * - * \param[in] r0 current kernel center location - * \param[in] one_over_h0 current kernel width - * \param[in] r data location - * \param[in] one_over_h kernel width at this data location - * \param[in] qindex quick indexing class - * \param[in] discontinuity_cutoff maximum size of value discrepancies to include in the - * reconstruction - * - * \return weight contribution to the current kernel - * - * \post the local reconstruction of the original data is returned. - */ -double kde::calc_spherical_weight(const std::array &r0, - const std::array &one_over_h0, - const std::array &r, - const std::array &one_over_h, - const quick_index &qindex, - const double &discontinuity_cutoff) const { - Require(one_over_h0[0] > 0.0); - Require(qindex.dim > 1 ? one_over_h0[1] > 0.0 : true); - Require(qindex.dim > 2 ? one_over_h0[2] > 0.0 : true); - Require(one_over_h[0] > 0.0); - Require(qindex.dim > 1 ? one_over_h[1] > 0.0 : true); - Require(qindex.dim > 2 ? one_over_h[2] > 0.0 : true); - - // largest active smoothing length - const auto r0_theta_phi = qindex.transform_r_theta(sphere_center, r0); - // if we are near the origin of the sphere, fall back to xyz reconstruction - if (r0_theta_phi[0] < sphere_min_radius || r0_theta_phi[0] > sphere_max_radius) - return calc_cartesian_weight(r0, one_over_h0, r, one_over_h, qindex, discontinuity_cutoff); - - const auto r_theta_phi = qindex.transform_r_theta(sphere_center, r); - const double radius = r0_theta_phi[0]; - double weight = 1.0; - for (size_t d = 0; d < qindex.dim; d++) { - const double arch_scale = d > 0 ? radius : 1.0; - const double u = (r0_theta_phi[d] - r_theta_phi[d]) * arch_scale * one_over_h0[d]; - const double scale = - fabs(one_over_h0[d] - one_over_h[d]) / std::max(one_over_h0[d], one_over_h[d]) > - discontinuity_cutoff - ? 0.0 - : 1.0; - // Apply Boundary Condition Weighting - double bc_weight = 1.0; - /* BC are a little tricky so I will implement them later - const bool low_reflect = reflect_boundary[d * 2]; - const bool high_reflect = reflect_boundary[d * 2 + 1]; - if (low_reflect) { - const double low_u = - ((r0_theta_phi[d] - qindex.bounding_box_min[d]) + (r[d] - qindex.bounding_box_min[d])) * - one_over_h0[d]; - bc_weight += epan_kernel(low_u); - } - if (high_reflect) { - const double high_u = - ((qindex.bounding_box_max[d] - r0[d]) + (qindex.bounding_box_max[d] - r[d])) * - one_over_h0[d]; - bc_weight += epan_kernel(high_u); - } - */ - weight *= scale * bc_weight * epan_kernel(u) * one_over_h0[d]; - } - Ensure(!(weight < 0.0)); - return weight; -} - //------------------------------------------------------------------------------------------------// /*! * \brief KDE reconstruction @@ -465,28 +403,9 @@ void kde::calc_win_min_max(const quick_index &qindex, const std::array 0.0); Require(dim > 1 ? one_over_bandwidth[1] > 0.0 : true); Require(dim > 2 ? one_over_bandwidth[2] > 0.0 : true); - if (use_spherical_reconstruction) { - const double dr = 1.0 / one_over_bandwidth[0]; - const double rmax = sqrt((sphere_center[0] - position[0]) * (sphere_center[0] - position[0]) + - (sphere_center[1] - position[1]) * (sphere_center[1] - position[1])) + - dr; - Check(rmax > 0.0); - const double dtheta = - std::min(1.0 / (one_over_bandwidth[1] * rmax), rtt_units::PI / 2.0 - 1e-12); - if (!(rmax < sphere_min_radius || rmax > sphere_max_radius)) { - // dtheta = arch_length_max/rmax - qindex.calc_wedge_xy_bounds(position, sphere_center, {dr, dtheta, 0.0}, win_min, win_max); - } else { - for (size_t d = 0; d < dim; d++) { - win_min[d] = position[d] - 1.0 / one_over_bandwidth[d]; - win_max[d] = position[d] + 1.0 / one_over_bandwidth[d]; - } - } - } else { - for (size_t d = 0; d < dim; d++) { - win_min[d] = position[d] - 1.0 / one_over_bandwidth[d]; - win_max[d] = position[d] + 1.0 / one_over_bandwidth[d]; - } + for (size_t d = 0; d < dim; d++) { + win_min[d] = position[d] - 1.0 / one_over_bandwidth[d]; + win_max[d] = position[d] + 1.0 / one_over_bandwidth[d]; } } diff --git a/src/kde/kde.hh b/src/kde/kde.hh index 24db717fd..dfef2274c 100644 --- a/src/kde/kde.hh +++ b/src/kde/kde.hh @@ -3,7 +3,7 @@ * \file kde/kde.hh * \author Mathew Cleveland * \brief Define class kernel density estimator class - * \note Copyright (C) 2020-2021 Triad National Security, LLC., All rights reserved. + * \note Copyright (C) 2021-2021 Triad National Security, LLC., All rights reserved. */ //------------------------------------------------------------------------------------------------// @@ -60,25 +60,6 @@ public: //! Move the solution back from log space inline double log_inv_transform(const double log_value, const double bias) const; - //! Calculate radius - inline double calc_radius(const std::array ¢er, - const std::array &location) const; - - //! Calculate Arch Length between two locations at a specific radius - inline double calc_arch_length(const std::array ¢er, const double radius, - const std::array &location_1, - const std::array &location_2) const; - //! Setup a spherical reconstruction with a - void set_sphere_center(const std::array &sph_center, const double min_radius, - const double max_radius) { - Insist(max_radius > min_radius, "Spherical KDE max radius must be larger then min radius"); - - use_spherical_reconstruction = true; - sphere_center = sph_center; - sphere_min_radius = min_radius; - sphere_max_radius = max_radius; - } - protected: // IMPLEMENTATION @@ -86,12 +67,7 @@ private: //! Private function to calculate kernel weight double calc_weight(const std::array &r0, const std::array &one_over_h0, const std::array &r, const std::array &one_over_h, - const quick_index &qindex, const double &discontinuity_cutoff) const { - return use_spherical_reconstruction - ? calc_spherical_weight(r0, one_over_h0, r, one_over_h, qindex, discontinuity_cutoff) - : calc_cartesian_weight(r0, one_over_h0, r, one_over_h, qindex, - discontinuity_cutoff); - } + const quick_index &qindex, const double &discontinuity_cutoff) const; void calc_win_min_max(const quick_index &qindex, const std::array &position, const std::array &one_over_bandwidth, std::array &min, @@ -112,11 +88,6 @@ private: // DATA //! reflecting boundary conditions [lower_x, upper_x, lower_y, upper_y, lower_z, upper_z] const std::array reflect_boundary; - //! Spherical Mesh Reconstruction Data - std::array sphere_center{0.0, 0.0, 0.0}; - double sphere_min_radius{0.0}; - double sphere_max_radius{0.0}; - bool use_spherical_reconstruction{false}; }; } // end namespace rtt_kde diff --git a/src/kde/kde.i.hh b/src/kde/kde.i.hh index a6c39a62e..ab9347304 100644 --- a/src/kde/kde.i.hh +++ b/src/kde/kde.i.hh @@ -4,7 +4,7 @@ * \author Mathew Cleveland * \date Nov. 10th 2020 * \brief Member definitions of class kde - * \note Copyright (C) 2021 Triad National Security, LLC., All rights reserved. + * \note Copyright (C) 2021-2021 Triad National Security, LLC., All rights reserved. */ //------------------------------------------------------------------------------------------------// @@ -71,71 +71,6 @@ inline double kde::log_inv_transform(const double log_value, const double bias) return exp(log_value) - bias; } -//! Lambda to calculate a vector -const auto calc_vec = [](const auto &v1, const auto &v2) { - Require(v1.size() == 3); - Require(v2.size() == 3); - return std::array{v2[0] - v1[0], v2[1] - v1[1], v2[2] - v1[2]}; -}; - -//! Lambda to calculate vector magnitude -const auto calc_mag = [](const auto &v) { - Require(v.size() == 3); - return sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]); -}; - -//! Lambda to calculate unit vector -const auto calc_unit_vec = [](const auto &v) { - Require(v.size() == 3); - const double mag = calc_mag(v); - return std::array{v[0] / mag, v[1] / mag, v[2] / mag}; -}; - -//------------------------------------------------------------------------------------------------// -/*! - * \brief - * Calculate the radius given a sphere center and the current location. - * - * - * \param[in] center the center location (x,y,z) or (r,z) of the sphere - * \param[in] location data location (x,y,z) or (r,z) - * - * \return radius from cell center - * - * Test of kde. - */ - -inline double kde::calc_radius(const std::array ¢er, - const std::array &location) const { - return calc_mag(calc_vec(center, location)); -} - -//------------------------------------------------------------------------------------------------// -/*! - * \brief Calculate the arch length between two points - * - * Calculate the arch length between two points (infinitely extended from sphere center) at a - * specified radius. - * - * - * \param[in] center the center location (x,y,z) or (r,z) of the sphere - * \param[in] radius from sphere center to calculate the arch length - * \param[in] location_1 data location (x,y,z) or (r,z) - * \param[in] location_2 data location (x,y,z) or (r,z) - * - * \return arch length - * - * Test of kde. - */ - -inline double kde::calc_arch_length(const std::array ¢er, const double radius, - const std::array &location_1, - const std::array &location_2) const { - const std::array v1{calc_unit_vec(calc_vec(center, location_1))}; - const std::array v2{calc_unit_vec(calc_vec(center, location_2))}; - const double cos_theta = v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2]; - return radius * acos(cos_theta); -} } // end namespace rtt_kde #endif // kde_kde_i_hh diff --git a/src/kde/quick_index.cc b/src/kde/quick_index.cc index 2c627a546..7c6b6f28f 100644 --- a/src/kde/quick_index.cc +++ b/src/kde/quick_index.cc @@ -3,17 +3,17 @@ * \file kde/quick_index.cc * \author Mathew Cleveland * \brief Explicitly defined quick_index functions. - * \note Copyright (C) 2021 Triad National Security, LLC., All rights reserved. + * \note Copyright (C) 2021-2021 Triad National Security, LLC., All rights reserved. */ //------------------------------------------------------------------------------------------------// #include "quick_index.hh" #include "ds++/dbc.hh" -#include #include #include namespace rtt_kde { + //------------------------------------------------------------------------------------------------// /*! * \brief quick_index constructor. @@ -33,9 +33,13 @@ namespace rtt_kde { */ quick_index::quick_index(const size_t dim_, const std::vector> &locations_, const double max_window_size_, const size_t bins_per_dimension_, - const bool domain_decomposed_) - : dim(dim_), domain_decomposed(domain_decomposed_), coarse_bin_resolution(bins_per_dimension_), - max_window_size(max_window_size_), locations(locations_), n_locations(locations_.size()) { + const bool domain_decomposed_, const bool spherical_, + const std::array &sphere_center_) + : dim(dim_), domain_decomposed(domain_decomposed_), spherical(spherical_), + sphere_center(sphere_center_), coarse_bin_resolution(bins_per_dimension_), + max_window_size(max_window_size_), + locations(spherical ? transform_spherical(dim_, sphere_center_, locations_) : locations_), + n_locations(locations_.size()) { Require(dim > 0); Require(coarse_bin_resolution > 0); @@ -62,15 +66,27 @@ quick_index::quick_index(const size_t dim_, const std::vector &window_min, for (size_t d = 0; d < dim; d++) { // because local bounds can extend beyond the mesh we need to force a // positive index if necessary + double wmin = window_min[d]; + if (spherical && d == 1 && window_min[d] < bounding_box_min[d]) + wmin = bounding_box_min[d]; // truncate to standard theta space + double wmax = window_max[d]; + if (spherical && d == 1 && window_max[d] > bounding_box_max[d]) + wmax = bounding_box_max[d]; // truncate to standard theta space index_min[d] = static_cast(std::floor(std::max( - crd * (window_min[d] - bounding_box_min[d]) / (bounding_box_max[d] - bounding_box_min[d]), - 0.0))); - index_max[d] = static_cast(std::floor(crd * (window_max[d] - bounding_box_min[d]) / + crd * (wmin - bounding_box_min[d]) / (bounding_box_max[d] - bounding_box_min[d]), 0.0))); + index_max[d] = static_cast(std::floor(crd * (wmax - bounding_box_min[d]) / (bounding_box_max[d] - bounding_box_min[d]))); // because local bounds can extend beyond the mesh we need to floor to // the max bin size @@ -436,22 +457,84 @@ quick_index::window_coarse_index_list(const std::array &window_min, } } } + + // Fill in the overflow around theta=0.0 + if (spherical && (window_min[1] < 0.0 || window_max[1] > 2.0 * rtt_units::PI)) { + // Only one bound of the window should every overshoot zero + Check(!(window_min[1] < 0.0 && window_max[1] > 2.0 * rtt_units::PI)); + size_t overlap_nbins = 1; + for (size_t d = 0; d < dim; d++) { + // because local bounds can extend beyond the mesh we need to force a + // positive index if necessary + double wmin = window_min[d]; + double wmax = window_max[d]; + if (spherical && d == 1 && window_min[d] < 0.0) { + // Capture the overshoot theta space + wmin = std::min(2.0 * rtt_units::PI - window_min[d], bounding_box_max[d]); + wmax = 2.0 * rtt_units::PI; + } + if (spherical && d == 1 && window_max[d] > 2.0 * rtt_units::PI) { + // Capture the overshoot theta space + wmin = 0.0; + wmax = std::max(2.0 * rtt_units::PI - window_max[d], bounding_box_min[d]); + } + // Truncate based on global space + if (spherical && d == 1 && wmin < bounding_box_min[d]) + wmin = bounding_box_min[d]; // truncate to standard theta space + if (spherical && d == 1 && wmax > bounding_box_max[d]) + wmax = bounding_box_max[d]; //trunacte to standard theta space + index_min[d] = static_cast(std::floor(std::max( + crd * (wmin - bounding_box_min[d]) / (bounding_box_max[d] - bounding_box_min[d]), 0.0))); + index_max[d] = static_cast(std::floor(std::max( + crd * (wmax - bounding_box_min[d]) / (bounding_box_max[d] - bounding_box_min[d]), 0.0))); + // because local bounds can extend beyond the mesh we need to floor to + // the max bin size + index_min[d] = std::min(index_min[d], coarse_bin_resolution - 1); + index_max[d] = std::min(index_max[d], coarse_bin_resolution - 1); + + // Use multiplicity to accumulate total bins; + if ((index_max[d] - index_min[d]) > 0) + overlap_nbins *= index_max[d] - index_min[d] + 1; + } + for (size_t k = index_min[2]; k <= index_max[2]; k++) { + for (size_t j = index_min[1]; j <= index_max[1]; j++) { + for (size_t i = index_min[0]; i <= index_max[0]; i++) { + size_t bin_index = + i + j * coarse_bin_resolution + k * coarse_bin_resolution * coarse_bin_resolution; + // make sure we don't duplicate a bin here + if (std::find(bin_list.begin(), bin_list.end(), bin_index) == bin_list.end()) { + bin_list.push_back(bin_index); + } + } + } + } + } + return bin_list; } //------------------------------------------------------------------------------------------------// // Lambda for getting the mapped window bin -auto get_window_bin = [](const auto dim, const auto &grid_bins, const auto &location, - const auto &window_min, const auto &window_max, +auto get_window_bin = [](auto spherical, const auto dim, const auto &grid_bins, + const auto &location, const auto &window_min, const auto &window_max, const auto &Remember(n_map_bins)) { // calculate local bin index bool valid = true; std::array bin_id{0, 0, 0}; - std::array bin_center{0, 0, 0}; + double distance_to_bin_center = 0.0; for (size_t d = 0; d < dim; d++) { Check((window_max[d] - window_min[d]) > 0.0); - const double bin_value = static_cast(grid_bins[d]) * (location[d] - window_min[d]) / - (window_max[d] - window_min[d]); + double loc = location[d]; + // transform location for zero theta overshoot + if (spherical && d == 1 && window_max[d] > 2.0 * rtt_units::PI && + location[d] < window_max[d] - 2.0 * rtt_units::PI) + loc += 2.0 * rtt_units::PI; + // transform location for zero theta overshoot + if (spherical && d == 1 && window_min[d] < 0 && + location[d] > 2.0 * rtt_units::PI - window_min[d]) + loc -= 2.0 * rtt_units::PI; + const double bin_value = + static_cast(grid_bins[d]) * (loc - window_min[d]) / (window_max[d] - window_min[d]); if (bin_value < 0.0 || bin_value > static_cast(grid_bins[d])) { valid = false; break; @@ -459,93 +542,36 @@ auto get_window_bin = [](const auto dim, const auto &grid_bins, const auto &loca bin_id[d] = static_cast(bin_value); // catch any values exactly on the edge of the top bin bin_id[d] = std::min(grid_bins[d] - 1, bin_id[d]); - bin_center[d] = + const double bin_center = window_min[d] + (static_cast(bin_id[d]) / static_cast(grid_bins[d]) + 0.5 / static_cast(grid_bins[d])) * (window_max[d] - window_min[d]); + // approximate in spherical geometry; + distance_to_bin_center += (bin_center - loc) * (bin_center - loc); } } + distance_to_bin_center = + rtt_dsxx::soft_equiv(distance_to_bin_center, 0.0) ? 0.0 : sqrt(distance_to_bin_center); const size_t local_window_bin = bin_id[0] + bin_id[1] * grid_bins[0] + bin_id[2] * grid_bins[0] * grid_bins[1]; Check(valid ? local_window_bin < n_map_bins : true); - return std::tuple>{valid, local_window_bin, bin_center}; -}; - -//------------------------------------------------------------------------------------------------// -// Lambda for getting the mapped window bin -auto get_sphere_window_bin = [](const auto &grid_bins, const auto &location, const auto &window_min, - const auto &window_max, const auto &Remember(n_map_bins), - const auto pi) { - // calculate local bin index - bool valid = true; - std::array bin_id{0, 0, 0}; - std::array bin_center{0, 0, 0}; - { - Check((window_max[0] - window_min[0]) > 0.0); - const double bin_value = static_cast(grid_bins[0]) * (location[0] - window_min[0]) / - (window_max[0] - window_min[0]); - if (bin_value < 0.0 || bin_value > static_cast(grid_bins[0])) { - valid = false; - } else { - bin_id[0] = static_cast(bin_value); - // catch any values exactly on the edge of the top bin - bin_id[0] = std::min(grid_bins[0] - 1, bin_id[0]); - bin_center[0] = - window_min[0] + (static_cast(bin_id[0]) / static_cast(grid_bins[0]) + - 0.5 / static_cast(grid_bins[0])) * - (window_max[0] - window_min[0]); - } - } - if (valid) { - // catch the window that wraps around the zero theta location - const double theta_location = - (window_max[1] - window_min[1]) > 0.0 - ? location[1] - : location[1] < window_max[1] ? location[1] + 2 * pi : location[1]; - const double theta_max = - (window_max[1] - window_min[1]) > 0.0 ? window_max[1] : 2 * pi + window_max[1]; - Check(!((theta_max - window_min[1]) < 0.0)); - const double bin_value = static_cast(grid_bins[1]) * (theta_location - window_min[1]) / - (theta_max - window_min[1]); - if (bin_value < 0.0 || bin_value > static_cast(grid_bins[1])) { - valid = false; - } else { - bin_id[1] = static_cast(bin_value); - // catch any values exactly on the edge of the top bin - bin_id[1] = std::min(grid_bins[1] - 1, bin_id[1]); - bin_center[1] = - window_min[1] + (static_cast(bin_id[1]) / static_cast(grid_bins[1]) + - 0.5 / static_cast(grid_bins[1])) * - (theta_max - window_min[1]); - bin_center[1] = bin_center[1] < 2.0 * pi ? bin_center[1] : bin_center[1] - 2.0 * pi; - } - } - - const size_t local_window_bin = - bin_id[0] + bin_id[1] * grid_bins[0] + bin_id[2] * grid_bins[0] * grid_bins[1]; - - Check(valid ? local_window_bin < n_map_bins : true); - - return std::tuple>{valid, local_window_bin, bin_center}; + return std::tuple{valid, local_window_bin, distance_to_bin_center}; }; //------------------------------------------------------------------------------------------------// // Lambda for mapping the data auto map_data = [](auto &bias_cell_count, auto &data_count, auto &grid_data, auto &min_distance, - const auto &dim, const auto &map_type, const auto &data, const auto &bin_center, - const auto &location, const auto &local_window_bin, const auto &data_bin) { + const auto &dim, const auto &map_type, const auto &data, + const auto &distance_to_bin_center, const auto &location, + const auto &local_window_bin, const auto &data_bin) { // regardless of map type if it is the first value to enter the bin it // gets set to that value if (data_count[local_window_bin] == 0) { bias_cell_count += 1.0; data_count[local_window_bin]++; - double distance = 0.0; - for (size_t d = 0; d < dim; d++) { - distance += (location[d] - bin_center[d]) * (location[d] - bin_center[d]); - } - min_distance[local_window_bin] = sqrt(distance); + min_distance[local_window_bin] = distance_to_bin_center; grid_data[local_window_bin] = data[data_bin]; } else if (map_type == "max") { if (data[data_bin] > grid_data[local_window_bin]) @@ -557,16 +583,11 @@ auto map_data = [](auto &bias_cell_count, auto &data_count, auto &grid_data, aut data_count[local_window_bin] += 1; grid_data[local_window_bin] += data[data_bin]; } else if (map_type == "nearest") { - double distance = 0.0; - for (size_t d = 0; d < dim; d++) { - distance += (location[d] - bin_center[d]) * (location[d] - bin_center[d]); - } - distance = sqrt(distance); - if (rtt_dsxx::soft_equiv(distance, min_distance[local_window_bin])) { + if (rtt_dsxx::soft_equiv(distance_to_bin_center, min_distance[local_window_bin])) { data_count[local_window_bin] += 1; grid_data[local_window_bin] += data[data_bin]; - } else if (distance < min_distance[local_window_bin]) { - min_distance[local_window_bin] = distance; + } else if (distance_to_bin_center < min_distance[local_window_bin]) { + min_distance[local_window_bin] = distance_to_bin_center; data_count[local_window_bin] = 1; grid_data[local_window_bin] = data[data_bin]; } // else exclude the far points. @@ -607,9 +628,9 @@ void quick_index::map_data_to_grid_window( Require(domain_decomposed ? (fabs(window_max[0] - window_min[0]) - max_window_size) / max_window_size < 1e-6 : true); - Require(domain_decomposed - ? (fabs(window_max[1] - window_min[1]) - max_window_size) / max_window_size < 1e-6 - : true); + Remember(double ymax = spherical ? std::min(rtt_units::PI / 2.0, max_window_size / window_max[0]) + : max_window_size); + Require(domain_decomposed ? (fabs(window_max[1] - window_min[1]) - ymax) / ymax < 1e-6 : true); Require(domain_decomposed ? (fabs(window_max[2] - window_min[2]) - max_window_size) / max_window_size < 1e-6 : true); @@ -674,9 +695,9 @@ void quick_index::map_data_to_grid_window( for (auto &l : mapItr->second) { bool valid; size_t local_window_bin; - std::array bin_center; - std::tie(valid, local_window_bin, bin_center) = - get_window_bin(dim, grid_bins, locations[l], window_min, window_max, n_map_bins); + double distance_to_bin_center; + std::tie(valid, local_window_bin, distance_to_bin_center) = get_window_bin( + spherical, dim, grid_bins, locations[l], window_min, window_max, n_map_bins); // If the bin is outside the window continue to the next poin if (!valid) @@ -684,7 +705,7 @@ void quick_index::map_data_to_grid_window( // lambda for mapping the data map_data(bias_cell_count, data_count, grid_data, min_distance, dim, map_type, local_data, - bin_center, locations[l], local_window_bin, l); + distance_to_bin_center, locations[l], local_window_bin, l); } // end local point loop } // if valid local bin loop @@ -696,9 +717,10 @@ void quick_index::map_data_to_grid_window( for (auto &g : gmapItr->second) { bool valid; size_t local_window_bin; - std::array bin_center; - std::tie(valid, local_window_bin, bin_center) = get_window_bin( - dim, grid_bins, local_ghost_locations[g], window_min, window_max, n_map_bins); + double distance_to_bin_center; + std::tie(valid, local_window_bin, distance_to_bin_center) = + get_window_bin(spherical, dim, grid_bins, local_ghost_locations[g], window_min, + window_max, n_map_bins); // If the bin is outside the window continue to the next poin if (!valid) @@ -706,7 +728,7 @@ void quick_index::map_data_to_grid_window( // lambda for mapping the data map_data(bias_cell_count, data_count, grid_data, min_distance, dim, map_type, ghost_data, - bin_center, local_ghost_locations[g], local_window_bin, g); + distance_to_bin_center, local_ghost_locations[g], local_window_bin, g); } // end ghost point loop } // if valid ghost bin } // if dd @@ -763,17 +785,14 @@ void quick_index::map_data_to_grid_window( // Lambda for mapping the vector data auto map_vector_data = [](auto &bias_cell_count, auto &data_count, auto &grid_data, auto &min_distance, const auto &dim, const auto &map_type, - const auto &data, const auto &bin_center, const auto &location, - const auto &local_window_bin, const auto &data_bin, const auto &vsize) { + const auto &data, const auto &distance_to_bin_center, + const auto &location, const auto &local_window_bin, const auto &data_bin, + const auto &vsize) { // regardless of map type if it is the first value to enter the bin it gets set to that value if (data_count[local_window_bin] == 0) { bias_cell_count += 1.0; data_count[local_window_bin]++; - double distance = 0.0; - for (size_t d = 0; d < dim; d++) { - distance += (location[d] - bin_center[d]) * (location[d] - bin_center[d]); - } - min_distance[local_window_bin] = sqrt(distance); + min_distance[local_window_bin] = distance_to_bin_center; for (size_t v = 0; v < vsize; v++) grid_data[v][local_window_bin] = data[v][data_bin]; } else if (map_type == "max") { @@ -789,17 +808,12 @@ auto map_vector_data = [](auto &bias_cell_count, auto &data_count, auto &grid_da for (size_t v = 0; v < vsize; v++) grid_data[v][local_window_bin] += data[v][data_bin]; } else if (map_type == "nearest") { - double distance = 0.0; - for (size_t d = 0; d < dim; d++) { - distance += (location[d] - bin_center[d]) * (location[d] - bin_center[d]); - } - distance = sqrt(distance); - if (rtt_dsxx::soft_equiv(distance, min_distance[local_window_bin])) { + if (rtt_dsxx::soft_equiv(distance_to_bin_center, min_distance[local_window_bin])) { data_count[local_window_bin] += 1; for (size_t v = 0; v < vsize; v++) grid_data[v][local_window_bin] += data[v][data_bin]; - } else if (distance < min_distance[local_window_bin]) { - min_distance[local_window_bin] = distance; + } else if (distance_to_bin_center < min_distance[local_window_bin]) { + min_distance[local_window_bin] = distance_to_bin_center; data_count[local_window_bin] = 1; for (size_t v = 0; v < vsize; v++) grid_data[v][local_window_bin] = data[v][data_bin]; @@ -845,9 +859,10 @@ void quick_index::map_data_to_grid_window(const std::vector> Require(domain_decomposed ? (fabs(window_max[0] - window_min[0]) - max_window_size) / max_window_size < 1e-6 : true); - Require(domain_decomposed - ? (fabs(window_max[1] - window_min[1]) - max_window_size) / max_window_size < 1e-6 - : true); + Remember(double ymax = spherical ? std::min(rtt_units::PI / 2.0, max_window_size / window_max[0]) + : max_window_size;) + Require(domain_decomposed ? (fabs(window_max[1] - window_min[1]) - ymax) / ymax < 1e-6 + : true); Require(domain_decomposed ? (fabs(window_max[2] - window_min[2]) - max_window_size) / max_window_size < 1e-6 : true); @@ -916,15 +931,16 @@ void quick_index::map_data_to_grid_window(const std::vector> for (auto &l : mapItr->second) { bool valid; size_t local_window_bin; - std::array bin_center; - std::tie(valid, local_window_bin, bin_center) = - get_window_bin(dim, grid_bins, locations[l], window_min, window_max, n_map_bins); + double distance_to_bin_center; + std::tie(valid, local_window_bin, distance_to_bin_center) = get_window_bin( + spherical, dim, grid_bins, locations[l], window_min, window_max, n_map_bins); // If the bin is outside the window continue to the next poin if (!valid) continue; Check(local_window_bin < n_map_bins); map_vector_data(bias_cell_count, data_count, grid_data, min_distance, dim, map_type, - local_data, bin_center, locations[l], local_window_bin, l, vsize); + local_data, distance_to_bin_center, locations[l], local_window_bin, l, + vsize); } // end local point loop } // if valid local bin loop if (domain_decomposed) { @@ -935,16 +951,17 @@ void quick_index::map_data_to_grid_window(const std::vector> for (auto &g : gmapItr->second) { bool valid; size_t local_window_bin; - std::array bin_center; - std::tie(valid, local_window_bin, bin_center) = get_window_bin( - dim, grid_bins, local_ghost_locations[g], window_min, window_max, n_map_bins); + double distance_to_bin_center; + std::tie(valid, local_window_bin, distance_to_bin_center) = + get_window_bin(spherical, dim, grid_bins, local_ghost_locations[g], window_min, + window_max, n_map_bins); // If the bin is outside the window continue to the next poin if (!valid) continue; map_vector_data(bias_cell_count, data_count, grid_data, min_distance, dim, map_type, - ghost_data, bin_center, local_ghost_locations[g], local_window_bin, g, - vsize); + ghost_data, distance_to_bin_center, local_ghost_locations[g], + local_window_bin, g, vsize); } // end ghost point loop } // if valid ghost bin } // if dd @@ -1013,531 +1030,24 @@ void quick_index::map_data_to_grid_window(const std::vector> //------------------------------------------------------------------------------------------------// /*! - * \brief Transform (x, y, z) position to (r, theta, phi) grid - * - * Calculate a relative r theta and phi coordinate relative to a sphere center location from a - * standard (x,y,z) or (r,z) coordinates - * - * \param[in] sphere_center center of sphere in (x,y,z) or (r,z) coordinates - * \param[in] location (x,y,z) or (r,z) location to transform to relative (r, theta, phi) space. - * - * \return relative r theta phi location - */ -std::array quick_index::transform_r_theta(const std::array &sphere_center, - const std::array &location) const { - Insist(dim == 2, "Transform_r_theta Only implemented in 2d"); - const std::array v{location[0] - sphere_center[0], location[1] - sphere_center[1], - 0.0}; - const double r = sqrt(v[0] * v[0] + v[1] * v[1]); - const double mag = sqrt(v[0] * v[0] + v[1] * v[1]); - double cos_theta = mag > 0.0 ? std::max(std::min(v[1] / mag, 1.0), -1.0) : 0.0; - return std::array{ - r, location[0] < sphere_center[0] ? 2.0 * pi - acos(cos_theta) : acos(cos_theta), 0.0}; -} - -//------------------------------------------------------------------------------------------------// -/*! - * \brief Calculate the bounding box of a wedge - * - * Calculates the (x,y,0.0) min and max bounds for a pre-defined wedge [origin (x,y,z), - * wedge_xyz_center (x,y,z), wedge_dr_dtheta (dr, dtheta, 0.0)]. This computes the bounding box for - * the truncated wedge (bounded by rmin and rmax). - * - * win_max - * ---------------(xmax,ymax,0.0) - * | x | __ - * | /|\ | - * | / | \ | dr - * | / * \ | -- - * | / | \ | dr - * win_min | / dt | dt \ |__ - * (xmin,ymin,0.0)--------------- - * - * '*' is the geometric center (not the centroid) - * 'x' is the wedge_origin (or center of the spherical grid) - * - * \param[in] wedge_xyz_center geometric center of the wedge in (x,y,z) or (r,z) coordinates - * \param[in] wedge_origin axisymetric origin (x,y,z) or (r,z) of the wedge. - * \param[in] wedge_dr_dtheta differential size of the wedge in each dimension. - * \param[in,out] win_min differential size of the wedge in each dimension. - * \param[in,out] win_max differential size of the wedge in each dimension. - * - */ -void quick_index::calc_wedge_xy_bounds(const std::array &wedge_xyz_center, - const std::array &wedge_origin, - const std::array &wedge_dr_dtheta, - std::array &win_min, - std::array &win_max) const { - Require(wedge_dr_dtheta[0] > 0.0); - Require(wedge_dr_dtheta[1] > 0.0); - // Some of the checks might not hold for large theta angles - Require(wedge_dr_dtheta[1] < pi / 2.0); - const auto r_theta = transform_r_theta(wedge_origin, wedge_xyz_center); - const double rmin = std::max(0.0, r_theta[0] - wedge_dr_dtheta[0]); - const double rmax = r_theta[0] + wedge_dr_dtheta[0]; - const double dtheta = wedge_dr_dtheta[1]; - const double theta_min = - dtheta < r_theta[1] ? r_theta[1] - dtheta : 2.0 * pi + r_theta[1] - dtheta; - const double theta_max = dtheta + r_theta[1]; - const double cos_theta = cos(r_theta[1]); - const double cos_theta_y_min = r_theta[1] < pi ? cos(theta_max) : cos(theta_min); - const double cos_theta_y_max = r_theta[1] < pi ? cos(theta_min) : cos(theta_max); - const double ymin = theta_max > pi && theta_min < pi - ? wedge_origin[1] - rmax - : cos_theta_y_min < 0.0 ? wedge_origin[1] + rmax * cos_theta_y_min - : wedge_origin[1] + rmin * cos_theta_y_min; - const double ymax = theta_max > 2.0 * pi || theta_min > theta_max - ? wedge_origin[1] + rmax - : cos_theta_y_max < 0.0 ? wedge_origin[1] + rmin * cos_theta_y_max - : wedge_origin[1] + rmax * cos_theta_y_max; - const double xmin_theta = cos_theta < 0 ? theta_max : theta_min; - const double xmin_r = xmin_theta < pi ? rmin : rmax; - const double xmax_theta = cos_theta < 0 ? theta_min : theta_max; - const double xmax_r = xmax_theta < pi ? rmax : rmin; - const double sign_min = xmin_theta < pi ? 1.0 : -1.0; - const double sign_max = xmax_theta < pi ? 1.0 : -1.0; - const double xmin = - theta_max > 3. / 2. * pi && theta_min < 3. / 2. * pi - ? wedge_origin[0] - rmax - : wedge_origin[0] + - sign_min * sqrt(xmin_r * xmin_r * (1.0 - cos(xmin_theta) * cos(xmin_theta))); - - const double xmax = - theta_max > pi / 2. && theta_min < pi / 2.0 - ? wedge_origin[0] + rmax - : wedge_origin[0] + - sign_max * sqrt(xmax_r * xmax_r * (1.0 - cos(xmax_theta) * cos(xmax_theta))); - win_min[0] = xmin; - win_min[1] = ymin; - win_max[0] = xmax; - win_max[1] = ymax; - Ensure(!(win_min[0] > win_max[0])); - Ensure(!(win_min[1] > win_max[1])); - return; -} - -//------------------------------------------------------------------------------------------------// -/*! - * \brief Map data to sphere grid window for vector data - * - * Maps local+ghost data to a fixed r-theta mesh grid based on a specified weighting type. This data - * can additionally be normalized and positively biased on the grid. + * \brief Calculate the orthogonal distance between two points * - * - * \param[in] local_data the local data on the processor to be mapped to the window - * \param[in] ghost_data the ghost data on the processor to be mapped to the window - * \param[in,out] grid_data the resulting data map - * \param[in] sphere_center the center location of the sphere mesh - * \param[in] wedge_window_center the geometric center (x,y,x) of the wedge window - * \param[in] wedge_dr_dtheta the differential size in each direction (dr, dtheta, 0.0) used to form - * the wedge - * \param[in] grid_bins number of equally spaced bins in each dir - * \param[in] map_type_in string indicating the mapping (max, min, ave) - * \param[in] normalize bool operator to specify if the data should be normalized to a pdf - * \param[in] bias bool operator to specify if the data should be moved to the positive domain space - */ -void quick_index::map_data_to_sphere_grid_window( - const std::vector &local_data, const std::vector &ghost_data, - std::vector &grid_data, const std::array &sphere_center, - const std::array &wedge_window_center, const std::array &wedge_dr_dtheta, - const std::array &grid_bins, const std::string &map_type_in, const bool normalize, - const bool bias) const { - Insist(dim > 1, "Sphere grid window is invalid in 1d geometry"); - const auto r_theta = transform_r_theta(sphere_center, wedge_window_center); - // Store some r-theta values - const std::array r_theta_phi_max{r_theta[0] + wedge_dr_dtheta[0], - r_theta[1] + wedge_dr_dtheta[1], 0.0}; - const std::array r_theta_phi_min{std::max(r_theta[0] - wedge_dr_dtheta[0], 0.0), - wedge_dr_dtheta[1] < r_theta[1] - ? r_theta[1] - wedge_dr_dtheta[1] - : r_theta[1] - wedge_dr_dtheta[1] + 2. * pi, - 0.0}; - Check(!(r_theta_phi_min[1] > 2. * pi)); - // setup the xy window_max_min - std::array window_max{0.0, 0.0, 0.0}; - std::array window_min{0.0, 0.0, 0.0}; - calc_wedge_xy_bounds(wedge_window_center, sphere_center, wedge_dr_dtheta, window_min, window_max); - - Require(local_data.size() == n_locations); - Require(!(window_max[0] < window_min[0])); - Require(!(window_max[1] < window_min[1])); - Require(!(window_max[2] < window_min[2])); - Require(domain_decomposed ? ghost_data.size() == local_ghost_buffer_size : true); - Require(domain_decomposed - ? (fabs(window_max[0] - window_min[0]) - max_window_size) / max_window_size < 1e-6 - : true); - Require(domain_decomposed - ? (fabs(window_max[1] - window_min[1]) - max_window_size) / max_window_size < 1e-6 - : true); - Require(domain_decomposed - ? (fabs(window_max[2] - window_min[2]) - max_window_size) / max_window_size < 1e-6 - : true); - - bool fill = false; - std::string map_type = map_type_in; - if (map_type_in == "max_fill") { - Insist((grid_bins[0] > 1 && grid_bins[1] <= 1 && grid_bins[2] <= 1) || - (grid_bins[1] > 1 && grid_bins[0] <= 1 && grid_bins[2] <= 1) || - (grid_bins[2] > 1 && grid_bins[0] <= 1 && grid_bins[1] <= 1), - "one of grid bins must be == 1, Grid must be 1D to use max_fill option"); - fill = true; - map_type = "max"; - } else if (map_type_in == "min_fill") { - Insist((grid_bins[0] > 1 && grid_bins[1] <= 1 && grid_bins[2] <= 1) || - (grid_bins[1] > 1 && grid_bins[0] <= 1 && grid_bins[2] <= 1) || - (grid_bins[2] > 1 && grid_bins[0] <= 1 && grid_bins[1] <= 1), - "one of grid bins must be == 1, Grid must be 1D to use min_fill option"); - fill = true; - map_type = "min"; - } else if (map_type_in == "ave_fill") { - Insist((grid_bins[0] > 1 && grid_bins[1] <= 1 && grid_bins[2] <= 1) || - (grid_bins[1] > 1 && grid_bins[0] <= 1 && grid_bins[2] <= 1) || - (grid_bins[2] > 1 && grid_bins[0] <= 1 && grid_bins[1] <= 1), - "one of grid bins must be == 1, Grid must be 1D to use ave_fill option"); - fill = true; - map_type = "ave"; - } else if (map_type_in == "nearest_fill") { - Insist((grid_bins[0] > 1 && grid_bins[1] <= 1 && grid_bins[2] <= 1) || - (grid_bins[1] > 1 && grid_bins[0] <= 1 && grid_bins[2] <= 1) || - (grid_bins[2] > 1 && grid_bins[0] <= 1 && grid_bins[1] <= 1), - "one of grid bins must be == 1, Grid must be 1D to use ave_fill option"); - fill = true; - map_type = "nearest"; - } - - for (size_t d = 0; d < dim; d++) - Insist(grid_bins[d] > 0, "Bin size must be greater then zero for each active dimension"); - - size_t n_map_bins = 1; - for (size_t d = 0; d < dim; d++) - n_map_bins *= grid_bins[d]; - - Insist(grid_data.size() == n_map_bins, - "grid_data must match the flatten grid_bin size for the active dimensions (in 3d " - "grid_data.size()==grib_bins[0]*grid_bins[1]*grid_bins[2])"); - - std::fill(grid_data.begin(), grid_data.end(), 0.0); - - // Grab the global bins that lie in this window - std::vector global_bins = window_coarse_index_list(window_min, window_max); - - std::vector data_count(n_map_bins, 0); - std::vector min_distance(n_map_bins, 0); - double bias_cell_count = 0.0; - // Loop over all possible bins - for (auto &cb : global_bins) { - // skip bins that aren't present in the map (can't use [] operator with constness) - // loop over the local data - auto mapItr = coarse_index_map.find(cb); - if (mapItr != coarse_index_map.end()) { - for (auto &l : mapItr->second) { - bool valid; - size_t local_window_bin; - std::array bin_center; - std::tie(valid, local_window_bin, bin_center) = - get_sphere_window_bin(grid_bins, transform_r_theta(sphere_center, locations[l]), - r_theta_phi_min, r_theta_phi_max, n_map_bins, pi); - - // If the bin is outside the window continue to the next point - if (!valid) - continue; - - // lambda for mapping the data - map_data(bias_cell_count, data_count, grid_data, min_distance, dim, map_type, local_data, - bin_center, locations[l], local_window_bin, l); - - } // end local point loop - } // if valid local bin loop - if (domain_decomposed) { - // loop over the ghost data - auto gmapItr = local_ghost_index_map.find(cb); - if (gmapItr != local_ghost_index_map.end()) { - // loop over ghost data - for (auto &g : gmapItr->second) { - bool valid; - size_t local_window_bin; - std::array bin_center; - std::tie(valid, local_window_bin, bin_center) = get_sphere_window_bin( - grid_bins, transform_r_theta(sphere_center, local_ghost_locations[g]), - r_theta_phi_min, r_theta_phi_max, n_map_bins, pi); - - // If the bin is outside the window continue to the next poin - if (!valid) - continue; - - // lambda for mapping the data - map_data(bias_cell_count, data_count, grid_data, min_distance, dim, map_type, ghost_data, - bin_center, local_ghost_locations[g], local_window_bin, g); - } // end ghost point loop - } // if valid ghost bin - } // if dd - } // end coarse bin loop - - if (map_type == "ave" || map_type == "nearest") { - for (size_t i = 0; i < n_map_bins; i++) { - if (data_count[i] > 0) { - grid_data[i] /= data_count[i]; - } - } - } - if (fill) { - double last_val = 0.0; - int last_data_count = 0; - for (size_t i = 0; i < n_map_bins; i++) { - if (data_count[i] > 0) { - last_val = grid_data[i]; - last_data_count = data_count[i]; - } else { - grid_data[i] = last_val; - data_count[i] = last_data_count; - } - } - } - - if (bias && normalize) { - // return a positive normalized distribution - const double bias_value = - fabs(std::min(0.0, *std::min_element(grid_data.begin(), grid_data.end()))); - const double sum = - std::accumulate(grid_data.begin(), grid_data.end(), 0.0) + bias_value * bias_cell_count; - // catch zero instance - const double scale = !rtt_dsxx::soft_equiv(sum, 0.0) ? 1.0 / sum : 1.0; - for (size_t i = 0; i < n_map_bins; i++) - grid_data[i] = (grid_data[i] + bias_value) * scale; - } else if (bias) { - // return a positive distribution - const double bias_value = - fabs(std::min(0.0, *std::min_element(grid_data.begin(), grid_data.end()))); - for (size_t i = 0; i < n_map_bins; i++) - grid_data[i] += bias_value; - } else if (normalize) { - // return a normalized distribution - const double sum = std::accumulate(grid_data.begin(), grid_data.end(), 0.0); - // catch zero instance - const double scale = !rtt_dsxx::soft_equiv(sum, 0.0) ? 1.0 / sum : 1.0; - for (size_t i = 0; i < n_map_bins; i++) - grid_data[i] *= scale; - } -} - -//------------------------------------------------------------------------------------------------// -/*! - * \brief map_data_to_sphere_grid_window for vector> data - * - * Maps local+ghost data to a fixed r-theta mesh grid based on a specified weighting type. This data - * can additionally be normalized and positively biased on the grid. + * Maps multiple local+ghost data vectors to a fixed mesh grid based on a specified weighting type. + * This data can additionally be normalized and positively biased on the grid. * * - * \param[in] local_data the local data on the processor to be mapped to the window - * \param[in] ghost_data the ghost data on the processor to be mapped to the window - * \param[in,out] grid_data the resulting data map - * \param[in] sphere_center the center location of the sphere mesh - * \param[in] wedge_window_center the geometric center (x,y,x) of the wedge window - * \param[in] wedge_dr_dtheta the differential size in each direction (dr, dtheta, 0.0) used to form - * the wedge - * \param[in] grid_bins number of equally spaced bins in each dir - * \param[in] map_type_in string indicating the mapping (max, min, ave) - * \param[in] normalize bool operator to specify if the data should be normalized to a pdf - * (independent of each data vector) - * \param[in] bias bool operator to specify if the data should be moved to the positive domain space - * (independent of each data vector) - * \return bin_list list of global bins requested for the current window. + * \param[in] r0 initial position + * \param[in] r final position + * \param[in] arch_radius this is used for spherical geometry to determine at what radial point the + * orthognal archlength is measured. This point must be bound by the initial and final point radius. + * \return orthognal distance between the initial and final point in each direction */ -void quick_index::map_data_to_sphere_grid_window( - const std::vector> &local_data, - const std::vector> &ghost_data, std::vector> &grid_data, - const std::array &sphere_center, const std::array &wedge_window_center, - const std::array &wedge_dr_dtheta, const std::array &grid_bins, - const std::string &map_type_in, const bool normalize, const bool bias) const { - Insist(dim > 1, "Sphere grid window is invalid in 1d geometry"); - const auto r_theta = transform_r_theta(sphere_center, wedge_window_center); - // Store some r-theta values - const std::array r_theta_phi_max{r_theta[0] + wedge_dr_dtheta[0], - r_theta[1] + wedge_dr_dtheta[1], 0.0}; - const std::array r_theta_phi_min{std::max(r_theta[0] - wedge_dr_dtheta[0], 0.0), - wedge_dr_dtheta[1] < r_theta[1] - ? r_theta[1] - wedge_dr_dtheta[1] - : r_theta[1] - wedge_dr_dtheta[1] + 2. * pi, - 0.0}; - Check(!(r_theta_phi_min[1] > 2. * pi)); - // setup the xy window_max_min - std::array window_max{0.0, 0.0, 0.0}; - std::array window_min{0.0, 0.0, 0.0}; - calc_wedge_xy_bounds(wedge_window_center, sphere_center, wedge_dr_dtheta, window_min, window_max); - - Require(domain_decomposed ? local_data.size() == ghost_data.size() : true); - Require(!(window_max[0] < window_min[0])); - Require(!(window_max[1] < window_min[1])); - Require(!(window_max[2] < window_min[2])); - Require(domain_decomposed - ? (fabs(window_max[0] - window_min[0]) - max_window_size) / max_window_size < 1e-6 - : true); - Require(domain_decomposed - ? (fabs(window_max[1] - window_min[1]) - max_window_size) / max_window_size < 1e-6 - : true); - Require(domain_decomposed - ? (fabs(window_max[2] - window_min[2]) - max_window_size) / max_window_size < 1e-6 - : true); - - bool fill = false; - std::string map_type = map_type_in; - if (map_type_in == "max_fill") { - Insist((grid_bins[0] > 1 && grid_bins[1] <= 1 && grid_bins[2] <= 1) || - (grid_bins[1] > 1 && grid_bins[0] <= 1 && grid_bins[2] <= 1) || - (grid_bins[2] > 1 && grid_bins[0] <= 1 && grid_bins[1] <= 1), - "one of grid bins must be == 1, Grid must be 1D to use max_fill option"); - fill = true; - map_type = "max"; - } else if (map_type_in == "min_fill") { - Insist((grid_bins[0] > 1 && grid_bins[1] <= 1 && grid_bins[2] <= 1) || - (grid_bins[1] > 1 && grid_bins[0] <= 1 && grid_bins[2] <= 1) || - (grid_bins[2] > 1 && grid_bins[0] <= 1 && grid_bins[1] <= 1), - "one of grid bins must be == 1, Grid must be 1D to use min_fill option"); - fill = true; - map_type = "min"; - } else if (map_type_in == "ave_fill") { - Insist((grid_bins[0] > 1 && grid_bins[1] <= 1 && grid_bins[2] <= 1) || - (grid_bins[1] > 1 && grid_bins[0] <= 1 && grid_bins[2] <= 1) || - (grid_bins[2] > 1 && grid_bins[0] <= 1 && grid_bins[1] <= 1), - "one of grid bins must be == 1, Grid must be 1D to use ave_fill option"); - fill = true; - map_type = "ave"; - } else if (map_type_in == "nearest_fill") { - Insist((grid_bins[0] > 1 && grid_bins[1] <= 1 && grid_bins[2] <= 1) || - (grid_bins[1] > 1 && grid_bins[0] <= 1 && grid_bins[2] <= 1) || - (grid_bins[2] > 1 && grid_bins[0] <= 1 && grid_bins[1] <= 1), - "one of grid bins must be == 1, Grid must be 1D to use ave_fill option"); - fill = true; - map_type = "nearest"; - } - - for (size_t d = 0; d < dim; d++) - Insist(grid_bins[d] > 0, "Bin size must be greater then zero for each active dimension"); - - const size_t vsize = local_data.size(); - // Grab the global bins that lie in this window - std::vector global_bins = window_coarse_index_list(window_min, window_max); - size_t n_map_bins = 1; - for (size_t d = 0; d < dim; d++) { - n_map_bins *= grid_bins[d]; - } - - for (size_t v = 0; v < vsize; v++) { - Insist(grid_data[v].size() == n_map_bins, - "grid_data[" + std::to_string(v) + - "] must match the flatten grid_bin size for the active dimensions (in 3d " - "grid_data.size()==grib_bins[0]*grid_bins[1]*grid_bins[2])"); - std::fill(grid_data[v].begin(), grid_data[v].end(), 0.0); - } - - // initialize grid data - std::vector data_count(n_map_bins, 0); - std::vector min_distance(n_map_bins, 0); - double bias_cell_count = 0.0; - // Loop over all possible bins - for (auto &cb : global_bins) { - // skip bins that aren't present in the map (can't use [] operator with constness) - // loop over the local data - auto mapItr = coarse_index_map.find(cb); - if (mapItr != coarse_index_map.end()) { - for (auto &l : mapItr->second) { - bool valid; - size_t local_window_bin; - std::array bin_center; - std::tie(valid, local_window_bin, bin_center) = - get_sphere_window_bin(grid_bins, transform_r_theta(sphere_center, locations[l]), - r_theta_phi_min, r_theta_phi_max, n_map_bins, pi); - // If the bin is outside the window continue to the next poin - if (!valid) - continue; - Check(local_window_bin < n_map_bins); - map_vector_data(bias_cell_count, data_count, grid_data, min_distance, dim, map_type, - local_data, bin_center, locations[l], local_window_bin, l, vsize); - } // end local point loop - } // if valid local bin loop - if (domain_decomposed) { - // loop over the ghost data - auto gmapItr = local_ghost_index_map.find(cb); - if (gmapItr != local_ghost_index_map.end()) { - // loop over ghost data - for (auto &g : gmapItr->second) { - bool valid; - size_t local_window_bin; - std::array bin_center; - std::tie(valid, local_window_bin, bin_center) = get_sphere_window_bin( - grid_bins, transform_r_theta(sphere_center, local_ghost_locations[g]), - r_theta_phi_min, r_theta_phi_max, n_map_bins, pi); - - // If the bin is outside the window continue to the next poin - if (!valid) - continue; - map_vector_data(bias_cell_count, data_count, grid_data, min_distance, dim, map_type, - ghost_data, bin_center, local_ghost_locations[g], local_window_bin, g, - vsize); - } // end ghost point loop - } // if valid ghost bin - } // if dd - } // end coarse bin loop - - if (map_type == "ave" || map_type == "nearest") { - for (size_t i = 0; i < n_map_bins; i++) { - for (size_t v = 0; v < vsize; v++) { - if (data_count[i] > 0) { - grid_data[v][i] /= data_count[i]; - } - } - } - } - if (fill) { - std::vector last_val(vsize, 0.0); - int last_data_count = 0; - for (size_t i = 0; i < n_map_bins; i++) { - for (size_t v = 0; v < vsize; v++) { - if (data_count[i] > 0) { - last_val[v] = grid_data[v][i]; - last_data_count = data_count[i]; - } else { - grid_data[v][i] = last_val[v]; - if (v == vsize - 1) - data_count[i] = last_data_count; - } - } - } - } - - if (bias && normalize) { - // return a positive normalized distribution - for (size_t v = 0; v < vsize; v++) { - const double bias_value = - fabs(std::min(0.0, *std::min_element(grid_data[v].begin(), grid_data[v].end()))); - const double sum = std::accumulate(grid_data[v].begin(), grid_data[v].end(), 0.0) + - bias_value * bias_cell_count; - // catch zero instance - const double scale = !rtt_dsxx::soft_equiv(sum, 0.0) ? 1.0 / sum : 1.0; - for (size_t i = 0; i < n_map_bins; i++) - if (data_count[i] > 0) - grid_data[v][i] = (grid_data[v][i] + bias_value) * scale; - } - } else if (bias) { - // return a positive distribution - for (size_t v = 0; v < vsize; v++) { - const double bias_value = - fabs(std::min(0.0, *std::min_element(grid_data[v].begin(), grid_data[v].end()))); - for (size_t i = 0; i < n_map_bins; i++) - if (data_count[i] > 0) - grid_data[v][i] += bias_value; - } - } else if (normalize) { - // return a normalized distribution - for (size_t v = 0; v < vsize; v++) { - const double sum = std::accumulate(grid_data[v].begin(), grid_data[v].end(), 0.0); - // catch zero instance - const double scale = !rtt_dsxx::soft_equiv(sum, 0.0) ? 1.0 / sum : 1.0; - for (size_t i = 0; i < n_map_bins; i++) - if (data_count[i] > 0) - grid_data[v][i] *= scale; - } - } +std::array quick_index::calc_orthogonal_distance(const std::array &r0, + const std::array &r, + const double arch_radius) const { + Require(spherical ? dim == 2 : true); + Require(spherical ? !(arch_radius < 0.0) : true); + return {r[0] - r0[0], spherical ? arch_radius * (r[1] - r0[1]) : r[1] - r0[1], r[2] - r0[2]}; } } // namespace rtt_kde diff --git a/src/kde/quick_index.hh b/src/kde/quick_index.hh index dc87e281f..93bd635dc 100644 --- a/src/kde/quick_index.hh +++ b/src/kde/quick_index.hh @@ -5,7 +5,7 @@ * \brief This class generates coarse spatial indexing to quickly access near-neighbor data. This * additionally provides simple interpolation schemes to map data to simple structured * meshes. - * \note Copyright (C) 2021 Triad National Security, LLC., All rights reserved. + * \note Copyright (C) 2021-2021 Triad National Security, LLC., All rights reserved. */ //------------------------------------------------------------------------------------------------// @@ -15,11 +15,41 @@ #include "c4/global.hh" #include "units/MathConstants.hh" #include +#include #include #include namespace rtt_kde { +//------------------------------------------------------------------------------------------------// +/*! + * \brief Transform location array (x, y, z) positions to (r, theta, phi) grid + * + * Calculate a relative r theta and phi coordinate relative to a sphere center location from a + * standard (x,y,z) or (r,z) coordinates + * + * \param[in] sphere_center center of sphere in (x,y,z) or (r,z) coordinates + * \param[in] locations (x,y,z) or (r,z) locations to transform to relative (r, theta, phi) space. + * + */ +inline std::vector> +transform_spherical(const size_t dim, const std::array &sphere_center, + const std::vector> &locations) { + Insist(dim == 2, "Transform_r_theta Only implemented in 2d"); + std::vector> r_theta_locations(locations); + for (auto &location : r_theta_locations) { + const std::array v{location[0] - sphere_center[0], location[1] - sphere_center[1], + 0.0}; + const double r = sqrt(v[0] * v[0] + v[1] * v[1]); + const double mag = sqrt(v[0] * v[0] + v[1] * v[1]); + double cos_theta = mag > 0.0 ? std::max(std::min(v[1] / mag, 1.0), -1.0) : 0.0; + location = std::array{ + r, location[0] < sphere_center[0] ? 2.0 * rtt_units::PI - acos(cos_theta) : acos(cos_theta), + 0.0}; + } + return r_theta_locations; +} + //================================================================================================// /*! * \brief quick_index @@ -31,10 +61,11 @@ namespace rtt_kde { class quick_index { public: - //! Default constructors. + //! cartsian constructor quick_index(const size_t dim, const std::vector> &locations, const double max_window_size, const size_t bins_per_dimension, - const bool domain_decomposed); + const bool domain_decomposed, const bool spherical = false, + const std::array &sphere_center = {0.0, 0.0, 0.0}); //! Collect Ghost Data void collect_ghost_data(const std::vector &local_data, @@ -70,40 +101,21 @@ public: const std::array &grid_bins, const std::string &map_type, const bool normalize, const bool bias) const; - //! Map local+ghost data to grid window - void map_data_to_sphere_grid_window( - const std::vector &local_data, const std::vector &ghost_data, - std::vector &grid_data, const std::array &sphere_center, - const std::array &window_min, const std::array &window_max, - const std::array &grid_bins, const std::string &map_type, const bool normalize, - const bool bias) const; + //! Calculate the orthogonal distance between to locations + std::array calc_orthogonal_distance(const std::array &r0, + const std::array &r, + const double arch_radius) const; - //! Map local+ghost data to grid window for multi-dimensional data - void map_data_to_sphere_grid_window(const std::vector> &local_data, - const std::vector> &ghost_data, - std::vector> &grid_data, - const std::array &sphere_center, - const std::array &wedge_window_center, - const std::array &wedge_dr_dtheta, - const std::array &grid_bins, - const std::string &map_type, const bool normalize, - const bool bias) const; - - std::array transform_r_theta(const std::array &sphere_center, - const std::array &location) const; - - void calc_wedge_xy_bounds(const std::array &wedge_xyz_center, - const std::array &wedge_origin, - const std::array &wedge_dr_dtheta, - std::array &win_min, std::array &win_max) const; // PUBLIC DATA // Quick index initialization data const size_t dim; const bool domain_decomposed; + const double spherical; + const std::array sphere_center; const size_t coarse_bin_resolution; const double max_window_size; - // keep a copy of the locations - const std::vector> locations; + // keep a mutable copy of the locations (can be change to spherical) + std::vector> locations; const size_t n_locations; // Global bounds @@ -133,7 +145,6 @@ private: std::map>> put_window_map; // max put buffer size; size_t max_put_buffer_size; - const double pi = rtt_units::PI; }; } // end namespace rtt_kde diff --git a/src/kde/test/tstkde.cc b/src/kde/test/tstkde.cc index 4484895c4..c546d7011 100644 --- a/src/kde/test/tstkde.cc +++ b/src/kde/test/tstkde.cc @@ -4,7 +4,7 @@ * \author Mathew Cleveland * \date Nov. 10th 2020 * \brief KDE function tests - * \note Copyright (C) 2020-2021 Triad National Security, LLC., All rights reserved. + * \note Copyright (C) 2021-2021 Triad National Security, LLC., All rights reserved. */ //------------------------------------------------------------------------------------------------// @@ -30,30 +30,15 @@ void test_replication(ParallelUnitTest &ut) { if (!rtt_dsxx::soft_equiv(value, 0.75)) ITFAILS; - // test some public sphere calculations - { - const std::array sphere_center{0.0, 0.0, 0.0}; - const std::array location{sqrt(2), sqrt(2), 0.0}; - const std::array location2{0.0, 2.0, 0.0}; - const double radius = 2.0; - const double small_radius = 1.0; - const double pi_over_4 = 0.78539816; - FAIL_IF_NOT(rtt_dsxx::soft_equiv(test_kde.calc_radius(sphere_center, location), 2.0)); - FAIL_IF_NOT( - rtt_dsxx::soft_equiv(test_kde.calc_arch_length(sphere_center, radius, location, location2), - 2.0 * pi_over_4, 1e-6)); - FAIL_IF_NOT(rtt_dsxx::soft_equiv( - test_kde.calc_arch_length(sphere_center, small_radius, location, location2), pi_over_4, - 1e-6)); - } - // spherical reconstruction { + const bool spherical = true; const std::array sphere_center{0.0, -1.0, 0.0}; const double max_radius = 1.0; const double min_radius = 0.0; kde sphere_kde; - sphere_kde.set_sphere_center(sphere_center, min_radius, max_radius); + // reflect on the theta boundary + kde theta_reflected_sphere_kde({false, false, true, true, false, false}); const std::array radial_edges{0.025, 0.05, 0.075, 0.1, 0.25, 0.5, 0.75, 1.0}; const std::array cosine_edges{-1.0, -0.75, -0.5, -0.25, 0.0, 0.25, 0.5, 0.75, 1.0}; const size_t data_size = radial_edges.size() * cosine_edges.size(); @@ -83,6 +68,7 @@ void test_replication(ParallelUnitTest &ut) { // zero reconstruction array { std::vector zero_data(data_size, 0.0); + // R-theta space bandwidths in spherical reconstruction std::vector> one_over_bandwidth_array( data_size, std::array{1.0, 1.0e12, 0.0}); const bool dd = false; @@ -90,7 +76,8 @@ void test_replication(ParallelUnitTest &ut) { const size_t n_coarse_bins = 5; const double max_window_size = 1.0; const size_t dim = 2; - quick_index qindex(dim, position_array, max_window_size, n_coarse_bins, dd); + quick_index qindex(dim, position_array, max_window_size, n_coarse_bins, dd, spherical, + sphere_center); std::vector smooth_result = sphere_kde.reconstruction(zero_data, one_over_bandwidth_array, qindex); @@ -126,7 +113,8 @@ void test_replication(ParallelUnitTest &ut) { const size_t n_coarse_bins = 5; const double max_window_size = 1.0; const size_t dim = 2; - quick_index qindex(dim, position_array, max_window_size, n_coarse_bins, dd); + quick_index qindex(dim, position_array, max_window_size, n_coarse_bins, dd, spherical, + sphere_center); std::vector smooth_result = sphere_kde.reconstruction(spoke_data, one_over_bandwidth_array, qindex); @@ -163,7 +151,8 @@ void test_replication(ParallelUnitTest &ut) { const size_t n_coarse_bins = 5; const double max_window_size = 1.0; const size_t dim = 2; - quick_index qindex(dim, position_array, max_window_size, n_coarse_bins, dd); + quick_index qindex(dim, position_array, max_window_size, n_coarse_bins, dd, spherical, + sphere_center); std::vector smooth_result = sphere_kde.reconstruction(shell_data, one_over_bandwidth_array, qindex); @@ -195,14 +184,14 @@ void test_replication(ParallelUnitTest &ut) { { std::vector spoke_smoothed_shells{ - 2.51488, 2.99904, 2.99904, 3.69002, 3.72457, 3.69002, 2.99904, 2.99904, 2.51488, - 2.51645, 3.00418, 3.00418, 3.72015, 3.7866, 3.72015, 3.00418, 3.00418, 2.51645, - 2.51803, 3.00928, 3.00928, 3.74919, 3.84522, 3.74919, 3.00928, 3.00928, 2.51803, - 2.51961, 3.01436, 3.01436, 3.77729, 3.90089, 3.77729, 3.01436, 3.01436, 2.51961, - 5.52169, 3.04531, 3.04531, 3.93334, 4.19165, 3.93334, 3.04531, 3.04531, 5.52169, - 5.55417, 6.52859, 6.95461, 4.19454, 4.61685, 4.19454, 6.95461, 6.52859, 5.55417, - 7.53548, 6.56107, 7, 4.58158, 5.14978, 4.58158, 7, 6.56107, 7.53548, - 7.56796, 8, 7.14194, 8, 6.33455, 8, 7.14194, 8, 7.56796}; + 3.7717, 3.7717, 3.7717, 3.7717, 3.7717, 3.7717, 3.7717, 3.7717, 3.7717, + 3.83452, 3.83452, 3.83452, 3.83452, 3.83452, 3.83452, 3.83452, 3.83452, 3.83452, + 3.89388, 3.89388, 3.89388, 3.89388, 3.89388, 3.89388, 3.89388, 3.89388, 3.89388, + 3.95025, 3.95025, 3.95025, 3.95025, 3.95025, 3.95025, 3.95025, 3.95025, 3.95025, + 4.24469, 4.24469, 4.24469, 4.24469, 4.24469, 4.24469, 4.24469, 4.24469, 4.24469, + 4.67528, 4.67528, 4.67528, 4.67528, 4.67528, 4.67528, 4.67528, 4.67528, 4.67528, + 5.21495, 5.21495, 5.21495, 5.21495, 5.21495, 5.21495, 5.21495, 5.21495, 5.21495, + 6.41471, 6.41472, 6.41472, 6.41472, 6.41471, 6.41471, 6.41471, 6.41471, 6.41471}; std::vector> one_over_bandwidth_array( data_size, std::array{1.0, 1.0e12, 0.0}); const bool dd = false; @@ -210,7 +199,8 @@ void test_replication(ParallelUnitTest &ut) { const size_t n_coarse_bins = 5; const double max_window_size = 1.0; const size_t dim = 2; - quick_index qindex(dim, position_array, max_window_size, n_coarse_bins, dd); + quick_index qindex(dim, position_array, max_window_size, n_coarse_bins, dd, spherical, + sphere_center); std::vector smooth_result = sphere_kde.reconstruction(shell_data, one_over_bandwidth_array, qindex); @@ -241,14 +231,14 @@ void test_replication(ParallelUnitTest &ut) { // shell smoothing on spoke array { std::vector shell_smoothed_spoke{ - 4.82322, 4.82519, 4.82608, 4.82681, 4.8275, 4.82819, 4.82892, 4.8298, 4.83177, - 4.81029, 4.81825, 4.8218, 4.82475, 4.8275, 4.83025, 4.8332, 4.83675, 4.8447, - 4.78839, 4.80659, 4.81462, 4.82129, 4.8275, 4.83371, 4.84037, 4.84841, 4.86661, - 4.75694, 4.79008, 4.8045, 4.81642, 4.8275, 4.83857, 4.8505, 4.86492, 4.89805, - 4.04388, 4.22611, 4.67522, 4.75503, 5, 4.97765, 5.14326, 6.35454, 5.61112, - 2.62832, 3.75795, 4.12091, 4.47922, 5, 6.9148, 6.66199, 7.35789, 7.02668, - 1.66976, 3.01482, 3.72878, 4.43378, 5, 8.6895, 8.00342, 6.64018, 7.98524, - 1, 2, 4.02682, 4.51075, 5, 5.14424, 5.62818, 8, 9}; + 2.06793, 3.61898, 4.65287, 4.65345, 5.17042, 5.68739, 5.68796, 5.91166, 6.47113, + 2.06722, 3.61803, 4.65136, 4.65367, 5.17042, 5.68716, 5.68948, 5.91951, 6.49562, + 2.06603, 3.61645, 4.64883, 4.65404, 5.17042, 5.68679, 5.69201, 5.93267, 6.53704, + 2.06437, 3.61422, 4.64527, 4.65457, 5.17042, 5.68627, 5.69557, 5.95127, 6.59632, + 2.04372, 3.58678, 4.60095, 4.66103, 5, 5.67981, 5.73989, 6.19173, 7.4492, + 1.95893, 3.47577, 4.41363, 4.68754, 5, 5.6533, 5.92721, 6.86506, 8.38191, + 1.76287, 3.22898, 3.99365, 4.74873, 5, 5.59211, 6.34718, 7.11186, 8.57797, + 1.36835, 2.8677, 3.80371, 4.44763, 5, 5.89321, 6.53713, 7.47314, 8.97248}; std::vector> one_over_bandwidth_array( data_size, std::array{1.0e12, 1.0, 0.0}); const bool dd = false; @@ -256,7 +246,8 @@ void test_replication(ParallelUnitTest &ut) { const size_t n_coarse_bins = 5; const double max_window_size = 1.0; const size_t dim = 2; - quick_index qindex(dim, position_array, max_window_size, n_coarse_bins, dd); + quick_index qindex(dim, position_array, max_window_size, n_coarse_bins, dd, spherical, + sphere_center); std::vector smooth_result = sphere_kde.reconstruction(spoke_data, one_over_bandwidth_array, qindex); @@ -283,6 +274,53 @@ void test_replication(ParallelUnitTest &ut) { std::accumulate(log_smooth_result.begin(), log_smooth_result.end(), 0.0))) ITFAILS; } + + // shell smoothing on spoke array with theta reflection + { + std::vector shell_smoothed_spoke{ + 2.06881, 3.62039, 4.65429, 4.65452, 5.17084, 5.68715, 5.68738, 5.91075, 6.4677, + 2.0703, 3.62288, 4.65595, 4.65685, 5.17084, 5.68482, 5.68573, 5.91436, 6.48016, + 2.07296, 3.62722, 4.65887, 4.66084, 5.17084, 5.68084, 5.6828, 5.92023, 6.50097, + 2.07714, 3.63377, 4.66332, 4.66664, 5.17084, 5.67503, 5.67835, 5.92809, 6.53015, + 2.03363, 3.51215, 4.56263, 4.63799, 5, 5.70368, 5.77904, 6.28319, 7.525, + 1.91731, 3.16957, 4.17923, 4.56555, 5, 5.77612, 6.16244, 7.1721, 8.42436, + 1.67727, 3.05587, 3.93229, 4.74886, 5, 5.59282, 6.40939, 7.2858, 8.6644, + 1.31387, 2.7775, 3.80402, 4.44799, 5, 5.89368, 6.53766, 7.56417, 9.02781}; + std::vector> one_over_bandwidth_array( + data_size, std::array{1.0e12, 1.0, 0.0}); + const bool dd = false; + // two bins per point + const size_t n_coarse_bins = 5; + const double max_window_size = 1.0; + const size_t dim = 2; + quick_index qindex(dim, position_array, max_window_size, n_coarse_bins, dd, spherical, + sphere_center); + + std::vector smooth_result = + theta_reflected_sphere_kde.reconstruction(spoke_data, one_over_bandwidth_array, qindex); + std::vector log_smooth_result = + theta_reflected_sphere_kde.reconstruction(spoke_data, one_over_bandwidth_array, qindex); + // Apply Conservation + sphere_kde.apply_conservation(spoke_data, smooth_result, qindex.domain_decomposed); + sphere_kde.apply_conservation(spoke_data, log_smooth_result, qindex.domain_decomposed); + + // Check smooth result + for (size_t i = 0; i < data_size; i++) { + if (!rtt_dsxx::soft_equiv(smooth_result[i], shell_smoothed_spoke[i], 1e-3)) + ITFAILS; + if (!rtt_dsxx::soft_equiv(log_smooth_result[i], shell_smoothed_spoke[i], 1e-3)) + ITFAILS; + } + + // Energy conservation + if (!rtt_dsxx::soft_equiv(std::accumulate(spoke_data.begin(), spoke_data.end(), 0.0), + std::accumulate(smooth_result.begin(), smooth_result.end(), 0.0))) + ITFAILS; + if (!rtt_dsxx::soft_equiv( + std::accumulate(spoke_data.begin(), spoke_data.end(), 0.0), + std::accumulate(log_smooth_result.begin(), log_smooth_result.end(), 0.0))) + ITFAILS; + } } // No mean reconstruction because of small basis functions @@ -1005,15 +1043,15 @@ void test_decomposition(ParallelUnitTest &ut) { // spherical reconstruction { + const bool spherical = true; const size_t local_size = 24; const std::array sphere_center{0.0, -1.0, 0.0}; const double max_radius = 1.0; const double min_radius = 0.0; const double shell_min_radius = 0.5; kde sphere_kde; - sphere_kde.set_sphere_center(sphere_center, min_radius, max_radius); - kde shell_kde; - shell_kde.set_sphere_center(sphere_center, shell_min_radius, max_radius); + // reflect on the theta boundary + kde theta_reflected_sphere_kde({false, false, true, true, false, false}); const std::array radial_edges{0.025, 0.05, 0.075, 0.1, 0.25, 0.5, 0.75, 1.0}; const std::array cosine_edges{-1.0, -0.75, -0.5, -0.25, 0.0, 0.25, 0.5, 0.75, 1.0}; const size_t data_size = radial_edges.size() * cosine_edges.size(); @@ -1039,35 +1077,47 @@ void test_decomposition(ParallelUnitTest &ut) { ri++; } std::vector spoke_smoothed_shells{ - 2.51488, 2.99904, 2.99904, 3.69002, 3.72457, 3.69002, 2.99904, 2.99904, 2.51488, - 2.51645, 3.00418, 3.00418, 3.72015, 3.7866, 3.72015, 3.00418, 3.00418, 2.51645, - 2.51803, 3.00928, 3.00928, 3.74919, 3.84522, 3.74919, 3.00928, 3.00928, 2.51803, - 2.51961, 3.01436, 3.01436, 3.77729, 3.90089, 3.77729, 3.01436, 3.01436, 2.51961, - 5.52169, 3.04531, 3.04531, 3.93334, 4.19165, 3.93334, 3.04531, 3.04531, 5.52169, - 5.55417, 6.52859, 6.95461, 4.19454, 4.61685, 4.19454, 6.95461, 6.52859, 5.55417, - 7.53548, 6.56107, 7, 4.58158, 5.14978, 4.58158, 7, 6.56107, 7.53548, - 7.56796, 8, 7.14194, 8, 6.33455, 8, 7.14194, 8, 7.56796}; + 3.7717, 3.7717, 3.7717, 3.7717, 3.7717, 3.7717, 3.7717, 3.7717, 3.7717, + 3.83452, 3.83452, 3.83452, 3.83452, 3.83452, 3.83452, 3.83452, 3.83452, 3.83452, + 3.89388, 3.89388, 3.89388, 3.89388, 3.89388, 3.89388, 3.89388, 3.89388, 3.89388, + 3.95025, 3.95025, 3.95025, 3.95025, 3.95025, 3.95025, 3.95025, 3.95025, 3.95025, + 4.24469, 4.24469, 4.24469, 4.24469, 4.24469, 4.24469, 4.24469, 4.24469, 4.24469, + 4.67528, 4.67528, 4.67528, 4.67528, 4.67528, 4.67528, 4.67528, 4.67528, 4.67528, + 5.21495, 5.21495, 5.21495, 5.21495, 5.21495, 5.21495, 5.21495, 5.21495, 5.21495, + 6.41471, 6.41472, 6.41472, 6.41472, 6.41471, 6.41471, 6.41471, 6.41471, 6.41471}; std::vector shell_smoothed_spoke{ - 4.82322, 4.82519, 4.82608, 4.82681, 4.8275, 4.82819, 4.82892, 4.8298, 4.83177, - 4.81029, 4.81825, 4.8218, 4.82475, 4.8275, 4.83025, 4.8332, 4.83675, 4.8447, - 4.78839, 4.80659, 4.81462, 4.82129, 4.8275, 4.83371, 4.84037, 4.84841, 4.86661, - 4.75694, 4.79008, 4.8045, 4.81642, 4.8275, 4.83857, 4.8505, 4.86492, 4.89805, - 4.04388, 4.22611, 4.67522, 4.75503, 5, 4.97765, 5.14326, 6.35454, 5.61112, - 2.62832, 3.75795, 4.12091, 4.47922, 5, 6.9148, 6.66199, 7.35789, 7.02668, - 1.66976, 3.01482, 3.72878, 4.43378, 5, 8.6895, 8.00342, 6.64018, 7.98524, - 1, 2, 4.02682, 4.51075, 5, 5.14424, 5.62818, 8, 9}; + 2.06793, 3.61898, 4.65287, 4.65345, 5.17042, 5.68739, 5.68796, 5.91166, 6.47113, + 2.06722, 3.61803, 4.65136, 4.65367, 5.17042, 5.68716, 5.68948, 5.91951, 6.49562, + 2.06603, 3.61645, 4.64883, 4.65404, 5.17042, 5.68679, 5.69201, 5.93267, 6.53704, + 2.06437, 3.61422, 4.64527, 4.65457, 5.17042, 5.68627, 5.69557, 5.95127, 6.59632, + 2.04372, 3.58678, 4.60095, 4.66103, 5, 5.67981, 5.73989, 6.19173, 7.4492, + 1.95893, 3.47577, 4.41363, 4.68754, 5, 5.6533, 5.92721, 6.86506, 8.38191, + 1.76287, 3.22898, 3.99365, 4.74873, 5, 5.59211, 6.34718, 7.11186, 8.57797, + 1.36835, 2.8677, 3.80371, 4.44763, 5, 5.89321, 6.53713, 7.47314, 8.97248}; + std::vector shell_smoothed_spoke_theta_ref{ + 2.06881, 3.62039, 4.65429, 4.65452, 5.17084, 5.68715, 5.68738, 5.91075, 6.4677, + 2.0703, 3.62288, 4.65595, 4.65685, 5.17084, 5.68482, 5.68573, 5.91436, 6.48016, + 2.07296, 3.62722, 4.65887, 4.66084, 5.17084, 5.68084, 5.6828, 5.92023, 6.50097, + 2.07714, 3.63377, 4.66332, 4.66664, 5.17084, 5.67503, 5.67835, 5.92809, 6.53015, + 2.03363, 3.51215, 4.56263, 4.63799, 5, 5.70368, 5.77904, 6.28319, 7.525, + 1.91731, 3.16957, 4.17923, 4.56555, 5, 5.77612, 6.16244, 7.1721, 8.42436, + 1.67727, 3.05587, 3.93229, 4.74886, 5, 5.59282, 6.40939, 7.2858, 8.6644, + 1.31387, 2.7775, 3.80402, 4.44799, 5, 5.89368, 6.53766, 7.56417, 9.02781}; std::vector dd_const_data(local_size, 1.0); std::vector dd_spoke_data(local_size); std::vector dd_shell_data(local_size); std::vector dd_spoke_smoothed_shells(local_size); std::vector dd_shell_smoothed_spoke(local_size); + std::vector dd_shell_smoothed_spoke_theta_ref(local_size); std::vector> dd_position_array(local_size); for (size_t i = 0; i < local_size; i++) { dd_spoke_data[i] = spoke_data[i + rtt_c4::node() * local_size]; dd_shell_data[i] = shell_data[i + rtt_c4::node() * local_size]; dd_spoke_smoothed_shells[i] = spoke_smoothed_shells[i + rtt_c4::node() * local_size]; dd_shell_smoothed_spoke[i] = shell_smoothed_spoke[i + rtt_c4::node() * local_size]; + dd_shell_smoothed_spoke_theta_ref[i] = + shell_smoothed_spoke_theta_ref[i + rtt_c4::node() * local_size]; dd_position_array[i] = position_array[i + rtt_c4::node() * local_size]; } @@ -1081,15 +1131,16 @@ void test_decomposition(ParallelUnitTest &ut) { const size_t n_coarse_bins = 5; const double max_window_size = 2.0; const size_t dim = 2; - quick_index qindex(dim, dd_position_array, max_window_size, n_coarse_bins, dd); + quick_index qindex(dim, dd_position_array, max_window_size, n_coarse_bins, dd, spherical, + sphere_center); std::vector smooth_result = - shell_kde.reconstruction(zero_data, one_over_bandwidth_array, qindex); + sphere_kde.reconstruction(zero_data, one_over_bandwidth_array, qindex); std::vector log_smooth_result = - shell_kde.log_reconstruction(zero_data, one_over_bandwidth_array, qindex); + sphere_kde.log_reconstruction(zero_data, one_over_bandwidth_array, qindex); // Apply Conservation - shell_kde.apply_conservation(zero_data, smooth_result, qindex.domain_decomposed); - shell_kde.apply_conservation(zero_data, log_smooth_result, qindex.domain_decomposed); + sphere_kde.apply_conservation(zero_data, smooth_result, qindex.domain_decomposed); + sphere_kde.apply_conservation(zero_data, log_smooth_result, qindex.domain_decomposed); // Check smooth result for (size_t i = 0; i < local_size; i++) { @@ -1121,7 +1172,8 @@ void test_decomposition(ParallelUnitTest &ut) { const size_t n_coarse_bins = 5; const double max_window_size = 2.0; const size_t dim = 2; - quick_index qindex(dim, dd_position_array, max_window_size, n_coarse_bins, dd); + quick_index qindex(dim, dd_position_array, max_window_size, n_coarse_bins, dd, spherical, + sphere_center); std::vector smooth_result = sphere_kde.reconstruction(dd_spoke_data, one_over_bandwidth_array, qindex); @@ -1163,7 +1215,8 @@ void test_decomposition(ParallelUnitTest &ut) { const size_t n_coarse_bins = 5; const double max_window_size = 2.0; const size_t dim = 2; - quick_index qindex(dim, dd_position_array, max_window_size, n_coarse_bins, dd); + quick_index qindex(dim, dd_position_array, max_window_size, n_coarse_bins, dd, spherical, + sphere_center); std::vector smooth_result = sphere_kde.reconstruction(dd_shell_data, one_over_bandwidth_array, qindex); @@ -1205,7 +1258,8 @@ void test_decomposition(ParallelUnitTest &ut) { const size_t n_coarse_bins = 5; const double max_window_size = 2.0; const size_t dim = 2; - quick_index qindex(dim, dd_position_array, max_window_size, n_coarse_bins, dd); + quick_index qindex(dim, dd_position_array, max_window_size, n_coarse_bins, dd, spherical, + sphere_center); std::vector smooth_result = sphere_kde.reconstruction(dd_shell_data, one_over_bandwidth_array, qindex); @@ -1247,7 +1301,8 @@ void test_decomposition(ParallelUnitTest &ut) { const size_t n_coarse_bins = 5; const double max_window_size = 1.0; const size_t dim = 2; - quick_index qindex(dim, dd_position_array, max_window_size, n_coarse_bins, dd); + quick_index qindex(dim, dd_position_array, max_window_size, n_coarse_bins, dd, spherical, + sphere_center); std::vector smooth_result = sphere_kde.reconstruction(dd_spoke_data, one_over_bandwidth_array, qindex); @@ -1279,6 +1334,49 @@ void test_decomposition(ParallelUnitTest &ut) { log_smooth_conservation)) ITFAILS; } + + // shell smoothing on spoke array with theta reflection + { + std::vector> one_over_bandwidth_array( + local_size, std::array{1.0e12, 1.0, 0.0}); + const bool dd = true; + // two bins per point + const size_t n_coarse_bins = 5; + const double max_window_size = 1.0; + const size_t dim = 2; + quick_index qindex(dim, dd_position_array, max_window_size, n_coarse_bins, dd, spherical, + sphere_center); + + std::vector smooth_result = theta_reflected_sphere_kde.reconstruction( + dd_spoke_data, one_over_bandwidth_array, qindex); + std::vector log_smooth_result = theta_reflected_sphere_kde.reconstruction( + dd_spoke_data, one_over_bandwidth_array, qindex); + // Apply Conservation + sphere_kde.apply_conservation(dd_spoke_data, smooth_result, qindex.domain_decomposed); + sphere_kde.apply_conservation(dd_spoke_data, log_smooth_result, qindex.domain_decomposed); + + // Check smooth result + for (size_t i = 0; i < local_size; i++) { + if (!rtt_dsxx::soft_equiv(smooth_result[i], dd_shell_smoothed_spoke_theta_ref[i], 1e-3)) + ITFAILS; + if (!rtt_dsxx::soft_equiv(log_smooth_result[i], dd_shell_smoothed_spoke_theta_ref[i], 1e-3)) + ITFAILS; + } + + double smooth_conservation = std::accumulate(smooth_result.begin(), smooth_result.end(), 0.0); + rtt_c4::global_sum(smooth_conservation); + double log_smooth_conservation = + std::accumulate(log_smooth_result.begin(), log_smooth_result.end(), 0.0); + rtt_c4::global_sum(log_smooth_conservation); + + // Energy conservation + if (!rtt_dsxx::soft_equiv(std::accumulate(spoke_data.begin(), spoke_data.end(), 0.0), + smooth_conservation)) + ITFAILS; + if (!rtt_dsxx::soft_equiv(std::accumulate(spoke_data.begin(), spoke_data.end(), 0.0), + log_smooth_conservation)) + ITFAILS; + } } int local_size = 3; diff --git a/src/kde/test/tstquick_index.cc b/src/kde/test/tstquick_index.cc index cc6cb1392..3360e23e8 100644 --- a/src/kde/test/tstquick_index.cc +++ b/src/kde/test/tstquick_index.cc @@ -4,7 +4,7 @@ * \author Mathew Cleveland * \date Aug. 10th 2021 * \brief quick_index testing function - * \note Copyright (C) 2021 Triad National Security, LLC., All rights reserved. + * \note Copyright (C) 2021-2021 Triad National Security, LLC., All rights reserved. */ //------------------------------------------------------------------------------------------------// @@ -85,6 +85,72 @@ void test_replication(ParallelUnitTest &ut) { } } +void test_replication_sphere(ParallelUnitTest &ut) { + { + std::vector data{0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1}; + std::vector> position_array(10, std::array{0.0, 0.0, 0.0}); + for (int i = 0; i < 10; i++) { + position_array[i][0] = i < 5 ? i % 5 : i % 5 + 0.5; + position_array[i][1] = i < 5 ? 0.5 : -0.5; + } + // in rep mode the max window size does nothing so set it large + const bool spherical = true; + const std::array sphere_center{2.0, 0.0, 0.0}; + const double max_window_size = 100.0; + const size_t bins_per_dim = 10UL; + const bool dd = false; + const size_t dim = 2; + quick_index qindex(dim, position_array, max_window_size, bins_per_dim, dd, spherical, + sphere_center); + // Check public data + //------------------------ + if (qindex.domain_decomposed) + ITFAILS; + if (qindex.coarse_bin_resolution != bins_per_dim) + ITFAILS; + if (!soft_equiv(qindex.max_window_size, max_window_size)) + ITFAILS; + // Check global bounding box + if (!soft_equiv(qindex.bounding_box_min[0], 0.5)) + ITFAILS; + if (!soft_equiv(qindex.bounding_box_min[1], 0.0)) + ITFAILS; + if (!soft_equiv(qindex.bounding_box_min[2], 0.0)) + ITFAILS; + if (!soft_equiv(qindex.bounding_box_max[0], 2.54951, 1e-4)) + ITFAILS; + if (!soft_equiv(qindex.bounding_box_max[1], 5.17604, 1e-4)) + ITFAILS; + if (!soft_equiv(qindex.bounding_box_max[2], 0.0)) + ITFAILS; + // Check local coarse_index map + // build up a global gold to check the map + std::map> gold_map; + gold_map[0] = {2}; + gold_map[23] = {3}; + gold_map[27] = {4}; + gold_map[35] = {8}; + gold_map[39] = {9}; + gold_map[41] = {7}; + gold_map[71] = {6}; + gold_map[85] = {5}; + gold_map[93] = {1}; + gold_map[97] = {0}; + if (gold_map.size() != qindex.coarse_index_map.size()) + ITFAILS; + for (auto &map : qindex.coarse_index_map) + for (size_t i = 0; i < map.second.size(); i++) + if (gold_map[map.first][i] != map.second[i]) + ITFAILS; + } + + if (ut.numFails == 0) { + PASSMSG("quick_index sphere checks pass"); + } else { + FAILMSG("quick_index sphere checks failed"); + } +} + void test_decomposition(ParallelUnitTest &ut) { if (rtt_c4::nodes() != 3) ITFAILS; @@ -1247,294 +1313,6 @@ void test_decomposition(ParallelUnitTest &ut) { ITFAILS; } - // check max sphere r window mapping (more bins then data) functions - { - // build a length=1.0 window around the first point on each node - const std::array center{dd_position_array[0][0], dd_position_array[0][1], 0.0}; - // wedge location +- 1 degree theta and 0.5 radius - const std::array dr_dtheta{0.4, 0.0174533, 0.0}; - const std::array bin_sizes{5, 1, 0}; - const bool normalize = false; - const bool bias = false; - const std::string map_type = "max"; - const std::array sphere_center{-1.0, 0.0, 0.0}; - std::vector sphere_window_data(5, 0.0); - qindex.map_data_to_sphere_grid_window(dd_data, ghost_data, sphere_window_data, sphere_center, - center, dr_dtheta, bin_sizes, map_type, normalize, - bias); - std::vector> sphere_window_3x_data(3, std::vector(5, 0.0)); - qindex.map_data_to_sphere_grid_window(dd_3x_data, ghost_3x_data, sphere_window_3x_data, - sphere_center, center, dr_dtheta, bin_sizes, map_type, - normalize, bias); - std::vector gold_window_data; - std::vector> gold_window_3x_data(3); - // different result then 1D because the 1.0 y offset of the data - if (rtt_c4::node() == 0) { - gold_window_data = {0.0, 0.0, 3.0, 0.0, 0.0}; - gold_window_3x_data[0] = {0.0, 0.0, 3.0, 0.0, 0.0}; - gold_window_3x_data[1] = {0.0, 0.0, 4.0, 0.0, 0.0}; - gold_window_3x_data[2] = {0.0, 0.0, -3.0, 0.0, 0.0}; - } else if (rtt_c4::node() == 1) { - gold_window_data = {0.0, 0.0, 6.0, 0.0, 0.0}; - gold_window_3x_data[0] = {0.0, 0.0, 6.0, 0.0, 0.0}; - gold_window_3x_data[1] = {0.0, 0.0, 7.0, 0.0, 0.0}; - gold_window_3x_data[2] = {0.0, 0.0, -6.0, 0.0, 0.0}; - } else { - gold_window_data = {0.0, 0.0, 9.0, 0.0, 0.0}; - gold_window_3x_data[0] = {0.0, 0.0, 9.0, 0.0, 0.0}; - gold_window_3x_data[1] = {0.0, 0.0, 10.0, 0.0, 0.0}; - gold_window_3x_data[2] = {0.0, 0.0, -9.0, 0.0, 0.0}; - } - - for (size_t i = 0; i < bin_sizes[0]; i++) - if (!rtt_dsxx::soft_equiv(sphere_window_data[i], gold_window_data[i])) - ITFAILS; - for (size_t v = 0; v < 3; v++) - for (size_t i = 0; i < bin_sizes[0]; i++) - if (!rtt_dsxx::soft_equiv(sphere_window_3x_data[v][i], gold_window_3x_data[v][i])) - ITFAILS; - } - - // check max sphere r window mapping (more bins then data) functions with bias - { - // build a length=1.0 window around the first point on each node - const std::array center{dd_position_array[0][0], dd_position_array[0][1], 0.0}; - // wedge location +- 1 degree theta and 0.5 radius - const std::array dr_dtheta{0.4, 0.0174533, 0.0}; - const std::array bin_sizes{5, 1, 0}; - const bool normalize = false; - const bool bias = true; - const std::string map_type = "max"; - const std::array sphere_center{-1.0, 0.0, 0.0}; - std::vector sphere_window_data(5, 0.0); - qindex.map_data_to_sphere_grid_window(dd_data, ghost_data, sphere_window_data, sphere_center, - center, dr_dtheta, bin_sizes, map_type, normalize, - bias); - std::vector> sphere_window_3x_data(3, std::vector(5, 0.0)); - qindex.map_data_to_sphere_grid_window(dd_3x_data, ghost_3x_data, sphere_window_3x_data, - sphere_center, center, dr_dtheta, bin_sizes, map_type, - normalize, bias); - std::vector gold_window_data; - std::vector> gold_window_3x_data(3); - // different result then 1D because the 1.0 y offset of the data - if (rtt_c4::node() == 0) { - gold_window_data = {0.0, 0.0, 3.0, 0.0, 0.0}; - gold_window_3x_data[0] = {0.0, 0.0, 3.0, 0.0, 0.0}; - gold_window_3x_data[1] = {0.0, 0.0, 4.0, 0.0, 0.0}; - gold_window_3x_data[2] = {0.0, 0.0, 0.0, 0.0, 0.0}; - } else if (rtt_c4::node() == 1) { - gold_window_data = {0.0, 0.0, 6.0, 0.0, 0.0}; - gold_window_3x_data[0] = {0.0, 0.0, 6.0, 0.0, 0.0}; - gold_window_3x_data[1] = {0.0, 0.0, 7.0, 0.0, 0.0}; - gold_window_3x_data[2] = {0.0, 0.0, 0.0, 0.0, 0.0}; - } else { - gold_window_data = {0.0, 0.0, 9.0, 0.0, 0.0}; - gold_window_3x_data[0] = {0.0, 0.0, 9.0, 0.0, 0.0}; - gold_window_3x_data[1] = {0.0, 0.0, 10.0, 0.0, 0.0}; - gold_window_3x_data[2] = {0.0, 0.0, 0.0, 0.0, 0.0}; - } - - for (size_t i = 0; i < bin_sizes[0]; i++) - if (!rtt_dsxx::soft_equiv(sphere_window_data[i], gold_window_data[i])) - ITFAILS; - for (size_t v = 0; v < 3; v++) - for (size_t i = 0; i < bin_sizes[0]; i++) - if (!rtt_dsxx::soft_equiv(sphere_window_3x_data[v][i], gold_window_3x_data[v][i])) - ITFAILS; - } - - // check max sphere r window mapping (more bins then data) normalized - { - // build a length=1.0 window around the first point on each node - const std::array center{dd_position_array[0][0], dd_position_array[0][1], 0.0}; - // wedge location +- 1 degree theta and 0.5 radius - const std::array dr_dtheta{0.4, 0.0174533, 0.0}; - const std::array bin_sizes{5, 1, 0}; - const bool normalize = true; - const bool bias = false; - const std::string map_type = "max"; - const std::array sphere_center{-1.0, 0.0, 0.0}; - std::vector sphere_window_data(5, 0.0); - qindex.map_data_to_sphere_grid_window(dd_data, ghost_data, sphere_window_data, sphere_center, - center, dr_dtheta, bin_sizes, map_type, normalize, - bias); - std::vector> sphere_window_3x_data(3, std::vector(5, 0.0)); - qindex.map_data_to_sphere_grid_window(dd_3x_data, ghost_3x_data, sphere_window_3x_data, - sphere_center, center, dr_dtheta, bin_sizes, map_type, - normalize, bias); - std::vector gold_window_data; - std::vector> gold_window_3x_data(3); - // different result then 1D because the 1.0 y offset of the data - if (rtt_c4::node() == 0) { - gold_window_data = {0.0, 0.0, 1.0, 0.0, 0.0}; - gold_window_3x_data[0] = {0.0, 0.0, 1.0, 0.0, 0.0}; - gold_window_3x_data[1] = {0.0, 0.0, 1.0, 0.0, 0.0}; - gold_window_3x_data[2] = {0.0, 0.0, 1.0, 0.0, 0.0}; - } else if (rtt_c4::node() == 1) { - gold_window_data = {0.0, 0.0, 1.0, 0.0, 0.0}; - gold_window_3x_data[0] = {0.0, 0.0, 1.0, 0.0, 0.0}; - gold_window_3x_data[1] = {0.0, 0.0, 1.0, 0.0, 0.0}; - gold_window_3x_data[2] = {0.0, 0.0, 1.0, 0.0, 0.0}; - } else { - gold_window_data = {0.0, 0.0, 1.0, 0.0, 0.0}; - gold_window_3x_data[0] = {0.0, 0.0, 1.0, 0.0, 0.0}; - gold_window_3x_data[1] = {0.0, 0.0, 1.0, 0.0, 0.0}; - gold_window_3x_data[2] = {0.0, 0.0, 1.0, 0.0, 0.0}; - } - - for (size_t i = 0; i < bin_sizes[0]; i++) - if (!rtt_dsxx::soft_equiv(sphere_window_data[i], gold_window_data[i])) - ITFAILS; - for (size_t v = 0; v < 3; v++) - for (size_t i = 0; i < bin_sizes[0]; i++) - if (!rtt_dsxx::soft_equiv(sphere_window_3x_data[v][i], gold_window_3x_data[v][i])) - ITFAILS; - } - - // check max sphere r window mapping (more bins then data) normalized bias fill - { - // build a length=1.0 window around the first point on each node - const std::array center{dd_position_array[0][0], dd_position_array[0][1], 0.0}; - // wedge location +- 1 degree theta and 0.5 radius - const std::array dr_dtheta{0.4, 0.0174533, 0.0}; - const std::array bin_sizes{5, 1, 0}; - const bool normalize = true; - const bool bias = true; - const std::string map_type = "min_fill"; - const std::array sphere_center{-1.0, 0.0, 0.0}; - std::vector sphere_window_data(5, 0.0); - qindex.map_data_to_sphere_grid_window(dd_data, ghost_data, sphere_window_data, sphere_center, - center, dr_dtheta, bin_sizes, map_type, normalize, - bias); - std::vector> sphere_window_3x_data(3, std::vector(5, 0.0)); - qindex.map_data_to_sphere_grid_window(dd_3x_data, ghost_3x_data, sphere_window_3x_data, - sphere_center, center, dr_dtheta, bin_sizes, map_type, - normalize, bias); - std::vector gold_window_data; - std::vector> gold_window_3x_data(3); - // different result then 1D because the 1.0 y offset of the data - if (rtt_c4::node() == 0) { - gold_window_data = {0.0, 0.0, 1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0}; - gold_window_3x_data[0] = {0.0, 0.0, 1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0}; - gold_window_3x_data[1] = {0.0, 0.0, 1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0}; - gold_window_3x_data[2] = {0.0, 0.0, 0.0, 0.0, 0.0}; - } else if (rtt_c4::node() == 1) { - gold_window_data = {0.0, 0.0, 1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0}; - gold_window_3x_data[0] = {0.0, 0.0, 1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0}; - gold_window_3x_data[1] = {0.0, 0.0, 1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0}; - gold_window_3x_data[2] = {0.0, 0.0, 0.0, 0.0, 0.0}; - } else { - gold_window_data = {0.0, 0.0, 1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0}; - gold_window_3x_data[0] = {0.0, 0.0, 1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0}; - gold_window_3x_data[1] = {0.0, 0.0, 1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0}; - gold_window_3x_data[2] = {0.0, 0.0, 0.0, 0.0, 0.0}; - } - - for (size_t i = 0; i < bin_sizes[0]; i++) - if (!rtt_dsxx::soft_equiv(sphere_window_data[i], gold_window_data[i])) - ITFAILS; - for (size_t v = 0; v < 3; v++) - for (size_t i = 0; i < bin_sizes[0]; i++) - if (!rtt_dsxx::soft_equiv(sphere_window_3x_data[v][i], gold_window_3x_data[v][i])) - ITFAILS; - } - - // check max sphere r window mapping (more bins then data) functions - { - // build a length=1.0 window around the first point on each node - const std::array center{dd_position_array[0][0], dd_position_array[0][1], 0.0}; - // wedge location +- 1 degree theta and 0.5 radius - const std::array dr_dtheta{0.4, 0.0174533, 0.0}; - const std::array bin_sizes{1, 5, 0}; - const bool normalize = false; - const bool bias = false; - const std::string map_type = "ave"; - const std::array sphere_center{-1.0, 0.0, 0.0}; - std::vector sphere_window_data(5, 0.0); - qindex.map_data_to_sphere_grid_window(dd_data, ghost_data, sphere_window_data, sphere_center, - center, dr_dtheta, bin_sizes, map_type, normalize, - bias); - std::vector> sphere_window_3x_data(3, std::vector(5, 0.0)); - qindex.map_data_to_sphere_grid_window(dd_3x_data, ghost_3x_data, sphere_window_3x_data, - sphere_center, center, dr_dtheta, bin_sizes, map_type, - normalize, bias); - std::vector gold_window_data; - std::vector> gold_window_3x_data(3); - // different result then 1D because the 1.0 y offset of the data - if (rtt_c4::node() == 0) { - gold_window_data = {0.0, 0.0, 3.0, 0.0, 0.0}; - gold_window_3x_data[0] = {0.0, 0.0, 3.0, 0.0, 0.0}; - gold_window_3x_data[1] = {0.0, 0.0, 4.0, 0.0, 0.0}; - gold_window_3x_data[2] = {0.0, 0.0, -3.0, 0.0, 0.0}; - } else if (rtt_c4::node() == 1) { - gold_window_data = {0.0, 0.0, 6.0, 0.0, 0.0}; - gold_window_3x_data[0] = {0.0, 0.0, 6.0, 0.0, 0.0}; - gold_window_3x_data[1] = {0.0, 0.0, 7.0, 0.0, 0.0}; - gold_window_3x_data[2] = {0.0, 0.0, -6.0, 0.0, 0.0}; - } else { - gold_window_data = {0.0, 0.0, 9.0, 0.0, 0.0}; - gold_window_3x_data[0] = {0.0, 0.0, 9.0, 0.0, 0.0}; - gold_window_3x_data[1] = {0.0, 0.0, 10.0, 0.0, 0.0}; - gold_window_3x_data[2] = {0.0, 0.0, -9.0, 0.0, 0.0}; - } - - for (size_t i = 0; i < bin_sizes[0]; i++) - if (!rtt_dsxx::soft_equiv(sphere_window_data[i], gold_window_data[i])) - ITFAILS; - for (size_t v = 0; v < 3; v++) - for (size_t i = 0; i < bin_sizes[0]; i++) - if (!rtt_dsxx::soft_equiv(sphere_window_3x_data[v][i], gold_window_3x_data[v][i])) - ITFAILS; - } - - // check max sphere r window mapping (more bins then data) functions - { - // build a length=1.0 window around the first point on each node - const std::array center{dd_position_array[0][0], dd_position_array[0][1], 0.0}; - // wedge location +- 1 degree theta and 0.5 radius - const std::array dr_dtheta{0.4, 0.0174533, 0.0}; - const std::array bin_sizes{1, 5, 0}; - const bool normalize = false; - const bool bias = false; - const std::string map_type = "nearest"; - const std::array sphere_center{-1.0, 0.0, 0.0}; - std::vector sphere_window_data(5, 0.0); - qindex.map_data_to_sphere_grid_window(dd_data, ghost_data, sphere_window_data, sphere_center, - center, dr_dtheta, bin_sizes, map_type, normalize, - bias); - std::vector> sphere_window_3x_data(3, std::vector(5, 0.0)); - qindex.map_data_to_sphere_grid_window(dd_3x_data, ghost_3x_data, sphere_window_3x_data, - sphere_center, center, dr_dtheta, bin_sizes, map_type, - normalize, bias); - std::vector gold_window_data; - std::vector> gold_window_3x_data(3); - // different result then 1D because the 1.0 y offset of the data - if (rtt_c4::node() == 0) { - gold_window_data = {0.0, 0.0, 3.0, 0.0, 0.0}; - gold_window_3x_data[0] = {0.0, 0.0, 3.0, 0.0, 0.0}; - gold_window_3x_data[1] = {0.0, 0.0, 4.0, 0.0, 0.0}; - gold_window_3x_data[2] = {0.0, 0.0, -3.0, 0.0, 0.0}; - } else if (rtt_c4::node() == 1) { - gold_window_data = {0.0, 0.0, 6.0, 0.0, 0.0}; - gold_window_3x_data[0] = {0.0, 0.0, 6.0, 0.0, 0.0}; - gold_window_3x_data[1] = {0.0, 0.0, 7.0, 0.0, 0.0}; - gold_window_3x_data[2] = {0.0, 0.0, -6.0, 0.0, 0.0}; - } else { - gold_window_data = {0.0, 0.0, 9.0, 0.0, 0.0}; - gold_window_3x_data[0] = {0.0, 0.0, 9.0, 0.0, 0.0}; - gold_window_3x_data[1] = {0.0, 0.0, 10.0, 0.0, 0.0}; - gold_window_3x_data[2] = {0.0, 0.0, -9.0, 0.0, 0.0}; - } - - for (size_t i = 0; i < bin_sizes[0]; i++) - if (!rtt_dsxx::soft_equiv(sphere_window_data[i], gold_window_data[i])) - ITFAILS; - for (size_t v = 0; v < 3; v++) - for (size_t i = 0; i < bin_sizes[0]; i++) - if (!rtt_dsxx::soft_equiv(sphere_window_3x_data[v][i], gold_window_3x_data[v][i])) - ITFAILS; - } - // check max_fill window mapping (more bins then data) functions { // build a length=1.0 window around the first point on each node @@ -1714,14 +1492,909 @@ void test_decomposition(ParallelUnitTest &ut) { } } -//------------------------------------------------------------------------------------------------// -int main(int argc, char *argv[]) { - ParallelUnitTest ut(argc, argv, release); - try { - // >>> UNIT TESTS - test_replication(ut); - if (nodes() == 3) - test_decomposition(ut); +void test_decomposition_sphere(ParallelUnitTest &ut) { + if (rtt_c4::nodes() != 3) + ITFAILS; + + //2d sphere + { + + const size_t local_size = 4; + const std::array sphere_center{0.0, 0.0, 0.0}; + const double max_radius = 1.0; + const double min_radius = 0.0; + const double shell_min_radius = 0.5; + const std::array radial_edges{0.5, 1.0}; + const std::array cosine_edges{-.99, 0, .99, -.99, 0, .99}; + const size_t data_size = radial_edges.size() * cosine_edges.size(); + std::vector> position_array(data_size, + std::array{0.0, 0.0, 0.0}); + + std::vector shell_data(data_size, 0.0); + std::vector spoke_data(data_size, 0.0); + size_t point_i = 0; + size_t mui = 0; + for (auto &mu : cosine_edges) { + size_t ri = 0; + for (auto &r : radial_edges) { + double sign = mui < 3 ? 1.0 : -1.0; + spoke_data[point_i] = static_cast(mui) + 1.0; + shell_data[point_i] = static_cast(ri) + 1.0; + double rel_y = r * mu; + position_array[point_i][0] = + rtt_dsxx::soft_equiv(r * r, rel_y * rel_y, 1e-6) ? 0.0 : sqrt(r * r - rel_y * rel_y); + position_array[point_i][0] *= sign; + position_array[point_i][1] = sphere_center[1] + rel_y; + point_i++; + ri++; + } + mui++; + } + + std::vector dd_data(local_size, 1.0); + std::vector> dd_3x_data(3, std::vector(local_size, 1.0)); + std::vector dd_spoke_data(local_size); + std::vector dd_shell_data(local_size); + std::vector> dd_position_array(local_size); + for (size_t i = 0; i < local_size; i++) { + dd_spoke_data[i] = spoke_data[i + rtt_c4::node() * local_size]; + dd_shell_data[i] = shell_data[i + rtt_c4::node() * local_size]; + dd_position_array[i] = position_array[i + rtt_c4::node() * local_size]; + } + + // in dd mode the max window size does nothing so set it large + const bool spherical = true; + const double max_window_size = 1.0; + const size_t bins_per_dim = 10UL; + const bool dd = true; + const size_t dim = 2; + quick_index qindex(dim, dd_position_array, max_window_size, bins_per_dim, dd, spherical, + sphere_center); + // Check the local state data + if (!qindex.domain_decomposed) + ITFAILS; + if (qindex.coarse_bin_resolution != bins_per_dim) + ITFAILS; + if (!soft_equiv(qindex.max_window_size, max_window_size)) + ITFAILS; + + // Check global bounding box + if (!soft_equiv(qindex.bounding_box_min[0], 0.5)) + ITFAILS; + if (!soft_equiv(qindex.bounding_box_min[1], 0.141539, 1e-4)) + ITFAILS; + if (!soft_equiv(qindex.bounding_box_min[2], 0.0)) + ITFAILS; + if (!soft_equiv(qindex.bounding_box_max[0], 1.0)) + ITFAILS; + if (!soft_equiv(qindex.bounding_box_max[1], 6.14165, 1e-4)) + ITFAILS; + if (!soft_equiv(qindex.bounding_box_max[2], 0.0)) + ITFAILS; + // Check local coarse_index map + // local indexing will not match the domain replicated case (different + // number of points per rank so different local indexing) + std::map> gold_map; + if (rtt_c4::node() == 0) { + gold_map[20] = {2}; + gold_map[29] = {3}; + gold_map[40] = {0}; + gold_map[49] = {1}; + } else if (rtt_c4::node() == 1) { + gold_map[0] = {0}; + gold_map[9] = {1}; + gold_map[50] = {2}; + gold_map[59] = {3}; + } else { + gold_map[70] = {0}; + gold_map[79] = {1}; + gold_map[90] = {2}; + gold_map[99] = {3}; + } + if (gold_map.size() != qindex.coarse_index_map.size()) + ITFAILS; + + for (auto &map : qindex.coarse_index_map) + for (size_t i = 0; i < map.second.size(); i++) + if (gold_map[map.first][i] != map.second[i]) + ITFAILS; + + // Check Domain Decomposed Data + // local bounding box extends beyond local data based on the window size + if (rtt_c4::node() == 0) { + if (!soft_equiv(qindex.local_bounding_box_min[0], 0.0)) + ITFAILS; + if (!soft_equiv(qindex.local_bounding_box_min[1], 1.23746, 1e-4)) + ITFAILS; + if (!soft_equiv(qindex.local_bounding_box_min[2], 0.0)) + ITFAILS; + if (!soft_equiv(qindex.local_bounding_box_max[0], 1.5)) + ITFAILS; + if (!soft_equiv(qindex.local_bounding_box_max[1], 3.33339, 1e-4)) + ITFAILS; + if (!soft_equiv(qindex.local_bounding_box_max[2], 0.0)) + ITFAILS; + } else if (rtt_c4::node() == 1) { + if (!soft_equiv(qindex.local_bounding_box_min[0], 0.0)) + ITFAILS; + // overlaps the theta=0=pi*2 boundary + if (!soft_equiv(qindex.local_bounding_box_min[1], -0.191794, 1e-4)) + ITFAILS; + if (!soft_equiv(qindex.local_bounding_box_min[2], 0.0)) + ITFAILS; + if (!soft_equiv(qindex.local_bounding_box_max[0], 1.5)) + ITFAILS; + if (!soft_equiv(qindex.local_bounding_box_max[1], 3.61647, 1e-4)) + ITFAILS; + if (!soft_equiv(qindex.local_bounding_box_max[2], 0.0)) + ITFAILS; + } else { + if (!soft_equiv(qindex.local_bounding_box_min[0], 0.0)) + ITFAILS; + if (!soft_equiv(qindex.local_bounding_box_min[1], 4.37906, 1e-4)) + ITFAILS; + if (!soft_equiv(qindex.local_bounding_box_min[2], 0.0)) + ITFAILS; + if (!soft_equiv(qindex.local_bounding_box_max[0], 1.5)) + ITFAILS; + // overlaps the theta=0=pi*2 boundary + if (!soft_equiv(qindex.local_bounding_box_max[1], 6.47498, 1e-4)) + ITFAILS; + if (!soft_equiv(qindex.local_bounding_box_max[2], 0.0)) + ITFAILS; + } + // global bins that span the local domains + std::vector gold_bins; + if (rtt_c4::node() == 0) { + gold_bins = {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, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59}; + } else if (rtt_c4::node() == 1) { + gold_bins = {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, 50, 51, 52, 53, + 54, 55, 56, 57, 58, 59, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99}; + } else { + gold_bins = {70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, + 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9}; + } + + if (gold_bins.size() != qindex.local_bins.size()) + ITFAILS; + for (size_t i = 0; i < qindex.local_bins.size(); i++) + if (gold_bins[i] != qindex.local_bins[i]) + ITFAILS; + + // local ghost index map (how to find general location of the ghost data) + std::map> gold_ghost_index_map; + if (rtt_c4::node() == 0) { + gold_ghost_index_map[50] = {0}; + gold_ghost_index_map[59] = {1}; + } else if (rtt_c4::node() == 1) { + gold_ghost_index_map[20] = {0}; + gold_ghost_index_map[29] = {1}; + gold_ghost_index_map[40] = {2}; + gold_ghost_index_map[49] = {3}; + gold_ghost_index_map[90] = {4}; + gold_ghost_index_map[99] = {5}; + } else { + gold_ghost_index_map[0] = {0}; + gold_ghost_index_map[9] = {1}; + } + + if (gold_ghost_index_map.size() != qindex.local_ghost_index_map.size()) + ITFAILS; + for (auto &map : qindex.local_ghost_index_map) { + if (gold_ghost_index_map[map.first].size() != map.second.size()) + ITFAILS; + for (size_t i = 0; i < map.second.size(); i++) { + if (map.second[i] != gold_ghost_index_map[map.first][i]) + ITFAILS; + } + } + + // Check the local ghost locations (this tangentially checks the private + // put_window_map which is used to build this local data). + std::vector> gold_ghost_locations; + if (rtt_c4::node() == 0) { + gold_ghost_locations = {{0.5, 3.28313, 0.0}, {1, 3.28313, 0.0}}; + } else if (rtt_c4::node() == 1) { + gold_ghost_locations = {{0.5, 1.5708, 0.0}, {1, 1.5708, 0.0}, {0.5, 3.00005, 0.0}, + {1, 3.00005, 0.0}, {0.5, 6.14165, 0.0}, {1, 6.14165, 0.0}}; + } else { + gold_ghost_locations = {{0.5, 0.141539, 0.0}, {1, 0.141539, 0.0}}; + } + + if (gold_ghost_locations.size() != qindex.local_ghost_locations.size()) + ITFAILS; + + for (size_t i = 0; i < qindex.local_ghost_locations.size(); i++) { + for (size_t d = 0; d < 2; d++) + if (!rtt_dsxx::soft_equiv(gold_ghost_locations[i][d], qindex.local_ghost_locations[i][d], + 1e-4)) + ITFAILS; + } + } + + //2d half sphere + { + + const size_t local_size = 24; + const std::array sphere_center{0.0, -1.0, 0.0}; + const double max_radius = 1.0; + const double min_radius = 0.0; + const double shell_min_radius = 0.5; + const std::array radial_edges{0.025, 0.05, 0.075, 0.1, 0.25, 0.5, 0.75, 1.0}; + const std::array cosine_edges{-1.0, -0.75, -0.5, -0.25, 0.0, 0.25, 0.5, 0.75, 1.0}; + const size_t data_size = radial_edges.size() * cosine_edges.size(); + std::vector> position_array(data_size, + std::array{0.0, 0.0, 0.0}); + + std::vector shell_data(data_size, 0.0); + std::vector spoke_data(data_size, 0.0); + size_t point_i = 0; + size_t ri = 0; + for (auto &r : radial_edges) { + size_t mui = 0; + for (auto &mu : cosine_edges) { + spoke_data[point_i] = static_cast(mui) + 1.0; + shell_data[point_i] = static_cast(ri) + 1.0; + double rel_y = r * mu; + position_array[point_i][0] = + rtt_dsxx::soft_equiv(r * r, rel_y * rel_y, 1e-6) ? 0.0 : sqrt(r * r - rel_y * rel_y); + position_array[point_i][1] = sphere_center[1] + rel_y; + point_i++; + mui++; + } + ri++; + } + + std::vector dd_data(local_size, 1.0); + std::vector> dd_2x_data(2, std::vector(local_size, 1.0)); + std::vector> dd_position_array(local_size); + for (size_t i = 0; i < local_size; i++) { + dd_data[i] = shell_data[i + rtt_c4::node() * local_size]; + dd_2x_data[0][i] = shell_data[i + rtt_c4::node() * local_size]; + dd_2x_data[1][i] = spoke_data[i + rtt_c4::node() * local_size]; + dd_position_array[i] = position_array[i + rtt_c4::node() * local_size]; + } + + // in dd mode the max window size does nothing so set it large + const bool spherical = true; + const double max_window_size = 1.0; + const size_t bins_per_dim = 10UL; + const bool dd = true; + const size_t dim = 2; + quick_index qindex(dim, dd_position_array, max_window_size, bins_per_dim, dd, spherical, + sphere_center); + // Check the local state data + if (!qindex.domain_decomposed) + ITFAILS; + if (qindex.coarse_bin_resolution != bins_per_dim) + ITFAILS; + if (!soft_equiv(qindex.max_window_size, max_window_size)) + ITFAILS; + + // Check global bounding box + if (!soft_equiv(qindex.bounding_box_min[0], 0.025)) + ITFAILS; + if (!soft_equiv(qindex.bounding_box_min[1], 0.0)) + ITFAILS; + if (!soft_equiv(qindex.bounding_box_min[2], 0.0)) + ITFAILS; + if (!soft_equiv(qindex.bounding_box_max[0], 1.0)) + ITFAILS; + if (!soft_equiv(qindex.bounding_box_max[1], rtt_units::PI, 1e-6)) + ITFAILS; + if (!soft_equiv(qindex.bounding_box_max[2], 0.0)) + ITFAILS; + // Check local coarse_index map + // local indexing will not match the domain replicated case (different + // number of points per rank so different local indexing) + std::map> gold_map; + if (rtt_c4::node() == 0) { + gold_map[0] = {8, 17}; + gold_map[20] = {7, 16}; + gold_map[30] = {6, 15}; + gold_map[40] = {5, 14, 23}; + gold_map[50] = {3, 4, 12, 13, 21, 22}; + gold_map[60] = {2, 11, 20}; + gold_map[70] = {1, 10, 19}; + gold_map[90] = {0, 9, 18}; + } else if (rtt_c4::node() == 1) { + gold_map[0] = {2, 11}; + gold_map[2] = {20}; + gold_map[20] = {1, 10}; + gold_map[22] = {19}; + gold_map[30] = {0, 9}; + gold_map[32] = {18}; + gold_map[40] = {8}; + gold_map[42] = {17}; + gold_map[50] = {6, 7}; + gold_map[52] = {15, 16}; + gold_map[60] = {5}; + gold_map[62] = {14}; + gold_map[64] = {23}; + gold_map[70] = {4}; + gold_map[72] = {13}; + gold_map[74] = {22}; + gold_map[90] = {3}; + gold_map[92] = {12}; + gold_map[94] = {21}; + } else { + gold_map[4] = {5}; + gold_map[7] = {14}; + gold_map[9] = {23}; + gold_map[24] = {4}; + gold_map[27] = {13}; + gold_map[29] = {22}; + gold_map[34] = {3}; + gold_map[37] = {12}; + gold_map[39] = {21}; + gold_map[44] = {2}; + gold_map[47] = {11}; + gold_map[49] = {20}; + gold_map[54] = {0, 1}; + gold_map[57] = {9, 10}; + gold_map[59] = {18, 19}; + gold_map[67] = {8}; + gold_map[69] = {17}; + gold_map[77] = {7}; + gold_map[79] = {16}; + gold_map[97] = {6}; + gold_map[99] = {15}; + } + if (gold_map.size() != qindex.coarse_index_map.size()) + ITFAILS; + for (auto &map : qindex.coarse_index_map) + for (size_t i = 0; i < map.second.size(); i++) + if (gold_map[map.first][i] != map.second[i]) + ITFAILS; + + // Check Domain Decomposed Data + // local bounding box extends beyond local data based on the window size + if (rtt_c4::node() == 0) { + if (!soft_equiv(qindex.local_bounding_box_min[0], 0.0)) + ITFAILS; + if (!soft_equiv(qindex.local_bounding_box_min[1], -0.5 / 0.575)) + ITFAILS; + if (!soft_equiv(qindex.local_bounding_box_min[2], 0.0)) + ITFAILS; + if (!soft_equiv(qindex.local_bounding_box_max[0], 0.575)) + ITFAILS; + if (!soft_equiv(qindex.local_bounding_box_max[1], rtt_units::PI + 0.5 / 0.575)) + ITFAILS; + if (!soft_equiv(qindex.local_bounding_box_max[2], 0.0)) + ITFAILS; + } else if (rtt_c4::node() == 1) { + if (!soft_equiv(qindex.local_bounding_box_min[0], 0.0)) + ITFAILS; + if (!soft_equiv(qindex.local_bounding_box_min[1], -0.5)) + ITFAILS; + if (!soft_equiv(qindex.local_bounding_box_min[2], 0.0)) + ITFAILS; + if (!soft_equiv(qindex.local_bounding_box_max[0], 1.0)) + ITFAILS; + if (!soft_equiv(qindex.local_bounding_box_max[1], rtt_units::PI + 0.5)) + ITFAILS; + if (!soft_equiv(qindex.local_bounding_box_max[2], 0.0)) + ITFAILS; + } else { + if (!soft_equiv(qindex.local_bounding_box_min[0], 0.0)) + ITFAILS; + if (!soft_equiv(qindex.local_bounding_box_min[1], -0.5 / 1.5)) + ITFAILS; + if (!soft_equiv(qindex.local_bounding_box_min[2], 0.0)) + ITFAILS; + if (!soft_equiv(qindex.local_bounding_box_max[0], 1.5)) + ITFAILS; + if (!soft_equiv(qindex.local_bounding_box_max[1], rtt_units::PI + 0.5 / 1.5)) + ITFAILS; + if (!soft_equiv(qindex.local_bounding_box_max[2], 0.0)) + ITFAILS; + } + // global bins that span the local domains + std::vector gold_bins; + if (rtt_c4::node() == 0) { + gold_bins = {0, 1, 2, 3, 4, 5, 10, 11, 12, 13, 14, 15, 20, 21, 22, 23, 24, 25, 30, 31, + 32, 33, 34, 35, 40, 41, 42, 43, 44, 45, 50, 51, 52, 53, 54, 55, 60, 61, 62, 63, + 64, 65, 70, 71, 72, 73, 74, 75, 80, 81, 82, 83, 84, 85, 90, 91, 92, 93, 94, 95}; + } else if (rtt_c4::node() == 1) { + gold_bins = {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, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, + 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, + 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99}; + } else { + gold_bins = {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, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, + 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, + 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99}; + } + + if (gold_bins.size() != qindex.local_bins.size()) + ITFAILS; + for (size_t i = 0; i < qindex.local_bins.size(); i++) + if (gold_bins[i] != qindex.local_bins[i]) + ITFAILS; + + // local ghost index map (how to find general location of the ghost data) + std::map> gold_ghost_index_map; + if (rtt_c4::node() == 0) { + gold_ghost_index_map[0] = {0, 1}; + gold_ghost_index_map[2] = {2}; + gold_ghost_index_map[4] = {24}; + gold_ghost_index_map[20] = {3, 4}; + gold_ghost_index_map[22] = {5}; + gold_ghost_index_map[24] = {25}; + gold_ghost_index_map[30] = {6, 7}; + gold_ghost_index_map[32] = {8}; + gold_ghost_index_map[34] = {26}; + gold_ghost_index_map[40] = {9}; + gold_ghost_index_map[42] = {10}; + gold_ghost_index_map[44] = {27}; + gold_ghost_index_map[50] = {11, 12}; + gold_ghost_index_map[52] = {13, 14}; + gold_ghost_index_map[54] = {28, 29}; + gold_ghost_index_map[60] = {15}; + gold_ghost_index_map[62] = {16}; + gold_ghost_index_map[64] = {17}; + gold_ghost_index_map[70] = {18}; + gold_ghost_index_map[72] = {19}; + gold_ghost_index_map[74] = {20}; + gold_ghost_index_map[90] = {21}; + gold_ghost_index_map[92] = {22}; + gold_ghost_index_map[94] = {23}; + } else if (rtt_c4::node() == 1) { + gold_ghost_index_map[0] = {0, 1}; + gold_ghost_index_map[4] = {24}; + gold_ghost_index_map[7] = {25}; + gold_ghost_index_map[9] = {26}; + gold_ghost_index_map[20] = {2, 3}; + gold_ghost_index_map[24] = {27}; + gold_ghost_index_map[27] = {28}; + gold_ghost_index_map[29] = {29}; + gold_ghost_index_map[30] = {4, 5}; + gold_ghost_index_map[34] = {30}; + gold_ghost_index_map[37] = {31}; + gold_ghost_index_map[39] = {32}; + gold_ghost_index_map[40] = {6, 7, 8}; + gold_ghost_index_map[44] = {33}; + gold_ghost_index_map[47] = {34}; + gold_ghost_index_map[49] = {35}; + gold_ghost_index_map[50] = {9, 10, 11, 12, 13, 14}; + gold_ghost_index_map[54] = {36, 37}; + gold_ghost_index_map[57] = {38, 39}; + gold_ghost_index_map[59] = {40, 41}; + gold_ghost_index_map[60] = {15, 16, 17}; + gold_ghost_index_map[67] = {42}; + gold_ghost_index_map[69] = {43}; + gold_ghost_index_map[70] = {18, 19, 20}; + gold_ghost_index_map[77] = {44}; + gold_ghost_index_map[79] = {45}; + gold_ghost_index_map[90] = {21, 22, 23}; + gold_ghost_index_map[97] = {46}; + gold_ghost_index_map[99] = {47}; + } else { + gold_ghost_index_map[0] = {0, 1, 24, 25}; + gold_ghost_index_map[2] = {26}; + gold_ghost_index_map[20] = {2, 3, 27, 28}; + gold_ghost_index_map[22] = {29}; + gold_ghost_index_map[30] = {4, 5, 30, 31}; + gold_ghost_index_map[32] = {32}; + gold_ghost_index_map[40] = {6, 7, 8, 33}; + gold_ghost_index_map[42] = {34}; + gold_ghost_index_map[50] = {9, 10, 11, 12, 13, 14, 35, 36}; + gold_ghost_index_map[52] = {37, 38}; + gold_ghost_index_map[60] = {15, 16, 17, 39}; + gold_ghost_index_map[62] = {40}; + gold_ghost_index_map[64] = {41}; + gold_ghost_index_map[70] = {18, 19, 20, 42}; + gold_ghost_index_map[72] = {43}; + gold_ghost_index_map[74] = {44}; + gold_ghost_index_map[90] = {21, 22, 23, 45}; + gold_ghost_index_map[92] = {46}; + gold_ghost_index_map[94] = {47}; + } + + if (gold_ghost_index_map.size() != qindex.local_ghost_index_map.size()) + ITFAILS; + for (auto &map : qindex.local_ghost_index_map) { + if (gold_ghost_index_map[map.first].size() != map.second.size()) + ITFAILS; + for (size_t i = 0; i < map.second.size(); i++) { + if (map.second[i] != gold_ghost_index_map[map.first][i]) + ITFAILS; + } + } + + // Check the local ghost locations (this tangentially checks the private + // put_window_map which is used to build this local data). + std::vector> gold_ghost_locations; + if (rtt_c4::node() == 0) { + gold_ghost_locations = { + {0.075, 0, 0.0}, {0.1, 0, 0.0}, {0.25, 0, 0.0}, {0.075, 0.722734, 0.0}, + {0.1, 0.722734, 0.0}, {0.25, 0.722734, 0.0}, {0.075, 1.0472, 0.0}, {0.1, 1.0472, 0.0}, + {0.25, 1.0472, 0.0}, {0.1, 1.31812, 0.0}, {0.25, 1.31812, 0.0}, {0.1, 1.82348, 0.0}, + {0.1, 1.5708, 0.0}, {0.25, 1.82348, 0.0}, {0.25, 1.5708, 0.0}, {0.1, 2.0944, 0.0}, + {0.25, 2.0944, 0.0}, {0.5, 2.0944, 0.0}, {0.1, 2.41886, 0.0}, {0.25, 2.41886, 0.0}, + {0.5, 2.41886, 0.0}, {0.1, 3.14159, 0.0}, {0.25, 3.14159, 0.0}, {0.5, 3.14159, 0.0}, + {0.5, 0, 0.0}, {0.5, 0.722734, 0.0}, {0.5, 1.0472, 0.0}, {0.5, 1.31812, 0.0}, + {0.5, 1.82348, 0.0}, {0.5, 1.5708, 0.0}}; + } else if (rtt_c4::node() == 1) { + gold_ghost_locations = {{0.025, 0, 0.0}, {0.05, 0, 0.0}, {0.025, 0.722734, 0.0}, + {0.05, 0.722734, 0.0}, {0.025, 1.0472, 0.0}, {0.05, 1.0472, 0.0}, + {0.025, 1.31812, 0.0}, {0.05, 1.31812, 0.0}, {0.075, 1.31812, 0.0}, + {0.025, 1.82348, 0.0}, {0.025, 1.5708, 0.0}, {0.05, 1.82348, 0.0}, + {0.05, 1.5708, 0.0}, {0.075, 1.82348, 0.0}, {0.075, 1.5708, 0.0}, + {0.025, 2.0944, 0.0}, {0.05, 2.0944, 0.0}, {0.075, 2.0944, 0.0}, + {0.025, 2.41886, 0.0}, {0.05, 2.41886, 0.0}, {0.075, 2.41886, 0.0}, + {0.025, 3.14159, 0.0}, {0.05, 3.14159, 0.0}, {0.075, 3.14159, 0.0}, + {0.5, 0, 0.0}, {0.75, 0, 0.0}, {1, 0, 0.0}, + {0.5, 0.722734, 0.0}, {0.75, 0.722734, 0.0}, {1, 0.722734, 0.0}, + {0.5, 1.0472, 0.0}, {0.75, 1.0472, 0.0}, {1, 1.0472, 0.0}, + {0.5, 1.31812, 0.0}, {0.75, 1.31812, 0.0}, {1, 1.31812, 0.0}, + {0.5, 1.82348, 0.0}, {0.5, 1.5708, 0.0}, {0.75, 1.82348, 0.0}, + {0.75, 1.5708, 0.0}, {1, 1.82348, 0.0}, {1, 1.5708, 0.0}, + {0.75, 2.0944, 0.0}, {1, 2.0944, 0.0}, {0.75, 2.41886, 0.0}, + {1, 2.41886, 0.0}, {0.75, 3.14159, 0.0}, {1, 3.14159, 0.0}}; + } else { + gold_ghost_locations = {{0.025, 0, 0.0}, {0.05, 0, 0.0}, {0.025, 0.722734, 0.0}, + {0.05, 0.722734, 0.0}, {0.025, 1.0472, 0.0}, {0.05, 1.0472, 0.0}, + {0.025, 1.31812, 0.0}, {0.05, 1.31812, 0.0}, {0.075, 1.31812, 0.0}, + {0.025, 1.82348, 0.0}, {0.025, 1.5708, 0.0}, {0.05, 1.82348, 0.0}, + {0.05, 1.5708, 0.0}, {0.075, 1.82348, 0.0}, {0.075, 1.5708, 0.0}, + {0.025, 2.0944, 0.0}, {0.05, 2.0944, 0.0}, {0.075, 2.0944, 0.0}, + {0.025, 2.41886, 0.0}, {0.05, 2.41886, 0.0}, {0.075, 2.41886, 0.0}, + {0.025, 3.14159, 0.0}, {0.05, 3.14159, 0.0}, {0.075, 3.14159, 0.0}, + {0.075, 0, 0.0}, {0.1, 0, 0.0}, {0.25, 0, 0.0}, + {0.075, 0.722734, 0.0}, {0.1, 0.722734, 0.0}, {0.25, 0.722734, 0.0}, + {0.075, 1.0472, 0.0}, {0.1, 1.0472, 0.0}, {0.25, 1.0472, 0.0}, + {0.1, 1.31812, 0.0}, {0.25, 1.31812, 0.0}, {0.1, 1.82348, 0.0}, + {0.1, 1.5708, 0.0}, {0.25, 1.82348, 0.0}, {0.25, 1.5708, 0.0}, + {0.1, 2.0944, 0.0}, {0.25, 2.0944, 0.0}, {0.5, 2.0944, 0.0}, + {0.1, 2.41886, 0.0}, {0.25, 2.41886, 0.0}, {0.5, 2.41886, 0.0}, + {0.1, 3.14159, 0.0}, {0.25, 3.14159, 0.0}, {0.5, 3.14159, 0.0}}; + } + if (gold_ghost_locations.size() != qindex.local_ghost_locations.size()) + ITFAILS; + for (size_t i = 0; i < qindex.local_ghost_locations.size(); i++) { + for (size_t d = 0; d < 2; d++) + if (!rtt_dsxx::soft_equiv(gold_ghost_locations[i][d], qindex.local_ghost_locations[i][d], + 1e-4)) + ITFAILS; + } + + // Check collect_ghost_data vector call + std::vector ghost_data(qindex.local_ghost_buffer_size, 0.0); + qindex.collect_ghost_data(dd_data, ghost_data); + std::vector> ghost_2x_data( + 2, std::vector(qindex.local_ghost_buffer_size, 0.0)); + qindex.collect_ghost_data(dd_2x_data, ghost_2x_data); + + std::vector gold_ghost_data; + std::vector> gold_2x_ghost_data(2); + if (rtt_c4::node() == 0) { + gold_ghost_data = {3, 4, 5, 3, 4, 5, 3, 4, 5, 4, 5, 4, 4, 5, 5, + 4, 5, 6, 4, 5, 6, 4, 5, 6, 6, 6, 6, 6, 6, 6}; + gold_2x_ghost_data[0] = {3, 4, 5, 3, 4, 5, 3, 4, 5, 4, 5, 4, 4, 5, 5, + 4, 5, 6, 4, 5, 6, 4, 5, 6, 6, 6, 6, 6, 6, 6}; + gold_2x_ghost_data[1] = {9, 9, 9, 8, 8, 8, 7, 7, 7, 6, 6, 4, 5, 4, 5, + 3, 3, 3, 2, 2, 2, 1, 1, 1, 9, 8, 7, 6, 4, 5}; + } else if (rtt_c4::node() == 1) { + gold_ghost_data = {1, 2, 1, 2, 1, 2, 1, 2, 3, 1, 1, 2, 2, 3, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, + 6, 7, 8, 6, 7, 8, 6, 7, 8, 6, 7, 8, 6, 6, 7, 7, 8, 8, 7, 8, 7, 8, 7, 8}; + gold_2x_ghost_data[0] = {1, 2, 1, 2, 1, 2, 1, 2, 3, 1, 1, 2, 2, 3, 3, 1, + 2, 3, 1, 2, 3, 1, 2, 3, 6, 7, 8, 6, 7, 8, 6, 7, + 8, 6, 7, 8, 6, 6, 7, 7, 8, 8, 7, 8, 7, 8, 7, 8}; + gold_2x_ghost_data[1] = {9, 9, 8, 8, 7, 7, 6, 6, 6, 4, 5, 4, 5, 4, 5, 3, + 3, 3, 2, 2, 2, 1, 1, 1, 9, 9, 9, 8, 8, 8, 7, 7, + 7, 6, 6, 6, 4, 5, 4, 5, 4, 5, 3, 3, 2, 2, 1, 1}; + } else { + gold_ghost_data = {1, 2, 1, 2, 1, 2, 1, 2, 3, 1, 1, 2, 2, 3, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, + 3, 4, 5, 3, 4, 5, 3, 4, 5, 4, 5, 4, 4, 5, 5, 4, 5, 6, 4, 5, 6, 4, 5, 6}; + gold_2x_ghost_data[0] = {1, 2, 1, 2, 1, 2, 1, 2, 3, 1, 1, 2, 2, 3, 3, 1, + 2, 3, 1, 2, 3, 1, 2, 3, 3, 4, 5, 3, 4, 5, 3, 4, + 5, 4, 5, 4, 4, 5, 5, 4, 5, 6, 4, 5, 6, 4, 5, 6}; + gold_2x_ghost_data[1] = {9, 9, 8, 8, 7, 7, 6, 6, 6, 4, 5, 4, 5, 4, 5, 3, + 3, 3, 2, 2, 2, 1, 1, 1, 9, 9, 9, 8, 8, 8, 7, 7, + 7, 6, 6, 4, 5, 4, 5, 3, 3, 3, 2, 2, 2, 1, 1, 1}; + } + + for (size_t i = 0; i < ghost_data.size(); i++) + if (!rtt_dsxx::soft_equiv(ghost_data[i], gold_ghost_data[i])) + ITFAILS; + for (size_t i = 0; i < ghost_2x_data[0].size(); i++) + if (!rtt_dsxx::soft_equiv(ghost_2x_data[0][i], gold_2x_ghost_data[0][i])) + ITFAILS; + for (size_t i = 0; i < ghost_2x_data[1].size(); i++) + if (!rtt_dsxx::soft_equiv(ghost_2x_data[1][i], gold_2x_ghost_data[1][i])) + ITFAILS; + + // check max sphere r window mapping (spoke shape more bins then data) functions + { + // wedge location +- 1 degree theta and 0.4 radius + const std::array dr_dtheta{0.4, 0.0174533, 0.0}; + // build a length=1.0 window around the first point on each node + const std::array center{qindex.locations[0]}; + const std::array win_min{center[0] - dr_dtheta[0], center[1] - dr_dtheta[1], 0.0}; + const std::array win_max{center[0] + dr_dtheta[0], center[1] + dr_dtheta[1], 0.0}; + const std::array bin_sizes{5, 1, 0}; + const bool normalize = false; + const bool bias = false; + const std::string map_type = "max"; + std::vector sphere_window_data(5, 0.0); + qindex.map_data_to_grid_window(dd_data, ghost_data, sphere_window_data, win_min, win_max, + bin_sizes, map_type, normalize, bias); + std::vector> sphere_window_2x_data(2, std::vector(5, 0.0)); + qindex.map_data_to_grid_window(dd_2x_data, ghost_2x_data, sphere_window_2x_data, win_min, + win_max, bin_sizes, map_type, normalize, bias); + std::vector gold_window_data; + std::vector> gold_window_2x_data(2); + // different result then 1D because the 1.0 y offset of the data + if (rtt_c4::node() == 0) { + gold_window_data = {0.0, 0.0, 4.0, 5.0, 0.0}; + gold_window_2x_data[0] = {0.0, 0.0, 4.0, 5.0, 0.0}; + gold_window_2x_data[1] = {0.0, 0.0, 1.0, 1.0, 0.0}; + } else if (rtt_c4::node() == 1) { + gold_window_data = {0.0, 0.0, 4.0, 5.0, 0.0}; + gold_window_2x_data[0] = {0.0, 0.0, 4.0, 5.0, 0.0}; + gold_window_2x_data[1] = {0.0, 0.0, 7.0, 7.0, 0.0}; + } else { + gold_window_data = {5.0, 0.0, 6.0, 0.0, 7.0}; + gold_window_2x_data[0] = {5.0, 0.0, 6.0, 0.0, 7.0}; + gold_window_2x_data[1] = {4.0, 0.0, 4.0, 0.0, 4.0}; + } + + for (size_t i = 0; i < bin_sizes[0]; i++) + if (!rtt_dsxx::soft_equiv(sphere_window_data[i], gold_window_data[i])) + ITFAILS; + for (size_t v = 0; v < 2; v++) + for (size_t i = 0; i < bin_sizes[0]; i++) + if (!rtt_dsxx::soft_equiv(sphere_window_2x_data[v][i], gold_window_2x_data[v][i])) + ITFAILS; + } + + // check nearest sphere r window mapping (spoke shape more bins then data) functions + { + // wedge location +- 30 degree theta and 0.4 radius + const std::array dr_dtheta{0.4, rtt_units::PI / 6.0, 0.0}; + // build a length=1.0 window around the first point on each node + const std::array center{qindex.locations[0]}; + const std::array win_min{center[0] - dr_dtheta[0], center[1] - dr_dtheta[1], 0.0}; + const std::array win_max{center[0] + dr_dtheta[0], center[1] + dr_dtheta[1], 0.0}; + const std::array bin_sizes{5, 1, 0}; + const bool normalize = false; + const bool bias = false; + const std::string map_type = "nearest"; + std::vector sphere_window_data(5, 0.0); + qindex.map_data_to_grid_window(dd_data, ghost_data, sphere_window_data, win_min, win_max, + bin_sizes, map_type, normalize, bias); + std::vector> sphere_window_2x_data(2, std::vector(5, 0.0)); + qindex.map_data_to_grid_window(dd_2x_data, ghost_2x_data, sphere_window_2x_data, win_min, + win_max, bin_sizes, map_type, normalize, bias); + std::vector gold_window_data; + std::vector> gold_window_2x_data(2); + // different result then 1D because the 1.0 y offset of the data + if (rtt_c4::node() == 0) { + gold_window_data = {0.0, 0.0, 1.0, 5.0, 0.0}; + gold_window_2x_data[0] = {0.0, 0.0, 1.0, 5.0, 0.0}; + gold_window_2x_data[1] = {0.0, 0.0, 1.0, 1.0, 0.0}; + } else if (rtt_c4::node() == 1) { + gold_window_data = {0.0, 0.0, 3.0, 5.0, 0.0}; + gold_window_2x_data[0] = {0.0, 0.0, 3.0, 5.0, 0.0}; + gold_window_2x_data[1] = {0.0, 0.0, 7.0, 7.0, 0.0}; + } else { + gold_window_data = {5.0, 0.0, 6.0, 0.0, 7.0}; + gold_window_2x_data[0] = {5.0, 0.0, 6.0, 0.0, 7.0}; + gold_window_2x_data[1] = {4.0, 0.0, 4.0, 0.0, 4.0}; + } + + for (size_t i = 0; i < bin_sizes[0]; i++) + if (!rtt_dsxx::soft_equiv(sphere_window_data[i], gold_window_data[i])) + ITFAILS; + for (size_t v = 0; v < 2; v++) + for (size_t i = 0; i < bin_sizes[0]; i++) + if (!rtt_dsxx::soft_equiv(sphere_window_2x_data[v][i], gold_window_2x_data[v][i])) + ITFAILS; + } + + // check max sphere r window mapping (shell shape more bins then data) functions + { + // wedge location +- 45 degree theta and 0.001 radius + const std::array dr_dtheta{0.001, rtt_units::PI / 4.0, 0.0}; + // build a length=1.0 window around the first point on each node + const std::array center{qindex.locations[0]}; + const std::array win_min{center[0] - dr_dtheta[0], center[1] - dr_dtheta[1], 0.0}; + const std::array win_max{center[0] + dr_dtheta[0], center[1] + dr_dtheta[1], 0.0}; + const std::array bin_sizes{1, 5, 0}; + const bool normalize = false; + const bool bias = false; + const std::string map_type = "max"; + std::vector sphere_window_data(5, 0.0); + qindex.map_data_to_grid_window(dd_data, ghost_data, sphere_window_data, win_min, win_max, + bin_sizes, map_type, normalize, bias); + std::vector> sphere_window_2x_data(2, std::vector(5, 0.0)); + qindex.map_data_to_grid_window(dd_2x_data, ghost_2x_data, sphere_window_2x_data, win_min, + win_max, bin_sizes, map_type, normalize, bias); + std::vector gold_window_data; + std::vector> gold_window_2x_data(2); + // different result then 1D because the 1.0 y offset of the data + if (rtt_c4::node() == 0) { + gold_window_data = {1.0, 0.0, 1.0, 0.0, 0.0}; + gold_window_2x_data[0] = {1.0, 0.0, 1.0, 0.0, 0.0}; + gold_window_2x_data[1] = {2.0, 0.0, 1.0, 0.0, 0.0}; + } else if (rtt_c4::node() == 1) { + gold_window_data = {0.0, 3.0, 3.0, 3.0, 3.0}; + gold_window_2x_data[0] = {0.0, 3.0, 3.0, 3.0, 3.0}; + gold_window_2x_data[1] = {0.0, 8.0, 7.0, 6.0, 5.0}; + } else { + gold_window_data = {6.0, 6.0, 6.0, 6.0, 6.0}; + gold_window_2x_data[0] = {6.0, 6.0, 6.0, 6.0, 6.0}; + gold_window_2x_data[1] = {7.0, 5.0, 4.0, 3.0, 2.0}; + } + + for (size_t i = 0; i < bin_sizes[1]; i++) + if (!rtt_dsxx::soft_equiv(sphere_window_data[i], gold_window_data[i])) + ITFAILS; + for (size_t v = 0; v < 2; v++) + for (size_t i = 0; i < bin_sizes[0]; i++) + if (!rtt_dsxx::soft_equiv(sphere_window_2x_data[v][i], gold_window_2x_data[v][i])) + ITFAILS; + } + + // check max sphere r window mapping (shell shape more bins then data) functions with fill + { + // wedge location +- 22.5 degree theta and 0.001 radius + const std::array dr_dtheta{0.001, rtt_units::PI / 4.0, 0.0}; + // build a length=1.0 window around the first point on each node + const std::array center{qindex.locations[0]}; + const std::array win_min{center[0] - dr_dtheta[0], center[1] - dr_dtheta[1], 0.0}; + const std::array win_max{center[0] + dr_dtheta[0], center[1] + dr_dtheta[1], 0.0}; + const std::array bin_sizes{1, 5, 0}; + const bool normalize = false; + const bool bias = false; + const std::string map_type = "max_fill"; + std::vector sphere_window_data(5, 0.0); + qindex.map_data_to_grid_window(dd_data, ghost_data, sphere_window_data, win_min, win_max, + bin_sizes, map_type, normalize, bias); + std::vector> sphere_window_2x_data(2, std::vector(5, 0.0)); + qindex.map_data_to_grid_window(dd_2x_data, ghost_2x_data, sphere_window_2x_data, win_min, + win_max, bin_sizes, map_type, normalize, bias); + std::vector gold_window_data; + std::vector> gold_window_2x_data(2); + // different result then 1D because the 1.0 y offset of the data + if (rtt_c4::node() == 0) { + gold_window_data = {1.0, 1.0, 1.0, 1.0, 1.0}; + gold_window_2x_data[0] = {1.0, 1.0, 1.0, 1.0, 1.0}; + gold_window_2x_data[1] = {2.0, 2.0, 1.0, 1.0, 1.0}; + } else if (rtt_c4::node() == 1) { + gold_window_data = {0.0, 3.0, 3.0, 3.0, 3.0}; + gold_window_2x_data[0] = {0.0, 3.0, 3.0, 3.0, 3.0}; + gold_window_2x_data[1] = {0.0, 8.0, 7.0, 6.0, 5.0}; + } else { + gold_window_data = {6.0, 6.0, 6.0, 6.0, 6.0}; + gold_window_2x_data[0] = {6.0, 6.0, 6.0, 6.0, 6.0}; + gold_window_2x_data[1] = {7.0, 5.0, 4.0, 3.0, 2.0}; + } + + for (size_t i = 0; i < bin_sizes[1]; i++) + if (!rtt_dsxx::soft_equiv(sphere_window_data[i], gold_window_data[i])) + ITFAILS; + for (size_t v = 0; v < 2; v++) + for (size_t i = 0; i < bin_sizes[0]; i++) + if (!rtt_dsxx::soft_equiv(sphere_window_2x_data[v][i], gold_window_2x_data[v][i])) + ITFAILS; + } + + // check min sphere r window mapping (shell shape more bins then data) functions with fill + { + // wedge location +- 22.5 degree theta and 0.001 radius + const std::array dr_dtheta{0.001, rtt_units::PI / 4.0, 0.0}; + // build a length=1.0 window around the first point on each node + const std::array center{qindex.locations[0]}; + const std::array win_min{center[0] - dr_dtheta[0], center[1] - dr_dtheta[1], 0.0}; + const std::array win_max{center[0] + dr_dtheta[0], center[1] + dr_dtheta[1], 0.0}; + const std::array bin_sizes{1, 5, 0}; + const bool normalize = false; + const bool bias = false; + const std::string map_type = "min_fill"; + std::vector sphere_window_data(5, 0.0); + qindex.map_data_to_grid_window(dd_data, ghost_data, sphere_window_data, win_min, win_max, + bin_sizes, map_type, normalize, bias); + std::vector> sphere_window_2x_data(2, std::vector(5, 0.0)); + qindex.map_data_to_grid_window(dd_2x_data, ghost_2x_data, sphere_window_2x_data, win_min, + win_max, bin_sizes, map_type, normalize, bias); + std::vector gold_window_data; + std::vector> gold_window_2x_data(2); + // different result then 1D because the 1.0 y offset of the data + if (rtt_c4::node() == 0) { + gold_window_data = {1.0, 1.0, 1.0, 1.0, 1.0}; + gold_window_2x_data[0] = {1.0, 1.0, 1.0, 1.0, 1.0}; + gold_window_2x_data[1] = {2.0, 2.0, 1.0, 1.0, 1.0}; + } else if (rtt_c4::node() == 1) { + gold_window_data = {0.0, 3.0, 3.0, 3.0, 3.0}; + gold_window_2x_data[0] = {0.0, 3.0, 3.0, 3.0, 3.0}; + gold_window_2x_data[1] = {0.0, 8.0, 7.0, 6.0, 5.0}; + } else { + gold_window_data = {6.0, 6.0, 6.0, 6.0, 6.0}; + gold_window_2x_data[0] = {6.0, 6.0, 6.0, 6.0, 6.0}; + gold_window_2x_data[1] = {6.0, 5.0, 4.0, 3.0, 2.0}; + } + + for (size_t i = 0; i < bin_sizes[1]; i++) + if (!rtt_dsxx::soft_equiv(sphere_window_data[i], gold_window_data[i])) + ITFAILS; + for (size_t v = 0; v < 2; v++) + for (size_t i = 0; i < bin_sizes[0]; i++) + if (!rtt_dsxx::soft_equiv(sphere_window_2x_data[v][i], gold_window_2x_data[v][i])) + ITFAILS; + } + + // check nearest sphere r window mapping (shell shape more bins then data) functions with fill + { + // wedge location +- 45 degree theta and 0.001 radius + const std::array dr_dtheta{0.001, rtt_units::PI / 4.0, 0.0}; + // build a length=1.0 window around the first point on each node + const std::array center{qindex.locations[0]}; + const std::array win_min{center[0] - dr_dtheta[0], center[1] - dr_dtheta[1], 0.0}; + const std::array win_max{center[0] + dr_dtheta[0], center[1] + dr_dtheta[1], 0.0}; + const std::array bin_sizes{1, 5, 0}; + const bool normalize = false; + const bool bias = false; + const std::string map_type = "nearest_fill"; + std::vector sphere_window_data(5, 0.0); + qindex.map_data_to_grid_window(dd_data, ghost_data, sphere_window_data, win_min, win_max, + bin_sizes, map_type, normalize, bias); + std::vector> sphere_window_2x_data(2, std::vector(5, 0.0)); + qindex.map_data_to_grid_window(dd_2x_data, ghost_2x_data, sphere_window_2x_data, win_min, + win_max, bin_sizes, map_type, normalize, bias); + std::vector gold_window_data; + std::vector> gold_window_2x_data(2); + // different result then 1D because the 1.0 y offset of the data + if (rtt_c4::node() == 0) { + gold_window_data = {1.0, 1.0, 1.0, 1.0, 1.0}; + gold_window_2x_data[0] = {1.0, 1.0, 1.0, 1.0, 1.0}; + gold_window_2x_data[1] = {2.0, 2.0, 1.0, 1.0, 1.0}; + } else if (rtt_c4::node() == 1) { + gold_window_data = {0.0, 3.0, 3.0, 3.0, 3.0}; + gold_window_2x_data[0] = {0.0, 3.0, 3.0, 3.0, 3.0}; + gold_window_2x_data[1] = {0.0, 8.0, 7.0, 6.0, 5.0}; + } else { + gold_window_data = {6.0, 6.0, 6.0, 6.0, 6.0}; + gold_window_2x_data[0] = {6.0, 6.0, 6.0, 6.0, 6.0}; + gold_window_2x_data[1] = {6.0, 5.0, 4.0, 3.0, 2.0}; + } + + for (size_t i = 0; i < bin_sizes[1]; i++) + if (!rtt_dsxx::soft_equiv(sphere_window_data[i], gold_window_data[i])) + ITFAILS; + for (size_t v = 0; v < 2; v++) + for (size_t i = 0; i < bin_sizes[0]; i++) + if (!rtt_dsxx::soft_equiv(sphere_window_2x_data[v][i], gold_window_2x_data[v][i])) + ITFAILS; + } + } + + if (ut.numFails == 0) { + PASSMSG("quick_index sphere DD checks pass"); + } else { + FAILMSG("quick_index sphere DD checks failed"); + } +} + +//------------------------------------------------------------------------------------------------// +int main(int argc, char *argv[]) { + ParallelUnitTest ut(argc, argv, release); + try { + // >>> UNIT TESTS + test_replication(ut); + test_replication_sphere(ut); + if (nodes() == 3) { + test_decomposition(ut); + test_decomposition_sphere(ut); + } } UT_EPILOG(ut); } From e06a99d2539d0f95138ce6e09faa728610367a2d Mon Sep 17 00:00:00 2001 From: Cleveland Date: Thu, 26 Aug 2021 15:52:18 -0600 Subject: [PATCH 2/4] regression cleanup --- src/kde/quick_index.cc | 30 ++++++++++++++---------------- src/kde/quick_index.hh | 3 ++- src/kde/test/tstkde.cc | 5 ----- src/kde/test/tstquick_index.cc | 16 ++++++++++------ 4 files changed, 26 insertions(+), 28 deletions(-) diff --git a/src/kde/quick_index.cc b/src/kde/quick_index.cc index 7c6b6f28f..903c4ff94 100644 --- a/src/kde/quick_index.cc +++ b/src/kde/quick_index.cc @@ -30,6 +30,8 @@ namespace rtt_kde { * \param[in] max_window_size_ maximum supported window size * \param[in] bins_per_dimension_ number of bins in each dimension * \param[in] domain_decomposed_ + * \param[in] spherical_ bool operator to enable spherical transform + * \param[in] sphere_center_ origin of spherical transform */ quick_index::quick_index(const size_t dim_, const std::vector> &locations_, const double max_window_size_, const size_t bins_per_dimension_, @@ -563,8 +565,7 @@ auto get_window_bin = [](auto spherical, const auto dim, const auto &grid_bins, //------------------------------------------------------------------------------------------------// // Lambda for mapping the data auto map_data = [](auto &bias_cell_count, auto &data_count, auto &grid_data, auto &min_distance, - const auto &dim, const auto &map_type, const auto &data, - const auto &distance_to_bin_center, const auto &location, + const auto &map_type, const auto &data, const auto &distance_to_bin_center, const auto &local_window_bin, const auto &data_bin) { // regardless of map type if it is the first value to enter the bin it // gets set to that value @@ -704,8 +705,8 @@ void quick_index::map_data_to_grid_window( continue; // lambda for mapping the data - map_data(bias_cell_count, data_count, grid_data, min_distance, dim, map_type, local_data, - distance_to_bin_center, locations[l], local_window_bin, l); + map_data(bias_cell_count, data_count, grid_data, min_distance, map_type, local_data, + distance_to_bin_center, local_window_bin, l); } // end local point loop } // if valid local bin loop @@ -727,8 +728,8 @@ void quick_index::map_data_to_grid_window( continue; // lambda for mapping the data - map_data(bias_cell_count, data_count, grid_data, min_distance, dim, map_type, ghost_data, - distance_to_bin_center, local_ghost_locations[g], local_window_bin, g); + map_data(bias_cell_count, data_count, grid_data, min_distance, map_type, ghost_data, + distance_to_bin_center, local_window_bin, g); } // end ghost point loop } // if valid ghost bin } // if dd @@ -784,10 +785,9 @@ void quick_index::map_data_to_grid_window( //------------------------------------------------------------------------------------------------// // Lambda for mapping the vector data auto map_vector_data = [](auto &bias_cell_count, auto &data_count, auto &grid_data, - auto &min_distance, const auto &dim, const auto &map_type, - const auto &data, const auto &distance_to_bin_center, - const auto &location, const auto &local_window_bin, const auto &data_bin, - const auto &vsize) { + auto &min_distance, const auto &map_type, const auto &data, + const auto &distance_to_bin_center, const auto &local_window_bin, + const auto &data_bin, const auto &vsize) { // regardless of map type if it is the first value to enter the bin it gets set to that value if (data_count[local_window_bin] == 0) { bias_cell_count += 1.0; @@ -938,9 +938,8 @@ void quick_index::map_data_to_grid_window(const std::vector> if (!valid) continue; Check(local_window_bin < n_map_bins); - map_vector_data(bias_cell_count, data_count, grid_data, min_distance, dim, map_type, - local_data, distance_to_bin_center, locations[l], local_window_bin, l, - vsize); + map_vector_data(bias_cell_count, data_count, grid_data, min_distance, map_type, local_data, + distance_to_bin_center, local_window_bin, l, vsize); } // end local point loop } // if valid local bin loop if (domain_decomposed) { @@ -959,9 +958,8 @@ void quick_index::map_data_to_grid_window(const std::vector> // If the bin is outside the window continue to the next poin if (!valid) continue; - map_vector_data(bias_cell_count, data_count, grid_data, min_distance, dim, map_type, - ghost_data, distance_to_bin_center, local_ghost_locations[g], - local_window_bin, g, vsize); + map_vector_data(bias_cell_count, data_count, grid_data, min_distance, map_type, + ghost_data, distance_to_bin_center, local_window_bin, g, vsize); } // end ghost point loop } // if valid ghost bin } // if dd diff --git a/src/kde/quick_index.hh b/src/kde/quick_index.hh index 93bd635dc..9208d42cf 100644 --- a/src/kde/quick_index.hh +++ b/src/kde/quick_index.hh @@ -28,6 +28,7 @@ namespace rtt_kde { * Calculate a relative r theta and phi coordinate relative to a sphere center location from a * standard (x,y,z) or (r,z) coordinates * + * \param[in] dim used to ensure it is only used in valid dimension ranges * \param[in] sphere_center center of sphere in (x,y,z) or (r,z) coordinates * \param[in] locations (x,y,z) or (r,z) locations to transform to relative (r, theta, phi) space. * @@ -110,7 +111,7 @@ public: // Quick index initialization data const size_t dim; const bool domain_decomposed; - const double spherical; + const bool spherical; const std::array sphere_center; const size_t coarse_bin_resolution; const double max_window_size; diff --git a/src/kde/test/tstkde.cc b/src/kde/test/tstkde.cc index c546d7011..3ddeac126 100644 --- a/src/kde/test/tstkde.cc +++ b/src/kde/test/tstkde.cc @@ -34,8 +34,6 @@ void test_replication(ParallelUnitTest &ut) { { const bool spherical = true; const std::array sphere_center{0.0, -1.0, 0.0}; - const double max_radius = 1.0; - const double min_radius = 0.0; kde sphere_kde; // reflect on the theta boundary kde theta_reflected_sphere_kde({false, false, true, true, false, false}); @@ -1046,9 +1044,6 @@ void test_decomposition(ParallelUnitTest &ut) { const bool spherical = true; const size_t local_size = 24; const std::array sphere_center{0.0, -1.0, 0.0}; - const double max_radius = 1.0; - const double min_radius = 0.0; - const double shell_min_radius = 0.5; kde sphere_kde; // reflect on the theta boundary kde theta_reflected_sphere_kde({false, false, true, true, false, false}); diff --git a/src/kde/test/tstquick_index.cc b/src/kde/test/tstquick_index.cc index 3360e23e8..2228e6ffd 100644 --- a/src/kde/test/tstquick_index.cc +++ b/src/kde/test/tstquick_index.cc @@ -76,6 +76,11 @@ void test_replication(ParallelUnitTest &ut) { for (size_t i = 0; i < map.second.size(); i++) if (gold_map[map.first][i] != map.second[i]) ITFAILS; + + // Check non-spherical orthogonal distance calculation + auto distance = qindex.calc_orthogonal_distance({-1, -1, -1}, {1, 1, 1}, 10.0); + for (auto &val : distance) + FAIL_IF_NOT(rtt_dsxx::soft_equiv(val, 2.0)); } if (ut.numFails == 0) { @@ -142,6 +147,11 @@ void test_replication_sphere(ParallelUnitTest &ut) { for (size_t i = 0; i < map.second.size(); i++) if (gold_map[map.first][i] != map.second[i]) ITFAILS; + + // Check non-spherical orthogonal distance calculation + auto distance = qindex.calc_orthogonal_distance({-1, 0.5, -1}, {1, 1, 1}, 4.0); + for (auto &val : distance) + FAIL_IF_NOT(rtt_dsxx::soft_equiv(val, 2.0)); } if (ut.numFails == 0) { @@ -1501,9 +1511,6 @@ void test_decomposition_sphere(ParallelUnitTest &ut) { const size_t local_size = 4; const std::array sphere_center{0.0, 0.0, 0.0}; - const double max_radius = 1.0; - const double min_radius = 0.0; - const double shell_min_radius = 0.5; const std::array radial_edges{0.5, 1.0}; const std::array cosine_edges{-.99, 0, .99, -.99, 0, .99}; const size_t data_size = radial_edges.size() * cosine_edges.size(); @@ -1721,9 +1728,6 @@ void test_decomposition_sphere(ParallelUnitTest &ut) { const size_t local_size = 24; const std::array sphere_center{0.0, -1.0, 0.0}; - const double max_radius = 1.0; - const double min_radius = 0.0; - const double shell_min_radius = 0.5; const std::array radial_edges{0.025, 0.05, 0.075, 0.1, 0.25, 0.5, 0.75, 1.0}; const std::array cosine_edges{-1.0, -0.75, -0.5, -0.25, 0.0, 0.25, 0.5, 0.75, 1.0}; const size_t data_size = radial_edges.size() * cosine_edges.size(); From 1dc4122eae4a546576fd8da2c94e423f94139cfd Mon Sep 17 00:00:00 2001 From: Cleveland Date: Thu, 26 Aug 2021 15:54:31 -0600 Subject: [PATCH 3/4] locations can remain non-mutable --- src/kde/quick_index.hh | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/kde/quick_index.hh b/src/kde/quick_index.hh index 9208d42cf..c39cc3144 100644 --- a/src/kde/quick_index.hh +++ b/src/kde/quick_index.hh @@ -115,8 +115,7 @@ public: const std::array sphere_center; const size_t coarse_bin_resolution; const double max_window_size; - // keep a mutable copy of the locations (can be change to spherical) - std::vector> locations; + const std::vector> locations; const size_t n_locations; // Global bounds From b41b77ca5ed8876c90ab23e6cf22dcd18006b160 Mon Sep 17 00:00:00 2001 From: Cleveland Date: Thu, 26 Aug 2021 18:29:09 -0600 Subject: [PATCH 4/4] fix overlap sign error --- src/kde/quick_index.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/kde/quick_index.cc b/src/kde/quick_index.cc index 903c4ff94..6bde89d45 100644 --- a/src/kde/quick_index.cc +++ b/src/kde/quick_index.cc @@ -529,11 +529,11 @@ auto get_window_bin = [](auto spherical, const auto dim, const auto &grid_bins, double loc = location[d]; // transform location for zero theta overshoot if (spherical && d == 1 && window_max[d] > 2.0 * rtt_units::PI && - location[d] < window_max[d] - 2.0 * rtt_units::PI) + location[d] < (window_max[d] - 2.0 * rtt_units::PI)) loc += 2.0 * rtt_units::PI; // transform location for zero theta overshoot if (spherical && d == 1 && window_min[d] < 0 && - location[d] > 2.0 * rtt_units::PI - window_min[d]) + location[d] > (2.0 * rtt_units::PI + window_min[d])) loc -= 2.0 * rtt_units::PI; const double bin_value = static_cast(grid_bins[d]) * (loc - window_min[d]) / (window_max[d] - window_min[d]);