Commit ceb70c31 authored by Giovanni La Mura's avatar Giovanni La Mura
Browse files

Set up ScattererConfiguration HDF5 output

parent 2f098737

src/include/Configuration.h

+9 −2
Original line number Diff line number Diff line
@@ -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;
@@ -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

src/libnptm/Configuration.cpp

+198 −115
Original line number Diff line number Diff line
@@ -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;
@@ -639,6 +533,161 @@ ScattererConfiguration* ScattererConfiguration::from_dedfb(string dedfb_file_nam 
  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;
@@ -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;

}