22 bool bothspaces(
char lhs,
char rhs) {
return (lhs == rhs) && (lhs ==
' '); }
37 if (str ==
"stride")
return xstride;
38 if (str ==
"span")
return xspan;
41 if (str ==
"nscales")
return xnscales;
78 const double value_avg,
const double value_stdev,
79 const bool endnode,
const nodeids& node_ancestors) {
112 graph.push_back(node);
136 std::cout <<
"debug flag1 create: open file\n";
145 std::cout <<
"debug flag2 create: read file\n";
152 while ( std::getline(source,line) ){
163 if ((line[0] ==
'%')&(line[1] ==
'%')) {
167 std::string::iterator line_end = std::unique(line.begin(), line.end(),
bothspaces);
168 line.erase(line_end, line.end());
170 std::string delimiter_fields =
" ";
171 std::string delimiter_values =
"=";
172 size_t pos_field = 0;
176 pos_field = line.find(delimiter_fields);
177 std::string field = line.substr(0, pos_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);
193 stride=atoi(value.c_str());
196 span=atoi(value.c_str());
211 std::cout <<
"Invalid token: " << token <<
"\n";
212 throw std::invalid_argument(
"Error: Graph file has an invalid header");
216 if (pos_field == std::string::npos)
220 line.erase(0, pos_field + delimiter_fields.length());
224 std::cout <<
"header: stride= " <<
stride <<
" span=" <<
span;
234 std::cout << line <<
"\n";
239 std::stringstream streamline(line);
245 if (nodeid != count_nodes)
246 throw std::invalid_argument(
"Error: Invalid file");
250 streamline >> time_idx;
254 streamline >> freq_idx;
258 streamline >> scale_idx;
262 streamline >> value_avg;
266 streamline >> value_stdev;
270 streamline >> endnode_flag;
275 bool read_ok =
static_cast<bool> (streamline >> node);
278 ancestors.insert(node);
282 add_node(nodeid, time_idx, freq_idx, scale_idx, value_avg, value_stdev, endnode_flag, ancestors);
299 std::cout <<
"stride=" <<
stride <<
" span=" <<
span;
303 for (
int node_id=0; node_id !=
graph.size(); ++node_id) {
304 std::cout << node_id <<
" (" <<
nodeid(node_id) <<
"): ";
307 std::cout <<
" [no ancestor] \n";
311 nodeids::iterator anc;
313 std::cout << *anc <<
" ";
329 std::set<int> marked_nodes;
345 marked_nodes.insert(
nodeid);
350 nodeids::iterator anc;
353 if ( marked_nodes.find(*anc) == marked_nodes.end() )
357 marked_nodes.insert(
nodeid);
367 std::vector<int>& best_ancestor,
const std::vector<double> weights) {
382 std::cout <<
"debug flag1 heaviest_path: looking at node " <<
nodeid <<
":";
388 if (ancestors.empty()) {
390 std::cout <<
"[no ancestors]\n";
394 best_ancestor[
nodeid] = -1;
399 nodeids::const_iterator first = ancestors.begin();
400 if (ancestors.size() == 1) {
402 std::cout <<
"[one ancestor] select " << *first <<
"\n";
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;
411 std::cout <<
"[many ancestors] find max ";
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) {
420 if (total_weight.at(*anc) > best_weight) {
421 best_weight = total_weight.at(*anc);
422 best_length = total_length.at(*anc);
428 std::cout <<
"select " << best_anc <<
"\n";
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;
445 std::cout <<
"data scalemin= " << log2(
get_scale(data.front())) <<
", graph scalemin= " <<
scalemin <<
"\n";
448 std::cerr <<
"Error: min scale of graph does not match that of data\n";
455 std::cout <<
"data scalemax= " << log2(
get_scale(data.back())) <<
", graph scalemax= " <<
scalemax <<
"\n";
458 std::cerr <<
"Error: max scale of graph does not match that of data\n";
464 std::cout <<
"data nscales= " << data.size() <<
", graph nscales= " <<
nscales <<
"\n";
467 std::cerr <<
"Error: scale axis of graph does not match that of data\n";
472 std::cerr <<
"Error: scale axis sampling is not contiguous\n";
477 std::cerr <<
"Error: sampling frequency of the graph does not match that of data\n";
487 size_t n = vec.size() / 2;
488 nth_element(vec.begin(), vec.begin()+
n, vec.end());
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){
520 throw std::invalid_argument(
"Error: graph and data are incompatible");
523 int strip_scalemax = 0;
525 strip_scalemax = strip_edges*(2*data.back()->resolution());
528 std::cout <<
"debug flag0 clustering: main loop\n";
529 std::cout <<
"strip_max = " << strip_scalemax <<
"\n";
533 std::vector < std::vector<double> > noise_levels;
534 for (
int k=0;
k<data.size();
k++) {
537 std::vector<double> non_zero_data;
541 non_zero_data.push_back(value);
543 if (non_zero_data.empty())
544 noise_level.push_back(0.0);
546 noise_level.push_back(
median(non_zero_data));
548 noise_levels.push_back(noise_level);
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) <<
" ";
566 std::vector<wavepath> selected_paths;
569 for (
int block_offset_scalemax = strip_scalemax; block_offset_scalemax <= end_scalemax; block_offset_scalemax +=
get_stride()) {
572 std::cout <<
"debug flag1 clustering: feed weight table\n";
576 std::vector<double> weights(
size(),0.0);
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";
589 throw std::string(
"Error: requested position exceed allocated limits");
593 weights[
nodeid]= pixel_data > energy_thresholds.at(
scale(
nodeid)) ? pixel_data : 0.0 ;
595 for (
int n=1;
n<=path_halfwidth; ++
n) {
597 weights[
nodeid]+= pixel_data > energy_thresholds.at(
scale(
nodeid)) ? pixel_data : 0.0 ;
599 weights[
nodeid]+= pixel_data > energy_thresholds.at(
scale(
nodeid)) ? pixel_data : 0.0 ;
605 std::cout <<
"debug flag2 clustering: dynamic programming\n";
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);
615 std::cout <<
"debug flag3 clustering: select and reconstruct best paths\n";
620 double best_weight=0;
621 for (
int node_id=0; node_id !=
size(); ++node_id) {
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;
633 if (best_node_id>0) {
637 int node=best_node_id;
642 std::cout <<
"path: " << node;
645 node = best_ancestors.at(node);
647 std::cout <<
"->" << node;
659 selected_paths.push_back(path);
663 std::cout << selected_paths.size() <<
" paths selected\n";
671 std::cout <<
"debug flag4 clustering: sort selected paths\n";
675 std::sort(selected_paths.begin(), selected_paths.end(),
compare_paths);
678 std::cout << selected_paths.size() <<
" selected paths\n";
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));
691 std::cout << clusters.size() <<
" selected clusters\n";
695 std::vector<cluster> univocal_clusters;
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";
711 return univocal_clusters;
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)
bool is_compatible(const std::vector< WSeries< double > * > &data)
std::vector< pixel > cluster
int num_of_time_bins(WSeries< double > *WS)
double median(std::vector< double > &vec)
std::vector< cluster > select_clusters_distinct_in_time(std::vector< cluster > &sorted_clusters, const double precision)
void heaviest_path(std::vector< double > &total_weight, std::vector< int > &total_length, std::vector< int > &best_ancestor, const std::vector< double > weights)
float get_map00(WSeries< double > *WS, int index)
void init(const wavenode node, const int offset, const int refscaleidx, const double path_weight)
wavenode get_node(int id)
bool is_topologically_sorted()
int num_of_freq_bins(WSeries< double > *WS)
bool bothspaces(char lhs, char rhs)
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)
nodeids get_ancestors(int id)
void create(const std::string &srcfile)
int get_scale(WSeries< double > *WS)
bool compare_paths(const wavepath &path1, const wavepath &path2)
void add_node(const wavenode node)
graph_header_code header_hash(std::string const &str)