Merge branch 'add-odf-reader' into 'master'

add odf reader and hdf5 writer

See merge request rnsd/sammy!48

sammy/src/io/CMakeLists.txt

@@ -2,20 +2,25 @@ INCLUDE(SammyPackageSetup)
 
 TRIBITS_PACKAGE(io)
 
+find_package(HDF5 REQUIRED COMPONENTS C CXX HL)
 
 configure_file(${CMAKE_CURRENT_SOURCE_DIR}/../sammy_conf.h.in sammy_conf.h @ONLY)
 
-INCLUDE_DIRECTORIES(${CMAKE_CURRENT_BINARY_DIR})
-
+INCLUDE_DIRECTORIES(${CMAKE_CURRENT_BINARY_DIR}
+                    ${HDF5_INCLUDE_DIRS}) # from find_package()
 
 # Set headers
 SET(HEADERS ExperimentalParametersIO.h
+            OdfIO.h
+            Hdf5IO.h
             interface/cpp/ExperimentalParametersIOInterface.h 
 )
 
 APPEND_SET(IO_SOURCES
 
         ExperimentalParametersIO.cpp
+        OdfIO.cpp
+        Hdf5IO.cpp
         interface/cpp/ExperimentalParametersIOInterface.cpp
         interface/fortran/ExperimentalParametersIO_M.f90   
         interface/fortran/ExperimentalParametersIO_I.f90      
@@ -34,3 +39,15 @@ INSTALL(FILES ${HEADERS} DESTINATION "src/io")
 TRIBITS_ADD_TEST_DIRECTORIES(tests)
 
 TRIBITS_PACKAGE_POSTPROCESS()
+
+target_link_libraries(SammyIOUtilsLib ${HDF5_LIBRARIES} ) # from find_package()
+
+# --------------------------------------------
+# Build separate executable to convert ODF
+# to HDF5 files, install to normal location
+# --------------------------------------------
+
+add_executable(convertodf convertodf.cpp)
+target_link_libraries(convertodf SammyIOUtilsLib)
+install(PROGRAMS ${CMAKE_CURRENT_BINARY_DIR}/convertodf  DESTINATION "bin/")
+

sammy/src/io/Hdf5IO.cpp

+#include "Hdf5IO.h"
+
+namespace sammy {
+    
+    Hdf5IO::Hdf5IO() {}
+
+    Hdf5IO::Hdf5IO(std::map<std::string,int> &header){
+        // copy function var directly to class header member
+        fillHDF5header(header);
+    }
+
+    void Hdf5IO::writeODFtoHDF5(std::map<std::string,int> odfhead, 
+                                std::vector<std::vector<double>> &data, 
+                                std::string hdf5Filename){
+
+        // ----------------------------------------------------------
+        // Set up data sections
+        // ----------------------------------------------------------
+        if( data.size() != 11 ){
+            std::cerr << "---------------------------------------" << std::endl;
+            std::cerr << "Reading of ODF files with < 11 sections"
+                      << "has not been implemented yet. " << std::endl;
+            std::cerr << "---------------------------------------" << std::endl;
+            throw std::runtime_error("Unrecognized ODF file data.");
+        }
+            
+        // attempt to make std::strings
+
+        // HDF5 requires either char * or char[]
+        std::vector<const char *> sect11 = {"energy","trans","dtrans",
+                              "cro","dcro","cs","zero","co","zero","idk","dcs"};
+        std::vector<const char *> sect9;
+
+        // HDF5 accepts set-size double arrays for datasets: copy over
+        const size_t ns = odfhead["numSection"];
+        const size_t nc = odfhead["numChan"];
+        double dset_data[ns][nc];
+        for( size_t i=0;i<data.size();++i ){
+            std::copy(data[i].begin(),data[i].end(), dset_data[i]);
+        }
+
+        // ----------------------------------------------------------
+        // prepare 3 datasets: data, data section names, and header
+        // ----------------------------------------------------------
+
+        // --- data space -----
+        hsize_t ddim[2];
+        ddim[0] = odfhead["numSection"];
+        ddim[1] = odfhead["numChan"];
+        int drank = sizeof(ddim) / sizeof(hsize_t);
+        std::string dataname("/data/sections");
+
+        // --- section names space -----
+        hsize_t sdim[1];
+        sdim[0] = odfhead["numSection"];
+        int srank = sizeof(sdim) / sizeof(hsize_t);
+        std::string sectionnames("/data/sectnames");
+
+        // --- header -----
+        fillHDF5header(odfhead);
+
+        hsize_t mdim[1];
+        mdim[0] = h5head.size();
+        int mrank = sizeof(mdim) / sizeof(hsize_t);
+
+        const std::string headername("/meta/header");
+        const std::string headKeyword("key");
+        const std::string headVal("value");
+
+        // --- define compound type for header struct -----
+        H5::CompType htype(sizeof(sammy::headKeypair));
+        htype.insertMember(headKeyword, HOFFSET(sammy::headKeypair, key   ),H5::StrType(H5::PredType::C_S1, H5T_VARIABLE));
+        htype.insertMember(headVal,     HOFFSET(sammy::headKeypair, value ),H5::PredType::NATIVE_INT);
+
+        // ----------------------------------------------------------
+        // open file and write 
+        // ----------------------------------------------------------
+        H5::DataSpace dspace(drank,ddim);
+        H5::DataSpace sspace(srank,sdim);
+        H5::DataSpace mspace(mrank,mdim);
+        H5::H5File  file = H5::H5File(hdf5Filename,H5F_ACC_TRUNC);
+        H5::Group   group = H5::Group( file.createGroup( "/data" ));
+        H5::Group   metaGroup = H5::Group( file.createGroup( "/meta" ));
+
+        H5::DataSet dataset = H5::DataSet(file.createDataSet(dataname,H5::PredType::NATIVE_DOUBLE,dspace));
+        H5::DataSet nameset = H5::DataSet(file.createDataSet(sectionnames,H5::StrType(H5::PredType::C_S1, H5T_VARIABLE),sspace));
+        H5::DataSet metaset = H5::DataSet(file.createDataSet(headername,htype,mspace));
+
+        dataset.write( dset_data,     H5::PredType::NATIVE_DOUBLE );
+        nameset.write( sect11.data(), H5::StrType(H5::PredType::C_S1, H5T_VARIABLE) );
+        metaset.write( h5head.data(), htype );
+    }
+
+    void Hdf5IO::fillHDF5header(std::map<std::string,int> &odfhead){
+        
+        for( auto & member : odfhead ){
+            headKeypair hd;
+            hd.key = member.first.c_str();
+            hd.value = odfhead[member.first];
+            h5head.push_back(hd);
+        }
+    }
+
+}
\ No newline at end of file

sammy/src/io/Hdf5IO.h

+#ifndef SAMMY_HDF5IO_H
+#define SAMMY_HDF5IO_H
+
+#include <vector>
+#include <fstream>
+#include <iostream>
+#include <cstring>
+#include <stdexcept>
+#include "OdfIO.h"
+#include "H5Cpp.h"
+
+namespace sammy{
+
+    /*****************************************************************
+     * This struct defines the keypair values that will be used to 
+     * compose the header for HDF5 metadata
+     ****************************************************************/
+    struct headKeypair{
+        const char * key; 
+        int value;
+    };
+
+    /*****************************************************************
+     * The Hdf5IO class is designed to house any HDF5 input/output 
+     * routines inlcuding functions to convert ORELA Data Format (ODF)
+     * files into HDF5 files.
+     ****************************************************************/
+    class Hdf5IO {
+    
+    public:
+        Hdf5IO();
+        Hdf5IO(std::map<std::string,int> &header);
+        virtual ~Hdf5IO() {}
+
+        /*************************************************************
+         * Write data from ODF file into HDF5
+         * 
+         * @param head the header information from the ODF file
+         * @param data the sections of data read from the ODF file
+         * @param hdf5Filename name of HDF5 file to be written
+         ************************************************************/
+        void writeODFtoHDF5(std::map<std::string,int> odfhead, 
+                            std::vector<std::vector<double>> &data,
+                            std::string hdf5Filename);
+        
+        /*************************************************************
+         * Move data into an HDF5 friendly data-structure: array of 
+         * keyword and value structs
+         * 
+         * @param head the header information from the ODF file
+         * @param h5head the array to be filled and written to HDF5
+         ************************************************************/
+        void fillHDF5header(std::map<std::string,int> &odfhead);
+
+    private:
+        std::vector<headKeypair> h5head;
+    };
+}
+
+#endif /* SAMMY_HDF5IO_H */
\ No newline at end of file

sammy/src/io/OdfIO.cpp

+#include "OdfIO.h"
+
+namespace sammy{
+    
+    OdfIO::OdfIO() {}
+
+    double OdfIO::readFloat(std::istream &fs, bool unix){
+    
+        unsigned char tmp1,tmp2,tmp3,tmp4;
+        int i=0;
+        float value;
+       
+        if(!unix){ // wacky VMS
+            fs.read( (char*)&tmp1, 1 ); // 1 byte each
+            fs.read( (char*)&tmp2, 1 );
+            fs.read( (char*)&tmp3, 1 );
+            fs.read( (char*)&tmp4, 1 );
+            if((tmp2 & 0x7f) == 0){
+                if(tmp1+tmp3+tmp4 != 0) {
+                    tmp1=tmp2=tmp3=0;
+                }   
+                tmp2=0;
+            }
+            else tmp2--;
+            
+            i=(tmp2 << 24) + (tmp1 << 16) + (tmp4 << 8) + tmp3;
+        }
+        else{
+            fs.read( (char*)&tmp1, 1 );
+            fs.read( (char*)&tmp2, 1 );
+            fs.read( (char*)&tmp3, 1 );
+            fs.read( (char*)&tmp4, 1 );
+            i=(tmp4 << 24) + (tmp3 << 16) + (tmp2 << 8 ) +tmp1;
+        }
+        // grab int pointer, cast to float, dereference
+        value = *reinterpret_cast<float *>( &i );
+        return (double)value;
+    }
+
+    int OdfIO::readInteger(std::istream &fs){
+        size_t intsize = sizeof(int);
+        int val;
+        fs.read( (char*)&val, intsize );
+        return val;
+    }
+
+    std::vector<double> OdfIO::getSection(int start, int readSize, int section, 
+                                          std::istream &in){
+        in.seekg(0,std::ios::beg);
+        std::vector<double> dataSec;
+        int block,skipped,channel,max;
+        size_t intsize = sizeof(int);
+
+        dataSec.resize(readSize);
+
+        block = start/125;
+        channel = start - block * 125;
+        block = head["ifb"] + block * head["numSection"];
+        block--; 
+        skipped=( block + (section-1) )*512;     // skip 512 bytes per block
+        skipped += 4;                            // account for control word
+        skipped += channel*4;                    // and add channels not to read
+        in.ignore(skipped);                      // jump skipped-many bytes, alt: in.seekg(skipped,std::ios::cur);
+
+        max = std::min(readSize,125-channel);
+
+        int read = 0;
+        for(int i=0;i<max;i++){
+            if(      head["mode"] == 3 ) dataSec[read] = readFloat(in,head["unix"]);
+            else if( head["mode"] == 1 ) dataSec[read] = readInteger(in);
+            read++;
+        }
+        in.ignore(intsize*2); // skip over two ints in the block which contain no data
+
+        while( read<readSize ){
+            skipped=(head["numSection"]-1)*512;  // skip over other sections
+            in.ignore(skipped);
+            in.ignore(intsize);          // discard control word
+
+            for( int i=0; i<125 && read<readSize; i++ ) {
+                if(     head["mode"] == 3) dataSec[read] = readFloat(in,head["unix"]);
+                else if(head["mode"] == 1) dataSec[read] = readInteger(in);
+                read++;
+            }
+            in.ignore(intsize*2); // skip over two ints in the block which contain no data
+        }
+
+        return dataSec;
+    }
+
+    std::vector<std::vector<double>> OdfIO::getAllSections(std::istream &in){
+        // --- go to start of file ---
+        in.seekg(0, std::ios::beg);
+
+        std::vector<std::vector<double>> data;
+        data.resize(head["numSection"]);
+
+        // --- grab each section ---
+        for( int section=1; section<=head["numSection"]; ++section ){
+            data[section-1].resize(head["numChan"]);
+            data[section-1] = getSection(0,head["numChan"],section,in);
+        }
+
+        return data;
+    }
+
+    void OdfIO::readHeader(std::istream &in, int unix){
+        size_t intsize = sizeof(int);
+        int bytes_per_block = 512;
+        int bytes_per_word = 4;
+        head["unix"] = unix;
+
+        // comments based on OPRODF manual: Jack Craven, report no. ORNL/CSD/TM-45, 1978
+        //  ** keep in mind that changes may have happened between that report and
+        //  ** the writing of the files!!! For example it seems that 36 bit ints are no 
+        //  ** longer used (32 bit now). 
+        in.read( (char*)&head["mode"], intsize ); // 0,1,3 for 18-bit ints, 36-bit ints, floating-point vals
+        in.read( (char*)&head["source"], intsize ); // 0,1,2 for SEL, CSISRS, or ENDF/B data
+        in.read( (char*)&head["irun"], intsize ); // run number
+        in.read( (char*)&head["ncblks"], intsize ); // starting block no. for the comment section
+        in.read( (char*)&head["ncwrds"], intsize ); // number of words in comment section
+        in.read( (char*)&head["nsblks"], intsize ); // starting block for scalar section
+        in.read( (char*)&head["nswrds"], intsize ); // number of words in scalar section
+        in.read( (char*)&head["ncstrt"], intsize ); // starting word number of DAQ scalar and counter section in scalar section
+        in.read( (char*)&head["ncntrs"], intsize ); // number of words in DAQ scalar and counter section
+        in.read( (char*)&head["nxstrt"], intsize ); // starting word number of DAQ variable section in scalar section
+        in.read( (char*)&head["nxwrds"], intsize ); // number of words in DAQ variable section
+        in.read( (char*)&head["npblks"], intsize ); // starting block number of parameter section
+        in.read( (char*)&head["npwrds"], intsize ); // number of words in parameter section
+        in.read( (char*)&head["ndtype"], intsize ); // 0,1 for data section described by table in parameter section, data section is dependent on dataset 1 (i.e. datasets 2-x correspond to energy in dataset 1)
+        in.read( (char*)&head["numSection"], intsize ); // number of sections in each dataset
+        in.read( (char*)&head["ifb"], intsize ); // 
+        in.read( (char*)&head["numChan"], intsize ); // number of channels in each section
+        in.read( (char*)&head["zan"], intsize ); // ENDF/B charge/mass id
+        in.read( (char*)&head["awr"], intsize ); // ENDF/B ratio of nuclear mass (iso,ele,molec,mix) to neutron
+        in.read( (char*)&head["mat"], intsize ); // ENDF/B MAT number
+        in.read( (char*)&head["mf"], intsize ); // ENDF/B file number
+        in.read( (char*)&head["mt"], intsize ); // ENDF/B reaction type
+        in.read( (char*)&head["varswt"], intsize ); // 0,1 for dataset 1 decreasing, and increasing in value (ignored if ndtype is 0)
+        in.read( (char*)&head["dctswt"], intsize ); // 0,1 for not dead-time corrected, and dead-time corrected
+        in.read( (char*)&head["strt"], intsize ); // starting word number of each dataset in data section (mode 0 only)
+        in.read( (char*)&head["ndwend"], intsize ); // ending word number of last word written in block NDBEND
+
+        
+        /* TODO: the following doesn't work properly yet, we ignore for now */
+        // ------------------------------------
+        // --- Read crunch table
+        // ------------------------------------
+        // energy index table. each word contains the largest energy for each n blocks in the data set. N=(NDWRDS/125)+1
+        ndtbl.resize(100);
+        for( int i=0;i<100;++i ) in.read( (char*)&ndtbl[i], intsize );
+        // ------------------------------------
+        // --- comments section ---
+        // ------------------------------------
+        int bytes_to_comments = (head["ncblks"]) * bytes_per_block;
+        std::vector<unsigned char> comments;
+        int bytes_in_comments = head["ncwrds"]*bytes_per_word;
+        comments.resize( bytes_in_comments );
+
+        in.seekg(bytes_to_comments,std::ios::beg);
+        in.ignore(bytes_per_word);
+
+        for( int i=0;i<bytes_in_comments;++i ){
+            in.read( (char*)&comments[i], 1 );
+            comments[i] = comments[i];
+        }
+        // ------------------------------------
+        // --- scalars section ---
+        // ------------------------------------
+        int bytes_to_scalars = (head["nsblks"]) * bytes_per_block;
+    }
+
+    void OdfIO::printHeader(){
+        for( auto member : head ){
+            std::cout << member.first << " = " << head[member.first] << std::endl;
+        }
+
+    }
+
+}
\ No newline at end of file

sammy/src/io/OdfIO.h

+#ifndef SAMMY_ODFIO_H
+#define SAMMY_ODFIO_H
+
+#include <vector>
+#include <fstream> // ifstream
+#include <iostream>
+#include <map>
+
+
+namespace sammy{
+
+    /*****************************************************************
+     * The OdfIO class houses ORELA Data Format (ODF) reading routines
+     * so that old data files can be read and used for SAMMY.
+     ****************************************************************/
+    class OdfIO {
+    
+    public:
+        OdfIO();
+        virtual ~OdfIO() {}
+
+        /*************************************************************
+         * Read a single floating point value, based on whether it is in
+         * VAX-VMS format (?) or unix, from an ODF file
+         * 
+         * @param fs a reference to the ODF file stream
+         * @param unix whether or not the format is unix-style
+         * @return val the stored value
+         ************************************************************/
+        double readFloat(std::istream &fs, bool unix);
+
+        /*************************************************************
+         * Read a single integer
+         * @param fs a reference to the ODF file stream
+         * @return val the stored value
+         ************************************************************/
+        int readInteger(std::istream &fs);
+
+        /*************************************************************
+         * Read one section of data. Sections are typcially energy, transmission,
+         * cross section, etc. 
+         * 
+         * @param start what channel to start in
+         * @param readSize how many channels to read
+         * @param section which section to read from
+         * @param head a struct containing the header info
+         * @param in the input data file stream
+         * @return dataSec all requested channels of requested section
+         ************************************************************/
+        std::vector<double> getSection(int start, int readSize, int section, 
+                                       std::istream &in);
+
+        /*************************************************************
+         * Read all sections of data into double-vector, grabbing all info.
+         * This calls getSection() for all sections and channels
+         * 
+         * @param head a struct containing the header info
+         * @param in the input data file stream
+         * @return data all channels and all sections of data
+         ************************************************************/
+        std::vector<std::vector<double>> getAllSections(std::istream &in);
+
+        /*************************************************************
+         * Read the ODF header 
+         * 
+         * @param fs a reference to the ODF file stream
+         * @param unix 1 or 0 for true or false unix format
+         ************************************************************/
+        void readHeader(std::istream &in, int unix);
+
+        const std::map<std::string, int> getHeader() const{
+            return head;
+        }
+
+        /*************************************************************
+         * Print the ODF header 
+         * 
+         * @param head an ODF-header struct
+         ************************************************************/
+        void printHeader();
+
+    private:
+        std::map<std::string, int> head;
+        std::vector<int> ndtbl;
+        
+    };
+}
+
+#endif /* SAMMY_ODFIO_H */
\ No newline at end of file

sammy/src/io/convertodf.cpp

+#include <iostream>
+#include "H5Cpp.h"
+#include "OdfIO.h"
+#include "Hdf5IO.h"
+
+int main(int argc, char const *argv[])
+{
+    // -------------------------------------------------------------------------
+    // --- Grabbing the input args
+    // -------------------------------------------------------------------------
+    if( argc < 3 ){
+        std::cerr << "Wrong number of arguments given." << std::endl
+                  << "Specify input ODF file and output HDF5 file." << std::endl;
+        return -1;
+    }
+    
+    std::string  odfFilename = argv[1];
+    std::string hdf5Filename = argv[2];
+
+    std::cout << "Reading from file: " << odfFilename << std::endl;
+    std::cout << "Writing to file:   " << hdf5Filename << std::endl;
+
+    // -------------------------------------------------------------------------
+    // --- Read ODF file
+    // -------------------------------------------------------------------------
+    sammy::OdfIO odf;
+    std::vector<std::vector<double>> data;
+
+    // --- Open file ---
+    std::ifstream in( odfFilename, std::ios::in | std::ios::binary );
+    // --- kill bad open  ---
+    if( !in.is_open() ){
+        std::cerr << "Error opening file: " << odfFilename << std::endl;
+        return 1;
+    }
+
+    int unix = false;
+    // --- read header ---
+    odf.readHeader(in,unix);
+
+    // --- read data sections ---
+    data = odf.getAllSections(in);
+
+    in.close();
+
+    // -------------------------------------------------------------------------
+    // --- Output HDF5 file
+    // -------------------------------------------------------------------------
+
+    sammy::Hdf5IO h5writer;
+    h5writer.writeODFtoHDF5(odf.getHeader(),data,hdf5Filename);
+
+    return 0;
+}
\ No newline at end of file

sammy/src/io/tests/CMakeLists.txt

 include(NemesisTest)
 
-ADD_NEMESIS_TEST(ExperimentalParametersIOTest.cpp    TIMEOUT 60000 NP 1 )
+ADD_NEMESIS_TEST(ExperimentalParametersIOTest.cpp TIMEOUT 60000 NP 1 )
+ADD_NEMESIS_TEST(OdfIOTest.cpp                    TIMEOUT 60000 NP 1 )

sammy/src/io/tests/OdfIOTest.cpp

+#include "Nemesis/gtest/Gtest_Functions.hh"
+#include "Nemesis/gtest/nemesis_gtest.hh"
+#include "../OdfIO.h"
+
+/*******************************************************************************
+ * Test whether we read floats with VAX format properly
+ ******************************************************************************/
+TEST(OdfIO,readFloatVax)
+{
+    sammy::OdfIO odf;
+    bool is_unix = false;
+    double readVal;
+    float * dummy;
+    double gold;
+    char byte1 = 0b00000001;
+    char byte2 = 0b00000010;
+    char byte3 = 0b00000011;
+    char byte4 = 0b00000100;
+    //                  |       |       |       |
+    int vax =  0b00000001000000010000010000000011; 
+    // byte order     2--,      1,      4,      3 
+
+    std::stringstream in;
+    in << byte1 << byte2 << byte3 << byte4;
+
+    readVal = odf.readFloat(in,is_unix);
+
+    dummy = (float *) &vax;
+    gold = *dummy;
+
+    ASSERT_EQ( readVal, gold );
+}
+
+/*******************************************************************************
+ * Test whether we read floats with UNIX format properly
+ ******************************************************************************/
+TEST(OdfIO,readFloatUnix)
+{
+    sammy::OdfIO odf;
+    bool is_unix = true;
+    double readVal;
+    float * dummy;
+    double gold;
+    char byte1 = 0b00000001;
+    char byte2 = 0b00000010;
+    char byte3 = 0b00000011;
+    char byte4 = 0b00000100;
+    //                   |       |       |       |
+    int unix =  0b00000100000000110000001000000001;
+
+    std::stringstream in;
+    in << byte1 << byte2 << byte3 << byte4;
+
+    readVal = odf.readFloat(in,is_unix);
+
+    dummy = (float *) &unix;
+    gold = *dummy;
+
+    ASSERT_EQ( readVal, gold );
+}
+
+/*******************************************************************************
+ * Test whether we read integers right
+ ******************************************************************************/
+TEST(OdfIO,readInteger)
+{
+    sammy::OdfIO odf;
+    int gold = 16909060;
+    char byte1 = 0b00000001;
+    char byte2 = 0b00000010;
+    char byte3 = 0b00000011;
+    char byte4 = 0b00000100;
+    int readVal;
+
+    std::stringstream in;
+    in << byte4 << byte3 << byte2 << byte1;
+
+    readVal = odf.readInteger(in);
+
+    ASSERT_EQ( readVal, gold );
+}
+