Merge pull request #294 from WireCell/apply-pointcloud

Apply pointcloud

.gitignore

@@ -78,3 +78,5 @@ configure
 
 *.gz
 styles
+*.heap
+*.prof

apps/apps/wcwires.cxx

@@ -12,23 +12,23 @@
 using namespace WireCell;
 using namespace WireCell::WireSchema;
 
-static void
-parse_param(std::string name,
-            const std::vector<std::string>& args,
-            std::map<std::string,std::string>& store)
-{
-    for (auto one : args) {
-        auto two = String::split(one, "=");
-        if (two.size() != 2) {
-            std::cerr
-                << name
-                << ": parameters are set as <name>=<value>, got "
-                << one << std::endl;
-            throw CLI::CallForHelp();
-        }
-        store[two[0]] = two[1];
-    }
-}
+// static void
+// parse_param(std::string name,
+//             const std::vector<std::string>& args,
+//             std::map<std::string,std::string>& store)
+// {
+//     for (auto one : args) {
+//         auto two = String::split(one, "=");
+//         if (two.size() != 2) {
+//             std::cerr
+//                 << name
+//                 << ": parameters are set as <name>=<value>, got "
+//                 << one << std::endl;
+//             throw CLI::CallForHelp();
+//         }
+//         store[two[0]] = two[1];
+//     }
+// }
 
 int main(int argc, char** argv)
 {

aux/docs/tensor-data-model.org

@@ -127,37 +127,100 @@ The *datatype* of *pcnamedset* indicates a tensor representing a ~std::map<std::
 
 - items :: an object representing the named point cloud set.  Each attribute name provides the name of the point cloud and the value provides the *datapath* to a *pcdataset*.
 
-** pctreenode
+** pctree
 
-The *datatype* of *pctreenode* indicates a tensor representing an $n$-ary tree node such as represented in C++ with ~WireCell::NaryTree::Node~ with a value type of ~WireCell::PointCloud::Tree::Points~.  The tensor has no array part and shall provide these metadata attributes:
-
-- parent :: a string which either holds a *datapath* to the parent *pctreenode* or is empty indicating the node is the root of the tree.
-
-- children :: a list of string with each element providing the *datapath* of a *pctreenode* constructed from a child node.  The order of this list reflects the sibling order.  The list may be empty if the node is a leaf.
-
-- pointclouds :: the "local" point cloud datasets held as a *pcnamedset*
+#+begin_note
+All node-local point clouds of a given name must be comprised of arrays of the same names and types.  An earlier data model called *pctreenode* did not carry this restriction but results in many small records which led to ruinously slow conversions.
+#+end_note
 
-The tree structure shall be definitely represented by the *datapath* values held in the *parent* and *children* attributes.  The tree structure may additionally be reflected in a hierarchical *datapath* naming convention.  If such a convention is used it shall reflect the descending order from root to leaf.  For example a tree with three nodes, one root and two children and the children have named point cloud sets each with a single point cloud named "3d" and each of those have arrays named "x", "y", "z", and "q":
+The *datatype* of *pctree* indicates a tensor representing an \(n\)-ary tree such as represented in C++ with ~WireCell::NaryTree::Node~ with a value type of ~WireCell::PointCloud::Tree::Points~.
+
+The *pctree* has an array part that represents the tree structure as a flattened *parentage map*.  
+The metadata of *pctree* has these attributes:
+
+- pointclouds :: the point cloud datasets of the tree held as a *pcnamedset*
+- lpcmaps :: the local point cloud mapping into the concatenated *pointclouds* as a *dataset*.
+
+#+begin_src dot :width 200pt :file tdm-pctree.svg :exports none :results none
+  graph {
+      node[shape=circle]
+      0 1 2 3 4 5
+      node[shape=plaintext,label=<
+           <TABLE BORDER="0" CELLPADDING="3" CELLBORDER="1" CELLSPACING="3">
+           <TR>
+           <TD>a </TD><TD>w </TD>
+           </TR>
+           <TR>
+           <TD rowspan="2">b </TD><TD>x </TD>
+           </TR>
+           <TR>
+           <TD>y </TD>
+           </TR>
+           <TR>
+           <TD>c </TD><TD>z </TD>
+           </TR>
+           </TABLE>>];
+      lpc0 lpc1 lpc2 lpc3 lpc4 lpc5
+
+      {rank=same 0 lpc0}
+      {rank=same 1 lpc1}
+      {rank=same 2 lpc2}
+      {rank=same 3 lpc3}
+      {rank=same 4 lpc4}
+      {rank=same 5 lpc5}
+
+      0 -- 1
+      0 -- 2
+      0 -- 3
+      3 -- 4
+      3 -- 5
+
+      0 -- lpc0
+      1 -- lpc1
+      2 -- lpc2
+      3 -- lpc3
+      4 -- lpc4
+      5 -- lpc5
+  }
+#+end_src
+
+
+To describe the *parentage map* array and the *pointclouds* and *lpcmaps* elements of a *pctree*
+consider the point cloud tree shown in figure [[fig:pctree-example]].  It consists of six nodes which are labeled in the order they are visited in a depth-first descent.  Each node has two local point clouds of zero or more points.  When represented as a ~NaryTree::Node~ with ~Tree::Points~ type, any PC with zero points are omitted but these empty point clouds may be conceptually included.  
+
+#+name: fig:pctree-example
+#+caption: Example point cloud tree with nodes labeled in order of depth-first descent.  Each node has a three local points clouds ("a", "b" and "c")
+#+attr_latex: :width 0.5\textwidth
+[[file:tdm-pctree.svg]]
+
+This tree will have a *parentage map* array as shown in table [[tab:parentage-example]].  Each element of the array represents a node as visited in depth-first descent order.  The value of the element is the index of the parent of the node.  To indicate the non-existent parent of a root node, the element value is set to the element index.
+This representation allows for a tree to have multiple roots (disjoint subtrees).
+
+#+name: tab:parentage-example
+#+caption: Example parentage map for the six-node example pctree.   
+#+ATTR_LATEX: :align |c|c|c|c|c|c|l|
+|---+---+---+---+---+---+--------------|
+| 0 | 1 | 2 | 3 | 4 | 5 | node index   |
+|---+---+---+---+---+---+--------------|
+| 0 | 0 | 0 | 0 | 3 | 3 | parent index |
+|---+---+---+---+---+---+--------------|
+
+If we assume that only node 0 has non-empty PC "a", only nodes 1, 2 and 3 have non-empty PC "b" and only nodes 4 and 5 have non-empty PC "c" the *lpcmaps* *dataset* will be as illustrated in table [[tab:lpcmaps-example]].  For simplicity, this example associates (non-empty) PCs with a "layer" in the tree.  However, the representation is not limited and can represent point clouds of a given name that span layers.
+
+#+name: tab:lpcmaps-example
+#+caption: Example local point clouds map for the six-node example pctree assuming a distribution of sizes of local point clouds across the nodes.
+#+ATTR_LATEX: :align |c|c|c|c|c|c|l|
+|-------+-------+-------+-------+-------+-------+--------------|
+|     0 |     1 |     2 |     3 |     4 |     5 | (node index) |
+|-------+-------+-------+-------+-------+-------+--------------|
+| $n_0$ |     0 |     0 |     0 |     0 |     0 | a (sizes)    |
+|     0 | $n_1$ | $n_2$ | $n_3$ |     0 |     0 | b (sizes)    |
+|     0 |     0 |     0 |     0 | $n_4$ | $n_5$ | c (sizes)    |
+|-------+-------+-------+-------+-------+-------+--------------|
+
+
+This example will have an entry "a" in *pointclouds* consisting of array "w" that has total size $n_0$, entry "b" with arrays "x" and "y" of total size $n_1 + n_2 + n_3$ and entry "c" with array "z" with total size "$n_4 + n_5$.  From these concatenated *pointclouds* and the *lpcmaps* and the depth-first descent ordering of the *parentage map* the individual local point clouds may be reconstructed.
 
-#+begin_example
-datatype             datapath
-----------           ---------
-pctreenode           root
-pctreenode           root/0
-pcnamedset           root/0/pointclouds
-pcdataset            root/0/pointclouds/namedpcs/3d
-pcarray              root/0/pointclouds/namedpcs/3d/arrays/q
-pcarray              root/0/pointclouds/namedpcs/3d/arrays/x
-pcarray              root/0/pointclouds/namedpcs/3d/arrays/y
-pcarray              root/0/pointclouds/namedpcs/3d/arrays/z
-pctreenode           root/1
-pcnamedset           root/1/pointclouds
-pcdataset            root/1/pointclouds/namedpcs/3d
-pcarray              root/1/pointclouds/namedpcs/3d/arrays/q
-pcarray              root/1/pointclouds/namedpcs/3d/arrays/x
-pcarray              root/1/pointclouds/namedpcs/3d/arrays/y
-pcarray              root/1/pointclouds/namedpcs/3d/arrays/z
-#+end_example
 
 ** trace
 

aux/inc/WireCellAux/ClusterArrays.h

@@ -4,7 +4,7 @@
 #include "ClusterHelpers.h"
 #include "WireCellIface/ICluster.h"
 #include "WireCellIface/ITensorSet.h"
-#include <boost/multi_array.hpp>
+#include "WireCellUtil/MultiArray.h"
 #include <boost/range/counting_range.hpp>
 
 #include <vector>

aux/inc/WireCellAux/ClusterHelpers.h

@@ -6,12 +6,13 @@
 #define WIRECELLAUX_CLUSTERJSONIFY
 
 #include "WireCellUtil/Configuration.h"
+#include "WireCellUtil/GraphTools.h"
 #include "WireCellIface/IFrame.h"
 #include "WireCellIface/ISlice.h"
 #include "WireCellIface/ICluster.h"
 #include "WireCellIface/IAnodePlane.h"
 
-#include <boost/graph/graphviz.hpp>
+#include "WireCellUtil/Graph.h"
 #include <sstream>
 #include <functional>
 #include <unordered_set>
@@ -142,6 +143,34 @@ namespace WireCell::Aux {
 
     /// Return counts indexed by node and/or edge code
     std::map<std::string, size_t> count(const cluster_graph_t& cgraph, bool nodes=true, bool edges=true);
+
+    /// function to merge  boost graph
+    template <typename GraphType>
+    GraphType merge_graphs(const std::vector<std::reference_wrapper<const GraphType>>& graphs)
+    {
+        GraphType merged_graph;
+        // merge the graphs
+        for (const auto& graphr : graphs) {
+            const auto& graph = graphr.get();
+            std::unordered_map<size_t, size_t> vertex_map;
+            // add the vertices
+            for (const auto& vtx : GraphTools::mir(boost::vertices(graph))) {
+                const auto& node = graph[vtx];
+                auto new_vtx = boost::add_vertex(node, merged_graph);
+                vertex_map[vtx] = new_vtx;
+            }
+            // add the edges
+            for (const auto& edg : GraphTools::mir(boost::edges(graph))) {
+                const auto& src = boost::source(edg, graph);
+                const auto& tgt = boost::target(edg, graph);
+                const auto& edge = graph[edg];
+                auto new_src = vertex_map[src];
+                auto new_tgt = vertex_map[tgt];
+                boost::add_edge(new_src, new_tgt, edge, merged_graph);
+            }
+        }
+        return merged_graph;
+    }
 }
 
 #endif

aux/inc/WireCellAux/DftTools.h

@@ -48,7 +48,7 @@
 
 #include "WireCellIface/IDFT.h"
 #include <vector>
-#include <Eigen/Core>
+#include "WireCellUtil/Eigen.h"
 
 namespace WireCell::Aux::DftTools {
 

aux/inc/WireCellAux/SimpleCluster.h

@@ -1,5 +1,5 @@
 #include "WireCellIface/ICluster.h"
-#include <boost/graph/copy.hpp>
+#include "WireCellUtil/Graph.h"
 namespace WireCell::Aux {
 
     class SimpleCluster : public ICluster {

aux/inc/WireCellAux/SimpleTensor.h

@@ -2,8 +2,7 @@
 #define WIRECELL_AUX_SIMPLETENSOR
 
 #include "WireCellIface/ITensor.h"
-
-#include <boost/multi_array.hpp>
+#include "WireCellUtil/MultiArray.h"
 #include <cstring>
 
 namespace WireCell::Aux {
@@ -12,13 +11,42 @@ namespace WireCell::Aux {
     class SimpleTensor : public WireCell::ITensor {
       public:
 
+        // Create a null tensor.
+        SimpleTensor()
+            : m_typeinfo(typeid(void))
+        {
+        }
         // Create a simple tensor with a null array and type, just metadata, if any
-        explicit SimpleTensor(const Configuration& md = Configuration())
-            : m_md(md)
-            , m_typeinfo(typeid(void))
+        explicit SimpleTensor(const Configuration& md)
+            : m_typeinfo(typeid(void))
+            , m_mdptr(std::make_unique<Configuration>(md))
         {
         }
 
+        // Create simple tensor, allocating space for data.  If data
+        // given it must have at least as many elements as implied by
+        // shape and that span will be copied into allocated memory.
+        // The SimpleTensor will allocate memory if a null data
+        // pointer is given but note the type is required so pass call like:
+        //   SimpleTensor(shape, (Type*)nullptr)
+        template <typename ElementType>
+        SimpleTensor(const shape_t& shape, const ElementType* data)
+            : m_typeinfo(typeid(ElementType))
+            , m_sizeof(sizeof(ElementType))
+            , m_shape(shape)
+        {
+            size_t nbytes = element_size();
+            for (const auto& s : m_shape) {
+                nbytes *= s;
+            }
+            if (data) {
+                const std::byte* bytes = reinterpret_cast<const std::byte*>(data);
+                m_store.assign(bytes, bytes+nbytes);
+            }
+            else {
+                m_store.resize(nbytes);
+            }
+        }
         // Create simple tensor, allocating space for data.  If data
         // given it must have at least as many elements as implied by
         // shape and that span will be copied into allocated memory.
@@ -28,11 +56,11 @@ namespace WireCell::Aux {
         template <typename ElementType>
         SimpleTensor(const shape_t& shape,
                      const ElementType* data,
-                     const Configuration& md = Configuration())
-            : m_shape(shape)
-            , m_md(md)
-            , m_typeinfo(typeid(ElementType))
+                     const Configuration& md)
+            : m_typeinfo(typeid(ElementType))
             , m_sizeof(sizeof(ElementType))
+            , m_shape(shape)
+            , m_mdptr(std::make_unique<Configuration>(md))
         {
             size_t nbytes = element_size();
             for (const auto& s : m_shape) {
@@ -59,7 +87,14 @@ namespace WireCell::Aux {
          md[1][2][3] = 42.0;
         */
         std::vector<std::byte>& store() { return m_store; }
-        Configuration& metadata() { return m_md; }
+        Configuration& metadata() {
+            if (!m_mdptr) {
+                // lazy construction.
+                m_mdptr = std::make_unique<Configuration>();
+            }
+            return *m_mdptr;
+        }
+
 
         // ITensor const interface.
         virtual const std::type_info& element_type() const { return m_typeinfo; }
@@ -70,14 +105,14 @@ namespace WireCell::Aux {
         virtual const std::byte* data() const { return m_store.data(); }
         virtual size_t size() const { return m_store.size(); }
 
-        virtual Configuration metadata() const { return m_md; }
+        virtual Configuration metadata() const { return *m_mdptr; }
 
       private:
-        std::vector<std::byte> m_store;
-        std::vector<size_t> m_shape;
-        Configuration m_md;
         const std::type_info& m_typeinfo;
         size_t m_sizeof{0};
+        std::vector<size_t> m_shape;
+        std::vector<std::byte> m_store;
+        std::unique_ptr<Configuration> m_mdptr; // avoid constructor if empty
     };
 
 }  // namespace WireCell::Aux

aux/inc/WireCellAux/TensorDMcluster.h

@@ -7,6 +7,7 @@
 #include "WireCellIface/ICluster.h"
 #include "WireCellIface/ITensor.h"
 #include "WireCellIface/IAnodePlane.h"
+#include "WireCellAux/TensorDMcommon.h"
 
 namespace WireCell::Aux::TensorDM {
 
@@ -39,6 +40,9 @@ namespace WireCell::Aux::TensorDM {
     ICluster::pointer as_cluster(const ITensor::vector& tens,
                                  const IAnodePlane::vector& anodes,
                                  const std::string& datapath="");
+    ICluster::pointer as_cluster(const TensorIndex& ti,
+                                 const IAnodePlane::vector& anodes,
+                                 const std::string& datapath="");
 }
 
 #endif