From 2389e7a06a75ff5247cac4ea6f2a4d593fc60052 Mon Sep 17 00:00:00 2001
From: Giovanni La Mura <giovanni.lamura@inaf.it>
Date: Thu, 11 Jan 2024 17:16:08 +0100
Subject: [PATCH] Implement HDF5 configuration binary output
---
src/cluster/Makefile | 11 +-
src/cluster/cluster.cpp | 5 +-
src/include/Configuration.h | 52 ++++-
src/include/file_io.h | 4 -
src/libnptm/Configuration.cpp | 361 ++++++++++++++++++++++------------
5 files changed, 295 insertions(+), 138 deletions(-)
diff --git a/src/cluster/Makefile b/src/cluster/Makefile
index ae2539bf..ad58bc12 100644
--- a/src/cluster/Makefile
+++ b/src/cluster/Makefile
@@ -3,8 +3,8 @@ FC=gfortran
FCFLAGS=-std=legacy -O3
LFLAGS=
CXX=g++
-CXXFLAGS=-O2 -ggdb -pg -coverage
-CXXLFLAGS=
+CXXFLAGS=-O3 -ggdb -pg -coverage
+CXXLFLAGS=-L/usr/lib64 -lhdf5_hl -lhdf5
all: clu edfb np_cluster
@@ -14,8 +14,8 @@ clu: clu.o
edfb: edfb.o
$(FC) $(FCFLAGS) -o $(BUILDDIR)/edfb $(BUILDDIR)/edfb.o $(LFLAGS)
-np_cluster: $(BUILDDIR)/np_cluster.o $(BUILDDIR)/Commons.o $(BUILDDIR)/Configuration.o $(BUILDDIR)/Parsers.o $(BUILDDIR)/sph_subs.o $(BUILDDIR)/clu_subs.o $(BUILDDIR)/cluster.o
- $(CXX) $(CXXFLAGS) $(CXXLFLAGS) -o $(BUILDDIR)/np_cluster $(BUILDDIR)/np_cluster.o $(BUILDDIR)/Commons.o $(BUILDDIR)/Configuration.o $(BUILDDIR)/Parsers.o $(BUILDDIR)/sph_subs.o $(BUILDDIR)/clu_subs.o $(BUILDDIR)/cluster.o
+np_cluster: $(BUILDDIR)/np_cluster.o $(BUILDDIR)/Commons.o $(BUILDDIR)/Configuration.o $(BUILDDIR)/file_io.o $(BUILDDIR)/Parsers.o $(BUILDDIR)/sph_subs.o $(BUILDDIR)/clu_subs.o $(BUILDDIR)/cluster.o
+ $(CXX) $(CXXFLAGS) $(CXXLFLAGS) -o $(BUILDDIR)/np_cluster $(BUILDDIR)/np_cluster.o $(BUILDDIR)/Commons.o $(BUILDDIR)/Configuration.o $(BUILDDIR)/file_io.o $(BUILDDIR)/Parsers.o $(BUILDDIR)/sph_subs.o $(BUILDDIR)/clu_subs.o $(BUILDDIR)/cluster.o
$(BUILDDIR)/np_cluster.o:
$(CXX) $(CXXFLAGS) -c np_cluster.cpp -o $(BUILDDIR)/np_cluster.o
@@ -26,6 +26,9 @@ $(BUILDDIR)/Commons.o:
$(BUILDDIR)/Configuration.o:
$(CXX) $(CXXFLAGS) -c ../libnptm/Configuration.cpp -o $(BUILDDIR)/Configuration.o
+$(BUILDDIR)/file_io.o:
+ $(CXX) $(CXXFLAGS) -c ../libnptm/file_io.cpp -o $(BUILDDIR)/file_io.o
+
$(BUILDDIR)/Parsers.o:
$(CXX) $(CXXFLAGS) -c ../libnptm/Parsers.cpp -o $(BUILDDIR)/Parsers.o
diff --git a/src/cluster/cluster.cpp b/src/cluster/cluster.cpp
index 559375ff..5c1f96e7 100644
--- a/src/cluster/cluster.cpp
+++ b/src/cluster/cluster.cpp
@@ -1,12 +1,12 @@
/*! \file cluster.cpp
*/
#include <cstdio>
+#include <complex>
+#include <exception>
#include <fstream>
#include <string>
-#include <complex>
#ifndef INCLUDE_CONFIGURATION_H_
-#include <exception>
#include "../include/Configuration.h"
#endif
@@ -42,6 +42,7 @@ void cluster(string config_file, string data_file, string output_path) {
}
sconf->write_formatted(output_path + "/c_OEDFB");
sconf->write_binary(output_path + "/c_TEDF");
+ sconf->write_binary(output_path + "/c_TEDF.hd5", "HDF5");
GeometryConfiguration *gconf = NULL;
try {
gconf = GeometryConfiguration::from_legacy(data_file);
diff --git a/src/include/Configuration.h b/src/include/Configuration.h
index a4cf7f44..1050a61a 100644
--- a/src/include/Configuration.h
+++ b/src/include/Configuration.h
@@ -218,6 +218,51 @@ protected:
//! \brief Flag to control whether to add an external layer.
bool use_external_sphere;
+ /*! \brief Build configuration from a HDF5 binary input file.
+ *
+ * This is the function called by the public method `from_binary()` in case of
+ * HDF5 mode selection. This function creates a configuration structure from
+ * a binary file written according to the HDF5 format standard.
+ *
+ * \param file_name: `string` Name of the binary configuration data file.
+ * \return config: `ScattererConfiguration*` Pointer to object containing the
+ * scatterer configuration data.
+ */
+ static ScattererConfiguration *from_hdf5(std::string file_name);
+
+ /*! \brief Build configuration from legacy binary input file.
+ *
+ * This is the function called by the public method `from_binary()` in case of
+ * legacy mode selection. This function creates a configuration structure from
+ * a binary file written according to the proprietary mode used by the original
+ * FORTRAN code.
+ *
+ * \param file_name: `string` Name of the binary configuration data file.
+ * \return config: `ScattererConfiguration*` Pointer to object containing the
+ * scatterer configuration data.
+ */
+ static ScattererConfiguration *from_legacy(std::string file_name);
+
+ /*! \brief Write the scatterer configuration data to HDF5 binary output.
+ *
+ * This function is invoked by the public method `write_binary()` with the
+ * "HDF5" format mode. It undertakes the task of writing the configuration
+ * information to a binary file using the standard HDF5 format.
+ *
+ * \param file_name: `string` Name of the binary configuration data file.
+ */
+ void write_hdf5(std::string file_name);
+
+ /*! \brief Write the scatterer configuration data to legacy binary output.
+ *
+ * This function is invoked by the public method `write_binary()` with the
+ * "LEGACY" format mode. It undertakes the task of writing the configuration
+ * information to a binary file using a proprietary format, as it was done
+ * originally in the FORTRAN code.
+ *
+ * \param file_name: `string` Name of the binary configuration data file.
+ */
+ void write_legacy(std::string file_name);
public:
/*! \brief Build a scatterer configuration structure.
*
@@ -306,14 +351,13 @@ public:
* be useful to save the configuration data. `ScattererConfiguration.write_binary()`
* performs the operation of saving the configuration in binary format. This function
* can work in legacy mode, to write backward compatible configuration files, as well
- * as by wrapping the data into common scientific formats (NB: this last option still
- * needs to be implemented).
+ * as by wrapping the data into common scientific formats.
*
* \param file_name: `string` Name of the file to be written.
- * \param mode: `string` Binary encoding. Can be one of "LEGACY", ... . Optional
+ * \param mode: `string` Binary encoding. Can be one of ["LEGACY", "HDF5"] . Optional
* (default is "LEGACY").
*/
- void write_binary(std::string file_name, std::string mode = "LEGACY");
+ void write_binary(std::string file_name, std::string mode="LEGACY");
/*! \brief Write the scatterer configuration data to formatted text output.
*
diff --git a/src/include/file_io.h b/src/include/file_io.h
index 9801804a..8c228098 100644
--- a/src/include/file_io.h
+++ b/src/include/file_io.h
@@ -120,10 +120,6 @@ class HDFFile {
*/
bool is_open() { return file_open_flag; }
- /* ! \brief Read data from attached file.
- */
- // herr_t read();
-
/*! \brief Write data to attached file.
*
* \param dataset_name: `string` Name of the dataset to write to.
diff --git a/src/libnptm/Configuration.cpp b/src/libnptm/Configuration.cpp
index 12524367..19a329b3 100644
--- a/src/libnptm/Configuration.cpp
+++ b/src/libnptm/Configuration.cpp
@@ -4,9 +4,11 @@
#include <cmath>
#include <complex>
#include <cstdio>
+#include <exception>
#include <fstream>
#include <regex>
#include <string>
+#include "hdf5.h"
#ifndef INCLUDE_LIST_H_
#include "../include/List.h"
@@ -20,6 +22,10 @@
#include "../include/Configuration.h"
#endif
+#ifndef INCLUDE_FILE_IO_H_
+#include "../include/file_io.h"
+#endif
+
using namespace std;
GeometryConfiguration::GeometryConfiguration(
@@ -203,8 +209,8 @@ ScattererConfiguration::ScattererConfiguration(
radii_of_spheres = ros_vector;
nshl_vec = nshl_vector;
rcf = rcf_vector;
- idfc = dielectric_func_type,
- dc0_matrix = dc_matrix;
+ idfc = dielectric_func_type;
+ dc0_matrix = dc_matrix;
use_external_sphere = is_external;
exdc = ex;
wp = w;
@@ -248,6 +254,24 @@ ScattererConfiguration::~ScattererConfiguration() {
}
ScattererConfiguration* ScattererConfiguration::from_binary(string file_name, string mode) {
+ ScattererConfiguration *conf = NULL;
+ if (mode.compare("LEGACY") == 0) {
+ conf = ScattererConfiguration::from_legacy(file_name);
+ } else if (mode.compare("HDF5") == 0) {
+ conf = ScattererConfiguration::from_hdf5(file_name);
+ } else {
+ string message = "Unknown format mode: \"" + mode + "\"";
+ throw UnrecognizedConfigurationException(message);
+ }
+ return conf;
+}
+
+ScattererConfiguration* ScattererConfiguration::from_hdf5(string file_name) {
+ ScattererConfiguration *conf = NULL;
+ return conf;
+}
+
+ScattererConfiguration* ScattererConfiguration::from_legacy(string file_name) {
int nsph;
int *iog;
double _exdc, _wp, _xip;
@@ -257,80 +281,78 @@ ScattererConfiguration* ScattererConfiguration::from_binary(string file_name, st
double *ros_vector;
double **rcf_vector;
complex<double> ***dc0m;
- if (mode.compare("LEGACY") == 0) { // Legacy mode was chosen.
- int max_ici = 0;
- fstream input;
- input.open(file_name.c_str(), ios::in | ios::binary);
- if (input.is_open()) {
- input.read(reinterpret_cast<char *>(&nsph), sizeof(int));
- iog = new int[nsph]();
- for (int i = 0; i < nsph; i++) {
- input.read(reinterpret_cast<char *>(&(iog[i])), sizeof(int));
+
+ int max_ici = 0;
+ fstream input;
+ input.open(file_name.c_str(), ios::in | ios::binary);
+ if (input.is_open()) {
+ input.read(reinterpret_cast<char *>(&nsph), sizeof(int));
+ iog = new int[nsph]();
+ for (int i = 0; i < nsph; i++) {
+ input.read(reinterpret_cast<char *>(&(iog[i])), sizeof(int));
+ }
+ input.read(reinterpret_cast<char *>(&_exdc), sizeof(double));
+ input.read(reinterpret_cast<char *>(&_wp), sizeof(double));
+ input.read(reinterpret_cast<char *>(&_xip), sizeof(double));
+ input.read(reinterpret_cast<char *>(&_idfc), sizeof(int));
+ input.read(reinterpret_cast<char *>(&nxi), sizeof(int));
+ try {
+ xi_vec = new double[nxi]();
+ } catch (const bad_alloc &ex) {
+ throw UnrecognizedConfigurationException("Wrong parameter set: invalid number of scales " + to_string(nxi));
+ }
+ 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];
+ 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 *>(&_exdc), sizeof(double));
- input.read(reinterpret_cast<char *>(&_wp), sizeof(double));
- input.read(reinterpret_cast<char *>(&_xip), sizeof(double));
- input.read(reinterpret_cast<char *>(&_idfc), sizeof(int));
- input.read(reinterpret_cast<char *>(&nxi), sizeof(int));
+ 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];
+ if (max_ici < (nsh + 1) / 2) max_ici = (nsh + 1) / 2;
try {
- xi_vec = new double[nxi]();
+ rcf_vector[i115 - 1] = new double[nsh]();
} catch (const bad_alloc &ex) {
- throw UnrecognizedConfigurationException("Wrong parameter set: invalid number of scales " + to_string(nxi));
+ throw UnrecognizedConfigurationException("Wrong parameter set: invalid number of layers " + to_string(nsh));
}
- for (int i = 0; i < nxi; i++) {
- input.read(reinterpret_cast<char *>(&(xi_vec[i])), sizeof(double));
+ for (int nsi = 0; nsi < nsh; nsi++) {
+ input.read(reinterpret_cast<char *>(&(rcf_vector[i115 - 1][nsi])), sizeof(double));
}
- nshl_vector = new int[nsph]();
- ros_vector = new double[nsph]();
- rcf_vector = new double*[nsph];
- 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];
- if (max_ici < (nsh + 1) / 2) max_ici = (nsh + 1) / 2;
- try {
- rcf_vector[i115 - 1] = new double[nsh]();
- } 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));
- }
- }
- dc0m = new complex<double>**[max_ici];
- for (int dim1 = 0; dim1 < max_ici; dim1++) {
- dc0m[dim1] = new complex<double>*[nsph];
- for (int dim2 = 0; dim2 < nsph; dim2++) {
- dc0m[dim1][dim2] = new complex<double>[nxi]();
- }
+ }
+ dc0m = new complex<double>**[max_ici];
+ int dim3 = (_idfc == 0) ? nxi : 1;
+ for (int dim1 = 0; dim1 < max_ici; dim1++) {
+ dc0m[dim1] = new complex<double>*[nsph];
+ for (int dim2 = 0; dim2 < nsph; dim2++) {
+ dc0m[dim1][dim2] = new complex<double>[dim3]();
}
- for (int jxi468 = 1; jxi468 <= nxi; jxi468++) {
- 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 ici = (nsh + 1) / 2; // QUESTION: is integer division really intended here?
- for (int i157 = 0; i157 < ici; i157++) {
- double dc0_real, dc0_img;
- input.read(reinterpret_cast<char *>(&dc0_real), sizeof(double));
- input.read(reinterpret_cast<char *>(&dc0_img), sizeof(double));
- dc0m[i157][i162 - 1][jxi468 - 1] = dc0_real + 1i * dc0_img;
- }
+ }
+ for (int jxi468 = 1; jxi468 <= nxi; jxi468++) {
+ 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 ici = (nsh + 1) / 2; // QUESTION: is integer division really intended here?
+ for (int i157 = 0; i157 < ici; i157++) {
+ double dc0_real, dc0_img;
+ input.read(reinterpret_cast<char *>(&dc0_real), sizeof(double));
+ input.read(reinterpret_cast<char *>(&dc0_img), sizeof(double));
+ dc0m[i157][i162 - 1][jxi468 - 1] = dc0_real + 1i * dc0_img;
}
}
- input.close();
- } else { // Opening of the input file did not succeed.
- OpenConfigurationFileException ex(file_name);
- throw ex;
}
- } else { // A different binary format was chosen.
- //TODO: this part is not yet implemented.
- // Functions to write optimized file formats may be invoked here.
+ input.close();
+ } else { // Opening of the input file did not succeed.
+ OpenConfigurationFileException ex(file_name);
+ throw ex;
}
+
ScattererConfiguration *conf = new ScattererConfiguration(
nsph,
xi_vec,
@@ -511,20 +533,10 @@ ScattererConfiguration* ScattererConfiguration::from_dedfb(string dedfb_file_nam
}
}
}
- //last_read_line++;
int *iog_vector = new int[nsph]();
double *ros_vector = new double[nsph]();
double **rcf_vector = new double*[nsph];
int *nshl_vector = new int[nsph]();
- /*for (int i = 0; i < nsph; i++) {
- string read_format = "";
- for (int j = 0; j < (i % 15); j++) read_format += " %*d";
- read_format += " %d";
- sscanf(file_lines[last_read_line].c_str(), read_format.c_str(), (iog_vector + i));
- if (i > 0 && i % 15 == 0) {
- last_read_line++;
- }
- }*/
int filled_iogs = 0;
re = regex("[0-9]+");
while (filled_iogs < nsph) {
@@ -561,8 +573,6 @@ ScattererConfiguration* ScattererConfiguration::from_dedfb(string dedfb_file_nam
rcf_vector[i113 - 1] = new double[nsh]();
for (int ns = 0; ns < nsh; ns++) {
double ns_rcf;
- //int ns_rcf_exp;
- //sscanf(file_lines[++last_read_line].c_str(), " %lf D%d", &ns_rcf, &ns_rcf_exp);
str_target = file_lines[++last_read_line];
regex_search(str_target, m, re);
str_number = m.str();
@@ -573,10 +583,11 @@ ScattererConfiguration* ScattererConfiguration::from_dedfb(string dedfb_file_nam
}
}
complex<double> ***dc0m = new complex<double>**[max_ici];
+ int dim3 = (_idfc == 0) ? nxi : 1;
for (int dim1 = 0; dim1 < max_ici; dim1++) {
dc0m[dim1] = new complex<double>*[nsph];
for (int dim2 = 0; dim2 < nsph; dim2++) {
- dc0m[dim1][dim2] = new complex<double>[nxi]();
+ dc0m[dim1][dim2] = new complex<double>[dim3]();
}
}
re = regex("-?[0-9]+\\.[0-9]+([eEdD][-+]?[0-9]+)?");
@@ -589,8 +600,6 @@ ScattererConfiguration* ScattererConfiguration::from_dedfb(string dedfb_file_nam
if (i162 == 1) ici = ici + ies;
for (int i157 = 0; i157 < ici; i157++) {
double dc0_real, dc0_img;
- //int dc0_real_exp, dc0_img_exp;
- //sscanf(file_lines[++last_read_line].c_str(), " (%lf D%d, %lf D%d)", &dc0_real, &dc0_real_exp, &dc0_img, &dc0_img_exp);
str_target = file_lines[++last_read_line];
regex_search(str_target, m, re);
string str_number = m.str();
@@ -603,7 +612,7 @@ ScattererConfiguration* ScattererConfiguration::from_dedfb(string dedfb_file_nam
str_number = regex_replace(str_number, regex("D"), "e");
str_number = regex_replace(str_number, regex("d"), "e");
dc0_img = stod(str_number);
- dc0m[i157][i162 - 1][jxi468 - 1] = std::complex<double>(dc0_real, dc0_img);
+ dc0m[i157][i162 - 1][jxi468 - 1] = complex<double>(dc0_real, dc0_img);
}
}
}
@@ -676,55 +685,159 @@ void ScattererConfiguration::print() {
}
void ScattererConfiguration::write_binary(string file_name, string mode) {
+ if (mode.compare("LEGACY") == 0) {
+ write_legacy(file_name);
+ } else if (mode.compare("HDF5") == 0) {
+ write_hdf5(file_name);
+ } else {
+ string message = "Unknown format mode: \"" + mode + "\"";
+ throw UnrecognizedConfigurationException(message);
+ }
+}
+
+void ScattererConfiguration::write_hdf5(string file_name) {
const double two_pi = acos(0.0) * 4.0;
const double evc = 6.5821188e-16;
+ int ies = (use_external_sphere)? 1 : 0;
int max_ici = 0;
- bool is_new_vector = false;
- if (mode.compare("LEGACY") == 0) { // Legacy mode was chosen.
- fstream output;
- int ies = (use_external_sphere)? 1 : 0;
- output.open(file_name.c_str(), ios::out | ios::binary);
- output.write(reinterpret_cast<char *>(&number_of_spheres), sizeof(int));
- for (int i = 0; i < number_of_spheres; i++)
- output.write(reinterpret_cast<char *>(&(iog_vec[i])), sizeof(int));
- output.write(reinterpret_cast<char *>(&exdc), sizeof(double));
- output.write(reinterpret_cast<char *>(&wp), sizeof(double));
- output.write(reinterpret_cast<char *>(&xip), sizeof(double));
- output.write(reinterpret_cast<char *>(&idfc), sizeof(int));
- output.write(reinterpret_cast<char *>(&number_of_scales), sizeof(int));
- for (int i = 0; i < number_of_scales; i++)
- 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];
- if (i115 == 1) nsh += ies;
- if (max_ici < (nsh + 1) / 2) max_ici = (nsh + 1) / 2;
- for (int nsi = 0; nsi < nsh; nsi++)
- output.write(reinterpret_cast<char *>(&(rcf[i115 - 1][nsi])), sizeof(double));
+ List<string> rec_name_list(1);
+ List<string> rec_type_list(1);
+ List<void *> rec_ptr_list(1);
+ string str_type, str_name;
+ rec_name_list.set(0, "NSPH");
+ rec_type_list.set(0, "INT32_(1)");
+ rec_ptr_list.set(0, &number_of_spheres);
+ rec_name_list.append("IOGVEC");
+ str_type = "INT32_(" + to_string(number_of_spheres) + ")";
+ rec_type_list.append(str_type);
+ rec_ptr_list.append(iog_vec);
+ rec_name_list.append("EXDC");
+ rec_type_list.append("FLOAT64_(1)");
+ rec_ptr_list.append(&exdc);
+ rec_name_list.append("WP");
+ rec_type_list.append("FLOAT64_(1)");
+ rec_ptr_list.append(&wp);
+ rec_name_list.append("XIP");
+ rec_type_list.append("FLOAT64_(1)");
+ rec_ptr_list.append(&xip);
+ rec_name_list.append("IDFC");
+ rec_type_list.append("INT32_(1)");
+ rec_ptr_list.append(&idfc);
+ rec_name_list.append("NXI");
+ rec_type_list.append("INT32_(1)");
+ rec_ptr_list.append(&number_of_scales);
+ rec_name_list.append("XIVEC");
+ str_type = "FLOAT64_(" + to_string(number_of_scales) + ")";
+ rec_type_list.append(str_type);
+ rec_ptr_list.append(scale_vec);
+ for (int i115 = 1; i115 <= number_of_spheres; i115++) {
+ if (iog_vec[i115 - 1] < i115) continue;
+ 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]));
+ 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];
+ if (i115 == 1) nsh += ies;
+ if (max_ici < (nsh + 1) / 2) max_ici = nsh + 1 / 2;
+ str_name = "RCF_" + to_string(i115);
+ 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]));
+ }
+
+ int dim3 = (idfc == 0) ? number_of_scales : 1;
+ int dc0m_size = 2 * dim3 * number_of_spheres * max_ici;
+ double *dc0m = new double[dc0m_size];
+ int dc0_index = 0;
+ 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 ici = (nsh + 1) / 2;
+ if (i162 == 1) ici = ici + ies;
+ for (int i157 = 0; i157 < ici; i157++) {
+ double dc0_real, dc0_imag;
+ dc0_real = dc0_matrix[i157][i162 - 1][jxi468 - 1].real();
+ dc0_imag = dc0_matrix[i157][i162 - 1][jxi468 - 1].imag();
+ dc0m[dc0_index++] = dc0_real;
+ dc0m[dc0_index++] = dc0_imag;
+ }
}
- 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 ici = (nsh + 1) / 2; // QUESTION: is integer division really intended here?
- if (i162 == 1) ici = ici + ies;
- for (int i157 = 0; i157 < ici; i157++) {
- double dc0_real, dc0_img;
- dc0_real = dc0_matrix[i157][i162 - 1][jxi468 - 1].real();
- dc0_img = dc0_matrix[i157][i162 - 1][jxi468 - 1].imag();
- // The FORTRAN code writes the complex numbers as a 16-byte long binary stream.
- // Here we assume that the 16 bytes are equally split in 8 bytes to represent the
- // real part and 8 bytes to represent the imaginary one.
- output.write(reinterpret_cast<char *>(&dc0_real), sizeof(double));
- output.write(reinterpret_cast<char *>(&dc0_img), sizeof(double));
- }
+ }
+ str_type = "FLOAT64_(" + to_string(dc0m_size) + ")";
+ rec_name_list.append("DC0M");
+ rec_type_list.append(str_type);
+ rec_ptr_list.append(dc0m);
+
+ string *rec_names = rec_name_list.to_array();
+ string *rec_types = rec_type_list.to_array();
+ void **rec_pointers = rec_ptr_list.to_array();
+ const int rec_num = rec_name_list.length();
+ FileSchema schema(rec_num, rec_types, rec_names);
+ HDFFile *hdf_file = HDFFile::from_schema(schema, "c_TEDF.hd5", H5F_ACC_TRUNC);
+ for (int ri = 0; ri < rec_num; ri++)
+ hdf_file->write(rec_names[ri], rec_types[ri], rec_pointers[ri]);
+ hdf_file->close();
+
+ // Clean memory
+ delete[] dc0m;
+ delete[] rec_names;
+ delete[] rec_types;
+ delete[] rec_pointers;
+ delete hdf_file;
+}
+
+void ScattererConfiguration::write_legacy(string file_name) {
+ const double two_pi = acos(0.0) * 4.0;
+ const double evc = 6.5821188e-16;
+ fstream output;
+ int ies = (use_external_sphere)? 1 : 0;
+ output.open(file_name.c_str(), ios::out | ios::binary);
+ output.write(reinterpret_cast<char *>(&number_of_spheres), sizeof(int));
+ for (int i = 0; i < number_of_spheres; i++)
+ output.write(reinterpret_cast<char *>(&(iog_vec[i])), sizeof(int));
+ output.write(reinterpret_cast<char *>(&exdc), sizeof(double));
+ output.write(reinterpret_cast<char *>(&wp), sizeof(double));
+ output.write(reinterpret_cast<char *>(&xip), sizeof(double));
+ output.write(reinterpret_cast<char *>(&idfc), sizeof(int));
+ output.write(reinterpret_cast<char *>(&number_of_scales), sizeof(int));
+ for (int i = 0; i < number_of_scales; i++)
+ 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];
+ if (i115 == 1) nsh += ies;
+ for (int nsi = 0; nsi < nsh; nsi++)
+ output.write(reinterpret_cast<char *>(&(rcf[i115 - 1][nsi])), sizeof(double));
+ }
+ 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 ici = (nsh + 1) / 2; // QUESTION: is integer division really intended here?
+ if (i162 == 1) ici = ici + ies;
+ for (int i157 = 0; i157 < ici; i157++) {
+ double dc0_real, dc0_img;
+ dc0_real = dc0_matrix[i157][i162 - 1][jxi468 - 1].real();
+ dc0_img = dc0_matrix[i157][i162 - 1][jxi468 - 1].imag();
+ // The FORTRAN code writes the complex numbers as a 16-byte long binary stream.
+ // Here we assume that the 16 bytes are equally split in 8 bytes to represent the
+ // real part and 8 bytes to represent the imaginary one.
+ output.write(reinterpret_cast<char *>(&dc0_real), sizeof(double));
+ output.write(reinterpret_cast<char *>(&dc0_img), sizeof(double));
}
}
- output.close();
}
+ output.close();
}
void ScattererConfiguration::write_formatted(string file_name) {
--
GitLab