Skip to content

Commit

Permalink
Fix two dumb bugs that my prior fixes added.
Browse files Browse the repository at this point in the history 
- The result of unnecessary resolution from point to major index was passed as a
point index.

- For some mysterious reason, I had the NaryTreeFacade::children() sorting on
child pointer value.

A few changes along the way:

- Make log pointer optional to the sanity() methods.

- Add sanity checking of the facade blob and scoped node ordering, now that they
are correctly, if accidentally/constructed to be the same.

- Rewrite get_closest_blob() to do slightly less work to get its job done.

- Add cluster length as most significant value to the mix of values when
comparing the "size" of two clusters.

- Rename test cases in clustering prototype test to better follow convention to
make the tests more user friendly for "-tc='clustering*'" type doctest executing.
  • Loading branch information
@brettviren
brettviren committed May 17, 2024
1 parent ddee41c commit 62240a3
Show file tree
Hide file tree
Showing 6 changed files with 82 additions and 76 deletions.
6 changes: 3 additions & 3 deletions img/inc/WireCellImg/PointCloudFacade.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -67,7 +67,7 @@ namespace WireCell::PointCloud::Facade {
size_t nbpoints() const;

// Check facade consistency
bool sanity(Log::logptr_t log) const;
bool sanity(Log::logptr_t log = nullptr) const;

private:

Expand Down Expand Up @@ -164,7 +164,7 @@ namespace WireCell::PointCloud::Facade {
using points_type = kd3d_t::points_type;
const points_type& points() const;

// Return charge-weighted average position of points of blobs withing distance of point.
// Return charge-weighted average position of points of blobs within distance of point.
geo_point_t calc_ave_pos(const geo_point_t& origin, const double dis) const;


Expand Down Expand Up @@ -259,7 +259,7 @@ namespace WireCell::PointCloud::Facade {
size_t hash() const;

// Check facade consistency between blob view and k-d tree view.
bool sanity(Log::logptr_t log) const;
bool sanity(Log::logptr_t log = nullptr) const;

private:
// start slice index (tick number) to blob facade pointer can be
Expand Down
118 changes: 63 additions & 55 deletions img/src/PointCloudFacade.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -135,7 +135,7 @@ size_t Blob::nbpoints() const
bool Blob::sanity(Log::logptr_t log) const
{
if (nbpoints() == (size_t)npoints()) return true;
log->debug("blob sanity: blob points mismatch: {}", *this);
if (log) log->debug("blob sanity: blob points mismatch: {}", *this);
return false;
}

Expand Down Expand Up @@ -424,9 +424,7 @@ Blob* Cluster::blob_with_point(size_t point_index)
const Blob* Cluster::blob_with_point(size_t point_index) const
{
const auto& sv = sv3d();
const auto& skd = sv.kd();
const size_t mai = skd.major_index(point_index);
const auto* node = sv.node_with_point(mai);
const auto* node = sv.node_with_point(point_index);
return node->value.facade<Blob>();
}

Expand All @@ -439,11 +437,10 @@ std::vector<const Blob*> Cluster::blobs_with_points(const kd_results_t& res) con
const size_t npts = res.size();
std::vector<const Blob*> ret(npts);
const auto& sv = sv3d();
const auto& skd = sv.kd();

for (size_t ind=0; ind<npts; ++ind) {
const size_t mai = skd.major_index(res[ind].first);
const auto* node = sv.node_with_point(mai);
const size_t point_index = res[ind].first;
const auto* node = sv.node_with_point(point_index);
ret[ind] = node->value.facade<Blob>();
}
return ret;
Expand All @@ -452,33 +449,28 @@ std::vector<const Blob*> Cluster::blobs_with_points(const kd_results_t& res) con

std::map<const Blob*, geo_point_t> Cluster::get_closest_blob(const geo_point_t& point, double radius) const
{
std::map<const Blob*, geo_point_t> blob_points; // return
std::map<const Blob*, double> blob_metrics; // temp

auto results = kd_radius(radius, point);
const auto& blobs = blobs_with_points(results);

const size_t npts = results.size();
for (size_t ind=0; ind<npts; ++ind) {
const auto& [point_index, metric] = results[ind];

const geo_point_t p1 = point3d(point_index);
const Blob* blob = blobs[ind];

auto it = blob_metrics.find(blob);
if (it == blob_metrics.end()) { // first time
blob_metrics[blob] = metric;
blob_points[blob] = p1;
struct Best { size_t point_index; double metric; };
std::unordered_map<size_t, Best> best_blob_point;

const auto& kd = kd3d();
auto results = kd.radius(radius*radius, point);
for (const auto& [point_index, metric] : results) {
const size_t major_index = kd.major_index(point_index);
auto it = best_blob_point.find(major_index);
if (it == best_blob_point.end()) { // first time seen
best_blob_point[major_index] = {point_index, metric};
continue;
}
else {
if (metric < it->second) { // check for yet closer
it->second = metric;
blob_points[blob] = p1;
}
if (metric < it->second.metric) {
it->second.point_index = point_index;
it->second.metric = metric;
}
}

return blob_points;
std::map<const Blob*, geo_point_t> ret;
for (const auto& [mi, bb] : best_blob_point) {
ret[blob_with_point(bb.point_index)] = point3d(bb.point_index);
}
return ret;
}

std::pair<geo_point_t, const Blob* > Cluster::get_closest_point_blob(const geo_point_t& point) const
Expand All @@ -499,7 +491,7 @@ geo_point_t Cluster::calc_ave_pos(const geo_point_t& origin, const double dis) c
double charge = 0;

auto blob_pts = get_closest_blob(origin, dis);
for (auto [blob, pt] : blob_pts) {
for (auto [blob, _] : blob_pts) {
double q = blob->charge();
if (q==0) q=1;
ret += blob->center_pos()*q;
Expand Down Expand Up @@ -770,29 +762,27 @@ std::string Facade::dump(const Facade::Grouping& grouping, int level)
}




bool Cluster::sanity(Log::logptr_t log) const
{
{
const auto* svptr = m_node->value.get_scoped(scope);
if (!svptr) {
log->debug("cluster sanity: note, not yet a scoped view {}", scope);
if (log) log->debug("cluster sanity: note, not yet a scoped view {}", scope);
}
}
if (!nchildren()) {
log->debug("cluster sanity: no children blobs");
if (log) log->debug("cluster sanity: no children blobs");
return false;
}

const auto& sv = m_node->value.scoped_view(scope);
const auto& snodes = sv.nodes();
if (snodes.empty()) {
log->debug("cluster sanity: no scoped nodes");
if (log) log->debug("cluster sanity: no scoped nodes");
return false;
}
if (sv.npoints() == 0) { // triggers a scoped view cache fill
log->debug("cluster sanity: no scoped points");
if (log) log->debug("cluster sanity: no scoped points");
return false;
}
// sv.force_invalid();
Expand All @@ -805,18 +795,35 @@ bool Cluster::sanity(Log::logptr_t log) const
}

if (skd.nblocks() != snodes.size()) {
log->debug("cluster sanity: k-d blocks={} scoped nodes={}", skd.nblocks(), fblobs.size());
if (log) log->debug("cluster sanity: k-d blocks={} scoped nodes={}", skd.nblocks(), fblobs.size());
return false;
}


/// Note, we do not expect k-d blocks order to be sync'ed to the Cluster
/// facade children blob order, but they should have equal number.
if (skd.nblocks() != fblobs.size()) {
log->debug("cluster sanity: k-d blocks={} cluster blobs={}", skd.nblocks(), fblobs.size());
if (log) log->debug("cluster sanity: k-d blocks={} cluster blobs={}", skd.nblocks(), fblobs.size());
return false;
}

for (size_t ind=0; ind<snodes.size(); ++ind) {
/// In general, the depth-first order of scoped view nodes is always the
/// same as k-d tree blocks but is not (again in general) expected to be
/// the same order as the children blobs of a parent cluster. After
/// all, a scoped view may span multiple essentially any subset of tree
/// nodes. However, in the special case of the "3d" SV and how the PC
/// tree is constructed, the depth-first and children blobs ordering
/// should be "accidentally" the same.
const auto* fblob = fblobs[ind];
const auto* sblob = snodes[ind]->value.facade<Blob>();
if (fblob != sblob) {
if (log) {
log->debug("cluster sanity: scoped node facade Blob differs from cluster child at {}", ind);
log->debug("cluster sanity: \tscoped blob: {}", *fblob);
log->debug("cluster sanity: \tfacade blob: {}", *sblob);
}
// return false;
}
}

const auto& majs = skd.major_indices();
const auto& mins = skd.minor_indices();
const size_t npts = skd.npoints();
Expand All @@ -833,7 +840,7 @@ bool Cluster::sanity(Log::logptr_t log) const
// scoped consistency
const node_t* tnode = sv.node_with_point(ind);
if (!tnode) {
log->debug("cluster sanity: scoped node facade not a Blob at majind={}", majind);
if (log) log->debug("cluster sanity: scoped node facade not a Blob at majind={}", majind);
return false;
}
const auto* tblob = tnode->value.facade<Blob>();
Expand All @@ -843,16 +850,17 @@ bool Cluster::sanity(Log::logptr_t log) const
}

if (minind >= spoints.size()) {
log->debug("cluster sanity: minind={} is beyond scoped blob npts={} majind={}, blob is: {}",
minind, spoints.size(), majind, *sblob);
if (log) log->debug("cluster sanity: minind={} is beyond scoped blob npts={} majind={}, blob is: {}",
minind, spoints.size(), majind, *sblob);
return false;
}
auto spt = spoints[minind];
if (spt != kdpt) {
log->debug("cluster sanity: scoped point mismatch at minind={} majind={} spt={} kdpt={}, blob is: {}",
minind, majind, spt, kdpt, *sblob);
if (log) log->debug("cluster sanity: scoped point mismatch at minind={} majind={} spt={} kdpt={}, blob is: {}",
minind, majind, spt, kdpt, *sblob);
return false;
}

}
return true;
}
Expand Down Expand Up @@ -968,13 +976,13 @@ void Facade::sort_blobs(std::vector<Blob*>& blobs)
bool Facade::cluster_less(const Cluster* a, const Cluster* b)
{
if (a==b) return false;
/// Keep this off for now to match the hacked March 24'th version
// {
// const double la = a->get_length();
// const double lb = b->get_length();
// if (la < lb) return true;
// if (lb < la) return false;
// }

{
const double la = a->get_length();
const double lb = b->get_length();
if (la < lb) return true;
if (lb < la) return false;
}
{
const int na = a->nchildren();
const int nb = b->nchildren();
Expand Down
16 changes: 13 additions & 3 deletions img/test/doctest_clustering_prototype.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -108,7 +108,7 @@ Points::node_ptr make_simple_pctree()
return root;
}

TEST_CASE("test PointTree API")
TEST_CASE("clustering point tree")
{
// this test does not touch the facades so needs to Grouping root
auto root = make_simple_pctree(); // a cluster node.
Expand Down Expand Up @@ -186,11 +186,10 @@ TEST_CASE("test PointTree API")
points[0][index], points[1][index], points[2][index], metric);
}
CHECK(rad.size() == 2);

}


TEST_CASE("test PointCloudFacade")
TEST_CASE("clustering facade")
{
Points::node_t root_node;
Grouping* grouping = root_node.value.facade<Grouping>();
Expand All @@ -201,7 +200,18 @@ TEST_CASE("test PointCloudFacade")
REQUIRE(pccptr->grouping() == grouping);
Cluster& pcc = *pccptr;

CHECK(pcc.sanity());

auto& blobs = pcc.children();

// (0.5 * 1 + 1.5 * 2) / 3 = 1.1666666666666665
debug("blob 0: q={}, r={}", blobs[0]->charge(), blobs[0]->center_x());
debug("blob 1: q={}, r={}", blobs[1]->charge(), blobs[1]->center_x());
double expect = 0;
expect += blobs[0]->charge() * blobs[0]->center_x();
expect += blobs[1]->charge() * blobs[1]->center_x();
expect /= blobs[0]->charge() + blobs[1]->charge();
debug("expect average pos {}", expect);
auto ave_pos = pcc.calc_ave_pos({1,0,0}, 1);
debug("ave_pos: {} | expecting (1.1666666666666665 0 0)", ave_pos);
auto l1 = fabs(ave_pos[0] - 1.1666666666666665) + fabs(ave_pos[1]) + fabs(ave_pos[2]);
Expand Down
1 change: 0 additions & 1 deletion util/inc/WireCellUtil/NaryTreeFacade.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -163,7 +163,6 @@ namespace WireCell::NaryTree {
}
m_children.push_back(child);
}
std::sort(m_children.begin(), m_children.end());
}
return m_children;
}
Expand Down
14 changes: 0 additions & 14 deletions util/inc/WireCellUtil/PointTree.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -171,20 +171,6 @@ namespace WireCell::PointCloud::Tree {
{
return const_cast<const ScopedView<ElementType>&>(
const_cast<self_t*>(this)->scoped_view(scope));

// using SV = ScopedView<ElementType>;
// auto const* sbptr = get_scoped(scope);
// if (sbptr) {
// auto const* svptr = dynamic_cast<const SV*>(sbptr);
// if (svptr) {
// return *svptr;
// }
// }
// auto uptr = std::make_unique<SV>(scope);
// auto& sv = *uptr;
// m_scoped[scope] = std::move(uptr);
// init(scope);
// return sv;
}
template<typename ElementType>
ScopedView<ElementType>& Points::scoped_view(const Scope& scope)
Expand Down
3 changes: 3 additions & 0 deletions util/test/doctest_nary_tree.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -345,6 +345,9 @@ TEST_CASE("nary tree depth limits")
}
CHECK(count == want[ind]);
}
for (const auto* child : root->children()[0]->children()[0]->children()) {
debug("gggc: {}", child->value.name);
}

}

Expand Down

0 comments on commit 62240a3

Please sign in to comment.