136 vector<art::Ptr<Hit>> hitlist;
137 vector<vector<size_t>> idsmap;
138 vector<NuGraphInput> graphinputs;
139 _loaderTool->loadData(e, hitlist, graphinputs, idsmap);
141 if (
debug) std::cout <<
"Hits size=" << hitlist.size() << std::endl;
142 if (hitlist.size() <
minHits) {
153 const vector<int32_t>* spids =
nullptr;
154 const vector<int32_t>* hitids_u =
nullptr;
155 const vector<int32_t>* hitids_v =
nullptr;
156 const vector<int32_t>* hitids_y =
nullptr;
157 const vector<int32_t>* hit_plane =
nullptr;
158 const vector<float>* hit_time =
nullptr;
159 const vector<int32_t>* hit_wire =
nullptr;
160 const vector<float>* hit_integral =
nullptr;
161 const vector<float>* hit_rms =
nullptr;
162 for (
const auto& gi : graphinputs) {
163 if (gi.input_name ==
"spacepoint_table_spacepoint_id")
164 spids = &gi.input_int32_vec;
165 else if (gi.input_name ==
"spacepoint_table_hit_id_u")
166 hitids_u = &gi.input_int32_vec;
167 else if (gi.input_name ==
"spacepoint_table_hit_id_v")
168 hitids_v = &gi.input_int32_vec;
169 else if (gi.input_name ==
"spacepoint_table_hit_id_y")
170 hitids_y = &gi.input_int32_vec;
171 else if (gi.input_name ==
"hit_table_local_plane")
172 hit_plane = &gi.input_int32_vec;
173 else if (gi.input_name ==
"hit_table_local_time")
174 hit_time = &gi.input_float_vec;
175 else if (gi.input_name ==
"hit_table_local_wire")
176 hit_wire = &gi.input_int32_vec;
177 else if (gi.input_name ==
"hit_table_integral")
178 hit_integral = &gi.input_float_vec;
179 else if (gi.input_name ==
"hit_table_rms")
180 hit_rms = &gi.input_float_vec;
184 vector<size_t> idsmapRev(hitlist.size(), hitlist.size());
185 for (
const auto& ipv : idsmap) {
186 for (
size_t ih = 0; ih < ipv.size(); ih++) {
187 idsmapRev[ipv[ih]] = ih;
196 if (this->n1 == other.n1 && this->n2 == other.n2)
204 auto start_preprocess1 = std::chrono::high_resolution_clock::now();
205 vector<vector<Edge>> edge2d(
planes.size(), vector<Edge>());
206 for (
size_t p = 0; p <
planes.size(); p++) {
207 vector<double> coords;
208 for (
size_t i = 0; i < hit_plane->size(); ++i) {
209 if (
size_t(hit_plane->at(i)) != p)
continue;
210 coords.push_back(hit_time->at(i) *
pos_norm[1]);
211 coords.push_back(hit_wire->at(i) *
pos_norm[0]);
213 if (
debug) std::cout <<
"Plane " << p <<
" has N hits=" << coords.size() / 2 << std::endl;
214 if (coords.size() / 2 < 3) {
continue; }
215 delaunator::Delaunator
d(coords);
216 if (
debug) std::cout <<
"Found N triangles=" <<
d.triangles.size() / 3 << std::endl;
217 for (std::size_t i = 0; i <
d.triangles.size(); i += 3) {
220 e.n1 =
d.triangles[i];
221 e.n2 =
d.triangles[i + 1];
222 edge2d[p].push_back(e);
223 e.n1 =
d.triangles[i + 1];
224 e.n2 =
d.triangles[i];
225 edge2d[p].push_back(e);
227 e.n1 =
d.triangles[i];
228 e.n2 =
d.triangles[i + 2];
229 edge2d[p].push_back(e);
230 e.n1 =
d.triangles[i + 2];
231 e.n2 =
d.triangles[i];
232 edge2d[p].push_back(e);
234 e.n1 =
d.triangles[i + 1];
235 e.n2 =
d.triangles[i + 2];
236 edge2d[p].push_back(e);
237 e.n1 =
d.triangles[i + 2];
238 e.n2 =
d.triangles[i + 1];
239 edge2d[p].push_back(e);
243 std::sort(edge2d[p].
begin(), edge2d[p].
end(), [](
const auto& i,
const auto& j) {
244 return (i.n1 != j.n1 ? i.n1 < j.n1 : i.n2 < j.n2);
247 for (
auto& e : edge2d[p]) {
248 std::cout <<
"sorted plane=" << p <<
" e1=" << e.n1 <<
" e2=" << e.n2 << std::endl;
251 edge2d[p].erase(std::unique(edge2d[p].
begin(), edge2d[p].
end()), edge2d[p].
end());
255 for (
size_t p = 0; p <
planes.size(); p++) {
256 for (
auto& e : edge2d[p]) {
257 std::cout <<
" plane=" << p <<
" e1=" << e.n1 <<
" e2=" << e.n2 << std::endl;
261 auto end_preprocess1 = std::chrono::high_resolution_clock::now();
262 std::chrono::duration<double> elapsed_preprocess1 = end_preprocess1 - start_preprocess1;
265 auto start_preprocess2 = std::chrono::high_resolution_clock::now();
266 vector<vector<Edge>> edge3d(
planes.size(), vector<Edge>());
267 for (
size_t i = 0; i < spids->size(); ++i) {
268 if (hitids_u->at(i) >= 0) {
270 e.n1 = idsmapRev[hitids_u->at(i)];
272 edge3d[0].push_back(e);
274 if (hitids_v->at(i) >= 0) {
276 e.n1 = idsmapRev[hitids_v->at(i)];
278 edge3d[1].push_back(e);
280 if (hitids_y->at(i) >= 0) {
282 e.n1 = idsmapRev[hitids_y->at(i)];
284 edge3d[2].push_back(e);
289 auto x = torch::Dict<std::string, torch::Tensor>();
290 auto batch = torch::Dict<std::string, torch::Tensor>();
291 for (
size_t p = 0; p <
planes.size(); p++) {
292 vector<float> nodeft;
293 for (
size_t i = 0; i < hit_plane->size(); ++i) {
294 if (
size_t(hit_plane->at(i)) != p)
continue;
295 nodeft.push_back((hit_wire->at(i) *
pos_norm[0] -
avgs[hit_plane->at(i)][0]) /
296 devs[hit_plane->at(i)][0]);
297 nodeft.push_back((hit_time->at(i) *
pos_norm[1] -
avgs[hit_plane->at(i)][1]) /
298 devs[hit_plane->at(i)][1]);
299 nodeft.push_back((hit_integral->at(i) -
avgs[hit_plane->at(i)][2]) /
300 devs[hit_plane->at(i)][2]);
301 nodeft.push_back((hit_rms->at(i) -
avgs[hit_plane->at(i)][3]) /
devs[hit_plane->at(i)][3]);
303 long int dim = nodeft.size() / 4;
304 torch::Tensor ix = torch::zeros({dim, 4}, torch::dtype(torch::kFloat32));
306 std::cout <<
"plane=" << p << std::endl;
307 std::cout << std::scientific;
308 for (
size_t n = 0;
n < nodeft.size();
n =
n + 4) {
309 std::cout << nodeft[
n] <<
" " << nodeft[
n + 1] <<
" " << nodeft[
n + 2] <<
" " 310 << nodeft[
n + 3] <<
" " << std::endl;
313 for (
size_t n = 0;
n < nodeft.size();
n =
n + 4) {
314 ix[
n / 4][0] = nodeft[
n];
315 ix[
n / 4][1] = nodeft[
n + 1];
316 ix[
n / 4][2] = nodeft[
n + 2];
317 ix[
n / 4][3] = nodeft[
n + 3];
320 torch::Tensor ib = torch::zeros({dim}, torch::dtype(torch::kInt64));
321 batch.insert(
planes[p], ib);
324 auto edge_index_plane = torch::Dict<std::string, torch::Tensor>();
325 for (
size_t p = 0; p <
planes.size(); p++) {
326 long int dim = edge2d[p].size();
327 torch::Tensor ix = torch::zeros({2, dim}, torch::dtype(torch::kInt64));
328 for (
size_t n = 0;
n < edge2d[p].size();
n++) {
329 ix[0][
n] = int(edge2d[p][
n].n1);
330 ix[1][
n] = int(edge2d[p][
n].n2);
332 edge_index_plane.insert(
planes[p], ix);
334 std::cout <<
"plane=" << p << std::endl;
335 std::cout <<
"2d edge size=" << edge2d[p].size() << std::endl;
336 for (
size_t n = 0;
n < edge2d[p].size();
n++) {
337 std::cout << edge2d[p][
n].n1 <<
" ";
339 std::cout << std::endl;
340 for (
size_t n = 0;
n < edge2d[p].size();
n++) {
341 std::cout << edge2d[p][
n].n2 <<
" ";
343 std::cout << std::endl;
347 auto edge_index_nexus = torch::Dict<std::string, torch::Tensor>();
348 for (
size_t p = 0; p <
planes.size(); p++) {
349 long int dim = edge3d[p].size();
350 torch::Tensor ix = torch::zeros({2, dim}, torch::dtype(torch::kInt64));
351 for (
size_t n = 0;
n < edge3d[p].size();
n++) {
352 ix[0][
n] = int(edge3d[p][
n].n1);
353 ix[1][
n] = int(edge3d[p][
n].n2);
355 edge_index_nexus.insert(
planes[p], ix);
357 std::cout <<
"plane=" << p << std::endl;
358 std::cout <<
"3d edge size=" << edge3d[p].size() << std::endl;
359 for (
size_t n = 0;
n < edge3d[p].size();
n++) {
360 std::cout << edge3d[p][
n].n1 <<
" ";
362 std::cout << std::endl;
363 for (
size_t n = 0;
n < edge3d[p].size();
n++) {
364 std::cout << edge3d[p][
n].n2 <<
" ";
366 std::cout << std::endl;
370 long int spdim = spids->size();
371 auto nexus =
torch::empty({spdim, 0}, torch::dtype(torch::kFloat32));
373 std::vector<torch::jit::IValue> inputs;
375 inputs.push_back(edge_index_plane);
376 inputs.push_back(edge_index_nexus);
377 inputs.push_back(nexus);
378 inputs.push_back(batch);
381 auto end_preprocess2 = std::chrono::high_resolution_clock::now();
382 std::chrono::duration<double> elapsed_preprocess2 = end_preprocess2 - start_preprocess2;
383 if (
debug) std::cout <<
"FORWARD!" << std::endl;
384 auto start = std::chrono::high_resolution_clock::now();
385 auto outputs =
model.forward(inputs).toGenericDict();
386 auto end = std::chrono::high_resolution_clock::now();
387 std::chrono::duration<double> elapsed =
end - start;
389 std::cout <<
"Time taken for inference: " 390 << elapsed_preprocess1.count() + elapsed_preprocess2.count() + elapsed.count()
391 <<
" seconds" << std::endl;
392 std::cout <<
"output =" << outputs << std::endl;
398 vector<NuGraphOutput> infer_output;
399 for (
const auto& elem1 : outputs) {
400 if (elem1.value().isTensor()) {
401 torch::Tensor tensor = elem1.value().toTensor();
402 std::vector<float> vec(tensor.data_ptr<
float>(), tensor.data_ptr<
float>() + tensor.numel());
403 infer_output.push_back(
NuGraphOutput(elem1.key().to<std::string>(), vec));
405 else if (elem1.value().isGenericDict()) {
406 for (
const auto& elem2 : elem1.value().toGenericDict()) {
407 torch::Tensor tensor = elem2.value().toTensor();
408 std::vector<float> vec(tensor.data_ptr<
float>(), tensor.data_ptr<
float>() + tensor.numel());
409 infer_output.push_back(
410 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 
