merged main into feature/polytrope

src/config/meson.build

@@ -0,0 +1,26 @@
+# Define the library
+config_sources = files(
+  'private/config.cpp',
+)
+
+config_headers = files(
+    'public/config.h'
+)
+
+# Define the libconfig library so it can be linked against by other parts of the build system
+libconfig = static_library('config',
+    config_sources, 
+    include_directories: include_directories('public'),
+    cpp_args: ['-fvisibility=default'],
+    dependencies: [yaml_cpp_dep],
+    install : true)
+
+config_dep = declare_dependency(
+    include_directories: include_directories('public'),
+    link_with: libconfig,
+    sources: config_sources,
+    dependencies: [yaml_cpp_dep],
+)
+
+# Make headers accessible
+install_headers(config_headers, subdir : '4DSSE/config')

src/config/private/config.cpp

@@ -0,0 +1,42 @@
+#include <string>
+#include <iostream>
+#include <fstream>
+#include <sstream>
+#include <vector>
+#include <map>
+
+#include "yaml-cpp/yaml.h"
+
+#include "config.h"
+
+Config::Config() {}
+
+Config::~Config() {}
+
+Config& Config::getInstance() {
+    static Config instance;
+    return instance;
+}
+
+bool Config::loadConfig(const std::string& configFile) {
+    configFilePath = configFile;
+    try {
+        yamlRoot = YAML::LoadFile(configFile);
+    } catch (YAML::BadFile& e) {
+        std::cerr << "Error: " << e.what() << std::endl;
+        return false;
+    }
+    return true;
+}
+
+bool Config::isKeyInCache(const std::string &key) {
+    return configMap.find(key) != configMap.end();
+}
+
+void Config::addToCache(const std::string &key, const YAML::Node &node) {
+    configMap[key] = node;
+}
+
+void Config::registerUnknownKey(const std::string &key) {
+    unknownKeys.push_back(key);
+}

src/config/public/config.h

@@ -0,0 +1,175 @@
+#ifndef CONFIG_H
+#define CONFIG_H
+
+#include <string>
+#include <iostream>
+#include <fstream>
+#include <sstream>
+#include <vector>
+#include <map>
+#include <algorithm>
+
+#include "yaml-cpp/yaml.h"
+
+/**
+ * @class Config
+ * @brief Singleton class to manage configuration settings loaded from a YAML file.
+ */
+class Config {
+private:
+   /**
+    * @brief Private constructor to prevent instantiation.
+    */
+   Config();
+
+   /**
+    * @brief Destructor.
+    */
+   ~Config();
+
+   YAML::Node yamlRoot; ///< Root node of the YAML configuration.
+   std::string configFilePath; ///< Path to the configuration file.
+   bool debug = false; ///< Flag to enable debug output.
+
+   std::map<std::string, YAML::Node> configMap; ///< Cache for the location of configuration settings.
+   std::vector<std::string> unknownKeys; ///< Cache for the existence of configuration settings.
+
+   /**
+    * @brief Get a value from the configuration cache.
+    * @tparam T Type of the value to retrieve.
+    * @param key Key of the configuration value.
+    * @param defaultValue Default value to return if the key does not exist.
+    * @return Configuration value of type T.
+    */
+   template <typename T>
+   T getFromCache(const std::string &key, T defaultValue) {
+      if (configMap.find(key) != configMap.end()) {
+         try {
+               return configMap[key].as<T>();
+         } catch (const YAML::Exception& e) {
+               return defaultValue;
+         }
+      }
+      return defaultValue;
+   }
+
+
+   /**
+    * @brief Check if a key exists in the configuration cache.
+    * @param key Key to check.
+    * @return True if the key exists in the cache, false otherwise.
+    */
+   bool isKeyInCache(const std::string &key);
+
+   /**
+    * @brief Add a key-value pair to the configuration cache.
+    * @param key Key of the configuration value.
+    * @param node YAML node containing the configuration value.
+    */
+   void addToCache(const std::string &key, const YAML::Node &node);
+
+   /**
+    * @brief Register a key as not found in the configuration.
+    * @param key Key that was not found.
+    */
+   void registerUnknownKey(const std::string &key);
+
+public:
+   /**
+    * @brief Get the singleton instance of the Config class.
+    * @return Reference to the Config instance.
+    */
+   static Config& getInstance();
+
+   Config (const Config&) = delete;
+   Config& operator= (const Config&) = delete;
+   Config (Config&&) = delete;
+   Config& operator= (Config&&) = delete;
+
+   void setDebug(bool debug) { this->debug = debug; }
+
+   /**
+    * @brief Load configuration from a YAML file.
+    * @param configFilePath Path to the YAML configuration file.
+    * @return True if the configuration was loaded successfully, false otherwise.
+    */
+   bool loadConfig(const std::string& configFilePath);
+
+   /**
+    * @brief Get the input table from the configuration.
+    * @return Input table as a string.
+    */
+   std::string getInputTable() const;
+
+   /**
+    * @brief Get a configuration value by key.
+    * @tparam T Type of the value to retrieve.
+    * @param key Key of the configuration value.
+    * @param defaultValue Default value to return if the key does not exist.
+    * @return Configuration value of type T.
+    * 
+    * @example
+    * @code
+    * Config& config = Config::getInstance();
+    * config.loadConfig("example.yaml");
+    * int maxIter = config.get<int>("opac:lowTemp:numeric:maxIter", 10);
+    */
+   template <typename T>
+   T get(const std::string &key, T defaultValue) {
+      // --- Check if the key has already been checked for existence 
+      if (std::find(unknownKeys.begin(), unknownKeys.end(), key) != unknownKeys.end()) {
+         return defaultValue; // If the key has already been added to the unknown cache do not traverse the YAML tree or hit the cache
+      }
+
+      // --- Check if the key is already in the cache (avoid traversing YAML nodes)
+      if (isKeyInCache(key)) {
+         return getFromCache<T>(key, defaultValue);
+      }
+      // --- If the key is not in the cache, check the YAML file
+      else {
+         YAML::Node node = YAML::Clone(yamlRoot);
+         std::istringstream keyStream(key);
+         std::string subKey;
+         while (std::getline(keyStream, subKey, ':')) {
+            if (!node[subKey]) {
+                  // Key does not exist
+                  registerUnknownKey(key);
+                  return defaultValue;
+            }
+            node = node[subKey]; // go deeper
+         }
+
+         try {
+            // Key exists and is of the requested type
+            addToCache(key, node);
+            return node.as<T>();
+         } catch (const YAML::Exception& e) {
+            // Key is not of the requested type
+            registerUnknownKey(key);
+            return defaultValue; // return default value if the key does not exist
+         }
+      }
+   }
+
+   /**
+    * @brief Print the configuration file path and the YAML root node.
+    * @param os Output stream.
+    * @param config Config object to print.
+    * @return Output stream.
+    */
+   friend std::ostream& operator<<(std::ostream& os, const Config& config) {
+      if (!config.debug) {
+         os << "Config file: " << config.configFilePath << std::endl;
+      } else{
+         // Print entire YAML file from root
+         os << "Config file: " << config.configFilePath << std::endl;
+         os << config.yamlRoot << std::endl;
+      }
+      return os;
+   }
+
+   // Setup gTest class as a friend
+   friend class configTestPrivateAccessor;
+};
+
+#endif

src/meson.build

@@ -6,4 +6,6 @@ subdir('dobj')
 subdir('const')
 subdir('opatIO')
 subdir('meshIO')
-subdir('poly')
+subdir('probe')
+subdir('poly')
+subdir('config')

src/opatIO/meson.build

@@ -12,7 +12,9 @@ libopatIO = library('opatIO',
     opatIO_sources, 
     include_directories: include_directories('public'),
     cpp_args: ['-fvisibility=default'],
-    install : true)
+    dependencies: [picosha2_dep],
+    install : true,
+    )
 
 # Make headers accessible
 install_headers(opatIO_headers, subdir : '4DSSE/opatIO')

src/opatIO/private/opatIO.cpp

@@ -10,6 +10,7 @@
 #include <cmath>
 #include <limits>
 #include <deque>
+#include "picosha2.h"
 
 // Function to check system endianness
 bool is_big_endian() {
@@ -75,13 +76,14 @@ void OpatIO::unload() {
 void OpatIO::readHeader(std::ifstream &file) {
     file.read(reinterpret_cast<char*>(&header), sizeof(Header));
     if (file.gcount() != sizeof(Header)) {
-        throw std::runtime_error("Error reading header from file: " + filename);
+        throw std::runtime_error("Error reading header from file");
     }
 
     if (is_big_endian()) {
         header.version = swap_bytes(header.version);
         header.numTables = swap_bytes(header.numTables);
         header.indexOffset = swap_bytes(header.indexOffset);
+        header.numIndex = swap_bytes(header.numIndex);
     }
 }
 
@@ -89,75 +91,86 @@ void OpatIO::readHeader(std::ifstream &file) {
 void OpatIO::readTableIndex(std::ifstream &file) {
     file.seekg(header.indexOffset, std::ios::beg);
     tableIndex.resize(header.numTables);
-    file.read(reinterpret_cast<char*>(tableIndex.data()), header.numTables * sizeof(TableIndex));
-    if (file.gcount() != static_cast<std::streamsize>(header.numTables * sizeof(TableIndex))) {
-        throw std::runtime_error("Error reading table index from file: " + filename);
+    long unsigned int indexReadBytes;
+    for (uint32_t i = 0; i < header.numTables; i++) {
+        indexReadBytes = 0;
+        // Read the index vector in based on numIndex
+        tableIndex.at(i).index.resize(header.numIndex);
+        file.read(reinterpret_cast<char*>(tableIndex.at(i).index.data()), header.numIndex * sizeof(double));
+        indexReadBytes += static_cast<int>(file.gcount());
+
+        // Read the start and end position of the table in 
+        file.read(reinterpret_cast<char*>(&tableIndex.at(i).byteStart), sizeof(uint64_t));
+        indexReadBytes += static_cast<int>(file.gcount());
+        file.read(reinterpret_cast<char*>(&tableIndex.at(i).byteEnd), sizeof(uint64_t));
+        indexReadBytes += static_cast<int>(file.gcount());
+
+        // Read the checksum in
+        file.read(tableIndex.at(i).sha256, 32);
+        indexReadBytes += static_cast<int>(file.gcount());
+
+        // validate that the size of read data is correct
+        if (indexReadBytes != header.numIndex * sizeof(double) + 32 + 2 * sizeof(uint64_t)) {
+            throw std::runtime_error("Error reading table index from file");
+        }
     }
 
-    buildTableIDToComposition();
+
+    buildTableIDToIndex();
 }
 
-void OpatIO::buildTableIDToComposition(){
-    tableIDToComposition.clear();
+void OpatIO::buildTableIDToIndex(){
+    tableIDToIndex.clear();
     int tableID = 0;
-    std::pair<double, double> comp;
-    for (const auto &index : tableIndex) {
-        comp.first = index.X;
-        comp.second = index.Z;
-        tableIDToComposition.emplace(tableID, comp);
+    std::vector<double> ind;
+    ind.resize(header.numIndex);
+    for (const auto &table : tableIndex) {
+        ind.clear();
+        for (const auto &index : table.index) {
+            ind.push_back(index);
+        }
+        tableIDToIndex.emplace(tableID, ind);
         tableID++;
     }
-    XZLookupEpsilon();
+    LookupEpsilon();
 }
 
-void OpatIO::XZLookupEpsilon() {
+void OpatIO::LookupEpsilon() {
     /* 
-     Get 10% of the minimum spacing between XZ values 
-     in the tableIDToComposition map. This can be used
+     Get 10% of the minimum spacing between index values 
+     in the tableIDToIndex map. This can be used
      to set the comparison distance when doing a reverse
-     lookup (composition -> tableID)
+     lookup (index -> tableID)
      */
-    std::vector<double> Xvalues, Zvalues;
-    double epsilonX, epsilonZ, xgap, zgap;
+    indexEpsilon.resize(header.numIndex);
 
-    // Start these out as larger than they will ever be
-    epsilonX = 1;
-    epsilonZ = 1;
-    for (const auto& pair : tableIDToComposition) {
-        Xvalues.push_back(pair.second.first);
-        Zvalues.push_back(pair.second.second);
-    }
-
-    // Sorting is required for this algorithm.
-    std::sort(Xvalues.begin(), Xvalues.end());
-    std::sort(Zvalues.begin(), Zvalues.end());
-
-    for (size_t i = 1; i < Xvalues.size(); ++i) {
-        xgap = Xvalues[i] - Xvalues[i - 1];
-        zgap = Zvalues[i] - Zvalues[i - 1];
-        if (xgap > 0 && xgap < epsilonX) {
-            epsilonX = xgap;
-        }
-        if (zgap > 0 && zgap < epsilonZ) {
-            epsilonZ = zgap;
+    double epsilon;
+    for (int i = 0; i < static_cast<int>(header.numIndex); i++) {
+        epsilon = std::numeric_limits<double>::max();
+        for (int j = 1; j < static_cast<int>(header.numTables); j++) {
+            epsilon = std::min(epsilon, std::fabs(tableIDToIndex.at(j).at(i) - tableIDToIndex.at(j-1).at(i)));
         }
+        indexEpsilon.at(i) = epsilon * 0.1;
     }
-    // 0.1 to extract 10% of min distance.
-    XZepsilon = {0.1*epsilonX, 0.1*epsilonZ};
 }
 
-int OpatIO::lookupTableID(double X, double Z){
-    bool XOkay;
-    bool ZOkay;
+int OpatIO::lookupTableID(std::vector<double> index) {
+    std::vector<bool> IndexOkay;
+    IndexOkay.resize(header.numIndex);
     int tableID = 0;
-    for (const auto &tableMap : tableIDToComposition){
-        XOkay = std::fabs(tableMap.second.first - X) < XZepsilon.first;
-        ZOkay = std::fabs(tableMap.second.second - Z) < XZepsilon.second;
-        if (XOkay and ZOkay){
+    for (const auto &tableMap : tableIDToIndex){
+        // Loop through all index values and check if they are within epsilon for that index
+        std::fill(IndexOkay.begin(), IndexOkay.end(), false);
+        for (long unsigned int i = 0; i < index.size(); i++) {
+            IndexOkay.at(i) = std::fabs(tableMap.second.at(i) - index.at(i)) < indexEpsilon.at(i);
+        }
+        // If all index values are within epsilon, return the table ID
+        if (std::all_of(IndexOkay.begin(), IndexOkay.end(), [](bool i){return i;})) {
             return tableID;
         }
         tableID++;
     }
+    // If no table is found, return -1 (sentinal value)
     return -1;
 }
 
@@ -201,10 +214,10 @@ uint16_t OpatIO::getOPATVersion() {
 }
 
 // Get a table for given X and Z
-OPATTable OpatIO::getTable(double X, double Z) {
-    int tableID = lookupTableID(X, Z);
+OPATTable OpatIO::getTable(std::vector<double> index) {
+    int tableID = lookupTableID(index);
     if (tableID == -1) {
-        throw std::out_of_range("X Z Pair Not found!");
+        throw std::out_of_range("Index Not found!");
     }
     try {
         return getTableFromQueue(tableID);
@@ -288,6 +301,7 @@ void OpatIO::printHeader() {
     std::cout << "Creation Date: " << header.creationDate << std::endl;
     std::cout << "Source Info: " << header.sourceInfo << std::endl;
     std::cout << "Comment: " << header.comment << std::endl;
+    std::cout << "Number of Indices: " << header.numIndex << std::endl;
 }
 
 // Print the table index
@@ -298,9 +312,11 @@ void OpatIO::printTableIndex() {
     }
 
     // Print table header
-    std::cout << std::left << std::setw(10) << "X"
-        << std::setw(10) << "Z"
-        << std::setw(15) << "Byte Start"
+    std::cout << std::left << std::setw(10);
+    for (int i = 0; i < header.numIndex; i++) {
+        std::cout << "Index " << i << std::setw(10);
+    }   
+    std::cout << std::setw(15) << "Byte Start"
         << std::setw(15) << "Byte End"
         << "Checksum (SHA-256)" << std::endl;
 
@@ -308,10 +324,11 @@ void OpatIO::printTableIndex() {
 
     // Print each entry in the table
     for (const auto &index : tableIndex) {
-        std::cout << std::fixed << std::setprecision(4)
-            << std::setw(10) << index.X
-            << std::setw(10) << index.Z
-            << std::setw(15) << index.byteStart
+        std::cout << std::fixed << std::setprecision(4) << std::setw(10);
+        for (int i = 0; i < header.numIndex; i++) {
+            std::cout << index.index[i] << std::setw(10);
+        }
+        std::cout << std::setw(5) << index.byteStart
             << std::setw(15) << index.byteEnd
             << std::hex; // Switch to hex mode for checksum
 
@@ -422,8 +439,8 @@ void OpatIO::printTable(OPATTable table, uint32_t truncateDigits) {
     std::cout << "]" << std::endl;
 }
 
-void OpatIO::printTable(double X, double Z, uint32_t truncateDigits) {
-    int tableID = lookupTableID(X, Z);
+void OpatIO::printTable(std::vector<double> index, uint32_t truncateDigits) {
+    int tableID = lookupTableID(index);
     OPATTable table = getTable(tableID);
     printTable(table, truncateDigits);
 }
@@ -438,23 +455,45 @@ Header OpatIO::getHeader() {
     return header;
 }
 
-// // Get the closest X tables
-// std::vector<OPATTable> OpatIO::getClosestXTables(double X, double ZExact, int numTables) {
-//     std::vector<OPATTable> closestTables;
-//     // Implement logic to find closest X tables
-//     return closestTables;
-// }
+// Get the size of the index vector used
+uint16_t OpatIO::getNumIndex() {
+    return header.numIndex;
+}
 
-// // Get the closest Z tables
-// std::vector<OPATTable> OpatIO::getClosestZTables(double XExact, double Z, int numTables) {
-//     std::vector<OPATTable> closestTables;
-//     // Implement logic to find closest Z tables
-//     return closestTables;
-// }
+TableIndex OpatIO::getTableIndex(std::vector<double> index) {
+    int tableID = lookupTableID(index);
+    return tableIndex.at(tableID);
+}
 
-// // Get the closest tables
-// std::vector<OPATTable> OpatIO::getClosestTables(double X, double Z, int numTables) {
-//     std::vector<OPATTable> closestTables;
-//     // Implement logic to find closest tables
-//     return closestTables;
-// }
+std::vector<unsigned char> OpatIO::computeChecksum(int tableID) {
+    OPATTable table = getTable(tableID);
+    std::vector<std::vector<double>> logKappa = table.logKappa;
+    std::vector<double> flatData(logKappa.size() * logKappa.size());
+    size_t offset = 0;
+    for (const auto& row : logKappa) {
+        for (const auto& val : row) {
+            flatData[offset++] = val;
+        }
+    }
+    std::vector<unsigned char> flatDataBytes(flatData.size() * sizeof(double));
+    std::memcpy(flatDataBytes.data(), flatData.data(), flatDataBytes.size());
+    std::vector<unsigned char> hash(picosha2::k_digest_size);
+    picosha2::hash256(flatDataBytes.begin(), flatDataBytes.end(), hash.begin(), hash.end());
+    return hash;
+}
+
+std::vector<unsigned char> OpatIO::computeChecksum(std::vector<double> index) {
+    int tableID = lookupTableID(index);
+    return computeChecksum(tableID);
+}
+
+bool OpatIO::validateAll() {
+    for (const auto &table : tableIDToIndex) {
+        std::vector<unsigned char> hash = computeChecksum(table.first);
+        std::vector<unsigned char> storedHash(tableIndex.at(table.first).sha256, tableIndex.at(table.first).sha256 + 32);
+        if (hash != storedHash) {
+            return false;
+        }
+    }
+    return true;
+}

src/opatIO/public/opatIO.h

@@ -22,7 +22,8 @@ struct Header {
     char creationDate[16];   ///< Creation date of the file
     char sourceInfo[64];     ///< Source information
     char comment[128];       ///< Comment section
-    char reserved[26];       ///< Reserved for future use
+    uint16_t numIndex;       ///< Size of index vector per table
+    char reserved[24];       ///< Reserved for future use
 };
 #pragma pack()
 
@@ -30,8 +31,7 @@ struct Header {
  * @brief Structure to hold the index information of a table in an OPAT file.
  */
 struct TableIndex {
-    double X;                ///< X composition value
-    double Z;                ///< Z composition value
+    std::vector<double> index; ///< Index vector for associated table
     uint64_t byteStart;      ///< Byte start position of the table
     uint64_t byteEnd;        ///< Byte end position of the table
     char sha256[32];         ///< SHA-256 hash of the table data
@@ -58,9 +58,9 @@ private:
     Header header; ///< Header information of the OPAT file
     std::vector<TableIndex> tableIndex; ///< Index information of the tables
     std::deque<std::pair<int, OPATTable>> tableQueue; ///< Queue to manage table caching
-    std::map<int, std::pair<double, double>> tableIDToComposition; ///< Map to store table ID to composition mapping
-    std::pair<double, double> XZepsilon; ///< Epsilon values for X and Z
-    int maxQDepth = 10; ///< Maximum depth of the table queue
+    std::map<int, std::vector<double>> tableIDToIndex; ///< Map to store table ID to indexing
+    std::vector<double> indexEpsilon; ///< Epsilon values for each index
+    int maxQDepth = 20; ///< Maximum depth of the table queue
     std::string filename; ///< Filename of the OPAT file
     bool loaded = false; ///< Flag to indicate if the file is loaded
 
@@ -115,14 +115,14 @@ private:
     void printTable(OPATTable table, uint32_t truncateDigits=5);
 
     /**
-     * @brief Lookup epsilon values for X and Z.
+     * @brief Lookup epsilon values for Index.
      */
-    void XZLookupEpsilon();
+    void LookupEpsilon();
 
     /**
-     * @brief Build the table ID to composition mapping.
+     * @brief Build the table ID to Index mapping.
      */
-    void buildTableIDToComposition();
+    void buildTableIDToIndex();
 
 public:
     /**
@@ -142,12 +142,12 @@ public:
     ~OpatIO();
 
     /**
-     * @brief Get a table by X and Z values.
-     * @param X The X composition value.
-     * @param Z The Z composition value.
+     * @brief Get a table by index vector
+     * @param index The index vector associated with the table to retrieve.
+     * @throw std::out_of_range if the index is not found.
      * @return The OPAT table.
      */
-    OPATTable getTable(double X, double Z);
+    OPATTable getTable(std::vector<double> index);
 
     /**
      * @brief Set the maximum depth of the table queue.
@@ -195,11 +195,10 @@ public:
 
     /**
      * @brief Print a table by X and Z values.
-     * @param X The X composition value.
-     * @param Z The Z composition value.
+     * @param index The index vector associated with the table to print.
      * @param truncateDigits Number of digits to truncate.
      */
-    void printTable(double X, double Z, uint32_t truncateDigits=5);
+    void printTable(std::vector<double> index, uint32_t truncateDigits=5);
 
     /**
      * @brief Get the table index.
@@ -213,60 +212,58 @@ public:
      */
     Header getHeader();
 
-    /**
-     * @brief Get the closest tables by X value.
-     * @param X The X composition value.
-     * @param ZExact The exact Z composition value.
-     * @param C The C composition value (default is 0).
-     * @param O The O composition value (default is 0).
-     * @param numTables The number of closest tables to retrieve (default is 1).
-     * @return A vector of OPAT tables.
-     */
-    std::vector<OPATTable> getClosestXTables(double X, double ZExact, double C=0, double O=0, int numTables=1);
-
-    /**
-     * @brief Get the closest tables by Z value.
-     * @param XExact The exact X composition value.
-     * @param Z The Z composition value.
-     * @param C The C composition value (default is 0).
-     * @param O The O composition value (default is 0).
-     * @param numTables The number of closest tables to retrieve (default is 1).
-     * @return A vector of OPAT tables.
-     */
-    std::vector<OPATTable> getClosestZTables(double XExact, double Z, double C=0, double O=0, int numTables=1);
-
-    /**
-     * @brief Get the closest tables by X and Z values.
-     * @param X The X composition value.
-     * @param Z The Z composition value.
-     * @param C The C composition value (default is 0).
-     * @param O The O composition value (default is 0).
-     * @param numTables The number of closest tables to retrieve (default is 1).
-     * @return A vector of OPAT tables.
-     */
-    std::vector<OPATTable> getClosestTables(double X, double Z, double C=0, double O=0, int numTables=1);
-
     /**
      * @brief Lookup the table ID by X and Z values.
-     * @param X The X composition value.
-     * @param Z The Z composition value.
-     * @return The table ID.
+     * @param index The index vector associated with the table to lookup.
+     * @return The table ID if index is found, otherwise -1.
      */
-    int lookupTableID(double X, double Z);
+    int lookupTableID(std::vector<double> index);
 
     /**
      * @brief Lookup the closest table ID by X and Z values.
-     * @param X The X composition value.
-     * @param Z The Z composition value.
+     * @param index The index vector associated with the table to lookup.
      * @return The closest table ID.
      */
-    int lookupClosestTableID(double X, double Z);
+    int lookupClosestTableID(std::vector<double> index);
 
     /**
      * @brief Get the version of the OPAT file format.
      * @return The version of the OPAT file format.
      */
     uint16_t getOPATVersion();
+
+    /**
+     * @brief Get size of the index vector per table in the OPAT file.
+     * @return The size of the index vector per table.
+     */
+    uint16_t getNumIndex();
+
+    /**
+     * @brief Get the index vector for a given table ID.
+     * @param index The index vector associated with the table to retrieve.
+     * @return The full TableIndex entry for the table
+     */
+    TableIndex getTableIndex(std::vector<double> index);
+
+    /**
+     * @brief Get the checksum (sha256) for a given table ID.
+     * @param tableID The ID of the table.
+     * @return The checksum vector for the table.
+     */
+    std::vector<unsigned char> computeChecksum(int tableID);
+
+    /**
+     * @brief Get the checksum (sha256) for a given index vector.
+     * @param index The index vector associated with the table to retrieve.
+     * @return The checksum vector for the table.
+     */
+    std::vector<unsigned char> computeChecksum(std::vector<double> index);
+
+    /**
+     * @brief Validate the checksum of all tables.
+     * @return True if all checksum are valid, false otherwise.
+     */
+    bool validateAll();
 };
 
 #endif

src/probe/meson.build

@@ -0,0 +1,22 @@
+# Define the library
+probe_sources = files(
+  'private/probe.cpp',
+)
+
+probe_headers = files(
+    'public/probe.h'
+)
+
+# Define the liblogger library so it can be linked against by other parts of the build system
+libprobe = static_library('probe',
+    probe_sources, 
+    include_directories: include_directories('public'),
+    cpp_args: ['-fvisibility=default'],
+    install : true,
+    dependencies: [config_dep, mfem_dep, quill_dep]
+    )
+
+probe_dep = declare_dependency(
+    include_directories: include_directories('public'),
+    link_with: libprobe,
+)

src/probe/private/probe.cpp

@@ -0,0 +1,99 @@
+#include "quill/Backend.h"
+#include "quill/Frontend.h"
+#include "quill/Logger.h"
+#include "quill/sinks/ConsoleSink.h"
+#include "quill/sinks/FileSink.h" 
+#include <stdexcept>              
+#include <string>
+#include <iostream>
+#include <chrono>
+
+#include "mfem.hpp"
+
+#include "config.h" 
+#include "probe.h" 
+
+
+namespace Probe {
+
+void pause() {
+  std::cout << "Execution paused. Please press enter to continue..."
+            << std::endl; // Use endl to flush
+  std::cin.get();
+}
+
+void wait(int seconds) {
+  std::this_thread::sleep_for(std::chrono::seconds(seconds));
+}
+
+void glVisView(mfem::GridFunction& u, mfem::Mesh& mesh,
+               const std::string& windowTitle) {
+  Config& config = Config::getInstance();
+  if (config.get<bool>("Probe:GLVis:Visualization", true)) {
+    std::string vishost = config.get<std::string>("Probe:GLVis:Host", "localhost");
+    int visport = config.get<int>("Probe:GLVis:Port", 19916); // Changed default port
+    mfem::socketstream sol_sock(vishost.c_str(), visport);
+    sol_sock.precision(8);
+    sol_sock << "solution\n" << mesh << u
+             << "window_title '" << windowTitle << "'\n" << std::flush; // Added title
+  }
+}
+
+LogManager::LogManager() {
+  Config& config = Config::getInstance();
+  quill::Backend::start();
+  auto CLILogger = quill::Frontend::create_or_get_logger(
+      "root",
+      quill::Frontend::create_or_get_sink<quill::ConsoleSink>("sink_id_1"));
+
+  newFileLogger(config.get<std::string>("Probe:LogManager:DefaultLogName", "4DSSE.log"), "log");
+  loggerMap.emplace("stdout", CLILogger);
+}
+
+LogManager::~LogManager() = default;
+
+quill::Logger* LogManager::getLogger(const std::string& loggerName) {
+  auto it = loggerMap.find(loggerName); // Find *once*
+  if (it == loggerMap.end()) {
+    throw std::runtime_error("Cannot find logger " + loggerName);
+  }
+  return it->second; // Return the raw pointer from the shared_ptr
+}
+
+std::vector<std::string> LogManager::getLoggerNames() {
+  std::vector<std::string> loggerNames;
+  loggerNames.reserve(loggerMap.size());
+  for (const auto& pair : loggerMap) { // Use range-based for loop and const auto&
+    loggerNames.push_back(pair.first);
+  }
+  return loggerNames;
+}
+
+std::vector<quill::Logger*> LogManager::getLoggers() {
+  std::vector<quill::Logger*> loggers;
+  loggers.reserve(loggerMap.size());
+  for (const auto& pair : loggerMap) {
+     loggers.push_back(pair.second); // Get the raw pointer
+  }
+  return loggers;
+}
+
+quill::Logger* LogManager::newFileLogger(const std::string& filename,
+                                        const std::string& loggerName) {
+  auto file_sink = quill::Frontend::create_or_get_sink<quill::FileSink>(
+      filename,
+      []() {
+        quill::FileSinkConfig cfg;
+        cfg.set_open_mode('w');
+        return cfg;
+      }(),
+      quill::FileEventNotifier{});
+  // Get the raw pointer.
+  quill::Logger* rawLogger = quill::Frontend::create_or_get_logger(loggerName, std::move(file_sink));
+
+  // Create a shared_ptr from the raw pointer.
+  loggerMap.emplace(loggerName, rawLogger);
+  return rawLogger;
+}
+
+} // namespace Probe

src/probe/public/probe.h

@@ -0,0 +1,99 @@
+//=== Probe.h ===
+#ifndef PROBE_H
+#define PROBE_H
+
+#include <string>
+#include <map>
+#include <vector>
+#include <iostream>
+
+#include "mfem.hpp"
+#include "quill/Logger.h"
+
+/**
+ * @brief The Probe namespace contains utility functions for debugging and logging.
+ */
+namespace Probe {
+  /**
+   * @brief Pause the execution and wait for user input.
+   */
+  void pause();
+
+  /**
+   * @brief Wait for a specified number of seconds.
+   * @param seconds The number of seconds to wait.
+   */
+  void wait(int seconds);
+
+  /**
+   * @brief Visualize a solution using GLVis.
+   * @param u The GridFunction to visualize.
+   * @param mesh The mesh associated with the GridFunction.
+   * @param windowTitle The title of the visualization window.
+   */
+  void glVisView(mfem::GridFunction& u, mfem::Mesh& mesh,
+                const std::string& windowTitle = "solution");
+
+  /**
+   * @brief Class to manage logging operations.
+   */
+  class LogManager {
+  private:
+      /**
+       * @brief Private constructor for singleton pattern.
+       */
+      LogManager();
+
+      /**
+       * @brief Destructor.
+       */
+      ~LogManager();
+
+      // Map to store pointers to quill loggers (raw pointers as quill deals with its own memory managment in a seperated, detatched, thread)
+      std::map<std::string, quill::Logger*> loggerMap;
+
+      // Prevent copying and assignment (Rule of Zero)
+      LogManager(const LogManager&) = delete;
+      LogManager& operator=(const LogManager&) = delete;
+
+  public:
+      /**
+       * @brief Get the singleton instance of LogManager.
+       * @return The singleton instance of LogManager.
+       */
+      static LogManager& getInstance() {
+          static LogManager instance;
+          return instance;
+      }
+
+      /**
+       * @brief Get a logger by name.
+       * @param loggerName The name of the logger.
+       * @return A pointer to the logger.
+       */
+      quill::Logger* getLogger(const std::string& loggerName);
+
+      /**
+       * @brief Get the names of all loggers.
+       * @return A vector of logger names.
+       */
+      std::vector<std::string> getLoggerNames();
+
+      /**
+       * @brief Get all loggers.
+       * @return A vector of pointers to the loggers.
+       */
+      std::vector<quill::Logger*> getLoggers();
+
+      /**
+       * @brief Create a new file logger.
+       * @param filename The name of the log file.
+       * @param loggerName The name of the logger.
+       * @return A pointer to the new logger.
+       */
+      quill::Logger* newFileLogger(const std::string& filename,
+                                  const std::string& loggerName);
+  };
+
+} // namespace Probe
+#endif