Logo coherent WaveBurst  
Library Reference Guide
Logo
wavepath.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 "wavepath.hh"
20 
21 ClassImp(wavepath)
22 
23 bool compare_paths(const wavepath& path1, const wavepath& path2) {
24  return path1.get_weight() > path2.get_weight();
25 }
26 
27 ////////////////////////////////////////////////////////////////////////////////
28 /* BEGIN_HTML
29 
30 <p>Wavegraph is a graph-based clustering algorithm for coherent WaveBurst (cWB).
31 Wavegraph generates structured clusters of pixels in the time-frequency-scale
32 energy maps computed by cWB. The structure of the clusters is dictated by an
33 auxiliary user-provided directed graph. Admissible clusters are one-dimensional
34 paths in this graph. Wavegraph finds paths of connected pixels that maximizes a
35 figure-of-merit related to the cumulated sum of the time-frequency-scale energy
36 along the path plus a regularization term.</p>
37 
38 <p>Wavegraph clustering analyses the data segment by successive blocks. The
39 block duration corresponds to that of the input graph. The algorithm collects
40 paths from the analysis of each blocks and then it sorts them and selects
41 distinct path that do not overlap significantly in time.</p>
42 
43 <p>A graph connects nodes together. A node is associated to a
44 (time-frequency-scale) pixel. The node object contains only geometrical
45 information. For instance it includes the list of nodes to which it is connected
46 (its ancestors). The pixel object contains physical information (center of time,
47 frequency and scale bin in physical units) but it does not include any
48 connectivity information.</p>
49 
50 <p>The wavepath class defines methods for paths in the graph. Wavepath are
51 ordered lists of nodes with some auxiliary information.</p>
52 
53 <p>The cluster object is also an ordered list of pixels.</p>
54 
55 <p>We use a fiducial scale as reference for time measurement. XXX TIME
56 CONVERSION FORMULA HERE! XXX</p>
57 
58 END_HTML */
59 ////////////////////////////////////////////////////////////////////////////////
60 
61 
62 void
63 wavepath::init(const wavenode node, const int path_offset, const int path_refscaleidx, const double path_weight) {
64  // Initialize a path and append a first node
65  //
66  // Input: node = node to append
67  // weight = initial path weight
68  // offset = time origin of the block (number of samples in the reference scale plane
69  // from the beginning of the data segment)
70  // refscale_idx = index of the reference scale
71  //
72 
73  weight=path_weight;
74  offset=path_offset;
75  refscaleidx=path_refscaleidx;
76  add_node(node);
77 
78 }
79 
80 void
82  // Append one node to a path
83  //
84  // Input: node = node to append
85 
86  path.push_back(node);
87 
88 }
89 
90 void
92  // Print path to standard output
93 
94  for (unsigned int rank = length(); rank-- > 0; ) {
95  std::cout << "node " << nodeid(rank);
96  std::cout << " (t=" << time(rank)+offset*pow(2,refscaleidx-scale(rank)) << ",f=" << freq(rank) << ",s=" << scale(rank) << ") ->\n";
97  }
98  std::cout << "weight=" << weight << "\n";
99 
100 }
101 
102 cluster
103 wavepath::convert_to_cluster(const int scalemin,const double deltaf_refscale, const int path_width) {
104  // Convert a path object to a cluster object
105  //
106  // Input: scalemin = value of the minimum scale
107  // deltaf_refscale = frequency resolution at reference scale in Hertz
108  // path_width = width of the path in time
109  //
110  // The time offset of the path/cluster is computed in the reference scale plane
111  //
112  // Output: cluster object
113 
114  cluster out;
115  for (unsigned int rank = length(); rank-- > 0; ) {
116  int refscale_to_thisscale = pow(2,refscaleidx-scale(rank));
117  double deltaf_thisscale = deltaf_refscale * refscale_to_thisscale;
118 
119  // init and add pixel in the center
120  pixel path_center = {
121  scale(rank), // scale index
122  time(rank)+offset*refscale_to_thisscale, // time index
123  freq(rank), // freq index
124  0, // log2(scale)
125  nodeid(rank), // node ID
126  0.0, // time
127  0.0}; // freq
128  path_center.log2scale=scalemin+path_center.scaleix; // log2(scale)
129  path_center.time=path_center.timeix/(2*deltaf_thisscale); // time (sec)
130  path_center.freq=path_center.freqix*deltaf_thisscale; // freq (Hz)
131 
132  out.push_back(path_center);
133 
134  // add pixels on the side if width > 1
135  // XXX The node ID of the side pixels is that of the path center and do not match with the graph. XXX
136  // XXX This will lead to issues in the chi^2 test. XXX
137  for (int n=1; n<path_width; ++n) {
138  pixel path_side=path_center;
139  path_side.timeix++;
140  path_side.time=path_side.timeix/(2*deltaf_thisscale);
141  out.push_back(path_side);
142 
143  path_side=path_center;
144  path_side.timeix--;
145  path_side.time=path_side.timeix/(2*deltaf_thisscale);
146  out.push_back(path_side);
147  }
148 
149  }
150  return out;
151 }
152 
153 class is_coincident_with {
154 public:
155  is_coincident_with(double ref_time, double precision): _ref_time(ref_time),_precision(precision){}
156 
157  bool operator()(cluster const & x) const {
158  return abs(_ref_time-(x.back()).time) <= _precision;
159  }
160 private:
161  double _ref_time;
162  double _precision;
163 };
164 
165 std::vector<cluster>
166 select_clusters_distinct_in_time(std::vector<cluster>& sorted_clusters,const double precision) {
167  // Select the loudest univocal clusters that do not overlap significantly with any others
168  //
169  // Input: sorted_clusters = list of clusters sorted by descending weights
170  // precision = duration of the coincidence window according to which two clusters
171  // are said to be overlapping
172  //
173  // Output: list of univocal clusters
174 
175  std::vector<cluster> selected;
176 
177 #ifdef NDEBUG
178  std::cout << "debug flag: select univocal clusters\n";
179  std::cout << "total number of clusters is " << sorted_clusters.size() << "\n";
180 #endif
181 
182  while (! sorted_clusters.empty()) {
183  cluster first=sorted_clusters.front();
184  selected.push_back(first);
185  double ref_time=(first.back()).time;
186  sorted_clusters.erase(std::remove_if(sorted_clusters.begin(),sorted_clusters.end(),is_coincident_with(ref_time,precision)),sorted_clusters.end());
187 #ifdef NDEBUG
188  std::cout << "remaining number of clusters is " << sorted_clusters.size() << "\n";
189 #endif
190  }
191 
192 #ifdef NDEBUG
193  std::cout << "retained " << selected.size() << " univocal clusters\n";
194 #endif
195 
196  return selected;
197 
198 }
double time
Definition: wavepath.hh:56
int refscaleidx
Definition: wavepath.hh:99
std::vector< pixel > cluster
Definition: wavepath.hh:61
int n
Definition: cwb_net.C:28
ofstream out
Definition: cwb_merge.C:214
cout<< endl;cout<< "ts size = "<< ts.size()<< " ts rate = "<< ts.rate()<< endl;tf.Forward(ts, wdm);int levels=tf.getLevel();cout<< "tf size = "<< tf.size()<< endl;double dF=tf.resolution();double dT=1./(2 *dF);cout<< "rate(hz) : "<< RATE<< "\ layers : "<< nLAYERS<< "\ dF(hz) : "<< dF<< "\ dT(ms) : "<< dT *1000.<< endl;int itime=TIME_PIXEL_INDEX;int ifreq=FREQ_PIXEL_INDEX;int index=(levels+1) *itime+ifreq;double time=itime *dT;double freq=(ifreq >0) ? ifreq *dF :dF/4;cout<< endl;cout<< "PIXEL TIME = "<< time<< " sec "<< endl;cout<< "PIXEL FREQ = "<< freq<< " Hz "<< endl;cout<< endl;wavearray< double > x
wavearray< double > rank
Definition: Regression_H1.C:80
std::vector< cluster > select_clusters_distinct_in_time(std::vector< cluster > &sorted_clusters, const double precision)
Definition: wavepath.cc:166
int length()
Definition: wavepath.hh:69
void init(const wavenode node, const int offset, const int refscaleidx, const double path_weight)
Definition: wavepath.cc:63
double precision
int scaleix
Definition: wavepath.hh:51
int log2scale
Definition: wavepath.hh:54
int nodeid(int id)
Definition: wavepath.hh:71
int freqix
Definition: wavepath.hh:53
int time(int id)
Definition: wavepath.hh:73
int offset
Definition: wavepath.hh:97
nodes path
Definition: wavepath.hh:96
double weight
Definition: wavepath.hh:100
int timeix
Definition: wavepath.hh:52
bool compare_paths(const wavepath &path1, const wavepath &path2)
Definition: wavepath.cc:23
void add_node(const wavenode node)
Definition: wavepath.cc:81
double freq
Definition: wavepath.hh:57
int scale(int id)
Definition: wavepath.hh:77
cluster convert_to_cluster(const int scalemin, const double fs, const int path_width)
Definition: wavepath.cc:103
int freq(int id)
Definition: wavepath.hh:75
void print()
Definition: wavepath.cc:91