Implement trapping configuration binary I/O and add it to library Makefile

src/include/tfrfme.h

@@ -6,6 +6,50 @@
 #ifndef INCLUDE_TFRFME_H_
 #define INCLUDE_TFRFME_H_
 
+/**
+ * \brief Exception for unrecognized parameters.
+ */
+class MatrixOutOfBoundsException: public std::exception {
+protected:
+  //! Description of the problem.
+  std::string message;
+public:
+  /**
+   * \brief Exception instance constructor.
+   *
+   * \param problem: `string` Description of the problem that occurred.
+   */
+  MatrixOutOfBoundsException(std::string problem) { message = problem; }
+  /**
+   * \brief Exception message.
+   */
+  virtual const char* what() const throw() {
+    return message.c_str();
+  }
+};
+
+/**
+ * \brief Exception for unrecognized parameters.
+ */
+class UnrecognizedParameterException: public std::exception {
+protected:
+  //! Description of the problem.
+  std::string message;
+public:
+  /**
+   * \brief Exception instance constructor.
+   *
+   * \param problem: `string` Description of the problem that occurred.
+   */
+  UnrecognizedParameterException(std::string problem) { message = problem; }
+  /**
+   * \brief Exception message.
+   */
+  virtual const char* what() const throw() {
+    return message.c_str();
+  }
+};
+
 /*! \brief Class to represent the trapping configuration.
  */
 class TFRFME {
@@ -51,28 +95,56 @@ class TFRFME {
   //! Vector of computed z positions
   double *zv;
   //! QUESTION: definition?
-  complex<double> **wsum;
+  std::complex<double> **wsum;
+
+  /*! \brief Load a configuration instance from a HDF5 binary file.
+   *
+   * \param file_name: `string` Name of the file to be loaded.
+   * \return instance: `TFRFME *` Pointer to a new trapping configuration
+   * instance.
+   */
+  static TFRFME *from_hdf5(std::string file_name);
+
+  /*! \brief Load a configuration instance from a legacy binary file.
+   *
+   * \param file_name: `string` Name of the file to be loaded.
+   * \return instance: `TFRFME *` Pointer to a new trapping configuration
+   * instance.
+   */
+  static TFRFME *from_legacy(std::string file_name);
+
+  /*! \brief Save a configuration instance to a HDF5 binary file.
+   *
+   * \param file_name: `string` Name of the file to be loaded.
+   */
+  void write_hdf5(std::string file_name);
+
+  /*! \brief Save a configuration instance to a legacy binary file.
+   *
+   * \param file_name: `string` Name of the file to be loaded.
+   */
+  void write_legacy(std::string file_name);
 
  public:
-  /*! \briesf Trapping configuration instance constructor.
+  /*! \brief Trapping configuration instance constructor.
    */
   TFRFME(
 	 int _lmode, int _lm, int _nkv, int _nxv, int _nyv, int _nzv
 );
   
-  /*! \briesf Trapping configuration instance destroyer.
+  /*! \brief Trapping configuration instance destroyer.
    */
   ~TFRFME();
 
-  /*! \briesf Load a trapping configuration instance from binary file.
+  /*! \brief Load a trapping configuration instance from binary file.
    *
    * \param file_name: `string` Name of the file.
-   * \param format: `string` Format of the file (can be either "HDF5"
+   * \param mode: `string` Format of the file (can be either "HDF5"
    * or "LGEACY". Default is "LEGACY").
    * \return instance: `TFRFME *` Pointer to anewly created configuration
    * instance.
    */
-  TFRFME* from_binary(std::string file_name, std::string format="LEGACY");
+  static TFRFME* from_binary(std::string file_name, std::string mode="LEGACY");
 
   /*! \brief Get an element from the WSUM matrix.
    *
@@ -89,6 +161,27 @@ class TFRFME {
    */
   double get_param(std::string param_name);
 
+  /*! \brief Get the X coordinate of the computed point at given index.
+   *
+   * \param index: `int` Index of the requested point.
+   * \return x: `double` X coordinate of the requested point.
+   */
+  double get_x(int index) { return xv[index]; }
+
+  /*! \brief Get the Y coordinate of the computed point at given index.
+   *
+   * \param index: `int` Index of the requested point.
+   * \return y: `double` Y coordinate of the requested point.
+   */
+  double get_y(int index) { return yv[index]; }
+
+  /*! \brief Get the Z coordinate of the computed point at given index.
+   *
+   * \param index: `int` Index of the requested point.
+   * \return z: `double` Z coordinate of the requested point.
+   */
+  double get_z(int index) { return zv[index]; }
+
   /*! \brief Set an element in the WSUM matrix.
    *
    * \param row: `int` Row index of the element to be read.
@@ -104,12 +197,33 @@ class TFRFME {
    */
   void set_param(std::string param_name, double value);
 
-  /*! \briesf Write a trapping configuration instance to binary file.
+  /*! \brief Set the X coordinate of the computed point at given index.
+   *
+   * \param index: `int` Index of the requested point.
+   * \param value: `double` X coordinate of the requested point.
+   */
+  void set_x(int index, double value) { xv[index] = value; }
+
+  /*! \brief Set the Y coordinate of the computed point at given index.
+   *
+   * \param index: `int` Index of the requested point.
+   * \param value: `double` Y coordinate of the requested point.
+   */
+  void set_y(int index, double value) { yv[index] = value; }
+
+  /*! \brief Set the Z coordinate of the computed point at given index.
+   *
+   * \param index: `int` Index of the requested point.
+   * \param value: `double` Z coordinate of the requested point.
+   */
+  void set_z(int index, double value) { zv[index] = value; }
+
+  /*! \brief Write a trapping configuration instance to binary file.
    *
    * \param file_name: `string` Name of the file.
-   * \param format: `string` Format of the file (can be either "HDF5"
+   * \param mode: `string` Format of the file (can be either "HDF5"
    * or "LGEACY". Default is "LEGACY").
    */
-  void write_binary(std::string file_name, std::string format="LEGACY");
+  void write_binary(std::string file_name, std::string mode="LEGACY");
 };
 #endif

src/libnptm/Makefile

@@ -19,9 +19,9 @@ endif
 include ../make.inc
 
 
-CXX_NPTM_OBJS=$(OBJDIR)/Commons.o $(OBJDIR)/Configuration.o $(OBJDIR)/file_io.o $(OBJDIR)/Parsers.o $(OBJDIR)/sph_subs.o $(OBJDIR)/clu_subs.o $(OBJDIR)/tra_subs.o $(OBJDIR)/TransitionMatrix.o
+CXX_NPTM_OBJS=$(OBJDIR)/Commons.o $(OBJDIR)/Configuration.o $(OBJDIR)/file_io.o $(OBJDIR)/Parsers.o $(OBJDIR)/sph_subs.o $(OBJDIR)/clu_subs.o $(OBJDIR)/tfrfme.o $(OBJDIR)/tra_subs.o $(OBJDIR)/TransitionMatrix.o
 
-CXX_NPTM_DYNOBJS=$(DYNOBJDIR)/Commons.o $(DYNOBJDIR)/Configuration.o $(DYNOBJDIR)/file_io.o $(DYNOBJDIR)/Parsers.o $(DYNOBJDIR)/sph_subs.o $(DYNOBJDIR)/clu_subs.o $(DYNOBJDIR)/tra_subs.o $(DYNOBJDIR)/TransitionMatrix.o
+CXX_NPTM_DYNOBJS=$(DYNOBJDIR)/Commons.o $(DYNOBJDIR)/Configuration.o $(DYNOBJDIR)/file_io.o $(DYNOBJDIR)/Parsers.o $(DYNOBJDIR)/sph_subs.o $(DYNOBJDIR)/clu_subs.o $(DYNOBJDIR)/tfrfme.o $(DYNOBJDIR)/tra_subs.o $(DYNOBJDIR)/TransitionMatrix.o
 
 
 all: $(BUILDDIR_NPTM)/libnptm.a $(BUILDDIR_NPTM)/libnptm.so

src/libnptm/TransitionMatrix.cpp

@@ -6,6 +6,7 @@
 #include <exception>
 #include <fstream>
 #include <hdf5.h>
+#include <string>
 
 #ifndef INCLUDE_LIST_H_
 #include "../include/List.h"

src/libnptm/tfrfme.cpp

@@ -3,11 +3,22 @@
  * \brief Implementation of the trapping configuration objects.
  */
 
-#include<complex>
-#include<string>
+#include <complex>
+#include <exception>
+#include <fstream>
+#include <hdf5.h>
+#include <string>
+
+#ifndef INCLUDE_LIST_H_
+#include "../include/List.h"
+#endif
 
 #ifndef INCLUDE_TFRFME_H_
-#include "../include/tfrfme.cpp"
+#include "../include/tfrfme.h"
+#endif
+
+#ifndef INCLUDE_FILE_IO_H_
+#include "../include/file_io.h"
 #endif
 
 using namespace std;
@@ -46,24 +57,292 @@ TFRFME::~TFRFME() {
   delete[] wsum;
 }
 
-TFRFME::TFRFME* from_binary(string file_name, string format) {
+TFRFME* TFRFME::from_binary(string file_name, string mode) {
+  TFRFME *instance = NULL;
+  if (mode.compare("LEGACY") == 0) {
+    instance = from_legacy(file_name);
+  } else if (mode.compare("HDF5") == 0) {
+    instance = from_hdf5(file_name);
+  } /* else {
+    string message = "Unknown format mode: \"" + mode + "\"";
+    throw UnrecognizedFormatException(message);
+    } */
+  return instance;
+}
+
+TFRFME* TFRFME::from_hdf5(string file_name) {
+  TFRFME *instance = 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 lmode, lm, nkv, nxv, nyv, nzv;
+    double vk, exri, an, ff, tra, spd, frsh, exril;
+    status = hdf_file->read("LMODE", "INT32", &lmode);
+    status = hdf_file->read("LM", "INT32", &lm);
+    status = hdf_file->read("NKV", "INT32", &nkv);
+    status = hdf_file->read("NXV", "INT32", &nxv);
+    status = hdf_file->read("NYV", "INT32", &nyv);
+    status = hdf_file->read("NZV", "INT32", &nzv);
+    status = hdf_file->read("VK", "FLOAT64", &vk);
+    status = hdf_file->read("EXRI", "FLOAT64", &exri);
+    status = hdf_file->read("AN", "FLOAT64", &an);
+    status = hdf_file->read("FF", "FLOAT64", &ff);
+    status = hdf_file->read("TRA", "FLOAT64", &tra);
+    status = hdf_file->read("SPD", "FLOAT64", &spd);
+    status = hdf_file->read("FRSH", "FLOAT64", &frsh);
+    status = hdf_file->read("EXRIL", "FLOAT64", &exril);
+    instance = new TFRFME(lmode, lm, nkv, nxv, nyv, nzv);
+    instance->set_param("vk", vk);
+    instance->set_param("exri", exri);
+    instance->set_param("an", an);
+    instance->set_param("ff", ff);
+    instance->set_param("tra", tra);
+    instance->set_param("spd", spd);
+    instance->set_param("frsh", frsh);
+    instance->set_param("exril", exril);
+    // DA QUI: aggiungere lettura di vettori e matrice.
+  }
+  return instance;
+}
+
+TFRFME* TFRFME::from_legacy(string file_name) {
+  fstream input;
   TFRFME *instance = NULL;
+  input.open(file_name.c_str(), ios::in | ios::binary);
+  if (input.is_open()) {
+    int lmode, lm, nkv, nxv, nyv, nzv;
+    double vk, exri, an, ff, tra, spd, frsh, exril;
+    double dval, rval, ival;
+    input.read(reinterpret_cast<char *>(&lmode), sizeof(int));
+    input.read(reinterpret_cast<char *>(&lm), sizeof(int));
+    input.read(reinterpret_cast<char *>(&nkv), sizeof(int));
+    input.read(reinterpret_cast<char *>(&nxv), sizeof(int));
+    input.read(reinterpret_cast<char *>(&nyv), sizeof(int));
+    input.read(reinterpret_cast<char *>(&nzv), sizeof(int));
+    input.read(reinterpret_cast<char *>(&vk), sizeof(double));
+    input.read(reinterpret_cast<char *>(&exri), sizeof(double));
+    input.read(reinterpret_cast<char *>(&an), sizeof(double));
+    input.read(reinterpret_cast<char *>(&ff), sizeof(double));
+    input.read(reinterpret_cast<char *>(&tra), sizeof(double));
+    input.read(reinterpret_cast<char *>(&spd), sizeof(double));
+    input.read(reinterpret_cast<char *>(&frsh), sizeof(double));
+    input.read(reinterpret_cast<char *>(&exril), sizeof(double));
+    instance = new TFRFME(lmode, lm, nkv, nxv, nyv, nzv);
+    instance->set_param("vk", vk);
+    instance->set_param("exri", exri);
+    instance->set_param("an", an);
+    instance->set_param("ff", ff);
+    instance->set_param("tra", tra);
+    instance->set_param("spd", spd);
+    instance->set_param("frsh", frsh);
+    instance->set_param("exril", exril);
+    for (int xi = 0; xi < nxv; xi++) {
+      input.read(reinterpret_cast<char *>(&dval), sizeof(double));
+      instance->set_x(xi, dval);
+    }
+    for (int yi = 0; yi < nyv; yi++) {
+      input.read(reinterpret_cast<char *>(&dval), sizeof(double));
+      instance->set_y(yi, dval);
+    }
+    for (int zi = 0; zi < nzv; zi++) {
+      input.read(reinterpret_cast<char *>(&dval), sizeof(double));
+      instance->set_z(zi, dval);
+    }
+    int _nlmmt = 2 * lm * (lm + 2);
+    int _nrvc = nxv * nyv * nzv;
+    for (int wi = 0; wi < _nlmmt; wi++) {
+      for (int wj = 0; wj < _nrvc; wj++) {
+	input.read(reinterpret_cast<char *>(&rval), sizeof(double));
+	input.read(reinterpret_cast<char *>(&ival), sizeof(double));
+	complex<double> value(rval, ival);
+	instance->set_matrix_element(wi, wj, value);
+      } // wj loop
+    } // wi loop
+  } else {
+    printf("ERROR: could not open input file \"%s\"\n", file_name.c_str());
+  }
   return instance;
 }
 
-TFRFME::double get_param(string param_name) {
+std::complex<double> TFRFME::get_matrix_element(int row, int col) {
+  return wsum[row][col];
+}
+
+double TFRFME::get_param(string param_name) {
   double value;
-  if (param_name.compare("exri") == 0) value = exri;
+  if (param_name.compare("vk") == 0) value = vk;
+  else if (param_name.compare("exri") == 0) value = exri;
   else if (param_name.compare("an") == 0) value = an;
   else if (param_name.compare("ff") == 0) value = ff;
   else if (param_name.compare("tra") == 0) value = tra;
   else if (param_name.compare("spd") == 0) value = spd;
   else if (param_name.compare("frsh") == 0) value = frsh;
   else if (param_name.compare("exril") == 0) value = exril;
-  //TODO: add else clause with exception.
+  else {
+    string message = "Unrecognized parameter name \"" + param_name + "\"";
+    UnrecognizedParameterException ex(message);
+    throw ex;
+  }
   return value;
 }
 
-TFRFME::void write_binary(string file_name, string format) {
-  return;
+void TFRFME::set_matrix_element(int row, int col, complex<double> value) {
+  wsum[row][col] = value;
+}
+
+void TFRFME::set_param(string param_name, double value) {
+  if (param_name.compare("vk") == 0) vk = value;
+  else if (param_name.compare("exri") == 0) exri = value;
+  else if (param_name.compare("an") == 0) an = value;
+  else if (param_name.compare("ff") == 0) ff = value;
+  else if (param_name.compare("tra") == 0) tra = value;
+  else if (param_name.compare("spd") == 0) spd = value;
+  else if (param_name.compare("frsh") == 0) frsh = value;
+  else if (param_name.compare("exril") == 0) exril = value;
+  else {
+    string message = "Unrecognized parameter name \"" + param_name + "\"";
+    UnrecognizedParameterException ex(message);
+    throw ex;
+  }
+}
+
+void TFRFME::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 UnrecognizedFormatException(message);
+    } */
+}
+
+void TFRFME::write_hdf5(string file_name) {
+  List<string> rec_name_list(1);
+  List<string> rec_type_list(1);
+  List<void *> rec_ptr_list(1);
+  string str_type;
+  rec_name_list.set(0, "LMODE");
+  rec_type_list.set(0, "INT32_(1)");
+  rec_ptr_list.set(0, &lmode);
+  rec_name_list.append("LM");
+  rec_type_list.append("INT32_(1)");
+  rec_ptr_list.append(&lm);
+  rec_name_list.append("NKV");
+  rec_type_list.append("INT32_(1)");
+  rec_ptr_list.append(&nkv);
+  rec_name_list.append("NXV");
+  rec_type_list.append("INT32_(1)");
+  rec_ptr_list.append(&nxv);
+  rec_name_list.append("NYV");
+  rec_type_list.append("INT32_(1)");
+  rec_ptr_list.append(&nyv);
+  rec_name_list.append("NZV");
+  rec_type_list.append("INT32_(1)");
+  rec_ptr_list.append(&nzv);
+  rec_name_list.append("VK");
+  rec_type_list.append("FLOAT64_(1)");
+  rec_ptr_list.append(&lm);
+  rec_name_list.append("EXRI");
+  rec_type_list.append("FLOAT64_(1)");
+  rec_ptr_list.append(&exri);
+  rec_name_list.append("AN");
+  rec_type_list.append("FLOAT64_(1)");
+  rec_ptr_list.append(&an);
+  rec_name_list.append("FF");
+  rec_type_list.append("FLOAT64_(1)");
+  rec_ptr_list.append(&ff);
+  rec_name_list.append("TRA");
+  rec_type_list.append("FLOAT64_(1)");
+  rec_ptr_list.append(&tra);
+  rec_name_list.append("SPD");
+  rec_type_list.append("FLOAT64_(1)");
+  rec_ptr_list.append(&spd);
+  rec_name_list.append("FRSH");
+  rec_type_list.append("FLOAT64_(1)");
+  rec_ptr_list.append(&frsh);
+  rec_name_list.append("EXRIL");
+  rec_type_list.append("FLOAT64_(1)");
+  rec_ptr_list.append(&exril);
+  str_type = "FLOAT64_(" + to_string(nxv) + ")";
+  rec_name_list.append("XVEC");
+  rec_type_list.append(str_type);
+  rec_ptr_list.append(xv);
+  str_type = "FLOAT64_(" + to_string(nyv) + ")";
+  rec_name_list.append("YVEC");
+  rec_type_list.append(str_type);
+  rec_ptr_list.append(yv);
+  str_type = "FLOAT64_(" + to_string(nzv) + ")";
+  rec_name_list.append("ZVEC");
+  rec_type_list.append(str_type);
+  rec_ptr_list.append(zv);
+  rec_name_list.append("WSUM");
+  str_type = "FLOAT64_(" + to_string(nlmmt) + "," + to_string(2 * nrvc) + ")";
+  rec_type_list.append(str_type);
+  // The (N x M) matrix of complex is converted to a (N x 2M) matrix of double,
+  // where REAL(E_i,j) -> E_i,(2 j) and IMAG(E_i,j) -> E_i,(2 j + 1)
+  int num_elements = 2 * nlmmt * nrvc;
+  double * ptr_elements = new double[num_elements]();
+  for (int wi = 0; wi < nlmmt; wi++) {
+    for (int wj = 0; wj < nrvc; wj++) {
+      int index = (2 * nrvc) * wi + 2 * wj;
+      ptr_elements[index] = wsum[wi][wj].real();
+      ptr_elements[index + 1] = wsum[wi][wj].imag();
+    } // wj loop
+  } // wi loop
+  rec_ptr_list.append(ptr_elements);
+
+  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, file_name, 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();
+
+  delete[] ptr_elements;
+  delete[] rec_names;
+  delete[] rec_types;
+  delete[] rec_pointers;
+  delete hdf_file;
+}
+
+void TFRFME::write_legacy(string file_name) {
+  fstream output;
+  output.open(file_name.c_str(), ios::out | ios::binary);
+  if (output.is_open()) {
+    output.write(reinterpret_cast<char *>(&lmode), sizeof(int));
+    output.write(reinterpret_cast<char *>(&lm), sizeof(int));
+    output.write(reinterpret_cast<char *>(&nkv), sizeof(int));
+    output.write(reinterpret_cast<char *>(&nxv), sizeof(int));
+    output.write(reinterpret_cast<char *>(&nyv), sizeof(int));
+    output.write(reinterpret_cast<char *>(&nzv), sizeof(int));
+    output.write(reinterpret_cast<char *>(&vk), sizeof(double));
+    output.write(reinterpret_cast<char *>(&exri), sizeof(double));
+    output.write(reinterpret_cast<char *>(&an), sizeof(double));
+    output.write(reinterpret_cast<char *>(&ff), sizeof(double));
+    output.write(reinterpret_cast<char *>(&tra), sizeof(double));
+    output.write(reinterpret_cast<char *>(&spd), sizeof(double));
+    output.write(reinterpret_cast<char *>(&frsh), sizeof(double));
+    output.write(reinterpret_cast<char *>(&exril), sizeof(double));
+    for (int xi = 0; xi < nxv; xi++)
+      output.write(reinterpret_cast<char *>(&(xv[xi])), sizeof(double));
+    for (int yi = 0; yi < nyv; yi++)
+      output.write(reinterpret_cast<char *>(&(yv[yi])), sizeof(double));
+    for (int zi = 0; zi < nzv; zi++)
+      output.write(reinterpret_cast<char *>(&(zv[zi])), sizeof(double));
+    for (int wi = 0; wi < nlmmt; wi++) {
+      for (int wj = 0; wj < nrvc; wj++) {
+	double rval = wsum[wi][wj].real();
+	double ival = wsum[wi][wj].imag();
+	output.write(reinterpret_cast<char *>(&rval), sizeof(double));
+	output.write(reinterpret_cast<char *>(&ival), sizeof(double));
+      } // wj loop
+    } // wi loop
+  } else { // Should never occur
+    printf("ERROR: could not open output file \"%s\"\n", file_name.c_str());
+  }
 }