Changeset 12727


Ignore:
Timestamp:
12/20/11 14:23:48 (13 years ago)
Author:
lyard
Message:
added minmax and nozero
File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/FACT++/src/fitsdump.cc

    r12726 r12727  
    5959    fits* fFile;
    6060    bool fDotsPlease;
     61    bool fNoZeroPlease;
    6162    string fFilename;
    6263
     
    99100    int doCurvesDisplay( const vector<string> &list, const string& tableName);
    100101#endif
     102    int doMinMaxPlease(const string& filename, const vector<string>& list, int precision);
    101103    int doStatsPlease(const string &filename, const vector<string>& list, int precision);
    102104    //    bool Plot(const vector<string> &list);
     
    113115//!        the ostream where to redirect the outputs
    114116//
    115 FitsDumper::FitsDumper() : fFile(0), fDotsPlease(false)
     117FitsDumper::FitsDumper() : fFile(0), fDotsPlease(false), fNoZeroPlease(false)
    116118{
    117119}
     
    453455            return -1;
    454456    }
    455 
     457#ifdef PLOTTING_PLEASE
    456458    if (conf.Get<bool>("stat") && conf.Get<bool>("graph"))
    457459    {
     
    459461        return -1;
    460462    }
    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    }
    462485    if (conf.Get<bool>("list"))
    463486        List();
     
    467490        if (!OpenTable(conf.Get<string>("tablename")))
    468491            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;
    469503    }
    470504
     
    532566{
    533567    //
     568}
     569
     570struct 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
     578int 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;
    534689}
    535690
     
    549704    out << "Max: " << reinterpret_cast<T*>(array)[numElems-1] << endl;
    550705    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;
    552707    else
    553708        out << "Med: " << reinterpret_cast<T*>(array)[numElems/2+1] << endl;
     
    9571112        ("header,h",    po_switch(),                "Dump header of given table")
    9581113        ("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")
    9591116#ifdef PLOTTING_PLEASE
    9601117        ("graph,g",     po_switch(),                "Plot the columns instead of dumping them")
Note: See TracChangeset for help on using the changeset viewer.