135 vector<art::Ptr<Hit>> hitlist;
136 vector<vector<size_t>> idsmap;
137 vector<NuGraphInput> graphinputs;
138 _loaderTool->loadData(e, hitlist, graphinputs, idsmap);
140 if (
debug) std::cout <<
"Hits size=" << hitlist.size() << std::endl;
141 if (hitlist.size() <
minHits) {
152 const vector<int32_t>* spids =
nullptr;
153 const vector<int32_t>* hitids_u =
nullptr;
154 const vector<int32_t>* hitids_v =
nullptr;
155 const vector<int32_t>* hitids_y =
nullptr;
156 const vector<int32_t>* hit_plane =
nullptr;
157 const vector<float>* hit_time =
nullptr;
158 const vector<int32_t>* hit_wire =
nullptr;
159 const vector<float>* hit_integral =
nullptr;
160 const vector<float>* hit_rms =
nullptr;
161 for (
const auto& gi : graphinputs) {
162 if (gi.input_name ==
"spacepoint_table_spacepoint_id")
163 spids = &gi.input_int32_vec;
164 else if (gi.input_name ==
"spacepoint_table_hit_id_u")
165 hitids_u = &gi.input_int32_vec;
166 else if (gi.input_name ==
"spacepoint_table_hit_id_v")
167 hitids_v = &gi.input_int32_vec;
168 else if (gi.input_name ==
"spacepoint_table_hit_id_y")
169 hitids_y = &gi.input_int32_vec;
170 else if (gi.input_name ==
"hit_table_local_plane")
171 hit_plane = &gi.input_int32_vec;
172 else if (gi.input_name ==
"hit_table_local_time")
173 hit_time = &gi.input_float_vec;
174 else if (gi.input_name ==
"hit_table_local_wire")
175 hit_wire = &gi.input_int32_vec;
176 else if (gi.input_name ==
"hit_table_integral")
177 hit_integral = &gi.input_float_vec;
178 else if (gi.input_name ==
"hit_table_rms")
179 hit_rms = &gi.input_float_vec;
183 vector<size_t> idsmapRev(hitlist.size(), hitlist.size());
184 for (
const auto& ipv : idsmap) {
185 for (
size_t ih = 0; ih < ipv.size(); ih++) {
186 idsmapRev[ipv[ih]] = ih;
195 if (this->n1 == other.n1 && this->n2 == other.n2)
203 auto start_preprocess1 = std::chrono::high_resolution_clock::now();
204 vector<vector<Edge>> edge2d(
planes.size(), vector<Edge>());
205 for (
size_t p = 0; p <
planes.size(); p++) {
206 vector<double> coords;
207 for (
size_t i = 0; i < hit_plane->size(); ++i) {
208 if (
size_t(hit_plane->at(i)) != p)
continue;
209 coords.push_back(hit_time->at(i) *
pos_norm[1]);
210 coords.push_back(hit_wire->at(i) *
pos_norm[0]);
212 if (
debug) std::cout <<
"Plane " << p <<
" has N hits=" << coords.size() / 2 << std::endl;
213 if (coords.size() / 2 < 3) {
continue; }
214 delaunator::Delaunator
d(coords);
215 if (
debug) std::cout <<
"Found N triangles=" <<
d.triangles.size() / 3 << std::endl;
216 for (std::size_t i = 0; i <
d.triangles.size(); i += 3) {
219 e.n1 =
d.triangles[i];
220 e.n2 =
d.triangles[i + 1];
221 edge2d[p].push_back(e);
222 e.n1 =
d.triangles[i + 1];
223 e.n2 =
d.triangles[i];
224 edge2d[p].push_back(e);
226 e.n1 =
d.triangles[i];
227 e.n2 =
d.triangles[i + 2];
228 edge2d[p].push_back(e);
229 e.n1 =
d.triangles[i + 2];
230 e.n2 =
d.triangles[i];
231 edge2d[p].push_back(e);
233 e.n1 =
d.triangles[i + 1];
234 e.n2 =
d.triangles[i + 2];
235 edge2d[p].push_back(e);
236 e.n1 =
d.triangles[i + 2];
237 e.n2 =
d.triangles[i + 1];
238 edge2d[p].push_back(e);
242 std::sort(edge2d[p].
begin(), edge2d[p].
end(), [](
const auto& i,
const auto& j) {
243 return (i.n1 != j.n1 ? i.n1 < j.n1 : i.n2 < j.n2);
246 for (
auto& e : edge2d[p]) {
247 std::cout <<
"sorted plane=" << p <<
" e1=" << e.n1 <<
" e2=" << e.n2 << std::endl;
250 edge2d[p].erase(std::unique(edge2d[p].
begin(), edge2d[p].
end()), edge2d[p].
end());
254 for (
size_t p = 0; p <
planes.size(); p++) {
255 for (
auto& e : edge2d[p]) {
256 std::cout <<
" plane=" << p <<
" e1=" << e.n1 <<
" e2=" << e.n2 << std::endl;
260 auto end_preprocess1 = std::chrono::high_resolution_clock::now();
261 std::chrono::duration<double> elapsed_preprocess1 = end_preprocess1 - start_preprocess1;
264 auto start_preprocess2 = std::chrono::high_resolution_clock::now();
265 vector<vector<Edge>> edge3d(
planes.size(), vector<Edge>());
266 for (
size_t i = 0; i < spids->size(); ++i) {
267 if (hitids_u->at(i) >= 0) {
269 e.n1 = idsmapRev[hitids_u->at(i)];
271 edge3d[0].push_back(e);
273 if (hitids_v->at(i) >= 0) {
275 e.n1 = idsmapRev[hitids_v->at(i)];
277 edge3d[1].push_back(e);
279 if (hitids_y->at(i) >= 0) {
281 e.n1 = idsmapRev[hitids_y->at(i)];
283 edge3d[2].push_back(e);
288 auto x = torch::Dict<std::string, torch::Tensor>();
289 auto batch = torch::Dict<std::string, torch::Tensor>();
290 for (
size_t p = 0; p <
planes.size(); p++) {
291 vector<float> nodeft;
292 for (
size_t i = 0; i < hit_plane->size(); ++i) {
293 if (
size_t(hit_plane->at(i)) != p)
continue;
294 nodeft.push_back((hit_wire->at(i) *
pos_norm[0] -
avgs[hit_plane->at(i)][0]) /
295 devs[hit_plane->at(i)][0]);
296 nodeft.push_back((hit_time->at(i) *
pos_norm[1] -
avgs[hit_plane->at(i)][1]) /
297 devs[hit_plane->at(i)][1]);
298 nodeft.push_back((hit_integral->at(i) -
avgs[hit_plane->at(i)][2]) /
299 devs[hit_plane->at(i)][2]);
300 nodeft.push_back((hit_rms->at(i) -
avgs[hit_plane->at(i)][3]) /
devs[hit_plane->at(i)][3]);
302 long int dim = nodeft.size() / 4;
303 torch::Tensor ix = torch::zeros({dim, 4}, torch::dtype(torch::kFloat32));
305 std::cout <<
"plane=" << p << std::endl;
306 std::cout << std::scientific;
307 for (
size_t n = 0;
n < nodeft.size();
n =
n + 4) {
308 std::cout << nodeft[
n] <<
" " << nodeft[
n + 1] <<
" " << nodeft[
n + 2] <<
" " 309 << nodeft[
n + 3] <<
" " << std::endl;
312 for (
size_t n = 0;
n < nodeft.size();
n =
n + 4) {
313 ix[
n / 4][0] = nodeft[
n];
314 ix[
n / 4][1] = nodeft[
n + 1];
315 ix[
n / 4][2] = nodeft[
n + 2];
316 ix[
n / 4][3] = nodeft[
n + 3];
319 torch::Tensor ib = torch::zeros({dim}, torch::dtype(torch::kInt64));
320 batch.insert(
planes[p], ib);
323 auto edge_index_plane = torch::Dict<std::string, torch::Tensor>();
324 for (
size_t p = 0; p <
planes.size(); p++) {
325 long int dim = edge2d[p].size();
326 torch::Tensor ix = torch::zeros({2, dim}, torch::dtype(torch::kInt64));
327 for (
size_t n = 0;
n < edge2d[p].size();
n++) {
328 ix[0][
n] = int(edge2d[p][
n].n1);
329 ix[1][
n] = int(edge2d[p][
n].n2);
331 edge_index_plane.insert(
planes[p], ix);
333 std::cout <<
"plane=" << p << std::endl;
334 std::cout <<
"2d edge size=" << edge2d[p].size() << std::endl;
335 for (
size_t n = 0;
n < edge2d[p].size();
n++) {
336 std::cout << edge2d[p][
n].n1 <<
" ";
338 std::cout << std::endl;
339 for (
size_t n = 0;
n < edge2d[p].size();
n++) {
340 std::cout << edge2d[p][
n].n2 <<
" ";
342 std::cout << std::endl;
346 auto edge_index_nexus = torch::Dict<std::string, torch::Tensor>();
347 for (
size_t p = 0; p <
planes.size(); p++) {
348 long int dim = edge3d[p].size();
349 torch::Tensor ix = torch::zeros({2, dim}, torch::dtype(torch::kInt64));
350 for (
size_t n = 0;
n < edge3d[p].size();
n++) {
351 ix[0][
n] = int(edge3d[p][
n].n1);
352 ix[1][
n] = int(edge3d[p][
n].n2);
354 edge_index_nexus.insert(
planes[p], ix);
356 std::cout <<
"plane=" << p << std::endl;
357 std::cout <<
"3d edge size=" << edge3d[p].size() << std::endl;
358 for (
size_t n = 0;
n < edge3d[p].size();
n++) {
359 std::cout << edge3d[p][
n].n1 <<
" ";
361 std::cout << std::endl;
362 for (
size_t n = 0;
n < edge3d[p].size();
n++) {
363 std::cout << edge3d[p][
n].n2 <<
" ";
365 std::cout << std::endl;
369 long int spdim = spids->size();
370 auto nexus =
torch::empty({spdim, 0}, torch::dtype(torch::kFloat32));
372 std::vector<torch::jit::IValue> inputs;
374 inputs.push_back(edge_index_plane);
375 inputs.push_back(edge_index_nexus);
376 inputs.push_back(nexus);
377 inputs.push_back(batch);
380 auto end_preprocess2 = std::chrono::high_resolution_clock::now();
381 std::chrono::duration<double> elapsed_preprocess2 = end_preprocess2 - start_preprocess2;
382 if (
debug) std::cout <<
"FORWARD!" << std::endl;
383 auto start = std::chrono::high_resolution_clock::now();
384 auto outputs =
model.forward(inputs).toGenericDict();
385 auto end = std::chrono::high_resolution_clock::now();
386 std::chrono::duration<double> elapsed =
end - start;
388 std::cout <<
"Time taken for inference: " 389 << elapsed_preprocess1.count() + elapsed_preprocess2.count() + elapsed.count()
390 <<
" seconds" << std::endl;
391 std::cout <<
"output =" << outputs << std::endl;
397 vector<NuGraphOutput> infer_output;
398 for (
const auto& elem1 : outputs) {
399 if (elem1.value().isTensor()) {
400 torch::Tensor tensor = elem1.value().toTensor();
401 std::vector<float> vec(tensor.data_ptr<
float>(), tensor.data_ptr<
float>() + tensor.numel());
402 infer_output.push_back(
NuGraphOutput(elem1.key().to<std::string>(), vec));
404 else if (elem1.value().isGenericDict()) {
405 for (
const auto& elem2 : elem1.value().toGenericDict()) {
406 torch::Tensor tensor = elem2.value().toTensor();
407 std::vector<float> vec(tensor.data_ptr<
float>(), tensor.data_ptr<
float>() + tensor.numel());
408 infer_output.push_back(
409 NuGraphOutput(elem1.key().to<std::string>() +
"_" + elem2.key().to<std::string>(), vec));
torch::jit::script::Module model
std::vector< std::unique_ptr< DecoderToolBase > > _decoderToolsVec
vector< std::string > planes
boost::graph_traits< ModuleGraph >::edge_descriptor Edge
decltype(auto) constexpr end(T &&obj)
ADL-aware version of std::end.
vector< vector< float > > devs
vector< vector< float > > avgs
decltype(auto) constexpr begin(T &&obj)
ADL-aware version of std::begin.
bool operator==(infinite_endcount_iterator< T > const &, count_iterator< T > const &)
decltype(auto) constexpr empty(T &&obj)
ADL-aware version of std::empty.
std::unique_ptr< LoaderToolBase > _loaderTool
1.8.11 
