diff --git a/aux/inc/WireCellAux/TensorDMcommon.h b/aux/inc/WireCellAux/TensorDMcommon.h index 53440c153..cd4b403e4 100644 --- a/aux/inc/WireCellAux/TensorDMcommon.h +++ b/aux/inc/WireCellAux/TensorDMcommon.h @@ -10,6 +10,11 @@ namespace WireCell::Aux::TensorDM { + using located_t = std::unordered_map; + + /// Index the tensors by their datapaths. + located_t index_datapaths(const ITensor::vector& tens); + /** Build metadata-only (array-less) tensor in the DM. */ ITensor::pointer make_metadata_tensor(const std::string& datatype, @@ -24,7 +29,8 @@ namespace WireCell::Aux::TensorDM { */ ITensor::pointer top_tensor(const ITensor::vector& tens, const std::string& datatype, - const std::string& datapath=""); + const std::string& datapath="", + const located_t& located = {}); /// Build a tensor set from set of tensors. The tensor data model @@ -35,10 +41,6 @@ namespace WireCell::Aux::TensorDM { const Configuration& tsetmd = Json::objectValue); - using located_t = std::map; - - /// Index the tensors by their datapaths. - located_t index_datapaths(const ITensor::vector& tens); /// Return first of type ITensor::pointer first_of(const ITensor::vector& tens, diff --git a/aux/inc/WireCellAux/TensorDMdataset.h b/aux/inc/WireCellAux/TensorDMdataset.h index c1b4eaf55..4d00821a2 100644 --- a/aux/inc/WireCellAux/TensorDMdataset.h +++ b/aux/inc/WireCellAux/TensorDMdataset.h @@ -5,6 +5,7 @@ #define WIRECELLAUX_TENSORDMDATASET #include "WireCellUtil/PointCloudDataset.h" +#include "WireCellAux/TensorDMcommon.h" #include "WireCellIface/ITensorSet.h" #include @@ -44,6 +45,7 @@ namespace WireCell::Aux::TensorDM { /// ValueError is thrown if tensors are badly formed. PointCloud::Dataset as_dataset(const ITensor::vector& tensors, const std::string& datapath="", + const located_t& located = {}, bool share=false); /// Convenience function calling above using the tensors from the @@ -51,6 +53,7 @@ namespace WireCell::Aux::TensorDM { /// ITensor with datatype "pcdataset" and not from ITensorSet. PointCloud::Dataset as_dataset(const ITensorSet::pointer& its, const std::string& datapath="", + const located_t& located = {}, bool share=false); diff --git a/aux/inc/WireCellAux/TensorDMpointtree.h b/aux/inc/WireCellAux/TensorDMpointtree.h index b7dd68efc..03d161828 100644 --- a/aux/inc/WireCellAux/TensorDMpointtree.h +++ b/aux/inc/WireCellAux/TensorDMpointtree.h @@ -64,7 +64,7 @@ namespace WireCell::Aux::TensorDM { empty in which case the first pcnamedset found is used. */ - named_pointclouds_t as_pcnamedset(const ITensor::vector& tens, const std::string& datapath = ""); + named_pointclouds_t as_pcnamedset(const ITensor::vector& tens, const std::string& datapath = "", const located_t& located = {}); diff --git a/aux/src/TensorDMcommon.cxx b/aux/src/TensorDMcommon.cxx index 63918f1be..e84f4107f 100644 --- a/aux/src/TensorDMcommon.cxx +++ b/aux/src/TensorDMcommon.cxx @@ -20,19 +20,21 @@ ITensor::pointer WireCell::Aux::TensorDM::make_metadata_tensor( ITensor::pointer WireCell::Aux::TensorDM::top_tensor(const ITensor::vector& tens, const std::string& datatype, - const std::string& datapath) + const std::string& datapath, + const located_t& located) { - for (const auto& iten : tens) { - const auto& tenmd = iten->metadata(); - const auto dtype = tenmd["datatype"].asString(); - const auto dpath = tenmd["datapath"].asString(); - if (dtype == datatype and (datapath.empty() or datapath == dpath)) { - return iten; - } + const auto it = located.find(datapath); + if (it == located.end()) { + raise("no array at datapath \"%s\"", datapath); } - raise("no array of datatype \"%s\" at datapath \"%s\"", - datatype, datapath); - return nullptr; + const auto& iten = it->second; + const auto& tenmd = iten->metadata(); + const auto dtype = tenmd["datatype"].asString(); + if (dtype != datatype) { + raise("array at datapath \"%s\" is of datatype \"%s\" not \"%s\"", + datapath, dtype, datatype); + } + return iten; } // /** Return a map from datapath to tensor. diff --git a/aux/src/TensorDMdataset.cxx b/aux/src/TensorDMdataset.cxx index 5d231bb54..a887732a2 100644 --- a/aux/src/TensorDMdataset.cxx +++ b/aux/src/TensorDMdataset.cxx @@ -106,27 +106,10 @@ PointCloud::Array WireCell::Aux::TensorDM::as_array(const ITensor::pointer& ten, PointCloud::Dataset WireCell::Aux::TensorDM::as_dataset(const ITensor::vector& tens, const std::string& datapath, + const located_t& located, bool share) { - // Index the tensors by datapath and find main "pcdataset" tensor. - ITensor::pointer top = nullptr; - std::unordered_map located; - for (const auto& iten : tens) { - const auto& tenmd = iten->metadata(); - const auto dtype = tenmd["datatype"].asString(); - const auto dpath = tenmd["datapath"].asString(); - if (!top and dtype == "pcdataset") { - if (datapath.empty() or datapath == dpath) { - top = iten; - } - continue; - } - if (dtype == "pcarray") { - located[dpath] = iten; - } - continue; - } - + ITensor::pointer top = top_tensor(tens, "pcdataset", datapath, located); if (!top) { THROW(ValueError() << errmsg{"no array of datatype \"pcdataset\""}); } @@ -137,7 +120,7 @@ WireCell::Aux::TensorDM::as_dataset(const ITensor::vector& tens, auto arrrefs = topmd["arrays"]; for (const auto& name : arrrefs.getMemberNames()) { const auto path = arrrefs[name].asString(); - auto it = located.find(path); + const auto it = located.find(path); if (it == located.end()) { THROW(ValueError() << errmsg{"no array \"" + name + "\" at path \"" + path + "\""}); } @@ -152,8 +135,9 @@ WireCell::Aux::TensorDM::as_dataset(const ITensor::vector& tens, PointCloud::Dataset WireCell::Aux::TensorDM::as_dataset(const ITensorSet::pointer& tens, const std::string& datapath, + const located_t& located, bool share) { auto sv = tens->tensors(); - return as_dataset(*(sv.get()), datapath, share); + return as_dataset(*(sv.get()), datapath, located, share); } diff --git a/aux/src/TensorDMframe.cxx b/aux/src/TensorDMframe.cxx index 07c78a1f3..d5ed17947 100644 --- a/aux/src/TensorDMframe.cxx +++ b/aux/src/TensorDMframe.cxx @@ -308,19 +308,9 @@ IFrame::pointer WireCell::Aux::TensorDM::as_frame(const ITensor::vector& tens, const std::string& datapath, std::function transform) { - ITensor::pointer ften=nullptr; - // map from datapath to tensor - std::unordered_map located; - for (const auto& iten : tens) { - auto dp = iten->metadata()["datapath"].asString(); - located[dp] = iten; - auto dt = iten->metadata()["datatype"].asString(); - if (!ften and dt == "frame") { - if (datapath.empty() or datapath == dp) { - ften = iten; - } - } - } + const auto& located = index_datapaths(tens); + ITensor::pointer ften=top_tensor(tens, "frame", datapath, located); + if (!ften) { THROW(ValueError() << errmsg{"no frame tensor found"}); } @@ -330,7 +320,7 @@ IFrame::pointer WireCell::Aux::TensorDM::as_frame(const ITensor::vector& tens, std::vector> straces; for (auto jtrpath : fmd["traces"]) { auto trpath = jtrpath.asString(); - auto trten = located[trpath]; + auto trten = located.at(trpath); auto trmd = trten->metadata(); int tbin = get(trmd, "tbin", 0); @@ -363,7 +353,7 @@ IFrame::pointer WireCell::Aux::TensorDM::as_frame(const ITensor::vector& tens, } auto tdpath = jtdpath.asString(); const bool share = false; // was true - PC::Dataset td = as_dataset(tens, tdpath, share); + PC::Dataset td = as_dataset(tens, tdpath, located, share); // dump_dataset(td, tdpath); auto tdmd = td.metadata(); diff --git a/aux/src/TensorDMpointgraph.cxx b/aux/src/TensorDMpointgraph.cxx index 36c84132e..470e5332c 100644 --- a/aux/src/TensorDMpointgraph.cxx +++ b/aux/src/TensorDMpointgraph.cxx @@ -13,31 +13,16 @@ WireCell::PointGraph WireCell::Aux::TensorDM::as_pointgraph(const ITensor::vector& tens, const std::string& datapath) { - std::unordered_map located; - ITensor::pointer top = nullptr; - for (const auto& iten : tens) { - const auto& tenmd = iten->metadata(); - const auto dtype = tenmd["datatype"].asString(); - const auto dpath = tenmd["datapath"].asString(); - if (!top and dtype == "pcgraph") { - if (datapath.empty() or datapath == dpath) { - top = iten; - } - continue; - } - if (dtype == "pcdataset") { - located[dpath] = iten; - } - continue; - } + const auto& located = index_datapaths(tens); + ITensor::pointer top = top_tensor(tens, "pcgraph", datapath, located); if (!top) { THROW(ValueError() << errmsg{"no array of datatype \"pcgraph\" at datapath \"" + datapath + "\"" }); } const auto& topmd = top->metadata(); - auto nodes = as_dataset(tens, topmd["nodes"].asString()); - auto edges = as_dataset(tens, topmd["edges"].asString()); + auto nodes = as_dataset(tens, topmd["nodes"].asString(), located); + auto edges = as_dataset(tens, topmd["edges"].asString(), located); return PointGraph(nodes, edges); } diff --git a/aux/src/TensorDMpointtree.cxx b/aux/src/TensorDMpointtree.cxx index b137eefde..e7719c502 100644 --- a/aux/src/TensorDMpointtree.cxx +++ b/aux/src/TensorDMpointtree.cxx @@ -1,6 +1,7 @@ #include "WireCellAux/TensorDMpointtree.h" #include "WireCellAux/TensorDMcommon.h" #include "WireCellAux/SimpleTensor.h" +#include using namespace WireCell; @@ -8,17 +9,17 @@ using namespace WireCell::Aux; using namespace WireCell::Aux::TensorDM; WireCell::Aux::TensorDM::named_pointclouds_t -WireCell::Aux::TensorDM::as_pcnamedset(const ITensor::vector& tens, const std::string& datapath) +WireCell::Aux::TensorDM::as_pcnamedset(const ITensor::vector& tens, const std::string& datapath, const located_t& located) { named_pointclouds_t ret; - auto top = top_tensor(tens, "pcnamedset", datapath); + auto top = top_tensor(tens, "pcnamedset", datapath, located); const auto& md = top->metadata(); auto items = md["items"]; for (const auto& name : items.getMemberNames()) { const auto path = items[name].asString(); - auto ds = as_dataset(tens, path); + auto ds = as_dataset(tens, path, located); ret.emplace(name, ds); } return ret; @@ -74,17 +75,17 @@ WireCell::Aux::TensorDM::as_pctree(const ITensor::vector& tens, std::unordered_map nodes_by_datapath; Points::node_t::owned_ptr root; - + const auto& located = index_datapaths(tens); std::function dochildren; dochildren = [&](const std::string& dpath) -> void { - auto top = top_tensor(tens, "pctreenode", dpath); + auto top = top_tensor(tens, "pctreenode", dpath, located); auto const& md = top->metadata(); named_pointclouds_t pcns; - if (! md["pointclouds"] ) { - pcns = as_pcnamedset(tens, md["pointclouds"].asString()); + if (! md["pointclouds"].asString().empty() ) { + pcns = as_pcnamedset(tens, md["pointclouds"].asString(), located); } std::string ppath = md["parent"].asString(); diff --git a/aux/test/doctest_tensordm_dataset.cxx b/aux/test/doctest_tensordm_dataset.cxx index 0311261b6..159e4377a 100644 --- a/aux/test/doctest_tensordm_dataset.cxx +++ b/aux/test/doctest_tensordm_dataset.cxx @@ -235,7 +235,8 @@ TEST_CASE("tensordm dataset dataset roundtrip") // back to dataset bool share = false; - Dataset d2 = TensorDM::as_dataset(itsp, datapath, share); + const auto& located = TensorDM::index_datapaths(itenvec); + Dataset d2 = TensorDM::as_dataset(itsp, datapath, located, share); CHECK(d2.size_major() == 3); diff --git a/cfg/pgrapher/experiment/dune-vd/dnnroi.jsonnet b/cfg/pgrapher/experiment/dune-vd/dnnroi.jsonnet index 61cc61e35..b38779795 100644 --- a/cfg/pgrapher/experiment/dune-vd/dnnroi.jsonnet +++ b/cfg/pgrapher/experiment/dune-vd/dnnroi.jsonnet @@ -80,7 +80,7 @@ function (anode, ts, prefix="dnnroi", output_scale=1.0) local retagger = pg.pnode({ type: 'Retagger', - name: 'dnnroi', + name: 'dnnroi%d' % apaid, data: { // Note: retagger keeps tag_rules an array to be like frame fanin/fanout. tag_rules: [{ diff --git a/img/inc/WireCellImg/MultiAlgBlobClustering.h b/img/inc/WireCellImg/MultiAlgBlobClustering.h new file mode 100644 index 000000000..c323c593d --- /dev/null +++ b/img/inc/WireCellImg/MultiAlgBlobClustering.h @@ -0,0 +1,44 @@ +#ifndef WIRECELLIMG_MULTIALGBLOBCLUSTERING +#define WIRECELLIMG_MULTIALGBLOBCLUSTERING + +#include "WireCellIface/ITensorSetFilter.h" +#include "WireCellIface/IConfigurable.h" +#include "WireCellAux/Logger.h" + +namespace WireCell::Img { + + class MultiAlgBlobClustering : public Aux::Logger, public ITensorSetFilter, public IConfigurable { + public: + MultiAlgBlobClustering(); + virtual ~MultiAlgBlobClustering() = default; + + virtual void configure(const WireCell::Configuration& cfg); + virtual WireCell::Configuration default_configuration() const; + + virtual bool operator()(const input_pointer& in, output_pointer& out); + + private: + // Count how many times we are called + size_t m_count{0}; + + /** Config: "inpath" + * + * The datapath for the input point graph data. This may be a + * regular expression which will be applied in a first-match + * basis against the input tensor datapaths. If the matched + * tensor is a pcdataset it is interpreted as providing the + * nodes dataset. Otherwise the matched tensor must be a + * pcgraph. + */ + std::string m_inpath{".*"}; + + /** Config: "outpath" + * + * The datapath for the resulting pcdataset. A "%d" will be + * interpolated with the ident number of the input tensor set. + */ + std::string m_outpath{""}; + }; +} + +#endif diff --git a/img/inc/WireCellImg/PointCloudFacade.h b/img/inc/WireCellImg/PointCloudFacade.h new file mode 100644 index 000000000..816b670a8 --- /dev/null +++ b/img/inc/WireCellImg/PointCloudFacade.h @@ -0,0 +1,42 @@ +/** + * + */ + +#ifndef WIRECELLIMG_POINTCLOUDFACADE +#define WIRECELLIMG_POINTCLOUDFACADE + +#include "WireCellUtil/PointCloudDataset.h" +#include "WireCellUtil/PointTree.h" +#include "WireCellUtil/Point.h" + +namespace WireCell::PointCloud { + using node_t = WireCell::PointCloud::Tree::Points::node_t; + using node_ptr = std::unique_ptr; + using Point = WireCell::Point; + + class Cluster { + public: + Cluster(const node_ptr& n) + : m_node(n.get()) + { + } + WireCell::PointCloud::Point calc_ave_pos(const Point& origin, const double dis) const; + + private: + node_t* m_node; /// do not own + }; + + class Blob { + public: + Blob(const node_ptr& n) + : m_node(n.get()) + { + } + double center_pos() const; + + private: + node_t* m_node; /// do not own + }; +} // namespace WireCell::PointCloud + +#endif \ No newline at end of file diff --git a/img/inc/WireCellImg/PointTreeBuilding.h b/img/inc/WireCellImg/PointTreeBuilding.h new file mode 100644 index 000000000..84afc5ef0 --- /dev/null +++ b/img/inc/WireCellImg/PointTreeBuilding.h @@ -0,0 +1,50 @@ +/** Sample blobs to make point cloud tree and output as tensors. + * Same as PointTreeBuilding but use ICluster as input. +*/ +#ifndef WIRECELL_IMG_POINTTREEBUILDING +#define WIRECELL_IMG_POINTTREEBUILDING + +#include "WireCellIface/IClusterTensorSet.h" +#include "WireCellIface/IBlobSampler.h" +#include "WireCellIface/IConfigurable.h" +#include "WireCellAux/Logger.h" + + +namespace WireCell::Img { + + class PointTreeBuilding : public Aux::Logger, public IClusterTensorSet, public IConfigurable + { + public: + PointTreeBuilding(); + virtual ~PointTreeBuilding(); + + // IConfigurable + virtual void configure(const WireCell::Configuration& cfg); + virtual WireCell::Configuration default_configuration() const; + + virtual bool operator()(const input_pointer& icluster, output_pointer& tensorset); + + private: + + /** Configuration: "samplers" + + An object with attributes providing names of + IBlobSamplers. The attribute names will be used to name + the point cloud produced by the samplers. + */ + std::map m_samplers; + + /** Config: "datapath" + + Set the datapath for the tensor representing the point + cloud tree. If a %d format code is found it wil be + interpolated with the IBlobSet::ident() value. + */ + std::string m_datapath = "pointtrees/%d"; + + size_t m_count{0}; + + }; +} + +#endif diff --git a/img/src/MultiAlgBlobClustering.cxx b/img/src/MultiAlgBlobClustering.cxx new file mode 100644 index 000000000..77445fe4b --- /dev/null +++ b/img/src/MultiAlgBlobClustering.cxx @@ -0,0 +1,70 @@ +#include "WireCellImg/MultiAlgBlobClustering.h" +#include "WireCellUtil/NamedFactory.h" +#include "WireCellAux/TensorDMpointtree.h" +#include "WireCellAux/TensorDMdataset.h" +#include "WireCellAux/TensorDMcommon.h" +#include "WireCellAux/SimpleTensorSet.h" + +WIRECELL_FACTORY(MultiAlgBlobClustering, WireCell::Img::MultiAlgBlobClustering, + WireCell::INamed, + WireCell::ITensorSetFilter, + WireCell::IConfigurable) + +using namespace WireCell; +using namespace WireCell::Img; +using namespace WireCell::Aux; +using namespace WireCell::Aux::TensorDM; + +MultiAlgBlobClustering::MultiAlgBlobClustering() + : Aux::Logger("MultiAlgBlobClustering", "img") +{ +} + +void MultiAlgBlobClustering::configure(const WireCell::Configuration& cfg) +{ + m_inpath = get(cfg, "inpath", m_inpath); + m_outpath = get(cfg, "outpath", m_outpath); +} + +WireCell::Configuration MultiAlgBlobClustering::default_configuration() const +{ + Configuration cfg; + cfg["inpath"] = m_inpath; + cfg["outpath"] = m_outpath; + return cfg; +} + + +bool MultiAlgBlobClustering::operator()(const input_pointer& ints, output_pointer& outts) +{ + outts = nullptr; + if (!ints) { + log->debug("EOS at call {}", m_count++); + return true; + } + + const int ident = ints->ident(); + std::string inpath = m_inpath; + if (inpath.find("%") != std::string::npos) { + inpath = String::format(inpath, ident); + } + + const auto& intens = *ints->tensors(); + log->debug("After merging, Got {} tensors", intens.size()); + const auto& root = as_pctree(intens, inpath); + if (!root) { + log->error("Failed to get point cloud tree from \"{}\"", inpath); + return false; + } + log->debug("Got pctree with {} children", root->children().size()); + + std::string outpath = m_outpath; + if (outpath.find("%") != std::string::npos) { + outpath = String::format(outpath, ident); + } + auto outtens = as_tensors(*root.get(), outpath); + log->debug("Made {} tensors", outtens.size()); + outts = as_tensorset(outtens, ints->ident()); + + return true; +} diff --git a/img/src/PointCloudFacade.cxx b/img/src/PointCloudFacade.cxx new file mode 100644 index 000000000..5de3bfb58 --- /dev/null +++ b/img/src/PointCloudFacade.cxx @@ -0,0 +1,49 @@ +#include "WireCellImg/PointCloudFacade.h" + +using namespace WireCell; +using namespace WireCell::PointCloud; +using namespace WireCell::PointCloud::Tree; // for "Points" node value type + +#include "WireCellUtil/Logging.h" +using spdlog::debug; + +WireCell::PointCloud::Point Cluster::calc_ave_pos(const Point& origin, const double dis) const { + spdlog::set_level(spdlog::level::debug); // Set global log level to debug + /// FIXME: there are many assumptions made, shoud we check these assumptions? + /// a bit worriying about the speed. + Scope scope = { "3d", {"x","y","z"} }; + const auto& sv = m_node->value.scoped_view(scope); + const auto& skd = sv.kd(); + auto rad = skd.radius(dis, origin); + const auto& snodes = sv.nodes(); + std::set maj_inds; + for (size_t pt_ind = 0; pt_indvalue.local_pcs(); + const auto& pc_scalar = lpcs.at("scalar"); + const auto charge = pc_scalar.get("charge")->elements()[0]; + const auto center_x = pc_scalar.get("center_x")->elements()[0]; + const auto center_y = pc_scalar.get("center_y")->elements()[0]; + const auto center_z = pc_scalar.get("center_z")->elements()[0]; + debug("charge {} center {{{} {} {}}}", charge, center_x, center_y, center_z); + Point inc(center_x, center_y, center_z); + inc = inc * charge; + ret += inc; + total_charge += charge; + } + if (total_charge != 0) { + ret = ret / total_charge; + } + debug("ret {{{} {} {}}}", ret.x(), ret.y(), ret.z()); + return ret; +} + +double Blob::center_pos() const { return 0; } \ No newline at end of file diff --git a/img/src/PointGraphProcessor.cxx b/img/src/PointGraphProcessor.cxx index 62cd9f547..59c2c120f 100644 --- a/img/src/PointGraphProcessor.cxx +++ b/img/src/PointGraphProcessor.cxx @@ -49,6 +49,7 @@ bool PointGraphProcessor::operator()(const input_pointer& ints, output_pointer& } const auto& tens = *ints->tensors(); + const auto& located = index_datapaths(tens); const auto pcgs = match_at(tens, m_inpath); if (pcgs.size() != 1) { outts = std::make_shared(ints->ident()); @@ -68,7 +69,7 @@ bool PointGraphProcessor::operator()(const input_pointer& ints, output_pointer& pg = as_pointgraph(tens, datapath); } else if (datatype == "pcdataset") { - auto nodes = as_dataset(tens, datapath); + auto nodes = as_dataset(tens, datapath, located); pg = PointGraph(nodes); } else { diff --git a/img/src/PointTreeBuilding.cxx b/img/src/PointTreeBuilding.cxx new file mode 100644 index 000000000..871f53d95 --- /dev/null +++ b/img/src/PointTreeBuilding.cxx @@ -0,0 +1,164 @@ +#include "WireCellImg/PointTreeBuilding.h" +#include "WireCellImg/Projection2D.h" +#include "WireCellUtil/PointTree.h" +#include "WireCellUtil/GraphTools.h" +#include "WireCellUtil/NamedFactory.h" + +#include "WireCellAux/ClusterHelpers.h" +#include "WireCellAux/TensorDMpointtree.h" +#include "WireCellAux/TensorDMcommon.h" + +WIRECELL_FACTORY(PointTreeBuilding, WireCell::Img::PointTreeBuilding, + WireCell::INamed, + WireCell::IClusterTensorSet, + WireCell::IConfigurable) + +using namespace WireCell; +using namespace WireCell::GraphTools; +using namespace WireCell::Img; +using WireCell::Img::Projection2D::get_geom_clusters; +using namespace WireCell::Aux; +using namespace WireCell::Aux::TensorDM; +using namespace WireCell::PointCloud::Tree; + +PointTreeBuilding::PointTreeBuilding() + : Aux::Logger("PointTreeBuilding", "img") +{ +} + + +PointTreeBuilding::~PointTreeBuilding() +{ +} + + +void PointTreeBuilding::configure(const WireCell::Configuration& cfg) +{ + m_datapath = get(cfg, "datapath", m_datapath); + auto samplers = cfg["samplers"]; + if (samplers.isNull()) { + raise("add at least one entry to the \"samplers\" configuration parameter"); + } + + for (auto name : samplers.getMemberNames()) { + auto tn = samplers[name].asString(); + if (tn.empty()) { + raise("empty type/name for sampler \"%s\"", name); + } + log->debug("point cloud \"{}\" will be made by sampler \"{}\"", + name, tn); + m_samplers[name] = Factory::find_tn(tn); + } + +} + + +WireCell::Configuration PointTreeBuilding::default_configuration() const +{ + Configuration cfg; + // eg: + // cfg["samplers"]["samples"] = "BlobSampler"; + cfg["datapath"] = m_datapath; + return cfg; +} + +namespace { + + std::string dump_ds(const WireCell::PointCloud::Dataset& ds) { + std::stringstream ss; + for (const auto& key : ds.keys()) {; + const auto& arr = ds.get(key); + ss << " {" << key << ":" << arr->dtype() << ":" << arr->shape()[0] << "} "; + // const auto& arr = ds.get(key)->elements(); + // for(auto elem : arr) { + // ss << elem << " "; + // } + } + return ss.str(); + } + // dump a NaryTree node + std::string dump_node(const WireCell::PointCloud::Tree::Points::node_t* node) + { + std::stringstream ss; + ss << "node: " << node; + if (node) { + const auto& lpcs = node->value.local_pcs(); + ss << " with " << lpcs.size() << " local pcs"; + for (const auto& [name, pc] : lpcs) { + ss << " " << name << ": " << dump_ds(pc); + } + } else { + ss << " null"; + } + return ss.str(); + } + // dump childs of a NaryTree node + std::string dump_children(const WireCell::PointCloud::Tree::Points::node_t* root) + { + std::stringstream ss; + ss << "NaryTree: " << root->children().size() << " children"; + const Points::node_ptr& first = root->children().front(); + ss << dump_node(first.get()); + return ss.str(); + } +} + +bool PointTreeBuilding::operator()(const input_pointer& icluster, output_pointer& tensorset) +{ + tensorset = nullptr; + + if (!icluster) { + log->debug("EOS at call {}", m_count++); + return true; + } + + const auto& gr = icluster->graph(); + log->debug("load cluster {} at call={}: {}", icluster->ident(), m_count, dumps(gr)); + + auto clusters = get_geom_clusters(gr); + log->debug("got {} clusters", clusters.size()); + size_t nblobs = 0; + Points::node_ptr root = std::make_unique(); + for (auto& [cluster_id, vdescs] : clusters) { + auto cnode = root->insert(std::move(std::make_unique())); + for (const auto& vdesc : vdescs) { + const char code = gr[vdesc].code(); + if (code != 'b') { + continue; + } + auto iblob = std::get(gr[vdesc].ptr); + named_pointclouds_t pcs; + for (auto& [name, sampler] : m_samplers) { + /// TODO: use nblobs or iblob->ident()? + pcs.emplace(name, sampler->sample_blob(iblob, nblobs)); + } + cnode->insert(Points(pcs)); + // log->debug("pcs {} cnode {}", pcs.size(), dump_children(cnode)); + // for (const auto& [name, pc] : pcs) { + // log->debug("{} -> keys {} size_major {}", name, pc.keys().size(), pc.size_major()); + // } + ++nblobs; + } + /// DEBUGONLY + // if (nblobs > 1000) { + // break; + // } + } + + const int ident = icluster->ident(); + std::string datapath = m_datapath; + if (datapath.find("%") != std::string::npos) { + datapath = String::format(datapath, ident); + } + auto tens = as_tensors(*root.get(), datapath); + log->debug("Made {} tensors", tens.size()); + tensorset = as_tensorset(tens, ident); + + log->debug("sampled {} blobs from set {} making {} tensors at call {}", + nblobs, ident, tens.size(), m_count++); + + return true; +} + + + diff --git a/img/test/doctest_clustering_prototype.cxx b/img/test/doctest_clustering_prototype.cxx index 277843488..f8facb8cd 100644 --- a/img/test/doctest_clustering_prototype.cxx +++ b/img/test/doctest_clustering_prototype.cxx @@ -1,10 +1,10 @@ #include "WireCellUtil/PointTree.h" #include "WireCellUtil/PointTesting.h" - #include "WireCellUtil/doctest.h" - #include "WireCellUtil/Logging.h" +#include "WireCellImg/PointCloudFacade.h" + #include using namespace WireCell; @@ -55,7 +55,13 @@ Points::node_ptr make_simple_pctree() /// here, we should be multiplying by some [length] unit. auto* n1 = root->insert(Points({ - {"center", make_janky_track(Ray(Point(0.5, 0, 0), Point(0.7, 0, 0)))}, + /// QUESTION: proper Array initiation? + {"scalar", Dataset({ + {"charge", Array({1.0})}, + {"center_x", Array({0.5})}, + {"center_y", Array({0.})}, + {"center_z", Array({0.})}, + })}, {"3d", make_janky_track(Ray(Point(0, 0, 0), Point(1, 0, 0)))} })); @@ -65,7 +71,12 @@ Points::node_ptr make_simple_pctree() // Ibid from a different track auto* n2 = root->insert(Points({ - {"center", make_janky_track(Ray(Point(1.5, 0, 0), Point(1.7, 0, 0)))}, + {"scalar", Dataset({ + {"charge", Array({2.0})}, + {"center_x", Array({1.5})}, + {"center_y", Array({0.})}, + {"center_z", Array({0.})}, + })}, {"3d", make_janky_track(Ray(Point(1, 0, 0), Point(2, 0, 0)))} })); @@ -78,7 +89,7 @@ Points::node_ptr make_simple_pctree() return root; } -TEST_CASE("PointCloudFacade test") +TEST_CASE("test PointTree API") { spdlog::set_level(spdlog::level::debug); // Set global log level to debug auto root = make_simple_pctree(); @@ -158,3 +169,12 @@ TEST_CASE("PointCloudFacade test") CHECK(rad.size() == 2); } + + +TEST_CASE("test PointCloudFacade") +{ + spdlog::set_level(spdlog::level::debug); // Set global log level to debug + auto root = make_simple_pctree(); + Cluster pcc(root); + auto ave_pos = pcc.calc_ave_pos({1,0,0}, 1); +} \ No newline at end of file diff --git a/sio/inc/WireCellSio/TensorFileSink.h b/sio/inc/WireCellSio/TensorFileSink.h index 075a59f3d..ac2fa3fe7 100644 --- a/sio/inc/WireCellSio/TensorFileSink.h +++ b/sio/inc/WireCellSio/TensorFileSink.h @@ -77,6 +77,9 @@ namespace WireCell::Sio { */ std::string m_prefix{""}; + // if true, disable writing to file which makes it faster to debug + bool m_dump_mode{false}; + private: using ostream_t = boost::iostreams::filtering_ostream; diff --git a/sio/src/TensorFileSink.cxx b/sio/src/TensorFileSink.cxx index aa2cd9492..053e3f6a0 100644 --- a/sio/src/TensorFileSink.cxx +++ b/sio/src/TensorFileSink.cxx @@ -42,6 +42,8 @@ void TensorFileSink::configure(const WireCell::Configuration& cfg) m_prefix = get(cfg, "prefix", m_prefix); log->debug("sink through {} filters to {} with prefix \"{}\"", m_out.size(), m_outname, m_prefix); + m_dump_mode = get(cfg, "dump_mode", m_dump_mode); + log->debug("dump_mode={}", m_dump_mode); } void TensorFileSink::finalize() @@ -82,6 +84,12 @@ bool TensorFileSink::operator()(const ITensorSet::pointer &in) return true; } + if(m_dump_mode) { + log->debug("dumping tensor set ident={} at call {}", + in->ident(), m_count); + return true; + } + const std::string pre = m_prefix + "tensor"; const std::string sident = std::to_string(in->ident()); auto tens = in->tensors();