Changeset 12727
- Timestamp:
- 12/20/11 14:23:48 (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/FACT++/src/fitsdump.cc
r12726 r12727 59 59 fits* fFile; 60 60 bool fDotsPlease; 61 bool fNoZeroPlease; 61 62 string fFilename; 62 63 … … 99 100 int doCurvesDisplay( const vector<string> &list, const string& tableName); 100 101 #endif 102 int doMinMaxPlease(const string& filename, const vector<string>& list, int precision); 101 103 int doStatsPlease(const string &filename, const vector<string>& list, int precision); 102 104 // bool Plot(const vector<string> &list); … … 113 115 //! the ostream where to redirect the outputs 114 116 // 115 FitsDumper::FitsDumper() : fFile(0), fDotsPlease(false) 117 FitsDumper::FitsDumper() : fFile(0), fDotsPlease(false), fNoZeroPlease(false) 116 118 { 117 119 } … … 453 455 return -1; 454 456 } 455 457 #ifdef PLOTTING_PLEASE 456 458 if (conf.Get<bool>("stat") && conf.Get<bool>("graph")) 457 459 { … … 459 461 return -1; 460 462 } 461 463 #endif 464 if (conf.Get<bool>("minmax") && conf.Get<bool>("stat")) 465 { 466 cout << "Invalid combination of options: cannot do stats and minmax. Aborting" << endl; 467 return -1; 468 } 469 #ifdef PLOTTING_PLEASE 470 if (conf.Get<bool>("minmax") && conf.Get<bool>("graph")) 471 { 472 cout << "Invalid combination of options: cannot graph minmax. Aborting" << endl; 473 return -1; 474 } 475 #endif 476 if (conf.Get<bool>("stat") && conf.Get<bool>("nozero")) 477 { 478 cout << "Invalid combination of options: nozero only works with minmax. Aborting" << endl; 479 return -1; 480 } 481 if (conf.Get<bool>("nozero")) 482 { 483 fNoZeroPlease = true; 484 } 462 485 if (conf.Get<bool>("list")) 463 486 List(); … … 467 490 if (!OpenTable(conf.Get<string>("tablename"))) 468 491 return -1; 492 } 493 494 if (conf.Get<bool>("minmax")) 495 { 496 if (!conf.Has("col")) 497 { 498 cout << "Please specify the columns that should be dumped as arguments. Aborting" << endl; 499 return 0; 500 } 501 doMinMaxPlease(conf.Get<string>("outfile"), conf.Get<vector<string>>("col"), conf.Get<int>("precision")); 502 return 0; 469 503 } 470 504 … … 532 566 { 533 567 // 568 } 569 570 struct minMaxStruct { 571 double min; 572 double max; 573 double average; 574 long numValues; 575 minMaxStruct() : min(1e10), max(-1e10), average(0), numValues(0) {} 576 }; 577 578 int FitsDumper::doMinMaxPlease(const string& filename, const vector<string>& list, int precision) 579 { 580 581 //first of all, let's separate the columns from their ranges and check that the requested columns are indeed part of the file 582 vector<pair<int, int> > ranges; 583 vector<string> listNamesOnly; 584 585 if (!separateColumnsFromRanges(list, ranges, listNamesOnly)) 586 { 587 cerr << "Something went wrong while extracting the columns names from parameters. Aborting" << endl; 588 return false; 589 } 590 591 ofstream out(filename=="-"?"/dev/stdout":filename); 592 if (!out) 593 { 594 cerr << "Cannot open file " << filename << ": " << strerror(errno) << endl; 595 return false; 596 } 597 598 out.precision(precision); 599 600 // Loop over all columns in our list of requested columns 601 vector<pair<char, char*> > columnsData; 602 vector<minMaxStruct> statData; 603 int numRows = fFile->GetInt("NAXIS2"); 604 auto rangesIt = ranges.begin(); 605 for (vector<string>::const_iterator it=listNamesOnly.begin(); it!=listNamesOnly.end(); it++, rangesIt++) 606 { 607 fits::Table::Column& cCol = fColMap.find(*it)->second; 608 columnsData.push_back(make_pair(cCol.type, new char[cCol.num*cCol.size])); 609 // minMaxStuct initData; 610 statData.push_back(minMaxStruct()); 611 fFile->SetPtrAddress(*it, columnsData[columnsData.size()-1].second); 612 } 613 614 int row = 0; 615 while (fFile->GetNextRow() && row < numRows) 616 { 617 rangesIt = ranges.begin(); 618 auto statsIt = statData.begin(); 619 for (auto it=columnsData.begin(); it != columnsData.end(); it++, rangesIt++, statsIt++) 620 { 621 double cValue = 0; 622 for (int i=rangesIt->first; i<rangesIt->second; i++) 623 { 624 switch (it->first) { 625 case 'L': 626 cValue = reinterpret_cast<bool*>(it->second)[i]; 627 break; 628 case 'B': 629 cValue = reinterpret_cast<bool*>(it->second)[i]; 630 break; 631 case 'I': 632 cValue = reinterpret_cast<int16_t*>(it->second)[i]; 633 break; 634 case 'J': 635 cValue = reinterpret_cast<int32_t*>(it->second)[i]; 636 break; 637 case 'K': 638 cValue = reinterpret_cast<int64_t*>(it->second)[i]; 639 break; 640 case 'E': 641 cValue = reinterpret_cast<float*>(it->second)[i]; 642 break; 643 case 'D': 644 cValue = reinterpret_cast<double*>(it->second)[i]; 645 break; 646 default: 647 ; 648 } 649 if (!fNoZeroPlease || cValue != 0) 650 { 651 statsIt->average += cValue; 652 if (cValue < statsIt->min) 653 statsIt->min = cValue; 654 if (cValue > statsIt->max) 655 statsIt->max = cValue; 656 statsIt->numValues++; 657 } 658 } 659 } 660 row++; 661 } 662 for (auto it = columnsData.begin(); it != columnsData.end(); it++) 663 delete[] it->second; 664 665 //okay. So now I've got ALL the data, loaded. 666 //let's do the summing and averaging in a safe way (i.e. avoid overflow of variables as much as possible) 667 rangesIt = ranges.begin(); 668 auto statsIt = statData.begin(); 669 670 auto nameIt = listNamesOnly.begin(); 671 for (auto it=columnsData.begin(); it != columnsData.end(); it++, rangesIt++, statsIt++, nameIt++) 672 { 673 int span = rangesIt->second - rangesIt->first; 674 cout << *nameIt << ": " << endl; 675 if (statsIt->numValues != 0) 676 { 677 statsIt->average /= statsIt->numValues; 678 out << "min: " << statsIt->min << endl; 679 out << "max: " << statsIt->max << endl; 680 out << "mea: " << statsIt->average << endl; 681 } 682 else 683 { 684 out << "min: 0" << endl << "max: 0" << endl << "mea: " << endl; 685 } 686 687 } 688 return true; 534 689 } 535 690 … … 549 704 out << "Max: " << reinterpret_cast<T*>(array)[numElems-1] << endl; 550 705 if (numElems%2 == 0) 551 out << "Med: " << (reinterpret_cast<T*>(array)[numElems/2] + reinterpret_cast< T*>(array)[numElems/2+1])/2.f << endl;706 out << "Med: " << (reinterpret_cast<T*>(array)[numElems/2] + reinterpret_cast<int16_t*>(array)[numElems/2+1])/2.f << endl; 552 707 else 553 708 out << "Med: " << reinterpret_cast<T*>(array)[numElems/2+1] << endl; … … 957 1112 ("header,h", po_switch(), "Dump header of given table") 958 1113 ("stat,s", po_switch(), "Perform statistics instead of dump") 1114 ("minmax,m", po_switch(), "Calculates min and max of data") 1115 ("nozero,z", po_switch(), "skip 0 values for stats") 959 1116 #ifdef PLOTTING_PLEASE 960 1117 ("graph,g", po_switch(), "Plot the columns instead of dumping them")
Note:
See TracChangeset
for help on using the changeset viewer.