Logo coherent WaveBurst  
Library Reference Guide
Logo
wavegraph.cc
Go to the documentation of this file.
1 /*
2 # Copyright (C) 2019 Eric Chassande-Mottin, Philippe Bacon, Gayathri V, Archana Pai
3 #
4 # This program is free software: you can redistribute it and/or modify
5 # it under the terms of the GNU General Public License as published by
6 # the Free Software Foundation, either version 3 of the License, or
7 # (at your option) any later version.
8 #
9 # This program is distributed in the hope that it will be useful,
10 # but WITHOUT ANY WARRANTY; without even the implied warranty of
11 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 # GNU General Public License for more details.
13 #
14 # You should have received a copy of the GNU General Public License
15 # along with this program. If not, see <https://www.gnu.org/licenses/>.
16 */
17 
18 
19 #include "wavegraph.hh"
20 
21 // Predicate function used to trim lines of text from double white spaces
22 bool bothspaces(char lhs, char rhs) { return (lhs == rhs) && (lhs == ' '); }
23 
24 // Symbols for admissible keys in the header of the graph file
32 };
33 
34 
35 graph_header_code header_hash (std::string const& str) {
36 // Hash function for the header of the graph file
37  if (str == "stride") return xstride;
38  if (str == "span") return xspan;
39  if (str == "scalemin") return xscalemin;
40  if (str == "scalemax") return xscalemax;
41  if (str == "nscales") return xnscales;
42  if (str == "samp_freq") return xsamp_freq;
43 }
44 
45 ClassImp(wavegraph)
46 
47 ////////////////////////////////////////////////////////////////////////////////
48 /* BEGIN_HTML
49 
50 <p>Wavegraph is a graph-based clustering algorithm for coherent WaveBurst (cWB).
51 Wavegraph generates structured clusters of pixels in the time-frequency-scale
52 energy maps computed by cWB. The structure of the clusters is dictated by an
53 auxiliary user-provided directed graph. Admissible clusters are one-dimensional
54 paths in this graph. Wavegraph finds paths of connected pixels that maximizes a
55 figure-of-merit related to the cumulated sum of the time-frequency-scale energy
56 along the path plus a regularization term.</p>
57 
58 <p>Wavegraph clustering analyses the data segment by successive blocks. The
59 block duration corresponds to that of the input graph. The algorithm collects
60 paths from the analysis of each blocks and then it sorts them and selects
61 distinct path that do not overlap significantly in time.</p>
62 
63 <p>A graph connects nodes together. A node is associated to a
64 (time-frequency-scale) pixel. The node object contains only geometrical
65 information. For instance it includes the list of nodes to which it is connected
66 (its ancestors). The pixel object contains physical information (center of time,
67 frequency and scale bin in physical units) but it does not include any
68 connectivity information.</p>
69 
70 <p>The wavegraph class includes methods for building and searching the graph.
71 The wavegraph object is an ordered lists of nodes with some auxiliary information.</p>
72 
73 END_HTML */
74 ////////////////////////////////////////////////////////////////////////////////
75 
76 void
77 wavegraph::add_node(const int nodeidx, const int timeidx, const int freqidx, const int scaleidx,
78  const double value_avg, const double value_stdev,
79  const bool endnode, const nodeids& node_ancestors) {
80  // Add a node to graph object
81  //
82  // Input: nodeidx = node index (should match the node rank in the graph -- this is for debugging purpose)
83  // timeidx = time index of the node
84  // freqidx = frequency index of the node
85  // scaleidx = scale index of the node
86  // value_avg = average value of the node
87  // value_stdev = value standard deviation of the node
88  // endnode = flag saying if the node is an endnode or not
89  // node_ancestors = list of indices of the node ancestors
90 
91  wavenode node = {
92  nodeidx, // node id in graph
93  timeidx, // time index
94  freqidx, // freq index
95  scaleidx, // scale index
96  value_avg, // node average value
97  value_stdev, // value standard deviation
98  endnode, // flag is true if node is end node
99  node_ancestors // ancestors of the node
100  };
101 
102  add_node(node);
103 
104 }
105 
106 void
108  // Add a node to graph object
109  //
110  // Input: node = node to append to graph
111 
112  graph.push_back(node);
113 
114 }
115 
116 void
117 wavegraph::create(const std::string& srcfile){
118 // Create the graph from the given source file
119 //
120 // Input: srcfile = file name
121 //
122 // The file format is as follows
123 //
124 // # Comments
125 // %% stride=XX span=XX scalemin=XX scalemax=XX nscales=XX samp_freq=XX
126 // 0 II TT FF SS E AA BB CC ... << first node
127 // 1 II TT FF SS E AA BB CC ... << second node
128 // ... etc ...
129 //
130 // II, TT, FF and SS: ID and time, freq and scale indices (integers)
131 // E: true if the node is an end node (boolean)
132 // AA, BB, CC, ...: node ancestors (there may be none)
133 //
134 
135 #ifdef NDEBUG
136  std::cout << "debug flag1 create: open file\n";
137 #endif
138 
139  std::ifstream source;
140  source.open(srcfile.c_str(), std::ifstream::in);
141 
142  clear();
143 
144 #ifdef NDEBUG
145  std::cout << "debug flag2 create: read file\n";
146 #endif
147 
148  std::string line;
149  int count_nodes=0;
150 
151  // read line
152  while ( std::getline(source,line) ){
153 
154  // ignore if empty line
155  if (line.empty())
156  continue;
157 
158  // ignore if commented
159  if (line[0] == '#')
160  continue;
161 
162  // parse header
163  if ((line[0] == '%')&(line[1] == '%')) {
164 
165  // trim leading characters and multiple spaces
166  line.erase(0,3);
167  std::string::iterator line_end = std::unique(line.begin(), line.end(), bothspaces);
168  line.erase(line_end, line.end());
169 
170  std::string delimiter_fields = " ";
171  std::string delimiter_values = "=";
172  size_t pos_field = 0;
173  while (true) {
174 
175  // extract field
176  pos_field = line.find(delimiter_fields);
177  std::string field = line.substr(0, pos_field);
178  if (field.empty())
179  continue;
180 
181  //std::cout << "line: " << line << "\n";
182  //std::cout << "field: " << field << "\n";
183 
184  // read token and values from field
185  size_t pos_value = field.find(delimiter_values);
186  std::string token = field.substr(0, pos_value);
187  field.erase(0, pos_value + delimiter_values.length());
188  std::string value = field.substr(0, pos_value);
189 
190  // set variables
191  switch (header_hash(token)) {
192  case xstride:
193  stride=atoi(value.c_str());
194  break;
195  case xspan:
196  span=atoi(value.c_str());
197  break;
198  case xscalemin:
199  scalemin=atoi(value.c_str());
200  break;
201  case xscalemax:
202  scalemax=atoi(value.c_str());
203  break;
204  case xnscales:
205  nscales=atoi(value.c_str());
206  break;
207  case xsamp_freq:
208  samp_freq=atoi(value.c_str());
209  break;
210  default:
211  std::cout << "Invalid token: " << token << "\n";
212  throw std::invalid_argument("Error: Graph file has an invalid header");
213  }
214 
215  // break if all done
216  if (pos_field == std::string::npos)
217  break;
218 
219  // erase processed field from header line
220  line.erase(0, pos_field + delimiter_fields.length());
221  }
222 
223 #ifdef NDEBUG
224  std::cout << "header: stride= " << stride << " span=" << span;
225  std::cout << " scalemin=" << scalemin << " scalemax=" << scalemax;
226  std::cout << " nscales=" << nscales << " samp_freq=" << samp_freq << "\n";
227 #endif
228 
229  // go to next line
230  continue;
231  }
232 
233 #ifdef NDEBUG
234  std::cout << line << "\n";
235 #endif
236 
237  // parse data
238  nodeids ancestors;
239  std::stringstream streamline(line);
240 
241  // read nodeid
242  int nodeid;
243  streamline >> nodeid;
244 
245  if (nodeid != count_nodes)
246  throw std::invalid_argument("Error: Invalid file");
247 
248  // read time index
249  int time_idx;
250  streamline >> time_idx;
251 
252  // read freq index
253  int freq_idx;
254  streamline >> freq_idx;
255 
256  // read scale index
257  int scale_idx;
258  streamline >> scale_idx;
259 
260  // read average graph value -- not used by wavegraph
261  double value_avg;
262  streamline >> value_avg;
263 
264  // read standard deviation of graph value -- not used by wavegraph
265  double value_stdev;
266  streamline >> value_stdev;
267 
268  // read end_node flag
269  bool endnode_flag;
270  streamline >> endnode_flag;
271 
272  // read node ancestors sequentially
273  while (true) {
274  int node;
275  bool read_ok = static_cast<bool> (streamline >> node);
276  if (!read_ok)
277  break;
278  ancestors.insert(node);
279  }
280 
281  // add node to graph
282  add_node(nodeid, time_idx, freq_idx, scale_idx, value_avg, value_stdev, endnode_flag, ancestors);
283  ancestors.clear();
284 
285  count_nodes++;
286  }
287 
288 #ifdef NDEBUG
289  print();
290 #endif
291 
292  source.close();
293 }
294 
295 void
297  // Print graph to standard output
298 
299  std::cout << "stride=" << stride << " span=" << span;
300  std::cout << " scalemin=" << scalemin << " scalemax=" << scalemax;
301  std::cout << " nscales=" << nscales << "samp_freq" << samp_freq << "\n";
302 
303  for (int node_id=0; node_id != graph.size(); ++node_id) {
304  std::cout << node_id << " (" << nodeid(node_id) << "): ";
305 
306  if (get_ancestors(node_id).empty()) {
307  std::cout << " [no ancestor] \n";
308  continue;
309  }
310 
311  nodeids::iterator anc;
312  for (anc = get_ancestors(node_id).begin(); anc != get_ancestors(node_id).end(); ++anc)
313  std::cout << *anc << " ";
314 
315  if (is_endnode(node_id))
316  std::cout << "*";
317 
318  std::cout << "\n";
319  }
320 }
321 
322 bool
324 // Test whether the node order defined by the node IDs is a topological order.
325 // The first nodes have no ancestor. The graph must be acyclic.
326 //
327 // Output: result of the test
328 
329  std::set<int> marked_nodes;
330 
331  bool result=true;
332 
333  // loop on nodes
334  for (int nodeid=0; nodeid != graph.size(); ++nodeid) {
335 
336  // std::set<node_type>::iterator node;
337  // std::cout << "looking to " << nodeid << ":\n";
338  // for (node = marked_nodes.begin(); node != marked_nodes.end(); ++node) {
339  // std::cout << *node << " ";
340  // }
341  // std::cout << "\n";
342 
343  // if node has no ancestor, mark the node and continue
344  if (get_ancestors(nodeid).empty()) {
345  marked_nodes.insert(nodeid);
346  continue;
347  }
348 
349  // loop on the node ancestors
350  nodeids::iterator anc;
351  for (anc = get_ancestors(nodeid).begin(); anc != get_ancestors(nodeid).end(); ++anc) {
352  // if one of the ancestors is NOT in marked_nodes, the graph is NOT topologically sorted
353  if ( marked_nodes.find(*anc) == marked_nodes.end() )
354  return false;
355  }
356  // mark current node
357  marked_nodes.insert(nodeid);
358 
359  }
360 
361  // the graph is topologically sorted
362  return true;
363 };
364 
365 void
366 wavegraph::heaviest_path(std::vector<double>& total_weight, std::vector<int>& total_length,
367  std::vector<int>& best_ancestor, const std::vector<double> weights) {
368 // Calculation of the "heaviest path" from the entry nodes (nodes with no
369 // ancestor) to each node of the graph, that has the largest total weight ( =
370 // sum of the weights of the node in the path). Uses dynamic programming.
371 //
372 // Input: weights = mapping from the node to its weight
373 // Output: total_weight = mapping from each node to the total weight of the heaviest path
374 // Output: total_length = mapping from each node to the total length of the heaviest path
375 // best_ancestor = mapping from each node to the node ancestor in the heaviest path
376 // (or -1 if none)
377 //
378 
379  for (int nodeid=0; nodeid != graph.size(); ++nodeid) {
380 
381 #ifdef NDEBUG
382  std::cout << "debug flag1 heaviest_path: looking at node " << nodeid << ":";
383 #endif
384 
385  nodeids ancestors=get_ancestors(nodeid);
386 
387  // case with no ancestors
388  if (ancestors.empty()) {
389 #ifdef NDEBUG
390  std::cout << "[no ancestors]\n";
391 #endif
392  total_weight[nodeid] = weights[nodeid];
393  total_length[nodeid] = 1;
394  best_ancestor[nodeid] = -1;
395  continue;
396  }
397 
398  // case with one ancestor
399  nodeids::const_iterator first = ancestors.begin();
400  if (ancestors.size() == 1) {
401 #ifdef NDEBUG
402  std::cout << "[one ancestor] select " << *first <<"\n";
403 #endif
404  total_weight.at(nodeid) = total_weight.at(*first)+weights.at(nodeid);
405  total_length.at(nodeid) = total_length.at(*first)+1;
406  best_ancestor.at(nodeid) = *first;
407  continue;
408  }
409 
410 #ifdef NDEBUG
411  std::cout << "[many ancestors] find max ";
412 #endif
413 
414  // case with many ancestors
415  double best_weight = total_weight.at(*first);
416  int best_length = total_length.at(*first);
417  int best_anc = *first;
418  for (nodeids::const_iterator anc = ancestors.begin(); anc != ancestors.end(); ++anc) {
419 
420  if (total_weight.at(*anc) > best_weight) {
421  best_weight = total_weight.at(*anc);
422  best_length = total_length.at(*anc);
423  best_anc = *anc;
424  }
425  }
426 
427 #ifdef NDEBUG
428  std::cout << "select " << best_anc << "\n";
429 #endif
430 
431  total_weight.at(nodeid) = best_weight+weights.at(nodeid);
432  total_length.at(nodeid) = best_length+1;
433  best_ancestor.at(nodeid) = best_anc;
434 
435  }
436 }
437 
438 bool wavegraph::is_compatible(const std::vector< WSeries<double>* >& data)
439 {
440 
441  // assert that the graph and cWB data cube have the same scale axis
442  if (log2(get_scale(data.front())) != scalemin) {
443 
444 #ifdef NDEBUG
445  std::cout << "data scalemin= " << log2(get_scale(data.front())) << ", graph scalemin= " << scalemin << "\n";
446 #endif
447 
448  std::cerr << "Error: min scale of graph does not match that of data\n";
449  return false;
450  }
451 
452  if (log2(get_scale(data.back())) != scalemax) {
453 
454 #ifdef NDEBUG
455  std::cout << "data scalemax= " << log2(get_scale(data.back())) << ", graph scalemax= " << scalemax << "\n";
456 #endif
457 
458  std::cerr << "Error: max scale of graph does not match that of data\n";
459  return false;
460  }
461 
462  if (data.size() != nscales) {
463 #ifdef NDEBUG
464  std::cout << "data nscales= " << data.size() << ", graph nscales= " << nscales << "\n";
465 #endif
466 
467  std::cerr << "Error: scale axis of graph does not match that of data\n";
468  return false;
469  }
470 
471  if (nscales != scalemax-scalemin+1) {
472  std::cerr << "Error: scale axis sampling is not contiguous\n";
473  return false;
474  }
475 
476  if (data.front()->rate() != samp_freq) {
477  std::cerr << "Error: sampling frequency of the graph does not match that of data\n";
478  return false;
479  }
480 
481  return true;
482 
483 }
484 
485 double median(std::vector<double> & vec)
486 {
487  size_t n = vec.size() / 2;
488  nth_element(vec.begin(), vec.begin()+n, vec.end());
489  return vec[n];
490 }
491 
492 std::vector<cluster>
493 wavegraph::clustering(const double threshold, const std::vector< WSeries<double>* >& data, const double strip_edges, const int path_halfwidth, const double penal_factor, const std::vector<double>& energy_thresholds){
494 // Compute clusters using wavegraph clustering. Main steps in pseudo-code:
495 //
496 // while data are available (avoid edges as required by strip_edges)
497 // fill the weight table of the graph with current data block
498 // compute the heaviest path for this block
499 // if its total weight > threshold
500 // trace back (the list of nodes of) this path and store
501 // endif
502 // go to next block
503 // endwhile
504 // sort paths by their total weight
505 // convert paths into clusters
506 // retain univocal clusters that do not overlap in time
507 //
508 // Input: threshold = selection threshold on path total weight
509 // data = list of Wilson Daubechies Meyer transforms computed by cWB
510 // Each WSeries in the list corresponds to a given scale (or "resolution").
511 // The list should go contiguously from the smallest to the largest scale.
512 // strip_edges = remove strip_edges seconds at the segment edges to prevent contamination
513 // by edge artifacts
514 // path_halfwidth = half width of the path (path weight is computed summing weights of +/- path_halfwidth pixels)
515 // penal_factor = leading factor of length penalization term
516 // energy_thresholds = vector of thresholds optionally applied to pixel energy
517 // one threshold per scale (no selection if threshold is 0)
518 
519  if (! is_compatible(data))
520  throw std::invalid_argument("Error: graph and data are incompatible");
521 
522  // compute initial and end offsets to strip the edges from analysis
523  int strip_scalemax = 0;
524  if (strip_edges > 0)
525  strip_scalemax = strip_edges*(2*data.back()->resolution());
526 
527 #ifdef NDEBUG
528  std::cout << "debug flag0 clustering: main loop\n";
529  std::cout << "strip_max = " << strip_scalemax << "\n";
530 #endif
531 
532  // compute median noise level for each frequency and for each scale
533  std::vector < std::vector<double> > noise_levels;
534  for (int k=0; k<data.size(); k++) {
535  std::vector<double> noise_level(num_of_freq_bins(data[k]));
536  for (int m=0; m<num_of_freq_bins(data[k]); m++) {
537  std::vector<double> non_zero_data;
538  for (int n=0; n<num_of_time_bins(data[k]); n++) {
539  double value = get_map00(data[k],n,m);
540  if (value > 0)
541  non_zero_data.push_back(value);
542  }
543  if (non_zero_data.empty())
544  noise_level.push_back(0.0);
545  else
546  noise_level.push_back(median(non_zero_data));
547  }
548  noise_levels.push_back(noise_level);
549  }
550 
551 #ifdef NDEBUG
552  std::cout << "noise level estimates:\n";
553  for (int k=0; k<noise_levels.size(); k++) {
554  std::cout << "scale=" << k << ": ";
555  for (int m=0; m<noise_levels[k].size(); m++) {
556  std::cout << noise_levels[k].at(m) << " ";
557  }
558  std::cout << "\n";
559  }
560 #endif
561 
562  //
563  // main loop: this loop runs by blocks of pixels in the max scale plane.
564  // The max scale is the last element in the data vector [hence data.back()]
565  //
566  std::vector<wavepath> selected_paths;
567 
568  int end_scalemax = num_of_time_bins(data.back())-get_span()-strip_scalemax;
569  for (int block_offset_scalemax = strip_scalemax; block_offset_scalemax <= end_scalemax; block_offset_scalemax += get_stride()) {
570 
571 #ifdef NDEBUG
572  std::cout << "debug flag1 clustering: feed weight table\n";
573 #endif
574 
575  // fill the weight table of the graph with current data block
576  std::vector<double> weights(size(),0.0);
577  for (int nodeid=0; nodeid != size(); ++nodeid) {
578  int node_offset = block_offset_scalemax * pow(2,scalemax-(scalemin+scale(nodeid)));
579  int time_index = node_offset + time(nodeid);
580 #ifdef NDEBUG
581  //std::cout << "debug flag2 clustering: accessing pixel t=" << node_offset+time(nodeid);
582  //std::cout << " f=" << freq(nodeid) << " scale=" << scale(nodeid) << " with block_offset=" << block_offset_scalemax << "\n";
583  if ((node_offset+time(nodeid) >= num_of_time_bins(data.at(scale(nodeid)))) | (freq(nodeid) >= num_of_freq_bins(data.at(scale(nodeid))))) {
584  std::cout << "debug flag2 clustering: accessing pixel t=" << node_offset+time(nodeid);
585  std::cout << " f=" << freq(nodeid) << " scale=" << scale(nodeid) << " with block_offset=" << block_offset_scalemax << "/" << end_scalemax << "\n";
586  std::cout << "at scale index= " << scale(nodeid);
587  std::cout << ": num of time bins= " << num_of_time_bins(data.at(scale(nodeid)));
588  std::cout << ", num of freq bins= " << num_of_freq_bins(data.at(scale(nodeid))) << "\n";
589  throw std::string("Error: requested position exceed allocated limits");
590  }
591 #endif
592  double pixel_data = get_map00(data.at(scale(nodeid)),time_index,freq(nodeid));
593  weights[nodeid]= pixel_data > energy_thresholds.at(scale(nodeid)) ? pixel_data : 0.0 ; // apply energy threshold
594  // if non zero path_width, sum over
595  for (int n=1; n<=path_halfwidth; ++n) {
596  pixel_data = get_map00(data.at(scale(nodeid)),time_index-n,freq(nodeid));
597  weights[nodeid]+= pixel_data > energy_thresholds.at(scale(nodeid)) ? pixel_data : 0.0 ; // apply energy threshold
598  pixel_data = get_map00(data.at(scale(nodeid)),time_index+n,freq(nodeid));
599  weights[nodeid]+= pixel_data > energy_thresholds.at(scale(nodeid)) ? pixel_data : 0.0 ; // apply energy threshold
600  }
601  weights[nodeid]-=penal_factor*(2*path_halfwidth+1)*noise_levels.at(scale(nodeid)).at(freq(nodeid));
602  }
603 
604 #ifdef NDEBUG
605  std::cout << "debug flag2 clustering: dynamic programming\n";
606 #endif
607 
608  // compute the heaviest path for this block
609  std::vector<double> total_weights(size());
610  std::vector<int> total_lengths(size());;
611  std::vector<int> best_ancestors(size());;
612  heaviest_path(total_weights,total_lengths,best_ancestors,weights);
613 
614 #ifdef NDEBUG
615  std::cout << "debug flag3 clustering: select and reconstruct best paths\n";
616 #endif
617 
618  // select path in block with largest total_weight
619  int best_node_id=-1;
620  double best_weight=0;
621  for (int node_id=0; node_id != size(); ++node_id) {
622 
623  // if its total weight > path_length * threshold
624  if (is_endnode(node_id) & (total_weights[node_id] >= total_lengths[node_id]*(2*path_halfwidth+1)*threshold)) {
625  if (total_weights[node_id] > best_weight) {
626  best_weight = total_weights[node_id];
627  best_node_id = node_id;
628  }
629  }
630  }
631 
632  // trace back (the list of nodes of) this path and store
633  if (best_node_id>0) {
634  wavepath path;
635  path.clear();
636 
637  int node=best_node_id;
638  // paths use the maximum scale as the reference scale
639  path.init(get_node(node),block_offset_scalemax,nscales-1,total_weights[node]);
640 
641 #ifdef NDEBUG
642  std::cout << "path: " << node;
643 #endif
644  while (true) {
645  node = best_ancestors.at(node);
646 #ifdef NDEBUG
647  std::cout << "->" << node;
648 #endif
649  if (node == -1)
650  break;
651  path.add_node(get_node(node));
652  }
653 #ifdef NDEBUG
654  std::cout << "\n";
655  path.print();
656 #endif
657 
658  // store the best path (largest total weight)
659  selected_paths.push_back(path);
660  }
661 
662 #ifdef NDEBUG
663  std::cout << selected_paths.size() << " paths selected\n";
664 #endif
665 
666  // go to next block
667  //block_offset_scalemax+=get_stride();
668  }
669 
670 #ifdef NDEBUG
671  std::cout << "debug flag4 clustering: sort selected paths\n";
672 #endif
673 
674  // sort paths by their total weight
675  std::sort(selected_paths.begin(), selected_paths.end(), compare_paths);
676 
677 #ifdef NDEBUG
678  std::cout << selected_paths.size() << " selected paths\n";
679  // wavepath best_path=selected_paths.front();
680  // std::cout << "best path is:\n";
681  // best_path.print();
682 #endif
683 
684  // convert paths to clusters
685  std::vector<cluster> clusters;
686  std::vector<wavepath>::const_iterator this_path;
687  for (this_path = selected_paths.begin(); this_path != selected_paths.end(); ++this_path)
688  clusters.push_back(((wavepath) *this_path).convert_to_cluster(log2(get_scale(data.front())),data.back()->resolution(),path_halfwidth));
689 
690 #ifdef NDEBUG
691  std::cout << clusters.size() << " selected clusters\n";
692 #endif
693 
694  // retain univocal clusters that are distinct in time
695  std::vector<cluster> univocal_clusters;
696  univocal_clusters=select_clusters_distinct_in_time(clusters,get_span()/(2*data.back()->resolution()));
697 
698 #ifdef NDEBUG
699  std::cout << univocal_clusters.size() << " univocal clusters\n";
700  if ( !univocal_clusters.empty() ) {
701  cluster best_cluster=univocal_clusters.front();
702  std::cout << "best cluster is:\n";
703  cluster::const_iterator this_pixel;
704  for (this_pixel=best_cluster.begin(); this_pixel != best_cluster.end(); ++this_pixel) {
705  std::cout << "timeix=" << ((pixel) *this_pixel).timeix << " freqix=" << ((pixel) *this_pixel).freqix << " scaleix=" << ((pixel) *this_pixel).scaleix << " : ";
706  std::cout << " time=" << ((pixel) *this_pixel).time << " freq=" << ((pixel) *this_pixel).freq << " log_2 scale=" << ((pixel) *this_pixel).log2scale << "\n";
707  }
708  }
709 #endif
710 
711  return univocal_clusters;
712 
713 }
std::vector< cluster > clustering(const double threshold, const std::vector< WSeries< double > * > &data, const double strip_edges, const int path_halfwidth, const double penal_factor, const std::vector< double > &energy_thresholds)
Definition: wavegraph.cc:493
par [0] value
bool is_compatible(const std::vector< WSeries< double > * > &data)
Definition: wavegraph.cc:438
std::vector< pixel > cluster
Definition: wavepath.hh:61
int num_of_time_bins(WSeries< double > *WS)
Definition: wavegraph.hh:51
int n
Definition: cwb_net.C:28
int time(int id)
Definition: wavegraph.hh:101
int nodeid(int id)
Definition: wavegraph.hh:99
double median(std::vector< double > &vec)
Definition: wavegraph.cc:485
int nscales
Definition: wavegraph.hh:129
int m
Definition: cwb_net.C:28
std::vector< cluster > select_clusters_distinct_in_time(std::vector< cluster > &sorted_clusters, const double precision)
Definition: wavepath.cc:166
void heaviest_path(std::vector< double > &total_weight, std::vector< int > &total_length, std::vector< int > &best_ancestor, const std::vector< double > weights)
Definition: wavegraph.cc:366
float get_map00(WSeries< double > *WS, int index)
Definition: wavegraph.hh:58
void init(const wavenode node, const int offset, const int refscaleidx, const double path_weight)
Definition: wavepath.cc:63
nodes graph
Definition: wavegraph.hh:124
graph_header_code
Definition: wavegraph.cc:25
wavenode get_node(int id)
Definition: wavegraph.hh:95
char str[1024]
int stride
Definition: wavegraph.hh:126
int size()
Definition: wavegraph.hh:93
int scalemax
Definition: wavegraph.hh:128
std::set< int > nodeids
Definition: wavepath.hh:37
bool is_topologically_sorted()
Definition: wavegraph.cc:323
int num_of_freq_bins(WSeries< double > *WS)
Definition: wavegraph.hh:48
bool bothspaces(char lhs, char rhs)
Definition: wavegraph.cc:22
void add_node(const int nodeidx, const int timeidx, const int freqidx, const int scaleidx, const double value_avg, const double value_stdev, const bool endnode, const nodeids &ancestors)
Definition: wavegraph.cc:77
int is_endnode(int id)
Definition: wavegraph.hh:107
nodeids get_ancestors(int id)
Definition: wavegraph.hh:97
int k
int scalemin
Definition: wavegraph.hh:127
int get_span()
Definition: wavegraph.hh:109
TObjArray * token
void create(const std::string &srcfile)
Definition: wavegraph.cc:117
int get_scale(WSeries< double > *WS)
Definition: wavegraph.hh:45
void clear()
Definition: wavepath.hh:84
ifstream in
int get_stride()
Definition: wavegraph.hh:111
bool compare_paths(const wavepath &path1, const wavepath &path2)
Definition: wavepath.cc:23
int samp_freq
Definition: wavegraph.hh:130
void add_node(const wavenode node)
Definition: wavepath.cc:81
void clear()
Definition: wavegraph.hh:83
int scale(int id)
Definition: wavegraph.hh:105
void print()
Definition: wavegraph.cc:296
char line[1024]
Definition: mdc.hh:231
graph_header_code header_hash(std::string const &str)
Definition: wavegraph.cc:35
int freq(int id)
Definition: wavegraph.hh:103
void print()
Definition: wavepath.cc:91