From ba5546534f6e6c83eeb863c702f3511302f97304 Mon Sep 17 00:00:00 2001
From: Yohai Meiron <sahmes@gmail.com>
Date: Fri, 27 Mar 2020 16:53:01 -0400
Subject: [PATCH] Continued working on HDF5 plugin

---
 io.cpp | 186 ++++++++++++++++++++++++++-------------------------------
 1 file changed, 84 insertions(+), 102 deletions(-)

diff --git a/io.cpp b/io.cpp
index 40f009c..b1286af 100644
--- a/io.cpp
+++ b/io.cpp
@@ -2,11 +2,13 @@
 #include <string>
 #include <vector>
 #include <algorithm>
-struct double3 {double x,y,z;};
+#include <stdexcept>
+struct double3 {double x,y,z;
+};
 
-void h5_write(const std::string file_name, const int step_num, const double t, const int N, const double *m, const double3 *x, const double3 *v, const double *pot, const double3 *acc, const double3 *jrk, const int write_mode=0, const bool use_double_precision=true)
+
+void h5_write(const std::string file_name, const int step_num, const int N, const double t, const double *m, const double3 *x, const double3 *v, const double *pot, const double3 *acc, const double3 *jrk, const int write_mode=0, const bool use_double_precision=true)
 {
-    char group_name[32], dataset_path[32];
     bool write_pot = (write_mode   )%2;
     bool write_acc = (write_mode>>1)%2;
     bool write_jrk = (write_mode>>2)%2;
@@ -19,6 +21,7 @@ void h5_write(const std::string file_name, const int step_num, const double t, c
     else h5_float_type = H5T_IEEE_F32LE;
 
     file_id = H5Fcreate(file_name.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
+    char group_name[32];
     sprintf(group_name, "/Step#%d", step_num);
     group_id = H5Gcreate2(file_id, group_name, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
     dataspace_id = H5Screate(H5S_SCALAR);
@@ -26,97 +29,89 @@ void h5_write(const std::string file_name, const int step_num, const double t, c
     H5Awrite(attribute_id, H5T_NATIVE_DOUBLE, &t);
     H5Sclose(dataspace_id);
 
-    dataspace_id = H5Screate_simple(1, dims, NULL);
-    sprintf(dataset_path, "%s/MASS", group_name);
-    dataset_id = H5Dcreate2(file_id, dataset_path, h5_float_type, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
-    H5Dwrite(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, m);
-    H5Dclose(dataset_id);
-    H5Sclose(dataspace_id);
-
-    dataspace_id = H5Screate_simple(2, dims, NULL);
-    sprintf(dataset_path, "%s/X", group_name);
-    dataset_id = H5Dcreate2(file_id, dataset_path, h5_float_type, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
-    H5Dwrite(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, x);
-    H5Dclose(dataset_id);
-    H5Sclose(dataspace_id);
-
-    dataspace_id = H5Screate_simple(2, dims, NULL);
-    sprintf(dataset_path, "%s/V", group_name);
-    dataset_id = H5Dcreate2(file_id, dataset_path, h5_float_type, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
-    H5Dwrite(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, v);
-    H5Dclose(dataset_id);
-    H5Sclose(dataspace_id);
-
-    if (write_pot) {
-        dataspace_id = H5Screate_simple(1, dims, NULL);
-        sprintf(dataset_path, "%s/POT", group_name);
-        dataset_id = H5Dcreate2(file_id, dataset_path, h5_float_type, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
-        H5Dwrite(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, pot);
+    auto write_dataset = [&](const char dataset_name[], int ndims, double *data) {
+        hid_t dataspace_id = H5Screate_simple(ndims, dims, NULL);
+        char dataset_path[32];
+        sprintf(dataset_path, "%s/%s", group_name, dataset_name);
+        hid_t dataset_id = H5Dcreate2(file_id, dataset_path, h5_float_type, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+        H5Dwrite(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, data);
         H5Dclose(dataset_id);
         H5Sclose(dataspace_id);
-    }
-    if (write_acc) {
-        dataspace_id = H5Screate_simple(2, dims, NULL);
-        sprintf(dataset_path, "%s/ACC", group_name);
-        dataset_id = H5Dcreate2(file_id, dataset_path, h5_float_type, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
-        H5Dwrite(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, acc);
-        H5Dclose(dataset_id);
-        H5Sclose(dataspace_id);
-    }
-    if (write_jrk) {
-        dataspace_id = H5Screate_simple(2, dims, NULL);
-        sprintf(dataset_path, "%s/JRK", group_name);
-        dataset_id = H5Dcreate2(file_id, dataset_path, h5_float_type, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
-        H5Dwrite(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, jrk);
-        H5Dclose(dataset_id);
-        H5Sclose(dataspace_id);
-    }
+    };
+
+    write_dataset("MASS", 1, (double*)m);
+    write_dataset("X",    2, (double*)x);
+    write_dataset("V",    2, (double*)v);
+    if (write_pot) write_dataset("POT", 1, (double*)pot);
+    if (write_acc) write_dataset("ACC", 2, (double*)acc);
+    if (write_jrk) write_dataset("JRK", 2, (double*)jrk);
+
     H5Gclose(group_id);
-
-    
-    group_id = H5Gcreate2(file_id, "/Step#22", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
-    H5Gclose(group_id);
-    
-    
     H5Fclose(file_id);
 }
 
 
-
-
-herr_t file_info(hid_t loc_id, const char *name, void *opdata)
+void h5_read(const std::string file_name, int *step_num, int *N, double *t, double m[], double3 x[], double3 v[])
 {
-    H5G_stat_t statbuf;
-
-    /*
-     * Get type of the object and display its name and type.
-     * The name of the object is passed to this function by 
-     * the Library. Some magic :-)
-     */
-    H5Gget_objinfo(loc_id, name, 0, &statbuf);
-    switch (statbuf.type) {
-    case H5G_GROUP: 
-         printf(" Object with name %s is a group \n", name);
-         break;
-    case H5G_DATASET: 
-         printf(" Object with name %s is a dataset \n", name);
-         break;
-    case H5G_TYPE: 
-         printf(" Object with name %s is a named datatype \n", name);
-         break;
-    default:
-         printf(" Unable to identify an object ");
+    // Open file and root group; count number of top level objects
+    hid_t file_id = H5Fopen(file_name.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
+    hid_t group_id = H5Gopen2(file_id, "/", H5P_DEFAULT);
+    H5G_info_t object_info;
+    H5Gget_info(group_id, &object_info);
+    // Iterate over objects and add the number of each "step" group into a vector
+    std::vector<int> step_num_arr;
+    for (int i=0; i<object_info.nlinks; i++) {
+        char name_cstr[256];
+        H5Gget_objname_by_idx(group_id, i, name_cstr, 256);
+        std::string name(name_cstr);
+        if (name.substr(0, 5) == "Step#") step_num_arr.push_back(std::stoi(name.substr(5, std::string::npos)));
     }
-    return 0;
- }
+    // Find the highest step in the file
+    *step_num = *max_element(step_num_arr.begin(), step_num_arr.end());
 
+    // Prepare to read the data sets dimensionality data.
+    char dataset_path[256];
+    hid_t dataset_id, dataspace_id;
+    int ndims;
+    hsize_t dims[2];
 
+    sprintf(dataset_path, "/Step#%d/MASS", *step_num);
+    dataset_id = H5Dopen2(file_id, dataset_path, H5P_DEFAULT);
+    dataspace_id = H5Dget_space(dataset_id);
+    ndims = H5Sget_simple_extent_ndims(dataspace_id);
+    if (ndims != 1)
+        throw std::runtime_error("Dataset MASS in Step#" + std::to_string(*step_num) + " of file " + file_name + " is " + std::to_string(ndims) + "-dimensional (expected 1-dimensional)");
+    H5Sget_simple_extent_dims(dataspace_id, dims, NULL);
+    H5Dread(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, m);
+    *N = dims[0];
+
+    sprintf(dataset_path, "/Step#%d/X", *step_num);
+    dataset_id = H5Dopen2(file_id, dataset_path, H5P_DEFAULT);
+    dataspace_id = H5Dget_space(dataset_id);
+    ndims = H5Sget_simple_extent_ndims(dataspace_id);
+    if (ndims != 2) throw std::runtime_error("Dataset X in Step#" + std::to_string(*step_num) + " of file " + file_name + " is " + std::to_string(ndims) + "-dimensional (expected 2-dimensional)");
+    H5Sget_simple_extent_dims(dataspace_id, dims, NULL);
+    if (dims[0] != *N) throw std::runtime_error("Expected " + std::to_string(*N) + "rows in Dataset X in Step#" + std::to_string(*step_num) + " of file " + file_name + " (based on MASS)");
+    if (dims[1] != 3) throw std::runtime_error("Dataset X in Step#" + std::to_string(*step_num) + " of file " + file_name + " has " + std::to_string(dims[1]) + " columns (expected 3 columns)");
+    H5Dread(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, (double*)x);
+
+    sprintf(dataset_path, "/Step#%d/V", *step_num);
+    dataset_id = H5Dopen2(file_id, dataset_path, H5P_DEFAULT);
+    dataspace_id = H5Dget_space(dataset_id);
+    ndims = H5Sget_simple_extent_ndims(dataspace_id);
+    if (ndims != 2) throw std::runtime_error("Dataset V in Step#" + std::to_string(*step_num) + " of file " + file_name + " is " + std::to_string(ndims) + "-dimensional (expected 2-dimensional)");
+    H5Sget_simple_extent_dims(dataspace_id, dims, NULL);
+    if (dims[0] != *N) throw std::runtime_error("Expected " + std::to_string(*N) + "rows in Dataset V in Step#" + std::to_string(*step_num) + " of file " + file_name + " (based on MASS)");
+    if (dims[1] != 3) throw std::runtime_error("Dataset V in Step#" + std::to_string(*step_num) + " of file " + file_name + " has " + std::to_string(dims[1]) + " columns (expected 3 columns)");
+    H5Dread(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, (double*)v);
+}
 
 
 int main()
 {
 //     double3 aaa;
 //     double pot;
+    std::string file_name = "new_file.h5";
     int N = 25;
     double3 x[N];
     for (int i=0; i<N; i++) {
@@ -124,31 +119,18 @@ int main()
         x[i].y=i+0.02;
         x[i].z=i+0.03;
     }
-    h5_write("new_file.h5", 34, 987.654, N, (double*)x, x, x, (double*)x, x, x, 7);
+    double m[N];
+    double3 v[N];
+    int step_num;
+    double t;
     
-    auto file_id = H5Fopen("new_file.h5", H5F_ACC_RDONLY, H5P_DEFAULT);
-    auto group_id = H5Gopen2(file_id, "/", H5P_DEFAULT);
-    H5G_info_t object_info;
-    H5Gget_info(group_id, &object_info);
-    std::vector<int> step_num;
-    for (int i=0; i<object_info.nlinks; i++) {
-        char name_cstr[256];
-        H5Gget_objname_by_idx(group_id, i, name_cstr, 256);
-        std::string name(name_cstr);
-        if (name.substr(0, 5) == "Step#") step_num.push_back(std::stoi(name.substr(5, std::string::npos)));
-    }
-    int step_num_max = *max_element(step_num.begin(), step_num.end());
+    h5_write("123.h5", 17, N, 123.987, (double*)x, x, x, (double*)x, x, x, 7, true);
+    h5_read("123.h5", &step_num, &N, &t, m, x, v);
+    printf("%d %d %f\n", step_num, N, t);
+#error still havent implemented reading the time attribute
     
-    
-    char dataset_path[256];
-    sprintf(dataset_path, "/Step#%d/X", step_num_max);
-    
-    
-    double dset_data[1024*1024];
-    
-    auto dataset_id = H5Dopen2(file_id, dataset_path, H5P_DEFAULT);
-    H5Dread(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, dset_data);
-    
-    
-//     printf("aaaa %d\n", step_num_max);
+//     h5_write(file_name, 34, 987.654, N, (double*)x, x, x, (double*)x, x, x, 7);
+
+//     h5_read(file_name, &N, (double*)x, x, x);
+//     printf("aaaa %d %d %d\n", ndims, dims[0], dims[1]);
 }