include/pops/amr/tagging/cluster.hpp Source FileΒΆ

adc_cpp: include/pops/amr/tagging/cluster.hpp Source File
adc_cpp 0.3.0
Model-free C++23 core for coupled hyperbolic-elliptic systems on adaptive (AMR) meshes, with MPI and GPU (Kokkos) backends
cluster.hpp
Go to the documentation of this file.
1
16
17#pragma once
18
19#include <pops/amr/tagging/tag_box.hpp>
20#include <pops/mesh/index/box2d.hpp>
21#include <pops/mesh/layout/box_array.hpp>
22
23#include <algorithm>
24#include <climits>
25#include <cstdlib>
26#include <vector>
27
28namespace pops {
29
31struct ClusterParams {
32 double min_efficiency = 0.7;
33 int min_box_size = 1;
34 int max_box_size = 32;
35};
36
37namespace detail {
38
40inline Box2D tag_bbox(const TagBox& tb, const Box2D& region) {
41 int lo0 = INT_MAX, lo1 = INT_MAX, hi0 = INT_MIN, hi1 = INT_MIN;
42 for (int j = region.lo[1]; j <= region.hi[1]; ++j)
43 for (int i = region.lo[0]; i <= region.hi[0]; ++i)
44 if (tb(i, j)) {
45 lo0 = std::min(lo0, i);
46 hi0 = std::max(hi0, i);
47 lo1 = std::min(lo1, j);
48 hi1 = std::max(hi1, j);
49 }
50 if (hi0 < lo0)
51 return Box2D{{0, 0}, {-1, -1}}; // empty
52 return Box2D{{lo0, lo1}, {hi0, hi1}};
53}
54
56inline long count_in(const TagBox& tb, const Box2D& r) {
57 long c = 0;
58 for (int j = r.lo[1]; j <= r.hi[1]; ++j)
59 for (int i = r.lo[0]; i <= r.hi[0]; ++i)
60 c += tb(i, j);
61 return c;
62}
63
65inline std::vector<long> signature(const TagBox& tb, const Box2D& r, int axis) {
66 const int len = (axis == 0) ? r.nx() : r.ny();
67 std::vector<long> s(len, 0);
68 for (int j = r.lo[1]; j <= r.hi[1]; ++j)
69 for (int i = r.lo[0]; i <= r.hi[0]; ++i)
70 if (tb(i, j))
71 s[(axis == 0) ? (i - r.lo[0]) : (j - r.lo[1])] += 1;
72 return s;
73}
74
76inline int best_hole(const std::vector<long>& s, int mb) {
77 const int len = static_cast<int>(s.size());
78 int best = -1, bestd = INT_MAX, c = len / 2;
79 for (int k = mb; k <= len - mb; ++k)
80 if (s[k] == 0) {
81 int d = std::abs(k - c);
82 if (d < bestd) {
83 bestd = d;
84 best = k;
85 }
86 }
87 return best;
88}
89
92inline int best_inflection(const std::vector<long>& s, int mb, long& score) {
93 const int len = static_cast<int>(s.size());
94 score = 0;
95 if (len < 3)
96 return -1;
97 std::vector<long> D(len, 0);
98 for (int k = 1; k < len - 1; ++k)
99 D[k] = s[k + 1] - 2 * s[k] + s[k - 1];
100 int best = -1;
101 const int klo = std::max(mb, 2), khi = std::min(len - mb, len - 2);
102 for (int k = klo; k <= khi; ++k) {
103 long delta = std::labs(D[k] - D[k - 1]);
104 if (delta > score) {
105 score = delta;
106 best = k;
107 }
108 }
109 return best;
110}
111
115inline void cluster_rec(const TagBox& tb, Box2D region, const ClusterParams& p,
116 std::vector<Box2D>& out) {
117 region = tag_bbox(tb, region);
118 if (region.empty())
119 return;
120
121 const long ntag = count_in(tb, region);
122 const double eff = double(ntag) / double(region.num_cells());
123 const int mb = std::max(1, p.min_box_size);
124 const bool sx = region.nx() >= 2 * mb;
125 const bool sy = region.ny() >= 2 * mb;
126 if (eff >= p.min_efficiency || (!sx && !sy)) {
127 out.push_back(region);
128 return;
129 }
130
131 const auto Sx = signature(tb, region, 0);
132 const auto Sy = signature(tb, region, 1);
133 const int hx = sx ? best_hole(Sx, mb) : -1;
134 const int hy = sy ? best_hole(Sy, mb) : -1;
135
136 int axis = -1, kcut = -1;
137 if (hx >= 0 && hy >= 0) {
138 axis = (region.nx() >= region.ny()) ? 0 : 1;
139 kcut = (axis == 0) ? hx : hy;
140 } else if (hx >= 0) {
141 axis = 0;
142 kcut = hx;
143 } else if (hy >= 0) {
144 axis = 1;
145 kcut = hy;
146 } else {
147 long scx = 0, scy = 0;
148 const int ix = sx ? best_inflection(Sx, mb, scx) : -1;
149 const int iy = sy ? best_inflection(Sy, mb, scy) : -1;
150 if (ix >= 0 && (scx >= scy || iy < 0)) {
151 axis = 0;
152 kcut = ix;
153 } else if (iy >= 0) {
154 axis = 1;
155 kcut = iy;
156 } else { // midpoint of the largest splittable dimension
157 if (sx && (region.nx() >= region.ny() || !sy)) {
158 axis = 0;
159 kcut = region.nx() / 2;
160 } else {
161 axis = 1;
162 kcut = region.ny() / 2;
163 }
164 }
165 }
166
167 Box2D left = region, right = region;
168 if (axis == 0) {
169 left.hi[0] = region.lo[0] + kcut - 1;
170 right.lo[0] = region.lo[0] + kcut;
171 } else {
172 left.hi[1] = region.lo[1] + kcut - 1;
173 right.lo[1] = region.lo[1] + kcut;
174 }
175 cluster_rec(tb, left, p, out);
176 cluster_rec(tb, right, p, out);
177}
178
179} // namespace detail
180
185inline std::vector<Box2D> berger_rigoutsos(const TagBox& tags, const ClusterParams& p = {}) {
186 std::vector<Box2D> raw;
187 detail::cluster_rec(tags, tags.box, p, raw);
188 std::vector<Box2D> result;
189 for (const auto& b : raw) {
190 BoxArray chopped = BoxArray::from_domain(b, p.max_box_size);
191 for (int k = 0; k < chopped.size(); ++k)
192 result.push_back(chopped[k]);
193 }
194 return result;
195}
196
197} // namespace pops
box2d.hpp
Box2D: the integer index space of a 2D cell-centered Cartesian grid.
box_array.hpp
BoxArray: the set of boxes tiling a level (disjoint, covering).
pops::BoxArray::from_domain
static BoxArray from_domain(const Box2D &domain, int max_grid_size)
Tile the domain into tiles of at most max_grid_size per direction, distributed evenly.
Definition box_array.hpp:30
pops::detail::cluster_rec
void cluster_rec(const TagBox &tb, Box2D region, const ClusterParams &p, std::vector< Box2D > &out)
Berger-Rigoutsos recursive core: trim, accept if efficient/not splittable, otherwise cut and recurse.
Definition cluster.hpp:115
pops::detail::tag_bbox
Box2D tag_bbox(const TagBox &tb, const Box2D &region)
Bounding box of the tagged cells in region; empty box (hi < lo) if none is tagged.
Definition cluster.hpp:40
pops::detail::best_inflection
int best_inflection(const std::vector< long > &s, int mb, long &score)
Inflection: index of the max |D[k] - D[k-1]| with D the discrete Laplacian of the signature,...
Definition cluster.hpp:92
pops::detail::best_hole
int best_hole(const std::vector< long > &s, int mb)
Interior hole (zero signature) closest to the center, in [mb, len-mb]; -1 if none.
Definition cluster.hpp:76
pops::detail::count_in
long count_in(const TagBox &tb, const Box2D &r)
Number of tagged cells in box r.
Definition cluster.hpp:56
pops::detail::signature
std::vector< long > signature(const TagBox &tb, const Box2D &r, int axis)
Signature of r along axis: number of tagged cells per column (axis 0) or per row (axis 1).
Definition cluster.hpp:65
pops
Definition amr_hierarchy.hpp:29
pops::berger_rigoutsos
std::vector< Box2D > berger_rigoutsos(const TagBox &tags, const ClusterParams &p={})
Cluster a TagBox into boxes covering the tagged cells (Berger-Rigoutsos), then final chop.
Definition cluster.hpp:185
pops::Box2D
2D integer index space, cell-centered.
Definition box2d.hpp:37
pops::Box2D::empty
bool empty() const
true if the box contains no cell (hi < lo in one direction).
Definition box2d.hpp:58
pops::Box2D::nx
POPS_HD int nx() const
Width (direction 0). POPS_HD (called from Geometry::dx() in a device kernel).
Definition box2d.hpp:50
pops::Box2D::hi
int hi[2]
Definition box2d.hpp:39
pops::Box2D::ny
POPS_HD int ny() const
Height (direction 1). POPS_HD (called from Geometry::dy() in a device kernel).
Definition box2d.hpp:52
pops::Box2D::lo
int lo[2]
Definition box2d.hpp:38
pops::Box2D::num_cells
std::int64_t num_cells() const
Total number of cells (nx*ny, floored at 0 per direction): 0 if the box is empty.
Definition box2d.hpp:54
pops::ClusterParams
Berger-Rigoutsos clustering parameters (configuration object).
Definition cluster.hpp:31
pops::ClusterParams::min_efficiency
double min_efficiency
efficiency threshold (tagged fraction) to accept a box.
Definition cluster.hpp:32
pops::ClusterParams::max_box_size
int max_box_size
max size of a box; accepted boxes are chopped to this size.
Definition cluster.hpp:34
pops::ClusterParams::min_box_size
int min_box_size
minimal size of a box; bounds the admissible cuts.
Definition cluster.hpp:33
pops::TagBox
Dense grid of 0/1 markers over a box, input to Berger-Rigoutsos clustering.
Definition tag_box.hpp:35
pops::TagBox::box
Box2D box
Definition tag_box.hpp:36
tag_box.hpp
TagBox: dense grid of markers (0/1) over a region, input to Berger-Rigoutsos clustering.
Generated by doxygen 1.9.8