add hdf5 writing to executable

sammy/src/io/CMakeLists.txt

@@ -2,15 +2,18 @@ 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_DIR})
 
 
 # Set headers
 SET(HEADERS ExperimentalParametersIO.h
             odfIO.h
+            hdf5IO.h
             interface/cpp/ExperimentalParametersIOInterface.h 
 )
 
@@ -18,6 +21,7 @@ APPEND_SET(IO_SOURCES
 
         ExperimentalParametersIO.cpp
         odfIO.cpp
+        hdf5IO.cpp
         interface/cpp/ExperimentalParametersIOInterface.cpp
         interface/fortran/ExperimentalParametersIO_M.f90   
         interface/fortran/ExperimentalParametersIO_I.f90      
@@ -36,3 +40,15 @@ INSTALL(FILES ${HEADERS} DESTINATION "src/io")
 TRIBITS_ADD_TEST_DIRECTORIES(tests)
 
 TRIBITS_PACKAGE_POSTPROCESS()
+
+
+# --------------------------------------------
+# Build separate executable to convert ODF
+# to HDF5 files, install to normal location
+# --------------------------------------------
+
+add_executable(readodf readodf.cpp)
+target_link_libraries(readodf PRIVATE HDF5::HDF5
+                                      SammyIOUtilsLib)
+install(PROGRAMS ${CMAKE_CURRENT_BINARY_DIR}/readodf  DESTINATION "bin/")
+

sammy/src/io/hdf5IO.cpp

+#include "hdf5IO.h"
+
+namespace sammy {
+    
+    hdf5IO::hdf5IO() {}
+
+    void hdf5IO::writeODFtoHDF5(){
+        std::cout << "hello world" << std::endl;
+    }
+
+}
\ No newline at end of file

sammy/src/io/hdf5IO.h

+#ifndef SAMMY_HDF5IO_H
+#define SAMMY_HDF5IO_H
+
+#include <vector>
+#include <fstream> // ifstream
+#include <iostream>
+#include "hdf5.h"
+
+namespace sammy{
+
+    class hdf5IO {
+    
+    public:
+        hdf5IO();
+        virtual ~hdf5IO() {}
+
+        /****************************************
+         * Write data from ODF file into HDF5
+         ***************************************/
+        void writeODFtoHDF5();
+        
+    };
+}
+
+#endif /* SAMMY_HDF5IO_H */
\ No newline at end of file

sammy/src/io/odfIO.cpp

@@ -16,7 +16,7 @@ namespace sammy{
             fs.read( (char*)&tmp2, 1 );
             fs.read( (char*)&tmp3, 1 );
             fs.read( (char*)&tmp4, 1 );
-            // std::cout << "1,2,3,4 = " << (int)tmp1<<","<< (int)tmp2<<","<< (int)tmp3<<","<< (int)tmp4<<","<<std::endl;
+            // to print: std::cout << "tmp1 = " << (int)tmp1 << std::endl;
             if((tmp2 & 0x7f) == 0){      
                 if(tmp1+tmp3+tmp4 != 0) {
                     tmp1=tmp2=tmp3=0;
@@ -39,8 +39,15 @@ namespace sammy{
         return value;
     }
 
+    int odfIO::readInteger(std::ifstream &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, 
-                                          header &hd, std::ifstream &in){
+                                          odfHeader &hd, std::ifstream &in){
         in.seekg(0,std::ios::beg);
         std::vector<double> dataSec;
         int block,skipped,channel,max;
@@ -62,7 +69,7 @@ namespace sammy{
         int read = 0;
         for(int i=0;i<max;i++){
             if(      hd.mode == 3 ) dataSec[read] = readFloat(in,hd.unix);
-            else if( hd.mode == 1 ) dataSec[read] = readFloat(in,hd.unix); // needs to read integer not float
+            else if( hd.mode == 1 ) dataSec[read] = readInteger(in);
             read++;
         }
         in.ignore(intsize*2); // skip over two ints in the block which contain no data
@@ -74,7 +81,7 @@ namespace sammy{
 
             for( int i=0; i<125 && read<readSize; i++ ) {
                 if(     hd.mode == 3) dataSec[read] = readFloat(in,hd.unix);
-                else if(hd.mode == 1) dataSec[read] = readFloat(in,hd.unix); // needs to read integer not float
+                else if(hd.mode == 1) dataSec[read] = readInteger(in);
                 read++;
             }
             in.ignore(intsize*2); // skip over two ints in the block which contain no data
@@ -83,11 +90,31 @@ namespace sammy{
         return dataSec;
     }
 
-    void odfIO::readHeader(std::ifstream &in, header &head){
+    std::vector<std::vector<double>> odfIO::getAllSections(odfHeader &head, 
+                                                           std::ifstream &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,head,in);
+        }
+
+        return data;
+    }
+
+    void odfIO::readHeader(std::ifstream &in, odfHeader &head){
         size_t intsize = sizeof(int);
         int bytes_per_block = 512;
         int bytes_per_word = 4;
         // 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. 
         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
@@ -142,7 +169,7 @@ namespace sammy{
         int bytes_to_scalars = (head.nsblks) * bytes_per_block;
     }
 
-    void odfIO::printHeader(header &head){
+    void odfIO::printHeader(odfHeader &head){
         std::cout << "mode       = " << head.mode       << std::endl       
                   << "source     = " << head.source     << std::endl     
                   << "irun       = " << head.irun       << std::endl       

sammy/src/io/odfIO.h

@@ -7,35 +7,43 @@
 
 namespace sammy{
 
+    /****************************************
+     * Object for ODF header info
+     ***************************************/
+    struct odfHeader {
+        // variables explained in readHeader()
+        int mode, source, irun;
+        int ncblks, ncwrds, nsblks, nswrds, ncstrt, ncntrs, nxstrt, nxwrds;
+        int npblks, npwrds, ndtype, numSection, ifb, numChan;
+        int zan, awr, mat, mf, mt; // ENDF identifiers
+        int varswt,dctswt,strt,ndwend;
+        bool unix;
+        std::vector<int> ndtbl;
+    };
+
     class odfIO {
     
     public:
         odfIO();
         virtual ~odfIO() {}
 
-        /****************************************
-         * Object for ODF header info
-         ***************************************/
-        struct header {
-            // variables explained in readHeader()
-            int mode, source, irun;
-            int ncblks, ncwrds, nsblks, nswrds, ncstrt, ncntrs, nxstrt, nxwrds;
-            int npblks, npwrds, ndtype, numSection, ifb, numChan;
-            int zan, awr, mat, mf, mt; // ENDF identifiers
-            int varswt,dctswt,strt,ndwend;
-            bool unix;
-            std::vector<int> ndtbl;
-        };
-
         /****************************************
          * 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::ifstream &fs, bool unix);
 
+        /****************************************
+         * Read a single integer
+         * @param fs a reference to the ODF file stream
+         * @return val the stored value
+         ***************************************/
+        int readInteger(std::ifstream &fs);
+
         /****************************************
          * Read one section of data. Sections are typcially energy, transmission,
          * cross section, etc. 
@@ -43,9 +51,23 @@ namespace sammy{
          * @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, 
-                                       header &hd, std::ifstream &in);
+                                       odfHeader &head, std::ifstream &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(odfHeader &head, 
+                                                        std::ifstream &in);
 
         /****************************************
          * Read the ODF header 
@@ -53,14 +75,14 @@ namespace sammy{
          * @param fs a reference to the ODF file stream
          * @param head an ODF-header struct
          ***************************************/
-        void readHeader(std::ifstream &in, header &head);
+        void readHeader(std::ifstream &in, odfHeader &head);
 
         /****************************************
          * Print the ODF header 
          * 
          * @param head an ODF-header struct
          ***************************************/
-        void printHeader(header &head);
+        void printHeader(odfHeader &head);
         
     };
 }

sammy/src/io/readodf.cpp

+#include <iostream>
+// #include "hdf5.h"
+#include "H5Cpp.h"
+#include "odfIO.h"
+
+#define MAX_NAME_LEN 32
+
+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;
+    sammy::odfHeader head;
+    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;
+    }
+
+    // --- read header ---
+    odf.readHeader(in, head);
+
+    // --- read data sections ---
+    data = odf.getAllSections(head, in);
+
+    in.close();
+
+    // -------------------------------------------------------------------------
+    // --- Output HDF5 file
+    // -------------------------------------------------------------------------
+    
+    char sectnames[11][MAX_NAME_LEN] = {"energy","trans","dtrans",
+                              "cro","dcro","cs","zero","co","zero","idk","dcs"};
+
+    // make dataset 
+    double dset_data[data.size()][data[0].size()];
+    for( size_t i=0;i<data.size();++i ){
+        std::copy(data[i].begin(),data[i].end(), dset_data[i]);
+    }
+
+    // write 2 datasets: data and data section names
+
+    // data space
+    hsize_t ddim[2];
+    ddim[0] = head.numSection;
+    ddim[1] = head.numChan;
+
+    int drank = sizeof(ddim) / sizeof(hsize_t);
+
+    std::string dataname("/data/sections");
+
+    // section names space
+    hsize_t sdim[1];
+    sdim[0] = head.numSection;
+
+    int srank = sizeof(sdim) / sizeof(hsize_t);
+
+    std::string sectionnames("/data/sectnames");
+
+
+    // open file and write
+    H5::DataSpace dspace(drank,ddim);
+    H5::DataSpace sspace(srank,sdim);
+    H5::H5File  *file = new H5::H5File(hdf5Filename,H5F_ACC_TRUNC);
+    H5::Group   *group = new H5::Group( file->createGroup( "/data" ));
+
+    H5::DataSet *dataset = new H5::DataSet(file->createDataSet(dataname,H5::PredType::NATIVE_DOUBLE,dspace));
+    H5::DataSet *nameset = new H5::DataSet(file->createDataSet(sectionnames,H5::StrType(H5::PredType::C_S1, MAX_NAME_LEN),sspace));
+
+    dataset->write(dset_data,H5::PredType::NATIVE_DOUBLE);
+    nameset->write(sectnames,H5::StrType(H5::PredType::C_S1, MAX_NAME_LEN));
+
+    delete file;
+    delete dataset;
+    delete nameset;
+
+    
+    /* old
+
+    // --- setup with HDF types ---
+    hid_t   file_id, dataspace_id, dataset_id, group_id, name_dataspace_id;
+    hid_t   name_dataset_id;
+    herr_t  status;
+    hsize_t dims[2];
+
+    //energy,trans,dtrans,cro,dcro,cs,zero,co,zero,idk,dcs
+
+    dims[0] = data.size();
+    dims[1] = data[0].size();
+    dataspace_id = H5Screate_simple(2, dims, NULL);
+
+    // --- open file ---
+    file_id = H5Fcreate(hdf5Filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
+
+    // --- make a dataset group ---
+    group_id = H5Gcreate2(file_id,"data",H5P_DEFAULT,H5P_DEFAULT,H5P_DEFAULT);
+
+    // --- create dataset ---
+    // H5T_IEEE_F64LE is an HDF5 predefined datatype for 64 bit floating point,
+    // little endian. Might replace with H5T_NATIVE_DOUBLE
+    dataset_id = H5Dcreate2(file_id, "/data/sections", H5T_IEEE_F64LE, dataspace_id, 
+                            H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+
+    double dset_data[data.size()][data[0].size()];
+    for( size_t i=0;i<data.size();++i ){
+        std::copy(data[i].begin(),data[i].end(), dset_data[i]);
+    }
+
+    // --- Write the dataset to file --- 
+    status = H5Dwrite(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT,
+                      dset_data);
+
+
+    // --- write column names for dataset ---
+    std::vector<std::string> colnames; 
+    colnames.resize(11);
+    colnames = {"energy","trans","dtrans","cro","dcro","cs",
+                                "zero","co","zero","idk","dcs"};
+    std::string name = "energy";
+    const char* cname = name.c_str();
+
+    hid_t datatype = H5Tcopy (H5T_C_S1);
+    size_t MAX_NAME_LENGTH = 32;
+    H5Tset_size (datatype, H5T_VARIABLE);
+    hsize_t name_dims[1];
+    name_dims[0] = 1;
+    name_dataspace_id = H5Screate_simple(1, name_dims, NULL);
+    name_dataset_id = H5Dcreate2(file_id, "/data/colnames", datatype, 
+                                 name_dataspace_id, H5P_DEFAULT, H5P_DEFAULT, 
+                                 H5P_DEFAULT);
+    for( auto name : colnames ){
+        const char* cname = name.c_str();
+        status = H5Dwrite(name_dataset_id, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT,
+                      &cname);
+    }
+
+
+    status = H5Fclose(file_id);
+
+    */
+
+    return 0;
+}
\ No newline at end of file