From ceb70c31cf11d85ec7f6c563e8473376b1489126 Mon Sep 17 00:00:00 2001
From: Giovanni La Mura <giovanni.lamura@inaf.it>
Date: Fri, 2 Feb 2024 17:16:47 +0100
Subject: [PATCH] Set up ScattererConfiguration HDF5 output
---
src/include/Configuration.h | 11 +-
src/libnptm/Configuration.cpp | 313 +++++++++++++++++++++-------------
2 files changed, 207 insertions(+), 117 deletions(-)
diff --git a/src/include/Configuration.h b/src/include/Configuration.h
index bdf3e170..208a69ba 100644
--- a/src/include/Configuration.h
+++ b/src/include/Configuration.h
@@ -213,7 +213,7 @@ class ScattererConfiguration {
friend void cluster(std::string, std::string, std::string);
friend void sphere(std::string, std::string, std::string);
protected:
- //! \brief Matrix of dielectric parameters with size [NON_TRANS_LAYERS x N_SPHERES x LAYERS].
+ //! \brief Matrix of dielectric parameters with size [NON_TRANS_LAYERS x N_SPHERES x N_SCALES].
std::complex<double> ***dc0_matrix;
//! \brief Vector of sphere radii expressed in m, with size [N_SPHERES].
double *radii_of_spheres;
@@ -325,7 +325,7 @@ public:
double exdc,
double wp,
double xip
- );
+ );
/*! \brief Destroy a scatterer configuration instance.
*/
@@ -394,6 +394,13 @@ public:
* \param file_name: `string` Name of the file to be written.
*/
void write_formatted(std::string file_name);
+
+ /*! \brief Test whether two instances of ScattererConfiguration are equal.
+ *
+ * \param other: `ScattererConfiguration &` Reference to the instance to be compared with.
+ * \return result: `bool` True, if the two instances are equal, false otherwise.
+ */
+ bool operator ==(ScattererConfiguration &other);
};
#endif
diff --git a/src/libnptm/Configuration.cpp b/src/libnptm/Configuration.cpp
index be5f3919..9ccb9ad8 100644
--- a/src/libnptm/Configuration.cpp
+++ b/src/libnptm/Configuration.cpp
@@ -39,7 +39,7 @@ GeometryConfiguration::GeometryConfiguration(
double in_ph_start, double in_ph_step, double in_ph_end,
double sc_ph_start, double sc_ph_step, double sc_ph_end,
int _jwtm
- ) {
+) {
number_of_spheres = _nsph;
l_max = _lm;
in_pol = _in_pol;
@@ -184,7 +184,7 @@ GeometryConfiguration* GeometryConfiguration::from_legacy(string file_name) {
in_ph_start, in_ph_step, in_ph_end,
sc_ph_start, sc_ph_step, sc_ph_end,
_jwtm
- );
+ );
return conf;
}
@@ -203,7 +203,7 @@ ScattererConfiguration::ScattererConfiguration(
double ex,
double w,
double x
- ) {
+) {
number_of_spheres = nsph;
number_of_scales = nxi;
reference_variable_name = variable_name;
@@ -268,112 +268,6 @@ ScattererConfiguration* ScattererConfiguration::from_binary(string file_name, st
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;
- int _idfc, nxi;
- int *nshl_vector;
- double *xi_vec;
- double *ros_vector;
- double **rcf_vector;
- complex<double> ***dc0m;
-
- 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 *>(&(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];
- 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;
- }
- }
- }
- input.close();
- } else { // Opening of the input file did not succeed.
- OpenConfigurationFileException ex(file_name);
- throw ex;
- }
-
- ScattererConfiguration *conf = new ScattererConfiguration(
- nsph,
- xi_vec,
- nxi,
- "XIV",
- iog,
- ros_vector,
- nshl_vector,
- rcf_vector,
- _idfc,
- dc0m,
- false,
- _exdc,
- _wp,
- _xip
- );
- return conf;
-}
-
ScattererConfiguration* ScattererConfiguration::from_dedfb(string dedfb_file_name) {
int num_lines = 0;
int last_read_line = 0;
@@ -634,11 +528,166 @@ ScattererConfiguration* ScattererConfiguration::from_dedfb(string dedfb_file_nam
_exdc,
_wp,
_xip
- );
+ );
delete[] file_lines;
return config;
}
+ScattererConfiguration* ScattererConfiguration::from_hdf5(string file_name) {
+ ScattererConfiguration *conf = NULL;
+ unsigned int flags = H5F_ACC_RDONLY;
+ HDFFile *hdf_file = new HDFFile(file_name, flags);
+ herr_t status = hdf_file->get_status();
+ if (status == 0) {
+ int nsph;
+ int *iog;
+ double _exdc, _wp, _xip;
+ int _idfc, nxi;
+ int *nshl_vector;
+ double *xi_vec;
+ double *ros_vector;
+ double **rcf_vector;
+ complex<double> ***dc0m;
+ status = hdf_file->read("NSPH", "INT32_(1)", &nsph);
+ status = hdf_file->read("EXDC", "FLOAT64_(1)", &_exdc);
+ status = hdf_file->read("WP", "FLOAT64_(1)", &_wp);
+ status = hdf_file->read("XIP", "FLOAT64_(1)", &_xip);
+ status = hdf_file->read("IDFC", "INT32_(1)", &_idfc);
+ status = hdf_file->read("NXI", "INT32_(1)", &nxi);
+ iog = new int[nsph];
+ string str_type = "INT32_(" + to_string(nsph) + ")";
+ status = hdf_file.read("IOGVEC", str_type, iog);
+ int configuration_count = 0;
+ for (int ci = 0; ci < nsph; ci++) {
+ if (iog[ci] < ci + 1) continue;
+ configuration_count++;
+ }
+ nshl_vector = new int[configuration_count];
+ ros_vector = new double[configuration_count];
+
+ status = hdf_file->close();
+ conf = new ScattererConfiguration(
+ nsph,
+ xi_vec,
+ nxi,
+ "XIV",
+ iog,
+ ros_vector,
+ nshl_vector,
+ rcf_vector,
+ _idfc,
+ dc0m,
+ false,
+ _exdc,
+ _wp,
+ _xip
+ );
+ }
+
+ return conf;
+}
+
+ScattererConfiguration* ScattererConfiguration::from_legacy(string file_name) {
+ int nsph;
+ int *iog;
+ double _exdc, _wp, _xip;
+ int _idfc, nxi;
+ int *nshl_vector;
+ double *xi_vec;
+ double *ros_vector;
+ double **rcf_vector;
+ complex<double> ***dc0m;
+
+ 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 *>(&(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];
+ 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;
+ }
+ }
+ }
+ input.close();
+ } else { // Opening of the input file did not succeed.
+ OpenConfigurationFileException ex(file_name);
+ throw ex;
+ }
+
+ ScattererConfiguration *conf = new ScattererConfiguration(
+ nsph,
+ xi_vec,
+ nxi,
+ "XIV",
+ iog,
+ ros_vector,
+ nshl_vector,
+ rcf_vector,
+ _idfc,
+ dc0m,
+ false,
+ _exdc,
+ _wp,
+ _xip
+ );
+ return conf;
+}
+
void ScattererConfiguration::print() {
int ies = (use_external_sphere)? 1 : 0;
int max_ici = 0;
@@ -880,7 +929,7 @@ void ScattererConfiguration::write_formatted(string file_name) {
wl_vec[i],
pu_vec[i],
ev_vec[i]
- );
+ );
}
break;
case 1:
@@ -900,7 +949,7 @@ void ScattererConfiguration::write_formatted(string file_name) {
pu_vec[i],
ev_vec[i],
xi_vec[i]
- );
+ );
}
break;
case 2:
@@ -920,7 +969,7 @@ void ScattererConfiguration::write_formatted(string file_name) {
pu_vec[i],
ev_vec[i],
xi_vec[i]
- );
+ );
}
break;
case 3:
@@ -940,7 +989,7 @@ void ScattererConfiguration::write_formatted(string file_name) {
wl_vec[i],
ev_vec[i],
xi_vec[i]
- );
+ );
}
break;
case 4:
@@ -960,7 +1009,7 @@ void ScattererConfiguration::write_formatted(string file_name) {
wl_vec[i],
pu_vec[i],
xi_vec[i]
- );
+ );
}
break;
default:
@@ -1016,3 +1065,37 @@ void ScattererConfiguration::write_formatted(string file_name) {
}
fclose(output);
}
+
+bool ScattererConfiguration::operator ==(ScattererConfiguration &other) {
+ if (reference_variable_name.compare(other.reference_variable_name) != 0) return false;
+ if (number_of_spheres != other.number_of_spheres) return false;
+ if (number_of_scales != other.number_of_scales) return false;
+ if (idfc != other.idfc) return false;
+ if (exdc != other.exdc) return false;
+ if (wp != other.wp) return false;
+ if (xip != other.xip) return false;
+ if (use_external_sphere != other.use_external_sphere) return false;
+ for (int svi = 0; svi < number_of_scales; svi++)
+ if (scale_vec[svi] != other.scale_vec[svi]) return false;
+ int dj_index = 0;
+ for (int dj = 0; dj < number_of_spheres; dj++) {
+ bool check matrixes = false;
+ if (iog_vec[dj] >= dj + 1) {
+ dj_index = iog_vec[dj] - 1;
+ check_matrixes = true;
+ }
+ if (iog_vec[dj] != other.iog_vec[dj]) return false;
+ if (radii_of_spheres[dj_index] != other.radii_of_spheres[dj_index]) return false;
+ int layers = nshl_vec[dj_index];
+ if (layers != other.nshl_vec[dj_index]) return false;
+ if (check_matrixes) {
+ for (int di = 0; di < layers; di++) {
+ if (rcf[dj_index][di] != other.rcf[dj_index][di]) return false;
+ for (int dk = 0; dk < number_of_scales; dk++) {
+ if (dc0_matrix[di][dj_index][dk] != other.dc0_matrix[di][dj_index][dk]) return false;
+ } // dk loop
+ } // di loop
+ }
+ } // dj loop
+ return true;
+}
--
GitLab