Implement complete configuration I/O to HDF5

src/libnptm/Configuration.cpp

@@ -538,6 +538,7 @@ ScattererConfiguration* ScattererConfiguration::from_hdf5(string file_name) {
   unsigned int flags = H5F_ACC_RDONLY;
   HDFFile *hdf_file = new HDFFile(file_name, flags);
   herr_t status = hdf_file->get_status();
+  string str_name, str_type;
   if (status == 0) {
     int nsph;
     int *iog;
@@ -562,10 +563,47 @@ ScattererConfiguration* ScattererConfiguration::from_hdf5(string file_name) {
       if (iog[ci] < ci + 1) continue;
       configuration_count++;
     }
-    nshl_vector = new int[configuration_count];
-    ros_vector = new double[configuration_count];
-    // DA QUI
+    nshl_vector = new int[configuration_count]();
+    ros_vector = new double[configuration_count]();
+    rcf_vector = new double*[configuration_count];
+    for (int i115 = 1; i115 <= nsph; i115++) {
+      if (iog[i115 - 1] < i115) continue;
+      str_name = "NSHL_" + to_string(iog[i115 - 1]);
+      str_type = "INT32_(1)";
+      status = hdf_file->read(str_name, str_type, (nshl_vector + i115 - 1));
+      str_name = "ROS_" + to_string(iog[i115 - 1]);
+      str_type = "FLOAT64_(1)";
+      status = hdf_file->read(str_name, str_type, (ros_vector + i115 - 1));
+      int nsh = nshl_vector[i115 - 1];
+      rcf_vector[i115 - 1] = new double[nsh]();
+      str_name = "RCF_" +  to_string(iog[i115 - 1]);
+      str_type = "FLOAT64_(" + to_string(nsh) + ")";
+      status = hdf_file->read(str_name, str_type, (rcf_vector[i115 - 1]));
+    }
+    xi_vec = new double[nxi];
+    str_name = "XIVEC";
+    str_type = "FLOAT64_(" + to_string(nxi) + ")";
+    status = hdf_file->read(str_name, str_type, xi_vec);
 
+    int dim3 = (_idfc == 0) ? nxi : 1;
+    int element_size = 2 * dim3 * nsph * configuration_count;
+    double *elements = new double[element_size]();
+    str_name = "DC0M";
+    str_type = "FLOAT64_(" + to_string(element_size) + ")";
+    status = hdf_file->read(str_name, str_type, elements);
+    dc0m = new complex<double>**[configuration_count];
+    for (int di = 0; di < configuration_count; di++) {
+      dc0m[di] = new complex<double>*[nsph];
+      for (int dj = 0; dj < nsph; dj++) {
+	dc0m[di][dj] = new complex<double>[dim3]();
+	for (int dk = 0; dk < dim3; dk++) {
+	  double rval = elements[(di * nsph) + (dj * dim3) + 2 * dk];
+	  double ival = elements[(di * nsph) + (dj * dim3) + 2 * dk + 1];
+	  dc0m[di][dj][dk] = complex<double>(rval, ival);
+	}
+      } // dj loop
+    } // di loop
+    delete[] elements;
     status = hdf_file->close();
     conf = new ScattererConfiguration(
 				      nsph,
@@ -608,6 +646,11 @@ ScattererConfiguration* ScattererConfiguration::from_legacy(string file_name) {
     for (int i = 0; i < nsph; i++) {
       input.read(reinterpret_cast<char *>(&(iog[i])), sizeof(int));
     }
+    int configuration_count = 0;
+    for (int i115 = 1; i115 <= nsph; i115++) {
+      if (iog[i115 - 1] < i115) continue;
+      configuration_count++;
+    }
     input.read(reinterpret_cast<char *>(&_exdc), sizeof(double));
     input.read(reinterpret_cast<char *>(&_wp), sizeof(double));
     input.read(reinterpret_cast<char *>(&_xip), sizeof(double));
@@ -621,25 +664,25 @@ ScattererConfiguration* ScattererConfiguration::from_legacy(string file_name) {
     for (int i = 0; i < nxi; i++) {
       input.read(reinterpret_cast<char *>(&(xi_vec[i])), sizeof(double));
     }
-    nshl_vector = new int[nsph]();
-    ros_vector = new double[nsph]();
-    rcf_vector = new double*[nsph];
+    nshl_vector = new int[configuration_count]();
+    ros_vector = new double[configuration_count]();
+    rcf_vector = new double*[configuration_count];
     for (int i115 = 1; i115 <= nsph; i115++) {
       if (iog[i115 - 1] < i115) {
 	rcf_vector[i115 - 1] = new double[1]();
 	continue;
       }
-      input.read(reinterpret_cast<char *>(&(nshl_vector[i115 - 1])), sizeof(int));
-      input.read(reinterpret_cast<char *>(&(ros_vector[i115 - 1])), sizeof(double));
-      int nsh = nshl_vector[i115 - 1];
+      input.read(reinterpret_cast<char *>(&(nshl_vector[i115 - 1])), sizeof(int)); // was not from IOG
+      input.read(reinterpret_cast<char *>(&(ros_vector[i115 - 1])), sizeof(double)); // was not from IOG
+      int nsh = nshl_vector[i115 - 1]; // was not from IOG
       if (max_ici < (nsh + 1) / 2) max_ici = (nsh + 1) / 2;
       try {
-	rcf_vector[i115 - 1] = new double[nsh]();
+	rcf_vector[i115 - 1] = new double[nsh](); // was not from IOG
       } catch (const bad_alloc &ex) {
 	throw UnrecognizedConfigurationException("Wrong parameter set: invalid number of layers " + to_string(nsh));
       }
       for (int nsi = 0; nsi < nsh; nsi++) {
-	input.read(reinterpret_cast<char *>(&(rcf_vector[i115 - 1][nsi])), sizeof(double));
+	input.read(reinterpret_cast<char *>(&(rcf_vector[i115 - 1][nsi])), sizeof(double)); // was not from IOG
       }
     }
     dc0m = new complex<double>**[max_ici];
@@ -654,7 +697,7 @@ ScattererConfiguration* ScattererConfiguration::from_legacy(string file_name) {
       if (_idfc != 0 && jxi468 > 1) continue;
       for (int i162 = 1; i162 <= nsph; i162++) {
 	if (iog[i162 - 1] < i162) continue;
-	int nsh = nshl_vector[i162 - 1];
+	int nsh = nshl_vector[iog[i162 - 1] - 1]; // was not from IOG
 	int ici = (nsh + 1) / 2; // QUESTION: is integer division really intended here?
 	for (int i157 = 0; i157 < ici; i157++) {
 	  double dc0_real, dc0_img;
@@ -787,19 +830,19 @@ void ScattererConfiguration::write_hdf5(string file_name) {
     str_name = "NSHL_" + to_string(i115);
     rec_name_list.append(str_name);
     rec_type_list.append("INT32_(1)");
-    rec_ptr_list.append(&(nshl_vec[i115 - 1]));
+    rec_ptr_list.append(&(nshl_vec[i115 - 1])); // was not from IOG
     str_name = "ROS_" + to_string(i115);
     rec_name_list.append(str_name);
     rec_type_list.append("FLOAT64_(1)");
-    rec_ptr_list.append(&(radii_of_spheres[i115 - 1]));
-    int nsh = nshl_vec[i115 - 1];
+    rec_ptr_list.append(&(radii_of_spheres[i115 - 1])); // was not from IOG
+    int nsh = nshl_vec[i115 - 1]; // was not from IOG
     if (i115 == 1) nsh += ies;
     if (max_ici < (nsh + 1) / 2) max_ici = nsh + 1 / 2;
-    str_name = "RCF_" + to_string(i115);
+    str_name = "RCF_" + to_string(i115); // was not from IOG
     str_type = "FLOAT64_(" + to_string(nsh) + ")";
     rec_name_list.append(str_name);
     rec_type_list.append(str_type);
-    rec_ptr_list.append(&(rcf[i115 - 1][0]));
+    rec_ptr_list.append(&(rcf[i115 - 1][0])); // was not from IOG
   }
 
   int dim3 = (idfc == 0) ? number_of_scales : 1;
@@ -810,7 +853,7 @@ void ScattererConfiguration::write_hdf5(string file_name) {
     if (idfc != 0 && jxi468 > 1) continue;
     for (int i162 = 1; i162 <= number_of_spheres; i162++) {
       if (iog_vec[i162 - 1] < i162) continue;
-      int nsh = nshl_vec[i162 - 1];
+      int nsh = nshl_vec[i162 - 1]; // was not from IOG
       int ici = (nsh + 1) / 2;
       if (i162 == 1) ici = ici + ies;
       for (int i157 = 0; i157 < ici; i157++) {
@@ -863,18 +906,18 @@ void ScattererConfiguration::write_legacy(string file_name) {
     output.write(reinterpret_cast<char *>(&(scale_vec[i])), sizeof(double));
   for (int i115 = 1; i115 <= number_of_spheres; i115++) {
     if (iog_vec[i115 - 1] < i115) continue;
-    output.write(reinterpret_cast<char *>(&(nshl_vec[i115 - 1])), sizeof(int));
-    output.write(reinterpret_cast<char *>(&(radii_of_spheres[i115 - 1])), sizeof(double));
-    int nsh = nshl_vec[i115 - 1];
+    output.write(reinterpret_cast<char *>(&(nshl_vec[i115 - 1])), sizeof(int)); // was not from IOG
+    output.write(reinterpret_cast<char *>(&(radii_of_spheres[i115 - 1])), sizeof(double)); // was not from IOG
+    int nsh = nshl_vec[i115 - 1]; // was not from IOG
     if (i115 == 1) nsh += ies;
     for (int nsi = 0; nsi < nsh; nsi++)
-      output.write(reinterpret_cast<char *>(&(rcf[i115 - 1][nsi])), sizeof(double));
+      output.write(reinterpret_cast<char *>(&(rcf[i115 - 1][nsi])), sizeof(double)); // was not from IOG
   }
   for (int jxi468 = 1; jxi468 <= number_of_scales; jxi468++) {
     if (idfc != 0 && jxi468 > 1) continue;
     for (int i162 = 1; i162 <= number_of_spheres; i162++) {
       if (iog_vec[i162 - 1] < i162) continue;
-      int nsh = nshl_vec[i162 - 1];
+      int nsh = nshl_vec[i162 - 1]; // was not from IOG
       int ici = (nsh + 1) / 2; // QUESTION: is integer division really intended here?
       if (i162 == 1) ici = ici + ies;
       for (int i157 = 0; i157 < ici; i157++) {