From 281e5c6cd19752d720048dee10041d4f118b6434 Mon Sep 17 00:00:00 2001
From: Giovanni La Mura <giovanni.lamura@inaf.it>
Date: Fri, 19 Apr 2024 15:38:34 +0200
Subject: [PATCH] Implement static output for single sphere TM

---
 src/include/TransitionMatrix.h   |  77 +++++++++++++++++++--
 src/libnptm/TransitionMatrix.cpp | 111 +++++++++++++++++++++++++++++--
 2 files changed, 178 insertions(+), 10 deletions(-)

diff --git a/src/include/TransitionMatrix.h b/src/include/TransitionMatrix.h
index 54903a95..0e51f1d8 100644
--- a/src/include/TransitionMatrix.h
+++ b/src/include/TransitionMatrix.h
@@ -60,8 +60,8 @@ class TransitionMatrix {
    *
    * This function takes care of the specific task of writing the transition
    * matrix memory data structure to a binary output file formatted according
-   * to the HDF5 standard. It is designed to work for the case of clusters of
-   * spheres.
+   * to the HDF5 standard without a pre-existing instance. It is designed to
+   * work for the case of a cluster of spheres.
    *
    * \param file_name: `string` Name of the binary configuration data file.
    * \param _nlemt: `np_int` Size of the matrix (2 * LE * (LE + 2)).
@@ -75,13 +75,33 @@ class TransitionMatrix {
 			 double _exri, dcomplex **_am0m
   );
   
+  /*! \brief Write transition matrix data to HDF5 binary output.
+   *
+   * This function takes care of the specific task of writing the transition
+   * matrix memory data structure to a binary output file formatted according
+   * to the HDF5 standard without a pre-existing instance. It is designed to
+   * work for the case of a single sphere.
+   *
+   * \param file_name: `string` Name of the binary configuration data file.
+   * \param _lm: `int` Maximum field expansion order.
+   * \param _vk: `double` Wave number in scale units.
+   * \param _exri: `double` External medium refractive index.
+   * \param _rmi: `complex double **`
+   * \param _rei: `complex double **`
+   * \param _sphere_radius: `double` Radius of the sphere.
+   */
+  static void write_hdf5(
+			 std::string file_name, int _lm, double _vk, double _exri,
+			 dcomplex **_rmi, dcomplex **_rei, double _sphere_radius
+  );
+
   /*! \brief Write the Transition Matrix to legacy binary output.
    *
    * \param file_name: `string` Name of the binary configuration data file.
    */
   void write_legacy(std::string file_name);
 
-  /*! \brief Write transition matrix data to HDF5 binary output.
+  /*! \brief Write transition matrix data to binary output using legacy format.
    *
    * This function takes care of the specific task of writing the transition
    * matrix memory data structure to a binary output file formatted according
@@ -100,6 +120,26 @@ class TransitionMatrix {
 			   double _exri, dcomplex **_am0m
   );
   
+  /*! \brief Write transition matrix data to binary output using legacy format.
+   *
+   * This function takes care of the specific task of writing the transition
+   * matrix memory data structure to a binary output file formatted according
+   * to the original code structure without a pre-existing instance. It is designed
+   * to work for the case of a single sphere.
+   *
+   * \param file_name: `string` Name of the binary configuration data file.
+   * \param _lm: `int` Maximum field expansion order.
+   * \param _vk: `double` Wave number in scale units.
+   * \param _exri: `double` External medium refractive index.
+   * \param _rmi: `complex double **`
+   * \param _rei: `complex double **`
+   * \param _sphere_radius: `double` Radius of the sphere.
+   */
+  static void write_legacy(
+			   std::string file_name, int _lm, double _vk, double _exri,
+			   dcomplex **_rmi, dcomplex **_rei, double _sphere_radius
+  );
+
  public:
   /*! \brief Default Transition Matrix instance constructor.
    *
@@ -181,7 +221,7 @@ class TransitionMatrix {
    */
   void write_binary(std::string file_name, std::string mode="LEGACY");
   
-  /*! \brief Write a Transition Matrix to a binary file without instanciating it.
+  /*! \brief Write a cluster Transition Matrix to a binary file without instanciating it.
    *
    * Transition Matrix data can take a large amount of memory. For such reason, attempts
    * to create TransitionMatrix instances only for writing purposes can create
@@ -208,6 +248,35 @@ class TransitionMatrix {
 			   double _exri, dcomplex **_am0m, std::string mode="LEGACY"
   );
   
+  /*! \brief Write a single sphere Transition Matrix to a binary file without instanciating it.
+   *
+   * Transition Matrix data can take a large amount of memory. For such reason, attempts
+   * to create TransitionMatrix instances only for writing purposes can create
+   * unnecessary resource consumption and computing time to duplicate the data into
+   * the output buffer. This function offers output to file as a static method. It
+   * takes the arguments of a constructor together with the usual arguments to specify
+   * the output file name and format, to write the required data directly to a file,
+   * without creating a new TransitionMatrix instance. The implementation works for
+   * TransitionMatrix objects built for the single sphere case. It belongs to the public
+   * class interface and it calls the proper versions of `write_legacy()` and `write_hdf5()`,
+   * depending on the requested output format.
+   *
+   * \param file_name: `string` Name of the file to be written.
+   * \param _lm: `int` Maximum field expansion order.
+   * \param _vk: `double` Wave number in scale units.
+   * \param _exri: `double` External medium refractive index.
+   * \param _rmi: `complex double **`
+   * \param _rei: `complex double **`
+   * \param _sphere_radius: `double` Radius of the sphere.
+   * \param mode: `string` Binary encoding. Can be one of ["LEGACY", "HDF5"] . Optional
+   * (default is "LEGACY").
+   */
+  static void write_binary(
+			   std::string file_name, int _lm, double _vk, double _exri,
+			   dcomplex **_rmi, dcomplex **_rei, double _sphere_radius,
+			   std::string mode="LEGACY"
+  );
+
   /*! \brief Test whether two instances of TransitionMatrix are equal.
    *
    * Transition matrices can be the result of a calculation or of a file input operation,
diff --git a/src/libnptm/TransitionMatrix.cpp b/src/libnptm/TransitionMatrix.cpp
index 7468ec8d..477264ef 100644
--- a/src/libnptm/TransitionMatrix.cpp
+++ b/src/libnptm/TransitionMatrix.cpp
@@ -94,7 +94,7 @@ TransitionMatrix::TransitionMatrix(
   }
 }
 
-TransitionMatrix* TransitionMatrix::from_binary(string file_name, string mode) {
+TransitionMatrix* TransitionMatrix::from_binary(std::string file_name, std::string mode) {
   TransitionMatrix *tm = NULL;
   if (mode.compare("LEGACY") == 0) {
     tm = TransitionMatrix::from_legacy(file_name);
@@ -107,7 +107,7 @@ TransitionMatrix* TransitionMatrix::from_binary(string file_name, string mode) {
   return tm;
 }
 
-TransitionMatrix* TransitionMatrix::from_hdf5(string file_name) {
+TransitionMatrix* TransitionMatrix::from_hdf5(std::string file_name) {
   TransitionMatrix *tm = NULL;
   unsigned int flags = H5F_ACC_RDONLY;
   HDFFile *hdf_file = new HDFFile(file_name, flags);
@@ -160,7 +160,7 @@ TransitionMatrix* TransitionMatrix::from_hdf5(string file_name) {
   return tm;
 }
 
-TransitionMatrix* TransitionMatrix::from_legacy(string file_name) {
+TransitionMatrix* TransitionMatrix::from_legacy(std::string file_name) {
   fstream ttms;
   TransitionMatrix *tm = NULL;
   ttms.open(file_name, ios::binary | ios::in);
@@ -207,7 +207,7 @@ TransitionMatrix* TransitionMatrix::from_legacy(string file_name) {
   return tm;
 }
 
-void TransitionMatrix::write_binary(string file_name, string mode) {
+void TransitionMatrix::write_binary(std::string file_name, std::string mode) {
   if (mode.compare("LEGACY") == 0) {
     write_legacy(file_name);
   } else if (mode.compare("HDF5") == 0) {
@@ -232,7 +232,22 @@ void TransitionMatrix::write_binary(
   }
 }
 
-void TransitionMatrix::write_hdf5(string file_name) {
+void TransitionMatrix::write_binary(
+				    std::string file_name, int _lm, double _vk, double _exri,
+				    dcomplex **_rmi, dcomplex **_rei, double _sphere_radius,
+				    std::string mode
+) {
+  if (mode.compare("LEGACY") == 0) {
+    write_legacy(file_name, _lm, _vk, _exri, _rmi, _rei, _sphere_radius);
+  } else if (mode.compare("HDF5") == 0) {
+    write_hdf5(file_name, _lm, _vk, _exri, _rmi, _rei, _sphere_radius);
+  } else {
+    string message = "Unknown format mode: \"" + mode + "\"";
+    throw UnrecognizedFormatException(message);
+  }
+}
+
+void TransitionMatrix::write_hdf5(std::string file_name) {
   if (is == 1 || is == 1111) {
     List<string> rec_name_list(1);
     List<string> rec_type_list(1);
@@ -328,7 +343,58 @@ void TransitionMatrix::write_hdf5(
   delete hdf_file;
 }
 
-void TransitionMatrix::write_legacy(string file_name) {
+void TransitionMatrix::write_hdf5(
+				  std::string file_name, int _lm, double _vk, double _exri,
+				  dcomplex **_rmi, dcomplex **_rei, double _sphere_radius
+) {
+  int is = 1111;
+  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, "IS");
+  rec_type_list.set(0, "INT32_(1)");
+  rec_ptr_list.set(0, &is);
+  rec_name_list.append("L_MAX");
+  rec_type_list.append("INT32_(1)");
+  rec_ptr_list.append(&_lm);
+  rec_name_list.append("VK");
+  rec_type_list.append("FLOAT64_(1)");
+  rec_ptr_list.append(&_vk);
+  rec_name_list.append("EXRI");
+  rec_type_list.append("FLOAT64_(1)");
+  rec_ptr_list.append(&_exri);
+  dcomplex *_elements = new dcomplex[2 * _lm]();
+  for (int ei = 0; ei < _lm; ei++) {
+    _elements[2 * ei] = -1.0 / _rmi[ei][0];
+    _elements[2 * ei + 1] = -1.0 / _rei[ei][0];
+  }
+  rec_name_list.append("ELEMENTS");
+  str_type = "COMPLEX128_(" + to_string(_lm) + "," + to_string(2) + ")";
+  rec_type_list.append(str_type);
+  rec_ptr_list.append(_elements);
+  rec_name_list.append("RADIUS");
+  rec_type_list.append("FLOAT64_(1)");
+  rec_ptr_list.append(&_sphere_radius);
+
+  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[] _elements;
+  delete[] rec_names;
+  delete[] rec_types;
+  delete[] rec_pointers;
+  delete hdf_file;
+}
+
+void TransitionMatrix::write_legacy(std::string file_name) {
   fstream ttms;
   if (is == 1111 || is == 1) {
     ttms.open(file_name, ios::binary | ios::out);
@@ -388,6 +454,39 @@ void TransitionMatrix::write_legacy(
   }
 }
 
+void TransitionMatrix::write_legacy(
+				    std::string file_name, int _lm, double _vk, double _exri,
+				    dcomplex **_rmi, dcomplex **_rei, double _sphere_radius
+) {
+  fstream ttms;
+  int is = 1111;
+  ttms.open(file_name, ios::binary | ios::out);
+  if (ttms.is_open()) {
+    ttms.write(reinterpret_cast<char *>(&is), sizeof(int));
+    ttms.write(reinterpret_cast<char *>(&_lm), sizeof(int));
+    ttms.write(reinterpret_cast<char *>(&_vk), sizeof(double));
+    ttms.write(reinterpret_cast<char *>(&_exri), sizeof(double));
+    double rval, ival;
+    dcomplex element;
+    for (int ei = 0; ei < _lm; ei++) {
+      element = -1.0 / _rmi[ei][0];
+      rval = real(element);
+      ival = imag(element);
+      ttms.write(reinterpret_cast<char *>(&rval), sizeof(double));
+      ttms.write(reinterpret_cast<char *>(&ival), sizeof(double));
+      element = -1.0 / _rei[ei][0];
+      rval = real(element);
+      ival = imag(element);
+      ttms.write(reinterpret_cast<char *>(&rval), sizeof(double));
+      ttms.write(reinterpret_cast<char *>(&ival), sizeof(double));
+    }
+    ttms.write(reinterpret_cast<char *>(&_sphere_radius), sizeof(double));
+    ttms.close();
+  } else { // Failed to open output file. Should never happen.
+    printf("ERROR: could not open Transition Matrix file for writing.\n");
+  }
+}
+
 bool TransitionMatrix::operator ==(TransitionMatrix &other) {
   if (is != other.is) {
     return false;
-- 
GitLab