Changeset 12850


Ignore:
Timestamp:
02/06/12 09:32:13 (13 years ago)
Author:
tbretz
Message:
Simplified all the loops looping over the columns to be displayed; do not handle the buffer memory but let fits handle it; unified the statistics output; added rms to statistics output; implemented the possibility to output ascii-columns ('A')
File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/FACT++/src/fitsdump.cc

    r12766 r12850  
    136136    {
    137137        case 'L': return "bool(8)";
     138        case 'A': return "char(8)";
    138139        case 'B': return "byte(8)";
    139140        case 'I': return "short(16)";
    … …  
    143144        case 'D': return "double(64)";
    144145    default:
    145         return "unknwown";
     146        return "unknown";
    146147    }
    147148}
    … …  
    151152    switch (type)
    152153    {
     154        case 'A': return sizeof(int8_t);
    153155        case 'L': return sizeof(uint8_t);
    154156        case 'B': return sizeof(int8_t);
    … …  
    396398    out << "#" << endl;
    397399
    398 
    399     // Loop over all columns in our list of requested columns
    400     vector<pair<char, char*> > columnsData;
    401 
    402     for (vector<string>::const_iterator it=listNamesOnly.begin(); it!=listNamesOnly.end(); it++)
    403     {
    404         fits::Table::Column& cCol = fColMap.find(*it)->second;
    405         columnsData.push_back(make_pair(cCol.type, new char[cCol.num*cCol.size]));
    406         fFile->SetPtrAddress(*it, columnsData[columnsData.size()-1].second);
    407     }
    408     int numRows = fFile->GetInt("NAXIS2");
    409     int row = 0;
    410     while (fFile->GetNextRow() && row < numRows)
    411     {
    412         row++;
     400    while (fFile->GetNextRow())
     401    {
     402        const size_t row = fFile->GetRow();
     403        if (row==fFile->GetNumRows())
     404            break;
     405
    413406        rangesIt = ranges.begin();
    414         for (auto it=columnsData.begin(); it != columnsData.end(); it++, rangesIt++)
     407
     408        for (vector<string>::const_iterator in=listNamesOnly.begin(); in!=listNamesOnly.end(); in++, rangesIt++)
    415409        {
     410            const void *ptr = fFile->SetPtrAddress(*in);
     411
     412            const fits::Table::Column &cCol = fColMap.find(*in)->second;
     413
     414            string msg;
    416415            for (int i=rangesIt->first; i<rangesIt->second; i++)
    417416            {
    418             switch (it->first) {
     417                switch (cCol.type)
     418                {
     419                case 'A':
     420                    msg += reinterpret_cast<const char*>(ptr)[i];
     421                    break;
     422                case 'B':
     423                    out << (unsigned int)reinterpret_cast<const unsigned char*>(ptr)[i] << " ";
     424                    break;
    419425                case 'L':
    420                         out << reinterpret_cast<bool*>(it->second)[i] << " ";
    421                         break;
    422                 case 'B':
    423                         out << reinterpret_cast<char*>(it->second)[i] << " ";
    424                         break;
     426                    out << reinterpret_cast<const bool*>(ptr)[i] << " ";
     427                    break;
    425428                case 'I':
    426                         out << reinterpret_cast<int16_t*>(it->second)[i] << " ";
    427                         break;
     429                    out << reinterpret_cast<const int16_t*>(ptr)[i] << " ";
     430                    break;
    428431                case 'J':
    429                         out << reinterpret_cast<int32_t*>(it->second)[i] << " ";
    430                         break;
     432                    out << reinterpret_cast<const int32_t*>(ptr)[i] << " ";
     433                    break;
    431434                case 'K':
    432                         out << reinterpret_cast<int64_t*>(it->second)[i] << " ";
    433                         break;
     435                    out << reinterpret_cast<const int64_t*>(ptr)[i] << " ";
     436                    break;
    434437                case 'E':
    435                         out << reinterpret_cast<float*>(it->second)[i] << " ";
    436                         break;
     438                    out << reinterpret_cast<const float*>(ptr)[i] << " ";
     439                    break;
    437440                case 'D':
    438                         out << reinterpret_cast<double*>(it->second)[i] << " ";
    439                         break;
     441                    out << reinterpret_cast<const double*>(ptr)[i] << " ";
     442                    break;
    440443                default:
    441444                    ;
     445                }
    442446            }
    443             }
     447
     448            if (cCol.type=='A')
     449                out << "'" << msg << "' ";
    444450        }
    445451        out << endl;
    446452    }
    447     for (auto it = columnsData.begin(); it != columnsData.end(); it++)
    448         delete[] it->second;
    449453    return true;
    450454}
    … …  
    591595    double max;
    592596    double average;
     597    double squared;
    593598    long numValues;
    594599    minMaxStruct() : min(1e10), max(-1e10), average(0), numValues(0) {}
    … …  
    599604    //first of all, let's separate the columns from their ranges and check that the requested columns are indeed part of the file
    600605    vector<pair<int, int> > ranges;
    601     vector<string> listNamesOnly;
    602 
    603     if (!separateColumnsFromRanges(list, ranges, listNamesOnly))
     606    vector<string> names;
     607
     608    if (!separateColumnsFromRanges(list, ranges, names))
    604609    {
    605610        cerr << "Something went wrong while extracting the columns names from parameters. Aborting" << endl;
    … …  
    616621    out.precision(precision);
    617622
     623    vector<minMaxStruct> statData(names.size());
     624
    618625    // Loop over all columns in our list of requested columns
     626    while (fFile->GetNextRow())
     627    {
     628        const size_t row = fFile->GetRow();
     629        if (row==fFile->GetNumRows())
     630            break;
     631
     632        auto rangesIt = ranges.begin();
     633        auto statsIt  = statData.begin();
     634
     635        for (auto in=names.begin(); in!=names.end(); in++, rangesIt++, statsIt++)
     636        {
     637            const void *ptr = fFile->SetPtrAddress(*in);
     638
     639            const fits::Table::Column &cCol = fColMap.find(*in)->second;
     640
     641            if (*in=="UnixTimeUTC" || *in=="PCTime")
     642            {
     643                const uint32_t *val = reinterpret_cast<const uint32_t*>(ptr);
     644                if (fNoZeroPlease && val[0]==0 && val[1]==0)
     645                    continue;
     646
     647                const double cValue = Time(val[0], val[1]).Mjd();
     648
     649                statsIt->average += cValue;
     650                statsIt->squared += cValue*cValue;
     651
     652                if (cValue < statsIt->min)
     653                    statsIt->min = cValue;
     654
     655                if (cValue > statsIt->max)
     656                    statsIt->max = cValue;
     657
     658                statsIt->numValues++;
     659
     660                continue;
     661            }
     662
     663            for (int i=rangesIt->first; i<rangesIt->second; i++)
     664            {
     665                double cValue = 0;
     666                switch (cCol.type)
     667                {
     668                case 'L':
     669                        cValue = reinterpret_cast<const bool*>(ptr)[i];
     670                        break;
     671                case 'B':
     672                        cValue = reinterpret_cast<const int8_t*>(ptr)[i];
     673                        break;
     674                case 'I':
     675                        cValue = reinterpret_cast<const int16_t*>(ptr)[i];
     676                        break;
     677                case 'J':
     678                        cValue = reinterpret_cast<const int32_t*>(ptr)[i];
     679                        break;
     680                case 'K':
     681                        cValue = reinterpret_cast<const int64_t*>(ptr)[i];
     682                        break;
     683                case 'E':
     684                        cValue = reinterpret_cast<const float*>(ptr)[i];
     685                        break;
     686                case 'D':
     687                        cValue = reinterpret_cast<const double*>(ptr)[i];
     688                        break;
     689                default:
     690                    ;
     691                }
     692
     693                if (fNoZeroPlease && cValue == 0)
     694                    continue;
     695
     696                statsIt->average += cValue;
     697                statsIt->squared += cValue*cValue;
     698
     699                if (cValue < statsIt->min)
     700                    statsIt->min = cValue;
     701
     702                if (cValue > statsIt->max)
     703                    statsIt->max = cValue;
     704
     705                statsIt->numValues++;
     706            }
     707        }
     708    }
     709
     710    // okay. So now I've got ALL the data, loaded.
     711    // let's do the summing and averaging in a safe way (i.e. avoid overflow
     712    // of variables as much as possible)
     713    auto statsIt = statData.begin();
     714
     715    for (auto it=names.begin(); it!=names.end(); it++, statsIt++)
     716    {
     717        cout << "[" << *it << "]" << endl;
     718        if (statsIt->numValues == 0)
     719        {
     720            out << "Min: -\nMax: -\nAvg: -\nRms: -" << endl;
     721            continue;
     722        }
     723
     724        double &avg = statsIt->average;
     725        double &rms = statsIt->squared;
     726        long   &num = statsIt->numValues;
     727
     728        avg /= num;
     729        rms  = sqrt(rms/num - avg*avg);
     730
     731        out << "Min: " << statsIt->min << '\n';
     732        out << "Max: " << statsIt->max << '\n';
     733        out << "Avg: " << avg << '\n';
     734        out << "Rms: " << rms << endl;
     735    }
     736
     737    /*
    619738    vector<pair<char, char*> > columnsData;
    620739    vector<minMaxStruct> statData;
    … …  
    644763            switch (it->first) {
    645764                case 'L':
    646                         cValue = reinterpret_cast<bool*>(it->second)[i];
     765                        cValue = reinterpret_cast<unsigned char*>(it->second)[i];
    647766                        break;
    648767                case 'B':
    … …  
    709828    for (auto it = columnsData.begin(); it != columnsData.end(); it++)
    710829        delete[] it->second;
    711 
    712     //okay. So now I've got ALL the data, loaded.
    713     //let's do the summing and averaging in a safe way (i.e. avoid overflow of variables as much as possible)
    714     rangesIt = ranges.begin();
    715     auto statsIt = statData.begin();
    716 
    717     auto nameIt = listNamesOnly.begin();
    718     for (auto it=columnsData.begin(); it != columnsData.end(); it++, rangesIt++, statsIt++, nameIt++)
    719     {
    720 //        int span = rangesIt->second - rangesIt->first;
    721         cout << *nameIt << ": " << endl;
    722         if (statsIt->numValues != 0)
    723         {
    724             statsIt->average /= statsIt->numValues;
    725             out << "min: " << statsIt->min << endl;
    726             out << "max: " << statsIt->max << endl;
    727             out << "mea: " << statsIt->average << endl;
    728         }
    729         else
    730         {
    731             out << "min: 0" << endl << "max: 0" << endl << "mea: " << endl;
    732         }
    733 
    734     }
     830*/
    735831    return true;
    736832}
    … …  
    854950void displayStats(char* array, int numElems, ofstream& out)
    855951{
    856         if (numElems == 0)
    857         {
    858                 out << "Min: 0" << endl;
    859                 out << "Max: 0" << endl;
    860                 out << "Med: 0" << endl;
    861                 out << "Mea: 0" << endl;
    862                 return;
    863         }
    864     sort(reinterpret_cast<T*>(array), &reinterpret_cast<T*>(array)[numElems]);
    865     out << "Min: " << reinterpret_cast<T*>(array)[0] << endl;
    866     out << "Max: " << reinterpret_cast<T*>(array)[numElems-1] << endl;
    867     if (numElems%2 == 0)
    868         out << "Med: " << (reinterpret_cast<T*>(array)[numElems/2] + reinterpret_cast<int16_t*>(array)[numElems/2+1])/2.f << endl;
    869     else
    870         out << "Med: " << reinterpret_cast<T*>(array)[numElems/2+1] << endl;
    871 
    872     double average = 0;
     952    if (numElems == 0)
     953    {
     954        out << "Min: -\nMax: -\nMed: -\nAvg: -\nRms: -" << endl;
     955        return;
     956    }
     957
     958    T *val = reinterpret_cast<T*>(array);
     959
     960    sort(val, val+numElems);
     961
     962    out << "Min: " << double(val[0]) << '\n';
     963    out << "Max: " << double(val[numElems-1]) << '\n';
     964
     965    if (numElems>2)
     966    {
     967        if (numElems%2 == 0)
     968            out << "Med: " << double(val[numElems/2] + val[numElems/2+1])/2 << '\n';
     969        else
     970            out << "Med: " << double(val[numElems/2+1]) << '\n';
     971    }
     972
     973    double avg = 0;
     974    double rms = 0;
    873975    for (int i=0;i<numElems;i++)
    874         average += reinterpret_cast<T*>(array)[i];
    875     out << "Mea: " << average/(double)numElems << endl;
     976    {
     977        avg += val[i];
     978        rms += val[i]*val[i];
     979    }
     980
     981    avg /= numElems;
     982    rms  = sqrt(rms/numElems - avg*avg);
     983
     984    out << "Avg: " << avg << '\n';
     985    out << "Rms: " << rms << endl;
    876986
    877987}
    … …  
    881991    //first of all, let's separate the columns from their ranges and check that the requested columns are indeed part of the file
    882992    vector<pair<int, int> > ranges;
    883     vector<string> listNamesOnly;
    884 
    885     if (!separateColumnsFromRanges(list, ranges, listNamesOnly))
     993    vector<string> names;
     994
     995    if (!separateColumnsFromRanges(list, ranges, names))
    886996    {
    887997        cerr << "Something went wrong while extracting the columns names from parameters. Aborting" << endl;
    … …  
    8991009
    9001010    // Loop over all columns in our list of requested columns
    901     vector<pair<char, char*> > columnsData;
    9021011    vector<char*> statData;
    903     int numRows = fFile->GetInt("NAXIS2");
     1012
    9041013    auto rangesIt = ranges.begin();
    905     for (vector<string>::const_iterator it=listNamesOnly.begin(); it!=listNamesOnly.end(); it++, rangesIt++)
    906     {
    907         fits::Table::Column& cCol = fColMap.find(*it)->second;
    908         columnsData.push_back(make_pair(cCol.type, new char[cCol.num*cCol.size]));
    909         statData.push_back(new char[(rangesIt->second - rangesIt->first)*cCol.size*numRows]);
    910         fFile->SetPtrAddress(*it, columnsData[columnsData.size()-1].second);
    911     }
    912 
    913     int row = 0;
    914     while (fFile->GetNextRow() && row < numRows)
    915     {
     1014    for (auto in=names.begin(); in!=names.end(); in++, rangesIt++)
     1015    {
     1016        cout << *in << endl;
     1017        const fits::Table::Column &cCol = fColMap.find(*in)->second;
     1018        statData.push_back(new char[(rangesIt->second - rangesIt->first)*cCol.size*fFile->GetNumRows()]);
     1019    }
     1020
     1021    while (fFile->GetNextRow())
     1022    {
     1023        const size_t row = fFile->GetRow();
     1024        if (row==fFile->GetNumRows())
     1025            break;
     1026
    9161027        rangesIt = ranges.begin();
    9171028        auto statsIt = statData.begin();
    918         for (auto it=columnsData.begin(); it != columnsData.end(); it++, rangesIt++, statsIt++)
     1029
     1030        for (auto in=names.begin(); in!=names.end(); in++, rangesIt++, statsIt++)
    9191031        {
    920             int span = rangesIt->second - rangesIt->first;
     1032            const void *ptr = fFile->SetPtrAddress(*in);
     1033
     1034            const fits::Table::Column &cCol = fColMap.find(*in)->second;
     1035
     1036            const int span = rangesIt->second - rangesIt->first;
     1037
    9211038            for (int i=rangesIt->first; i<rangesIt->second; i++)
    9221039            {
    923             switch (it->first) {
     1040                switch (cCol.type)
     1041                {
    9241042                case 'L':
    925                         reinterpret_cast<bool*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<bool*>(it->second)[i];
     1043                        reinterpret_cast<bool*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<const bool*>(ptr)[i];
    9261044                        break;
    9271045                case 'B':
    928                         reinterpret_cast<char*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<char*>(it->second)[i];
     1046                        reinterpret_cast<char*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<const char*>(ptr)[i];
    9291047                        break;
    9301048                case 'I':
    931                         reinterpret_cast<int16_t*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<int16_t*>(it->second)[i];
     1049                        reinterpret_cast<int16_t*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<const int16_t*>(ptr)[i];
    9321050                        break;
    9331051                case 'J':
    934                         reinterpret_cast<int32_t*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<int32_t*>(it->second)[i];
     1052                        reinterpret_cast<int32_t*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<const int32_t*>(ptr)[i];
    9351053                        break;
    9361054                case 'K':
    937                         reinterpret_cast<int64_t*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<int64_t*>(it->second)[i];
     1055                        reinterpret_cast<int64_t*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<const int64_t*>(ptr)[i];
    9381056                        break;
    9391057                case 'E':
    940                         reinterpret_cast<float*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<float*>(it->second)[i];
     1058                        reinterpret_cast<float*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<const float*>(ptr)[i];
    9411059                        break;
    9421060                case 'D':
    943                         reinterpret_cast<double*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<double*>(it->second)[i];
     1061                        reinterpret_cast<double*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<const double*>(ptr)[i];
    9441062                        break;
    9451063                default:
    9461064                    ;
     1065                }
    9471066            }
    948             }
    949         }
    950         row++;
    951     }
    952     for (auto it = columnsData.begin(); it != columnsData.end(); it++)
    953         delete[] it->second;
     1067        }
     1068    }
    9541069
    9551070    //okay. So now I've got ALL the data, loaded.
    … …  
    9581073    auto statsIt = statData.begin();
    9591074
    960     auto nameIt = list.begin();
    961     for (auto it=columnsData.begin(); it != columnsData.end(); it++, rangesIt++, statsIt++, nameIt++)
    962     {
    963         int span = rangesIt->second - rangesIt->first;
    964         out << *nameIt << ": " << endl;
    965         switch (it->first) {
     1075    for (auto it=list.begin(); it!=list.end(); it++, rangesIt++, statsIt++)
     1076    {
     1077        const fits::Table::Column &cCol = fColMap.find(*it)->second;
     1078
     1079        const int span = rangesIt->second - rangesIt->first;
     1080
     1081        out << "[" << *it << "]" << endl;
     1082        switch (cCol.type)
     1083        {
    9661084            case 'L':
    967                     displayStats<bool>(*statsIt, numRows*span, out);
     1085                    displayStats<bool>(*statsIt, fFile->GetNumRows()*span, out);
    9681086                    break;
    9691087            case 'B':
    970                     displayStats<char>(*statsIt, numRows*span, out);
     1088                    displayStats<char>(*statsIt, fFile->GetNumRows()*span, out);
    9711089                    break;
    9721090            case 'I':
    973                     displayStats<int16_t>(*statsIt, numRows*span, out);
     1091                    displayStats<int16_t>(*statsIt, fFile->GetNumRows()*span, out);
    9741092                    break;
    9751093            case 'J':
    976                     displayStats<int32_t>(*statsIt, numRows*span, out);
     1094                    displayStats<int32_t>(*statsIt, fFile->GetNumRows()*span, out);
    9771095                    break;
    9781096            case 'K':
    979                     displayStats<int64_t>(*statsIt, numRows*span, out);
     1097                    displayStats<int64_t>(*statsIt, fFile->GetNumRows()*span, out);
    9801098                    break;
    9811099            case 'E':
    982                     displayStats<float>(*statsIt, numRows*span, out);
     1100                    displayStats<float>(*statsIt, fFile->GetNumRows()*span, out);
    9831101                    break;
    9841102            case 'D':
    985                     displayStats<double>(*statsIt, numRows*span, out);
     1103                    displayStats<double>(*statsIt, fFile->GetNumRows()*span, out);
    9861104                    break;
    9871105            default:
    … …  
    9891107        }
    9901108    }
     1109
    9911110    return true;
    9921111}
Note: See TracChangeset for help on using the changeset viewer.