Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Simplify spherical quick_index and kde reconstruction #1123

Merged
merged 4 commits into from
Aug 27, 2021
Merged

Simplify spherical quick_index and kde reconstruction #1123

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
35c3ff9
rely on quick_index for spherical reconstruction
Aug 26, 2021
e06a99d
regression cleanup
Aug 26, 2021
1dc4122
locations can remain non-mutable
Aug 26, 2021
b41b77c
fix overlap sign error
Aug 27, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
135 changes: 27 additions & 108 deletions src/kde/kde.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
*/
//------------------------------------------------------------------------------------------------//

Expand All @@ -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
Expand All @@ -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<double, 3> &r0,
const std::array<double, 3> &one_over_h0,
const std::array<double, 3> &r,
const std::array<double, 3> &one_over_h,
const quick_index &qindex,
const double &discontinuity_cutoff) const {
double kde::calc_weight(const std::array<double, 3> &r0, const std::array<double, 3> &one_over_h0,
const std::array<double, 3> &r, const std::array<double, 3> &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<double, 3> distance = qindex.calc_orthogonal_distance(r0, r, arch_radius);
std::array<double, 3> low_reflect_r0_distance =
qindex.calc_orthogonal_distance(qindex.bounding_box_min, r0, arch_radius);
std::array<double, 3> low_reflect_r_distance =
qindex.calc_orthogonal_distance(qindex.bounding_box_min, r, arch_radius);
std::array<double, 3> high_reflect_r0_distance =
qindex.calc_orthogonal_distance(r0, qindex.bounding_box_max, arch_radius);
std::array<double, 3> 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
Expand All @@ -71,14 +83,12 @@ double kde::calc_cartesian_weight(const std::array<double, 3> &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];
Expand All @@ -87,78 +97,6 @@ double kde::calc_cartesian_weight(const std::array<double, 3> &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<double, 3> &r0,
const std::array<double, 3> &one_over_h0,
const std::array<double, 3> &r,
const std::array<double, 3> &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
Expand Down Expand Up @@ -465,28 +403,9 @@ void kde::calc_win_min_max(const quick_index &qindex, const std::array<double, 3
Require(one_over_bandwidth[0] > 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];
}
}

Expand Down
33 changes: 2 additions & 31 deletions src/kde/kde.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
*/
//------------------------------------------------------------------------------------------------//

Expand Down Expand Up @@ -60,38 +60,14 @@ 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<double, 3> &center,
const std::array<double, 3> &location) const;

//! Calculate Arch Length between two locations at a specific radius
inline double calc_arch_length(const std::array<double, 3> &center, const double radius,
const std::array<double, 3> &location_1,
const std::array<double, 3> &location_2) const;
//! Setup a spherical reconstruction with a
void set_sphere_center(const std::array<double, 3> &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

private:
//! Private function to calculate kernel weight
double calc_weight(const std::array<double, 3> &r0, const std::array<double, 3> &one_over_h0,
const std::array<double, 3> &r, const std::array<double, 3> &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<double, 3> &position,
const std::array<double, 3> &one_over_bandwidth, std::array<double, 3> &min,
Expand All @@ -112,11 +88,6 @@ private:
// DATA
//! reflecting boundary conditions [lower_x, upper_x, lower_y, upper_y, lower_z, upper_z]
const std::array<bool, 6> reflect_boundary;
//! Spherical Mesh Reconstruction Data
std::array<double, 3> 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
Expand Down
67 changes: 1 addition & 66 deletions src/kde/kde.i.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
*/
//------------------------------------------------------------------------------------------------//

Expand Down Expand Up @@ -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<double, 3>{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<double, 3>{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<double, 3> &center,
const std::array<double, 3> &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<double, 3> &center, const double radius,
const std::array<double, 3> &location_1,
const std::array<double, 3> &location_2) const {
const std::array<double, 3> v1{calc_unit_vec(calc_vec(center, location_1))};
const std::array<double, 3> 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
Expand Down
Loading