From b4659d20d6a50bb871c57e72cd891eb291d454c1 Mon Sep 17 00:00:00 2001
From: Giovanni La Mura <giovanni.lamura@inaf.it>
Date: Tue, 13 Feb 2024 15:31:33 +0100
Subject: [PATCH] Correct vector to matrix mapping when reading configuration
data from HDF5
---
src/libnptm/Configuration.cpp | 24 +++++++++++++++---------
1 file changed, 15 insertions(+), 9 deletions(-)
diff --git a/src/libnptm/Configuration.cpp b/src/libnptm/Configuration.cpp
index 7eb58c4c..572b2d12 100644
--- a/src/libnptm/Configuration.cpp
+++ b/src/libnptm/Configuration.cpp
@@ -475,9 +475,10 @@ ScattererConfiguration* ScattererConfiguration::from_dedfb(string dedfb_file_nam
complex<double> ***dc0m = new complex<double>**[configurations];
complex<double> *dc0 = new complex<double>[configurations]();
+ int dim3 = (_idfc == 0) ? nxi : 1;
for (int di = 0; di < configurations; di++) {
dc0m[di] = new complex<double>*[nsph];
- for (int dj = 0; dj < nsph; dj++) dc0m[di][dj] = new complex<double>[nxi]();
+ for (int dj = 0; dj < nsph; dj++) dc0m[di][dj] = new complex<double>[dim3]();
} // di loop
for (int jxi468 = 1; jxi468 <= nxi; jxi468++) {
if (_idfc != 0 && jxi468 > 1) continue; // jxi468 loop
@@ -574,7 +575,7 @@ ScattererConfiguration* ScattererConfiguration::from_hdf5(string file_name) {
str_type = "FLOAT64_(" + to_string(nsh) + ")";
status = hdf_file->read(str_name, str_type, (rcf_vector[i115 - 1]));
}
- xi_vec = new double[nxi];
+ xi_vec = new double[nxi]();
str_name = "XIVEC";
str_type = "FLOAT64_(" + to_string(nxi) + ")";
status = hdf_file->read(str_name, str_type, xi_vec);
@@ -586,13 +587,15 @@ ScattererConfiguration* ScattererConfiguration::from_hdf5(string file_name) {
str_type = "FLOAT64_(" + to_string(element_size) + ")";
status = hdf_file->read(str_name, str_type, elements);
dc0m = new complex<double>**[configuration_count];
+ int dc_index = 0;
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];
+ double rval = elements[2 * dc_index];
+ double ival = elements[2 * dc_index + 1];
+ dc_index++;
dc0m[di][dj][dk] = complex<double>(rval, ival);
}
} // dj loop
@@ -665,9 +668,10 @@ ScattererConfiguration* ScattererConfiguration::from_legacy(string file_name) {
input.read(reinterpret_cast<char *>(&(_rcf_vec[i115 - 1][nsi])), sizeof(double));
}
_dc0m = new complex<double>**[configurations];
+ int dim3 = (_idfc == 0) ? _nxi : 1;
for (int di = 0; di < configurations; di++) {
_dc0m[di] = new complex<double>*[_nsph];
- for (int dj = 0; dj < _nsph; dj++) _dc0m[di][dj] = new complex<double>[_nxi]();
+ for (int dj = 0; dj < _nsph; dj++) _dc0m[di][dj] = new complex<double>[dim3]();
} // di loop
for (int jxi468 = 1; jxi468 <= _nxi; jxi468++) {
if (_idfc != 0 && jxi468 > 1) continue;
@@ -822,7 +826,7 @@ void ScattererConfiguration::write_hdf5(string file_name) {
int dim3 = (idfc == 0) ? number_of_scales : 1;
int dc0m_size = 2 * dim3 * number_of_spheres * configurations;
- double *dc0m = new double[dc0m_size];
+ 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;
@@ -835,8 +839,9 @@ void ScattererConfiguration::write_hdf5(string file_name) {
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;
+ dc0m[2 * dc0_index] = dc0_real;
+ dc0m[2 * dc0_index + 1] = dc0_imag;
+ dc0_index++;
}
}
}
@@ -1109,6 +1114,7 @@ bool ScattererConfiguration::operator ==(ScattererConfiguration &other) {
return false;
}
int dj_index = 0;
+ int dim3 = (idfc == 0) ? number_of_scales : 1;
for (int dj = 0; dj < number_of_spheres; dj++) {
bool check_matrixes = false;
if (iog_vec[dj] >= dj + 1) {
@@ -1130,7 +1136,7 @@ bool ScattererConfiguration::operator ==(ScattererConfiguration &other) {
if (rcf[dj_index][di] != other.rcf[dj_index][di]) {
return false;
}
- for (int dk = 0; dk < number_of_scales; dk++) {
+ for (int dk = 0; dk < dim3; dk++) {
if (dc0_matrix[di][dj_index][dk] != other.dc0_matrix[di][dj_index][dk]) {
return false;
}
--
GitLab