int main	(	int	argc,
		const char **	argv
	)

Definition at line 1505 of file fitsCompressor.cc.

References Checksum::add(), CompressedFitsWriter::addColumn(), ofits::AddColumn(), BlockHeader(), buffer, izstream::close(), CompressedFitsWriter::close(), ofits::close(), count, Configuration::DoParse(), first, Configuration::Get(), ofits::GetChecksumData(), fits::GetColumns(), fits::GetKeys(), fits::GetNextRow(), zfits::GetNumRows(), fits::GetSortedColumns(), fits::GetStr(), fits::GetUInt(), Configuration::Has(), fits::HasKey(), i, fits::IsCompressedFITS(), zfits::IsFileOk(), FITS::kFactHuffman16, FITS::kFactRaw, FITS::kFactSmoothing, FITS::kOrderByRow, ofits::open(), CompressedFitsWriter::open(), printHelp(), printUsage(), second, CompressedFitsWriter::setDrsCalib(), CompressedFitsWriter::setHeaderKey(), ofits::SetKeyComment(), CompressedFitsWriter::setNumWorkingThreads(), Configuration::SetPrintUsage(), fits::SetPtrAddress(), setupConfiguration(), str, type, Checksum::val(), Configuration::Vec(), CompressedFitsWriter::writeBinaryRow(), ofits::WriteRow(), and ofits::WriteTableHeader().

{
     Configuration conf(argv[0]);
     conf.SetPrintUsage(printUsage);
     setupConfiguration(conf);

     if (!conf.DoParse(argc, argv, printHelp))
         return -1;

     //initialize the file names to nothing.
     string fileNameIn = "";
     string fileNameOut = "";
     string drsFileName = "";
     uint32_t numRowsPerTile = 100;
     bool displayText=true;

    //parse configuration
    if (conf.Get<bool>("quiet")) displayText = false;
    const vector<string> inputFileNameVec = conf.Vec<string>("inputFile");
    if (inputFileNameVec.size() != 1)
    {
       cout << "Error: ";
       if (inputFileNameVec.size() == 0) cout << "no";
       else cout << inputFileNameVec.size();
       cout << " input file(s) given. Expected one. Aborting. Input:" << endl;;
       for (unsigned int i=0;i<inputFileNameVec.size(); i++)
           cout << inputFileNameVec[i] << endl;
       return -1;
    }

    //Assign the input filename
    fileNameIn = inputFileNameVec[0];

    //Check if we have a drs calib too
    if (conf.Has("drs")) drsFileName = conf.Get<string>("drs");

    //Should we verify the data ?
    bool verifyDataPlease = false;
    if (conf.Has("verify")) verifyDataPlease = conf.Get<bool>("verify");


    //should we use a specific output filename ?
    if (conf.Has("output"))
        fileNameOut = conf.Get<string>("output");
    else
    {
        size_t pos = fileNameIn.find(".fits");
        if (pos == string::npos)
        {
            cout << "ERROR: input file does not seems ot be fits. Aborting." << endl;
            return -1;
        }
        fileNameOut = fileNameIn + ".fz";
    }


    //should we use specific compression on some columns ?
    const vector<string> columnsCompression = conf.Vec<string>("compression");

    //split up values between column names and compression scheme
    vector<std::pair<string, string>> compressions;
    for (unsigned int i=0;i<columnsCompression.size();i++)
    {
        size_t pos = columnsCompression[i].find_first_of("=");
        if (pos == string::npos)
        {
            cout << "ERROR: Something wrong occured while parsing " << columnsCompression[i] << ". Aborting." << endl;
            return -1;
        }
        string comp = columnsCompression[i].substr(pos+1);
        if (comp != "UNCOMPRESSED" && comp != "AMPLITUDE" && comp != "HUFFMAN" &&
            comp != "SMOOTHMAN" && comp != "INT_WAVELET")
        {
            cout << "Unkown compression scheme requested (" << comp << "). Aborting." << endl;
            return -1;
        }
        compressions.push_back(make_pair(columnsCompression[i].substr(0, pos), comp));
    }

    //How many rows per tile should we use ?
    if (conf.Has("rowPerTile")) numRowsPerTile = conf.Get<unsigned int>("rowPerTile");

    /************************************************************************************
     *  Done reading configuration. Open relevant files
     ************************************************************************************/

    //Open input's fits file
    factfits inFile(fileNameIn);

    if (inFile.IsCompressedFITS())
    {
        cout << "ERROR: input file is already a compressed fits. Cannot be compressed again: Aborting." << endl;
        return -1;
    }

    //decide how many tiles should be put in the compressed file
    uint32_t originalNumRows = inFile.GetNumRows();
    uint32_t numTiles = (originalNumRows%numRowsPerTile) ? (originalNumRows/numRowsPerTile)+1 : originalNumRows/numRowsPerTile;
    CompressedFitsWriter outFile(numTiles, numRowsPerTile);

    //should we use a specific number of threads for compressing ?
    unsigned int numThreads = 1;
    if (conf.Has("threads"))
    {
        numThreads = conf.Get<unsigned int>("threads");
        outFile.setNumWorkingThreads(numThreads);
    }

    if (!outFile.open(fileNameOut))
    {
        cout << "Error: could not open " << fileNameOut << " for writing" << endl;
        return -1;
    }

    //Because the file to open MUST be given by the constructor, I must use a pointer instead
    factfits* drsFile = NULL;
    //try to open the Drs file. If any.
    if (drsFileName != "")
    {
        try
        {
            drsFile = new factfits(drsFileName);
        }
        catch (...)
        {
            cout << "Error: could not open " << drsFileName << " for calibration" << endl;
            return -1;
        }
    }

    if (displayText)
    {
        cout << endl;
        cout << "**********************" << endl;
        cout << "Will compress from    : " << fileNameIn << endl;
        cout << "to                    : " << fileNameOut << endl;
        if (drsFileName != "")
            cout << "while calibrating with: " << drsFileName << endl;
        cout << "Compression will use  : " << numThreads << " worker threads" << endl;
        cout << "Data will be verified : ";
        if (verifyDataPlease)
            cout << "yes" << endl;
        else
            cout << "no (WARNING !)" << endl;
        cout << "**********************" << endl;
        cout << endl;
    }

    /************************************************************************************
     *  Done opening input files. Allocate memory and configure output file
     ************************************************************************************/

    //allocate the buffer for temporary storage of each read/written row
    uint32_t rowWidth = inFile.GetUInt("NAXIS1");
    char* buffer = new char[rowWidth + 12];
    memset(buffer, 0, 4);
    buffer = buffer+4;

    //get the source columns
    const fits::Table::Columns& columns = inFile.GetColumns();
    const fits::Table::SortedColumns& sortedColumns = inFile.GetSortedColumns();
    if (displayText)
        cout << "Input file has " << columns.size() << " columns and " << inFile.GetNumRows() << " rows" << endl;

    //Add columns.
    uint32_t totalRowWidth = 0;
    for (uint32_t i=0;i<sortedColumns.size(); i++)
    {
        //get column name
        ostringstream str;
        str << "TTYPE" << i+1;
        string colName = inFile.GetStr(str.str());
        if (displayText)
        {
            cout << "Column " << colName;
            for (uint32_t j=colName.size();j<21;j++)
                cout << " ";
            cout << " -> ";
        }

        //get header structures
        BlockHeader rawHeader;
        BlockHeader smoothmanHeader(0, FITS::kOrderByRow, 2);
        vector<uint16_t> rawProcessings(1);
        rawProcessings[0] = FITS::kFactRaw;
        vector<uint16_t> smoothmanProcessings(2);
        smoothmanProcessings[0] = FITS::kFactSmoothing;
        smoothmanProcessings[1] = FITS::kFactHuffman16;
//        smoothmanProcessings[2] = FACT_RAW;

        totalRowWidth += sortedColumns[i].bytes;

        //first lets see if we do have an explicit request
        bool explicitRequest = false;
        for (unsigned int j=0;j<compressions.size();j++)
        {
            if (compressions[j].first == colName)
            {
                explicitRequest = true;
                if (displayText) cout << compressions[j].second << endl;
                if (compressions[j].second == "UNCOMPRESSED")
                    outFile.addColumn(CompressedFitsFile::ColumnEntry(colName, sortedColumns[i].type, sortedColumns[i].num, rawHeader, rawProcessings));
                if (compressions[j].second == "SMOOTHMAN")
                    outFile.addColumn(CompressedFitsFile::ColumnEntry(colName, sortedColumns[i].type,  sortedColumns[i].num, smoothmanHeader, smoothmanProcessings));
                break;
            }
        }

        if (explicitRequest) continue;

        if (colName != "Data")
        {
            if (displayText) cout << "UNCOMPRESSED" << endl;
            outFile.addColumn(CompressedFitsFile::ColumnEntry(colName, sortedColumns[i].type,  sortedColumns[i].num, rawHeader, rawProcessings));
        }
        else
        {
            if (displayText) cout << "SMOOTHMAN" << endl;
            outFile.addColumn(CompressedFitsFile::ColumnEntry(colName, sortedColumns[i].type,  sortedColumns[i].num, smoothmanHeader, smoothmanProcessings));
        }
    }

    //translate original header entries to their Z-version
    const fits::Table::Keys& header = inFile.GetKeys();
    for (fits::Table::Keys::const_iterator i=header.begin(); i!= header.end(); i++)
    {
        string k = i->first;//header[i].key();
        if (k == "XTENSION" || k == "BITPIX" || k == "PCOUNT" || k == "GCOUNT" || k == "TFIELDS")
            continue;
        if (k == "CHECKSUM")
        {
            outFile.setHeaderKey(CompressedFitsFile::HeaderEntry("ZCHKSUM", i->second.value, i->second.comment));
            continue;
        }
        if (k == "DATASUM")
        {
            outFile.setHeaderKey(CompressedFitsFile::HeaderEntry("ZDTASUM", i->second.value, i->second.comment));
            continue;
        }
        k = k.substr(0,5);
        if (k == "TTYPE" || k == "TFORM")
        {
            string tmpKey = i->second.fitsString;
            tmpKey[0] = 'Z';
            outFile.setHeaderKey(tmpKey);
            continue;
        }
        if (k == "NAXIS")
            continue;
        outFile.setHeaderKey(i->second.fitsString);
//        cout << i->first << endl;
    }

    outFile.setHeaderKey(CompressedFitsFile::HeaderEntry("RAWSUM", "         0", "Checksum of raw littlen endian data"));

    //deal with the DRS calib
    int16_t* drsCalib16 = NULL;

    //load the drs calib. data
    int32_t startCellOffset = -1;
    if (drsFileName != "")
    {
        drsCalib16 = new int16_t[1440*1024];
        float* drsCalibFloat = NULL;
        try
        {
            drsCalibFloat = reinterpret_cast<float*>(drsFile->SetPtrAddress("BaselineMean"));
        }
        catch (...)
        {
            cout << "Could not find column BaselineMean in drs calibration file " << drsFileName << ". Aborting" << endl;
            return -1;
        }

        //read the calibration and calculate its integer value
        drsFile->GetNextRow();
        for (uint32_t i=0;i<1440*1024;i++)
            drsCalib16[i] = (int16_t)(drsCalibFloat[i]*4096.f/2000.f);


        //get the start cells offsets
        for (fits::Table::Columns::const_iterator it=columns.begin(); it!= columns.end(); it++)
            if (it->first == "StartCellData")
            {
                startCellOffset = it->second.offset;
                if (it->second.type != 'I')
                {
                    cout << "Wrong type for the StartCellData Column: " << it->second.type << " instead of I expected"<< endl;
                    return -1;
                }
            }

        if (startCellOffset == -1)
        {
            cout << "WARNING: Could not find StartCellData in input file " << fileNameIn << ". Doing it uncalibrated"<< endl;
        }
        else
        {
            //assign it to the ouput file
            outFile.setDrsCalib(drsCalib16);
        }
    }

    /************************************************************************************
     *  Done configuring compression. Do the real job now !
     ************************************************************************************/

    if (displayText) cout << "Converting file..." << endl;

    int numSlices = -1;
    int32_t dataOffset = -1;

    //Get the pointer to the column that must be drs-calibrated
    for (fits::Table::Columns::const_iterator it=columns.begin(); it!= columns.end(); it++)
        if (it->first == "Data")
        {
            numSlices = it->second.num;
            dataOffset = it->second.offset;
        }
    if (numSlices % 1440 != 0)
    {
        cout << "seems like the number of samples is not a multiple of 1440. Aborting." << endl;
        return -1;
    }
    if (dataOffset == -1)
    {
        cout << "Could not find the column Data in the input file. Aborting." << endl;
        return -1;
    }

    numSlices /= 1440;

    //set pointers to the readout data to later be able to gather it to "buffer".
    vector<void*>   readPointers;
    vector<int32_t> readOffsets;
    vector<int32_t> readElemSize;
    vector<int32_t> readNumElems;
    for (fits::Table::Columns::const_iterator it=columns.begin(); it!= columns.end(); it++)
    {
        readPointers.push_back(inFile.SetPtrAddress(it->first));
        readOffsets.push_back(it->second.offset);
        readElemSize.push_back(it->second.size);
        readNumElems.push_back(it->second.num);
    }

    Checksum rawsum;
    //Convert each row one after the other
    ostringstream wrongChannels;
    map<int, int> wrongChannelsMap;
    map<int, int> wrongChannelsLastValues;
    for (uint32_t i=0;i<1440;i++)
    {
        wrongChannelsMap[i] = 0;
        wrongChannelsLastValues[i] = 0;
    }
    for (uint32_t i=0;i<inFile.GetNumRows();i++)
    {
        if (displayText) cout << "\r Row " << i+1 << flush;
        if (!inFile.GetNextRow())
        {
            cout << "ERROR: file has less rows than advertized. aborting" << endl;
            exit(0);
        }
        //copy from inFile internal structures to buffer
        int32_t count=0;
        for (fits::Table::Columns::const_iterator it=columns.begin(); it!= columns.end();it++)
        {
            memcpy(&buffer[readOffsets[count]], readPointers[count], readElemSize[count]*readNumElems[count]);
            count++;
        }

        char* checkSumPointer = buffer;
        const int chkOffset = (i*totalRowWidth)%4;
        checkSumPointer -= chkOffset;
        uint32_t sizeToChecksum = totalRowWidth + chkOffset;
        if (sizeToChecksum%4 != 0)
        {
            int32_t extraBytes = 4 - (sizeToChecksum%4);
            memset(checkSumPointer+sizeToChecksum, 0, extraBytes);
            sizeToChecksum += extraBytes;
        }
        rawsum.add(checkSumPointer, sizeToChecksum, false);

        if (startCellOffset != -1)
        {//drs calibrate this data
            for (int j=0;j<1440;j++)
            {
                const int thisStartCell = reinterpret_cast<int16_t*>(&buffer[startCellOffset])[j];
                if (thisStartCell > 1023)
                {
                    wrongChannelsMap[j]++;
                    wrongChannelsLastValues[j] = thisStartCell;
                    wrongChannels << j;
                }
                if (thisStartCell < 0)  continue;
                for (int k=0;k<numSlices;k++)
                    reinterpret_cast<int16_t*>(&buffer[dataOffset])[numSlices*j + k] -= drsCalib16[1024*j + (thisStartCell+k)%1024];
            }
        }
        outFile.writeBinaryRow(buffer);
    };

    if (wrongChannels.str() != "")
    {
        cout << "ERROR: Wrong channels: ";
        for (uint32_t i=0;i<1440;i++)
        {
            if (wrongChannelsMap[i] != 0)
                cout << i << "(" << wrongChannelsMap[i] << "|" << wrongChannelsLastValues[i] << ") ";
        }
        cout << endl;
        exit(-1);
    }
    ostringstream strSum;
    strSum << rawsum.val();
    outFile.setHeaderKey(CompressedFitsFile::HeaderEntry("RAWSUM", strSum.str(), "Checksum of raw littlen endian data"));

    //Get table name for later use in case the compressed file is to be verified
    string tableName = inFile.GetStr("EXTNAME");

    if (displayText) cout << endl << "Done. Flushing output file..." << endl;

    inFile.close();
    if (!outFile.close())
    {
        cout << "Something went wrong while writing the catalog: negative index" << endl;
        return false;
    }

    if (drsCalib16 != NULL)
        delete[] drsCalib16;

    if (displayText) cout << "Done." << endl;

    /************************************************************************************
     *  Actual job done. Should we verify what we did ?
     ************************************************************************************/

    if (verifyDataPlease)
    {
        if (displayText) cout << "Now verify data..." << endl;
    }
    else
        return 0;

    //get a compressed reader
//TEMP try to copy the file too
//    string copyName("/gpfs/scratch/fact/etienne/copyFile.fz");
//    string copyName(fileNameOut+".copy");
    string copyName("");
    factfits verifyFile(fileNameOut, copyName, tableName, false);

    //and the header of the compressed file
    const fits::Table::Keys& header2 = verifyFile.GetKeys();

    //get a non-compressed writer
    ofits reconstructedFile;

    //figure out its name: /dev/null unless otherwise specified
    string reconstructedName = fileNameOut+".recons";
    reconstructedName = "/dev/null";
    reconstructedFile.open(reconstructedName.c_str(), false);

    //reconstruct the original columns from the compressed file.
    string origChecksumStr;
    string origDatasum;

    //reset tablename value so that it is re-read from compressed table's header
    tableName = "";

    /************************************************************************************
     *  Reconstruction setup done. Rebuild original header
     ************************************************************************************/

    //re-tranlate the keys
    for (fits::Table::Keys::const_iterator it=header2.begin(); it!= header2.end(); it++)
    {
        string k = it->first;
        if (k == "XTENSION" || k == "BITPIX"  || k == "PCOUNT"   || k == "GCOUNT" ||
            k == "TFIELDS"  || k == "ZTABLE"  || k == "ZNAXIS1"  || k == "ZNAXIS2" ||
            k == "ZHEAPPTR" || k == "ZPCOUNT" || k == "ZTILELEN" || k == "THEAP" ||
            k == "CHECKSUM" || k == "DATASUM" || k == "FCTCPVER" || k == "ZHEAP")
        {
            continue;
        }

        if (k == "ZCHKSUM")
        {
            reconstructedFile.SetKeyComment("CHECKSUM", it->second.comment);
            origChecksumStr = it->second.value;
            continue;
        }
        if (k == "RAWSUM")
        {
            continue;
        }

        if (k == "ZDTASUM")
        {
            reconstructedFile.SetKeyComment("DATASUM",  it->second.comment);
            origDatasum = it->second.value;
            continue;
        }

        if (k == "EXTNAME")
        {
            tableName = it->second.value;
        }

        k = k.substr(0,5);

        if (k == "TTYPE")
        {//we have an original column name here.
         //manually deal with these in order to preserve the ordering (easier than re-constructing yet another list on the fly)
            continue;
        }

        if (k == "TFORM" || k == "NAXIS" || k == "ZCTYP" )
        {
            continue;
        }

        if (k == "ZFORM" || k == "ZTYPE")
        {
            string tmpKey = it->second.fitsString;
            tmpKey[0] = 'T';
            reconstructedFile.SetKeyFromFitsString(tmpKey);
            continue;
        }

        reconstructedFile.SetKeyFromFitsString(it->second.fitsString);
    }

    if (tableName == "")
    {
        cout << "Error: table name from file " << fileNameOut << " could not be found. Aborting" << endl;
        return -1;
    }

    //Restore the original columns
    for (uint32_t numCol=1; numCol<10000; numCol++)
    {
        ostringstream str;
        str << numCol;
        if (!verifyFile.HasKey("TTYPE"+str.str())) break;

        string ttype    = verifyFile.GetStr("TTYPE"+str.str());
        string tform    = verifyFile.GetStr("ZFORM"+str.str());
        char   type     = tform[tform.size()-1];
        string number   = tform.substr(0, tform.size()-1);
        int    numElems = atoi(number.c_str());

        if (number == "") numElems=1;

        reconstructedFile.AddColumn(numElems, type, ttype, "", "", false);
    }

    reconstructedFile.WriteTableHeader(tableName.c_str());

    /************************************************************************************
     *  Original header restored. Do the data
     ************************************************************************************/

    //set pointers to the readout data to later be able to gather it to "buffer".
    readPointers.clear();
    readOffsets.clear();
    readElemSize.clear();
    readNumElems.clear();
    for (fits::Table::Columns::const_iterator it=columns.begin(); it!= columns.end(); it++)
    {
        readPointers.push_back(verifyFile.SetPtrAddress(it->first));
        readOffsets.push_back(it->second.offset);
        readElemSize.push_back(it->second.size);
        readNumElems.push_back(it->second.num);
    }

    //do the actual reconstruction work
    uint32_t i=1;
    while (i<=verifyFile.GetNumRows() && verifyFile.GetNextRow())
    {
        int count=0;
        for (fits::Table::Columns::const_iterator it=columns.begin(); it!= columns.end();it++)
        {
            memcpy(&buffer[readOffsets[count]], readPointers[count], readElemSize[count]*readNumElems[count]);
            count++;
        }
        if (displayText) cout << "\r Row " << i << flush;
        reconstructedFile.WriteRow(buffer, rowWidth);
        i++;
    }

    if (displayText) cout << endl;

    //close reconstruction input and output
//    Do NOT close the verify file, otherwise data cannot be flushed to copy file
//    verifyFile.close();
    if (!verifyFile.IsFileOk())
        cout << "ERROR: file checksums seems wrong" << endl;

    if (!reconstructedFile.close())
    {
        cout << "ERROR: disk probably full..." << endl;
        return -1;
    }

    //get original and reconstructed checksum and datasum
    std::pair<string, int> origChecksum = make_pair(origChecksumStr, atoi(origDatasum.c_str()));
    std::pair<string, int> newChecksum = reconstructedFile.GetChecksumData();

    //verify that no mistake was made
    if (origChecksum.second != newChecksum.second)
    {
        cout << "ERROR: datasums are NOT identical: " << (uint32_t)(origChecksum.second) << " vs " << (uint32_t)(newChecksum.second) << endl;
        return -1;
    }
    if (origChecksum.first != newChecksum.first)
    {
        cout << "WARNING: checksums are NOT Identical: " << origChecksum.first << " vs " << newChecksum.first << endl;
    }
    else
    {
        if (true) cout << "Ok" << endl;
    }

    buffer = buffer-4;
    delete[] buffer;
    return 0;
}

Here is the call graph for this function: