Fix #211.

The "measure" is found by applying the connected components algorithm
to a special per-slice, per plane graph holding nodes representing
blobs and channels.  The nodes of this graph hold the vertex
descriptor of the original cluster graph node.

The bug was due to a channel being added to this special graph each
time it was seen across the blobs.  Thus, by construction, the
connected components (ie, measures) always had exactly one blob.

The fix was to only add any given channel exactly once.

The following test is extended to check this bug:

```
$ bats img/test/test-wct-uboone-img.bats
test-wct-uboone-img.bats
 ✓ create graph
 ✓ check log files
 ✓ inspect blobs
 ✓ dump blobs
 ✓ at least one multi-blob measure

5 tests, 0 failures
```

img/inc/WireCellImg/BlobGrouping.h

@@ -1,15 +1,14 @@
 /** BlobGrouping takes in a channel-level cluster and produces another
- * with channels merged measure ('m' nodes)
+ * with channels merged into measure ('m' nodes)
  *
- * The input cluster must have (b,w), (c,w) and (s,b) edges and may have (b,b) edges.
+ * The input cluster must have (b-w), (c-w) and (b-s) edges and may have (b-b) edges.
  *
- * The output cluster will have (m, b) and (s,b) edges and if the
- * input has (b,b) edges, they are preserved.
+ * The output cluster will add (b-m) and (c-m) edges.
  *
  * The created m-node IMeasures will have sequential ident() unique to
  * the context of the cluster.
  *
- * Grouping is done in the "coarse grained" strategy.
+ * Grouping is done in the "coarse grained" strategy.  (see raygrid.pdf).
  *
  * See manual for more info.
  */

img/src/BlobGrouping.cxx

@@ -22,7 +22,9 @@ namespace bcdesc {
     struct nprop_t { cluster_vertex_t desc; };
     struct eprop_t { };
 
-    using graph_t = boost::adjacency_list<boost::setS, boost::vecS, boost::undirectedS,
+    using graph_t = boost::adjacency_list<boost::setS, 
+                                          boost::vecS, 
+                                          boost::undirectedS,
                                           nprop_t, eprop_t>;
 
     using vdesc_t = boost::graph_traits<graph_t>::vertex_descriptor;
@@ -48,18 +50,21 @@ void doit(cluster_graph_t& cgraph)
 
         // Recieve b and c reached from this s per plane
         std::vector<bcdesc::graph_t> bcs(3); // fixme: hard-code 3 planes
+        std::vector< std::unordered_map<int, bcdesc::vdesc_t> > uniq_chans(3);
 
         for (auto bvtx : mir(boost::adjacent_vertices(svtx, cgraph))) {
             const auto& bnode = cgraph[bvtx];
             const char bcode = bnode.code();
             if (bcode != 'b') {
                 continue;
             }
-
+
+            const auto iblob = std::get<IBlob::pointer>(bnode.ptr);
+
             // add this blob to each plane graph
             std::vector<bcdesc::vdesc_t> bidx;
             for (auto& bc : bcs) {
-                bidx.push_back(boost::add_vertex({bvtx}, bc));
+                bidx.push_back(boost::add_vertex({ bvtx }, bc));
             }
 
             // wires
@@ -70,6 +75,8 @@ void doit(cluster_graph_t& cgraph)
                     continue;
                 }
 
+                const auto iwire = std::get<IWire::pointer>(wnode.ptr);
+
                 // channels
                 for (auto cvtx : mir(boost::adjacent_vertices(wvtx, cgraph))) {
                     const auto& cnode = cgraph[cvtx];
@@ -81,41 +88,54 @@ void doit(cluster_graph_t& cgraph)
                     auto ich = std::get<IChannel::pointer>(cnode.ptr);
                     const auto wpid = ich->planeid();
                     const auto pind = wpid.index();
+                    const auto cident = ich->ident();
+
                     auto& bc = bcs[pind];
 
-                    const auto cidx = boost::add_vertex({ cvtx }, bc);
-                    boost::add_edge(bidx[pind], cidx, bc);
+                    bcdesc::vdesc_t cidx;
+                    auto& uniq_chan = uniq_chans[pind];
+                    auto ucit = uniq_chan.find(cident);
+                    if (ucit == uniq_chan.end()) {
+                        cidx = boost::add_vertex({ cvtx }, bc);
+                        uniq_chan[cident] = cidx;
+                    }
+                    else {
+                        cidx = ucit->second;
+                    }
 
+                    boost::add_edge(bidx[pind], cidx, bc);
                 } // channels
             } // wires
         } // blobs
 
         // now have a per-slice, per-plane graph of blob-channel vertices
-        for (const auto& bc : bcs) {
-            // every cc's int locally names a measure
-            std::unordered_map<bcdesc::vdesc_t, int> cc;
-            boost::connected_components(bc, boost::make_assoc_property_map(cc));
-
-            std::unordered_map<int, std::pair<Aux::SimpleMeasure*,cluster_vertex_t>> measures;
+        for (size_t iplane = 0; iplane < 3; ++iplane) {
+            const auto& bc = bcs[iplane];
 
-            for (const auto& [idx, num] : cc) {
-                Aux::SimpleMeasure* sm=nullptr;
-                cluster_vertex_t mvtx;
+            size_t nverts = boost::num_vertices(bc);
 
-                const auto it = measures.find(num);
-                if (it == measures.end()) {
-                    sm = new Aux::SimpleMeasure(tot_meas++);
-                    mvtx = boost::add_vertex(IMeasure::pointer(sm), cgraph);
-                    measures[num] = std::make_pair(sm, mvtx);
-                }
-                else {
-                    sm = it->second.first;
-                    mvtx = it->second.second;
+            // every cc's int locally names a measure
+            // std::unordered_map<bcdesc::vdesc_t, int> cc;
+            // auto ncomponents = boost::connected_components(bc, boost::make_assoc_property_map(cc));
+            std::vector<int> cc(nverts);
+            auto ncomponents = boost::connected_components(bc, &cc[0]);
+
+            std::vector<Aux::SimpleMeasure*> measures(ncomponents, nullptr);
+            std::vector<cluster_vertex_t> mvtxs(ncomponents);
+
+            for (size_t idx = 0; idx < nverts; ++idx) {
+            // for (const auto& [idx, ncomp] : cc) {
+                size_t ncomp = cc[idx];
+                Aux::SimpleMeasure* sm=measures[ncomp];
+                cluster_vertex_t mvtx=mvtxs[ncomp];
+                if (!sm) {
+                    sm = measures[ncomp] = new Aux::SimpleMeasure(tot_meas++);
+                    mvtx = mvtxs[ncomp] = boost::add_vertex(IMeasure::pointer(sm), cgraph);
                 }
 
                 // edge to measure for either blob or channel
                 const auto vtx = bc[idx].desc;
-                boost::add_edge(mvtx, vtx, cgraph);
+                boost::add_edge(vtx, mvtx, cgraph);
 
                 // if channel, accumulate signal
                 const auto& node = cgraph[vtx];
@@ -127,8 +147,7 @@ void doit(cluster_graph_t& cgraph)
                     sm->wpid = wpid;
                 }
             }
-        }
-
+        } // planes
     } // slices
 }
 

img/test/test-wct-uboone-img.bats

@@ -124,8 +124,8 @@ setup_file () {
         logs[$fmt]="$log"
 
         local dat="${base}.tar.gz"
-        echo "$wcimg inspect -o $log $dat"
-        run $wcimg inspect -o "$log" "$dat"
+        echo "$wcimg inspect --verbose -o $log $dat"
+        run $wcimg inspect --verbose -o "$log" "$dat"
         echo "$output"
         [[ "$status" -eq 0 ]]
         [[ -s "$log" ]]
@@ -152,8 +152,17 @@ setup_file () {
         [[ -s "${base}.dump" ]] 
     done
 
-    diff -u <(head -20 $(base_name json).dump) <(head -20 $(base_name numpy).dump)
-    local delta="$(diff -u <(head -20 $(base_name json).dump) <(head -20 $(base_name numpy).dump))"
+    run diff -u $(base_name json).dump $(base_name numpy).dump
+    echo "$output"
+    [[ "$status" -eq 0 ]]
     [[ -z "$delta" ]]
 }
 
+@test "at least one multi-blob measure" {
+    local fname="$(base_name numpy).inspect"
+    echo $fname
+
+    local got=$(grep 'nn for m' "${fname}" | grep -v 'b=1\b')
+    echo "$got"
+    [[ -n "$got" ]]
+}

img/test/test-wct-uboone-img.jsonnet

@@ -184,11 +184,13 @@ local img = {
             edges=[pg.edge(cs0,lcbr), pg.edge(lcbr,cs1)],
             name="chargesolving-" + aname),
         local solver = cs0,
-        ret: pg.intern(
+        local full = pg.intern(
             innodes=[bc], outnodes=[solver], centernodes=[bg],
             edges=[pg.edge(bc,bg), pg.edge(bg,solver)],
             name="solving-" + aname),
-        // ret: bc,
+        local nosolve = pg.pipeline([bc,bg]),
+        ret: full,
+
     }.ret,
 
     dump(outfile, fmt="json") :: {