Merge pull request #13 from tboudreaux/feature/opatOPALII

Updated OPAT file format and handlers to work for OPAL type I or type II opacity files

meson.build

@@ -5,6 +5,7 @@ add_project_arguments('-fvisibility=default', language: 'cpp')
 
 # Build external dependencies first so that all the embedded resources are available to the other targets
 subdir('build-config')
+subdir('subprojects/PicoSHA2')
 
 # Build the main project
 subdir('src')

mk

@@ -2,9 +2,9 @@
 
 # Check for the --noTest flag
 if [[ "$1" == "--noTest" ]]; then
-  meson setup build -Dbuild_tests=false
+  meson setup build -Dbuild_tests=false --buildtype=release
 else
-  meson setup build -Dbuild_tests=true
+  meson setup build -Dbuild_tests=true --buildtype=debug
 fi
 
 # Compile the project

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));
+    long unsigned int indexReadBytes;
-    if (file.gcount() != static_cast<std::streamsize>(header.numTables * sizeof(TableIndex))) {
+    for (uint32_t i = 0; i < header.numTables; i++) {
-        throw std::runtime_error("Error reading table index from file: " + filename);
+        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(){
+void OpatIO::buildTableIDToIndex(){
-    tableIDToComposition.clear();
+    tableIDToIndex.clear();
     int tableID = 0;
-    std::pair<double, double> comp;
+    std::vector<double> ind;
-    for (const auto &index : tableIndex) {
+    ind.resize(header.numIndex);
-        comp.first = index.X;
+    for (const auto &table : tableIndex) {
-        comp.second = index.Z;
+        ind.clear();
-        tableIDToComposition.emplace(tableID, comp);
+        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 
+     Get 10% of the minimum spacing between index values 
-     in the tableIDToComposition map. This can be used
+     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;
+    indexEpsilon.resize(header.numIndex);
-    double epsilonX, epsilonZ, xgap, zgap;
 
-    // Start these out as larger than they will ever be
+    double epsilon;
-    epsilonX = 1;
+    for (int i = 0; i < static_cast<int>(header.numIndex); i++) {
-    epsilonZ = 1;
+        epsilon = std::numeric_limits<double>::max();
-    for (const auto& pair : tableIDToComposition) {
+        for (int j = 1; j < static_cast<int>(header.numTables); j++) {
-        Xvalues.push_back(pair.second.first);
+            epsilon = std::min(epsilon, std::fabs(tableIDToIndex.at(j).at(i) - tableIDToIndex.at(j-1).at(i)));
-        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;
         }
+        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){
+int OpatIO::lookupTableID(std::vector<double> index) {
-    bool XOkay;
+    std::vector<bool> IndexOkay;
-    bool ZOkay;
+    IndexOkay.resize(header.numIndex);
     int tableID = 0;
-    for (const auto &tableMap : tableIDToComposition){
+    for (const auto &tableMap : tableIDToIndex){
-        XOkay = std::fabs(tableMap.second.first - X) < XZepsilon.first;
+        // Loop through all index values and check if they are within epsilon for that index
-        ZOkay = std::fabs(tableMap.second.second - Z) < XZepsilon.second;
+        std::fill(IndexOkay.begin(), IndexOkay.end(), false);
-        if (XOkay and ZOkay){
+        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) {
+OPATTable OpatIO::getTable(std::vector<double> index) {
-    int tableID = lookupTableID(X, Z);
+    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::cout << std::left << std::setw(10);
-        << std::setw(10) << "Z"
+    for (int i = 0; i < header.numIndex; i++) {
-        << std::setw(15) << "Byte Start"
+        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::cout << std::fixed << std::setprecision(4) << std::setw(10);
-            << std::setw(10) << index.X
+        for (int i = 0; i < header.numIndex; i++) {
-            << std::setw(10) << index.Z
+            std::cout << index.index[i] << std::setw(10);
-            << std::setw(15) << index.byteStart
+        }
+        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) {
+void OpatIO::printTable(std::vector<double> index, uint32_t truncateDigits) {
-    int tableID = lookupTableID(X, Z);
+    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
+// Get the size of the index vector used
-// std::vector<OPATTable> OpatIO::getClosestXTables(double X, double ZExact, int numTables) {
+uint16_t OpatIO::getNumIndex() {
-//     std::vector<OPATTable> closestTables;
+    return header.numIndex;
-//     // Implement logic to find closest X tables
+}
-//     return closestTables;
-// }
 
-// // Get the closest Z tables
+TableIndex OpatIO::getTableIndex(std::vector<double> index) {
-// std::vector<OPATTable> OpatIO::getClosestZTables(double XExact, double Z, int numTables) {
+    int tableID = lookupTableID(index);
-//     std::vector<OPATTable> closestTables;
+    return tableIndex.at(tableID);
-//     // Implement logic to find closest Z tables
+}
-//     return closestTables;
-// }
 
-// // Get the closest tables
+std::vector<unsigned char> OpatIO::computeChecksum(int tableID) {
-// std::vector<OPATTable> OpatIO::getClosestTables(double X, double Z, int numTables) {
+    OPATTable table = getTable(tableID);
-//     std::vector<OPATTable> closestTables;
+    std::vector<std::vector<double>> logKappa = table.logKappa;
-//     // Implement logic to find closest tables
+    std::vector<double> flatData(logKappa.size() * logKappa.size());
-//     return closestTables;
+    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
+    std::vector<double> index; ///< Index vector for associated table
-    double Z;                ///< Z composition value
     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::map<int, std::vector<double>> tableIDToIndex; ///< Map to store table ID to indexing
-    std::pair<double, double> XZepsilon; ///< Epsilon values for X and Z
+    std::vector<double> indexEpsilon; ///< Epsilon values for each index
-    int maxQDepth = 10; ///< Maximum depth of the table queue
+    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.
+     * @brief Get a table by index vector
-     * @param X The X composition value.
+     * @param index The index vector associated with the table to retrieve.
-     * @param Z The Z composition value.
+     * @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 index The index vector associated with the table to print.
-     * @param Z The Z composition value.
      * @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 index The index vector associated with the table to lookup.
-     * @param Z The Z composition value.
+     * @return The table ID if index is found, otherwise -1.
-     * @return The table ID.
      */
-    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 index The index vector associated with the table to lookup.
-     * @param Z The Z composition value.
      * @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

subprojects/PicoSHA2/LICENSE

@@ -0,0 +1,21 @@
+# MIT License
+
+Copyright (c) 2017 okdshin
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.

subprojects/PicoSHA2/README.md

@@ -0,0 +1,138 @@
+# PicoSHA2 - a C++ SHA256 hash generator
+
+Copyright © 2017 okdshin
+
+## Introduction
+
+PicoSHA2 is a tiny SHA256 hash generator for C++ with following properties:
+
+- header-file only
+- no external dependencies (only uses standard C++ libraries)
+- STL-friendly
+- licensed under MIT License
+
+## Generating SHA256 hash and hash hex string
+
+```c++
+// any STL sequantial container (vector, list, dequeue...)
+std::string src_str = "The quick brown fox jumps over the lazy dog";
+
+std::vector<unsigned char> hash(picosha2::k_digest_size);
+picosha2::hash256(src_str.begin(), src_str.end(), hash.begin(), hash.end());
+
+std::string hex_str = picosha2::bytes_to_hex_string(hash.begin(), hash.end());
+```
+
+## Generating SHA256 hash and hash hex string from byte stream
+
+```c++
+picosha2::hash256_one_by_one hasher;
+...
+hasher.process(block.begin(), block.end());
+...
+hasher.finish();
+
+std::vector<unsigned char> hash(picosha2::k_digest_size);
+hasher.get_hash_bytes(hash.begin(), hash.end());
+
+std::string hex_str = picosha2::get_hash_hex_string(hasher);
+```
+
+The file `example/interactive_hasher.cpp` has more detailed information.
+
+## Generating SHA256 hash from a binary file
+
+```c++
+std::ifstream f("file.txt", std::ios::binary);
+std::vector<unsigned char> s(picosha2::k_digest_size);
+picosha2::hash256(f, s.begin(), s.end());
+```
+
+This `hash256` may use less memory than reading whole of the file.
+
+## Generating SHA256 hash hex string from std::string
+
+```c++
+std::string src_str = "The quick brown fox jumps over the lazy dog";
+std::string hash_hex_str;
+picosha2::hash256_hex_string(src_str, hash_hex_str);
+std::cout << hash_hex_str << std::endl;
+//this output is "d7a8fbb307d7809469ca9abcb0082e4f8d5651e46d3cdb762d02d0bf37c9e592"
+```
+
+```c++
+std::string src_str = "The quick brown fox jumps over the lazy dog";
+std::string hash_hex_str = picosha2::hash256_hex_string(src_str);
+std::cout << hash_hex_str << std::endl;
+//this output is "d7a8fbb307d7809469ca9abcb0082e4f8d5651e46d3cdb762d02d0bf37c9e592"
+```
+
+```c++
+std::string src_str = "The quick brown fox jumps over the lazy dog.";//add '.'
+std::string hash_hex_str = picosha2::hash256_hex_string(src_str.begin(), src_str.end());
+std::cout << hash_hex_str << std::endl;
+//this output is "ef537f25c895bfa782526529a9b63d97aa631564d5d789c2b765448c8635fb6c"
+```
+
+## Generating SHA256 hash hex string from byte sequence
+
+```c++
+std::vector<unsigned char> src_vect(...);
+std::string hash_hex_str;
+picosha2::hash256_hex_string(src_vect, hash_hex_str);
+```
+
+```c++
+std::vector<unsigned char> src_vect(...);
+std::string hash_hex_str = picosha2::hash256_hex_string(src_vect);
+```
+
+```c++
+unsigned char src_c_array[picosha2::k_digest_size] = {...};
+std::string hash_hex_str;
+picosha2::hash256_hex_string(src_c_array, src_c_array+picosha2::k_digest_size, hash_hex_str);
+```
+
+```c++
+unsigned char src_c_array[picosha2::k_digest_size] = {...};
+std::string hash_hex_str = picosha2::hash256_hex_string(src_c_array, src_c_array+picosha2::k_digest_size);
+```
+
+
+## Generating SHA256 hash byte sequence from STL sequential container
+
+```c++
+//any STL sequantial container (vector, list, dequeue...)
+std::string src_str = "The quick brown fox jumps over the lazy dog";
+
+//any STL sequantial containers (vector, list, dequeue...)
+std::vector<unsigned char> hash(picosha2::k_digest_size);
+
+// in: container, out: container
+picosha2::hash256(src_str, hash);
+```
+
+```c++
+//any STL sequantial container (vector, list, dequeue...)
+std::string src_str = "The quick brown fox jumps over the lazy dog";
+
+//any STL sequantial containers (vector, list, dequeue...)
+std::vector<unsigned char> hash(picosha2::k_digest_size);
+
+// in: iterator pair, out: contaner
+picosha2::hash256(src_str.begin(), src_str.end(), hash);
+```
+
+```c++
+std::string src_str = "The quick brown fox jumps over the lazy dog";
+unsigned char hash_byte_c_array[picosha2::k_digest_size];
+// in: container, out: iterator(pointer) pair
+picosha2::hash256(src_str, hash_byte_c_array, hash_byte_c_array+picosha2::k_digest_size);
+```
+
+```c++
+std::string src_str = "The quick brown fox jumps over the lazy dog";
+std::vector<unsigned char> hash(picosha2::k_digest_size);
+// in: iterator pair, out: iterator pair
+picosha2::hash256(src_str.begin(), src_str.end(), hash.begin(), hash.end());
+```

subprojects/PicoSHA2/meson.build

@@ -0,0 +1,4 @@
+
+picosha2_dep = declare_dependency(
+  include_directories: include_directories('public')
+)

subprojects/PicoSHA2/public/picosha2.h

@@ -0,0 +1,377 @@
+/*
+The MIT License (MIT)
+
+Copyright (C) 2017 okdshin
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
+*/
+#ifndef PICOSHA2_H
+#define PICOSHA2_H
+// picosha2:20140213
+
+#ifndef PICOSHA2_BUFFER_SIZE_FOR_INPUT_ITERATOR
+#define PICOSHA2_BUFFER_SIZE_FOR_INPUT_ITERATOR \
+    1048576  //=1024*1024: default is 1MB memory
+#endif
+
+#include <algorithm>
+#include <cassert>
+#include <iterator>
+#include <sstream>
+#include <vector>
+#include <fstream>
+namespace picosha2 {
+typedef unsigned long word_t;
+typedef unsigned char byte_t;
+
+static const size_t k_digest_size = 32;
+
+namespace detail {
+inline byte_t mask_8bit(byte_t x) { return x & 0xff; }
+
+inline word_t mask_32bit(word_t x) { return x & 0xffffffff; }
+
+const word_t add_constant[64] = {
+    0x428a2f98, 0x71374491, 0xb5c0fbcf, 0xe9b5dba5, 0x3956c25b, 0x59f111f1,
+    0x923f82a4, 0xab1c5ed5, 0xd807aa98, 0x12835b01, 0x243185be, 0x550c7dc3,
+    0x72be5d74, 0x80deb1fe, 0x9bdc06a7, 0xc19bf174, 0xe49b69c1, 0xefbe4786,
+    0x0fc19dc6, 0x240ca1cc, 0x2de92c6f, 0x4a7484aa, 0x5cb0a9dc, 0x76f988da,
+    0x983e5152, 0xa831c66d, 0xb00327c8, 0xbf597fc7, 0xc6e00bf3, 0xd5a79147,
+    0x06ca6351, 0x14292967, 0x27b70a85, 0x2e1b2138, 0x4d2c6dfc, 0x53380d13,
+    0x650a7354, 0x766a0abb, 0x81c2c92e, 0x92722c85, 0xa2bfe8a1, 0xa81a664b,
+    0xc24b8b70, 0xc76c51a3, 0xd192e819, 0xd6990624, 0xf40e3585, 0x106aa070,
+    0x19a4c116, 0x1e376c08, 0x2748774c, 0x34b0bcb5, 0x391c0cb3, 0x4ed8aa4a,
+    0x5b9cca4f, 0x682e6ff3, 0x748f82ee, 0x78a5636f, 0x84c87814, 0x8cc70208,
+    0x90befffa, 0xa4506ceb, 0xbef9a3f7, 0xc67178f2};
+
+const word_t initial_message_digest[8] = {0x6a09e667, 0xbb67ae85, 0x3c6ef372,
+                                          0xa54ff53a, 0x510e527f, 0x9b05688c,
+                                          0x1f83d9ab, 0x5be0cd19};
+
+inline word_t ch(word_t x, word_t y, word_t z) { return (x & y) ^ ((~x) & z); }
+
+inline word_t maj(word_t x, word_t y, word_t z) {
+    return (x & y) ^ (x & z) ^ (y & z);
+}
+
+inline word_t rotr(word_t x, std::size_t n) {
+    assert(n < 32);
+    return mask_32bit((x >> n) | (x << (32 - n)));
+}
+
+inline word_t bsig0(word_t x) { return rotr(x, 2) ^ rotr(x, 13) ^ rotr(x, 22); }
+
+inline word_t bsig1(word_t x) { return rotr(x, 6) ^ rotr(x, 11) ^ rotr(x, 25); }
+
+inline word_t shr(word_t x, std::size_t n) {
+    assert(n < 32);
+    return x >> n;
+}
+
+inline word_t ssig0(word_t x) { return rotr(x, 7) ^ rotr(x, 18) ^ shr(x, 3); }
+
+inline word_t ssig1(word_t x) { return rotr(x, 17) ^ rotr(x, 19) ^ shr(x, 10); }
+
+template <typename RaIter1, typename RaIter2>
+void hash256_block(RaIter1 message_digest, RaIter2 first, RaIter2 last) {
+    assert(first + 64 == last);
+    static_cast<void>(last);  // for avoiding unused-variable warning
+    word_t w[64];
+    std::fill(w, w + 64, word_t(0));
+    for (std::size_t i = 0; i < 16; ++i) {
+        w[i] = (static_cast<word_t>(mask_8bit(*(first + i * 4))) << 24) |
+               (static_cast<word_t>(mask_8bit(*(first + i * 4 + 1))) << 16) |
+               (static_cast<word_t>(mask_8bit(*(first + i * 4 + 2))) << 8) |
+               (static_cast<word_t>(mask_8bit(*(first + i * 4 + 3))));
+    }
+    for (std::size_t i = 16; i < 64; ++i) {
+        w[i] = mask_32bit(ssig1(w[i - 2]) + w[i - 7] + ssig0(w[i - 15]) +
+                          w[i - 16]);
+    }
+
+    word_t a = *message_digest;
+    word_t b = *(message_digest + 1);
+    word_t c = *(message_digest + 2);
+    word_t d = *(message_digest + 3);
+    word_t e = *(message_digest + 4);
+    word_t f = *(message_digest + 5);
+    word_t g = *(message_digest + 6);
+    word_t h = *(message_digest + 7);
+
+    for (std::size_t i = 0; i < 64; ++i) {
+        word_t temp1 = h + bsig1(e) + ch(e, f, g) + add_constant[i] + w[i];
+        word_t temp2 = bsig0(a) + maj(a, b, c);
+        h = g;
+        g = f;
+        f = e;
+        e = mask_32bit(d + temp1);
+        d = c;
+        c = b;
+        b = a;
+        a = mask_32bit(temp1 + temp2);
+    }
+    *message_digest += a;
+    *(message_digest + 1) += b;
+    *(message_digest + 2) += c;
+    *(message_digest + 3) += d;
+    *(message_digest + 4) += e;
+    *(message_digest + 5) += f;
+    *(message_digest + 6) += g;
+    *(message_digest + 7) += h;
+    for (std::size_t i = 0; i < 8; ++i) {
+        *(message_digest + i) = mask_32bit(*(message_digest + i));
+    }
+}
+
+}  // namespace detail
+
+template <typename InIter>
+void output_hex(InIter first, InIter last, std::ostream& os) {
+    os.setf(std::ios::hex, std::ios::basefield);
+    while (first != last) {
+        os.width(2);
+        os.fill('0');
+        os << static_cast<unsigned int>(*first);
+        ++first;
+    }
+    os.setf(std::ios::dec, std::ios::basefield);
+}
+
+template <typename InIter>
+void bytes_to_hex_string(InIter first, InIter last, std::string& hex_str) {
+    std::ostringstream oss;
+    output_hex(first, last, oss);
+    hex_str.assign(oss.str());
+}
+
+template <typename InContainer>
+void bytes_to_hex_string(const InContainer& bytes, std::string& hex_str) {
+    bytes_to_hex_string(bytes.begin(), bytes.end(), hex_str);
+}
+
+template <typename InIter>
+std::string bytes_to_hex_string(InIter first, InIter last) {
+    std::string hex_str;
+    bytes_to_hex_string(first, last, hex_str);
+    return hex_str;
+}
+
+template <typename InContainer>
+std::string bytes_to_hex_string(const InContainer& bytes) {
+    std::string hex_str;
+    bytes_to_hex_string(bytes, hex_str);
+    return hex_str;
+}
+
+class hash256_one_by_one {
+   public:
+    hash256_one_by_one() { init(); }
+
+    void init() {
+        buffer_.clear();
+        std::fill(data_length_digits_, data_length_digits_ + 4, word_t(0));
+        std::copy(detail::initial_message_digest,
+                  detail::initial_message_digest + 8, h_);
+    }
+
+    template <typename RaIter>
+    void process(RaIter first, RaIter last) {
+        add_to_data_length(static_cast<word_t>(std::distance(first, last)));
+        std::copy(first, last, std::back_inserter(buffer_));
+        std::size_t i = 0;
+        for (; i + 64 <= buffer_.size(); i += 64) {
+            detail::hash256_block(h_, buffer_.begin() + i,
+                                  buffer_.begin() + i + 64);
+        }
+        buffer_.erase(buffer_.begin(), buffer_.begin() + i);
+    }
+
+    void finish() {
+        byte_t temp[64];
+        std::fill(temp, temp + 64, byte_t(0));
+        std::size_t remains = buffer_.size();
+        std::copy(buffer_.begin(), buffer_.end(), temp);
+        temp[remains] = 0x80;
+
+        if (remains > 55) {
+            std::fill(temp + remains + 1, temp + 64, byte_t(0));
+            detail::hash256_block(h_, temp, temp + 64);
+            std::fill(temp, temp + 64 - 4, byte_t(0));
+        } else {
+            std::fill(temp + remains + 1, temp + 64 - 4, byte_t(0));
+        }
+
+        write_data_bit_length(&(temp[56]));
+        detail::hash256_block(h_, temp, temp + 64);
+    }
+
+    template <typename OutIter>
+    void get_hash_bytes(OutIter first, OutIter last) const {
+        for (const word_t* iter = h_; iter != h_ + 8; ++iter) {
+            for (std::size_t i = 0; i < 4 && first != last; ++i) {
+                *(first++) = detail::mask_8bit(
+                    static_cast<byte_t>((*iter >> (24 - 8 * i))));
+            }
+        }
+    }
+
+   private:
+    void add_to_data_length(word_t n) {
+        word_t carry = 0;
+        data_length_digits_[0] += n;
+        for (std::size_t i = 0; i < 4; ++i) {
+            data_length_digits_[i] += carry;
+            if (data_length_digits_[i] >= 65536u) {
+                carry = data_length_digits_[i] >> 16;
+                data_length_digits_[i] &= 65535u;
+            } else {
+                break;
+            }
+        }
+    }
+    void write_data_bit_length(byte_t* begin) {
+        word_t data_bit_length_digits[4];
+        std::copy(data_length_digits_, data_length_digits_ + 4,
+                  data_bit_length_digits);
+
+        // convert byte length to bit length (multiply 8 or shift 3 times left)
+        word_t carry = 0;
+        for (std::size_t i = 0; i < 4; ++i) {
+            word_t before_val = data_bit_length_digits[i];
+            data_bit_length_digits[i] <<= 3;
+            data_bit_length_digits[i] |= carry;
+            data_bit_length_digits[i] &= 65535u;
+            carry = (before_val >> (16 - 3)) & 65535u;
+        }
+
+        // write data_bit_length
+        for (int i = 3; i >= 0; --i) {
+            (*begin++) = static_cast<byte_t>(data_bit_length_digits[i] >> 8);
+            (*begin++) = static_cast<byte_t>(data_bit_length_digits[i]);
+        }
+    }
+    std::vector<byte_t> buffer_;
+    word_t data_length_digits_[4];  // as 64bit integer (16bit x 4 integer)
+    word_t h_[8];
+};
+
+inline void get_hash_hex_string(const hash256_one_by_one& hasher,
+                                std::string& hex_str) {
+    byte_t hash[k_digest_size];
+    hasher.get_hash_bytes(hash, hash + k_digest_size);
+    return bytes_to_hex_string(hash, hash + k_digest_size, hex_str);
+}
+
+inline std::string get_hash_hex_string(const hash256_one_by_one& hasher) {
+    std::string hex_str;
+    get_hash_hex_string(hasher, hex_str);
+    return hex_str;
+}
+
+namespace impl {
+template <typename RaIter, typename OutIter>
+void hash256_impl(RaIter first, RaIter last, OutIter first2, OutIter last2, int,
+                  std::random_access_iterator_tag) {
+    hash256_one_by_one hasher;
+    // hasher.init();
+    hasher.process(first, last);
+    hasher.finish();
+    hasher.get_hash_bytes(first2, last2);
+}
+
+template <typename InputIter, typename OutIter>
+void hash256_impl(InputIter first, InputIter last, OutIter first2,
+                  OutIter last2, int buffer_size, std::input_iterator_tag) {
+    std::vector<byte_t> buffer(buffer_size);
+    hash256_one_by_one hasher;
+    // hasher.init();
+    while (first != last) {
+        int size = buffer_size;
+        for (int i = 0; i != buffer_size; ++i, ++first) {
+            if (first == last) {
+                size = i;
+                break;
+            }
+            buffer[i] = *first;
+        }
+        hasher.process(buffer.begin(), buffer.begin() + size);
+    }
+    hasher.finish();
+    hasher.get_hash_bytes(first2, last2);
+}
+}
+
+template <typename InIter, typename OutIter>
+void hash256(InIter first, InIter last, OutIter first2, OutIter last2,
+             int buffer_size = PICOSHA2_BUFFER_SIZE_FOR_INPUT_ITERATOR) {
+    picosha2::impl::hash256_impl(
+        first, last, first2, last2, buffer_size,
+        typename std::iterator_traits<InIter>::iterator_category());
+}
+
+template <typename InIter, typename OutContainer>
+void hash256(InIter first, InIter last, OutContainer& dst) {
+    hash256(first, last, dst.begin(), dst.end());
+}
+
+template <typename InContainer, typename OutIter>
+void hash256(const InContainer& src, OutIter first, OutIter last) {
+    hash256(src.begin(), src.end(), first, last);
+}
+
+template <typename InContainer, typename OutContainer>
+void hash256(const InContainer& src, OutContainer& dst) {
+    hash256(src.begin(), src.end(), dst.begin(), dst.end());
+}
+
+template <typename InIter>
+void hash256_hex_string(InIter first, InIter last, std::string& hex_str) {
+    byte_t hashed[k_digest_size];
+    hash256(first, last, hashed, hashed + k_digest_size);
+    std::ostringstream oss;
+    output_hex(hashed, hashed + k_digest_size, oss);
+    hex_str.assign(oss.str());
+}
+
+template <typename InIter>
+std::string hash256_hex_string(InIter first, InIter last) {
+    std::string hex_str;
+    hash256_hex_string(first, last, hex_str);
+    return hex_str;
+}
+
+inline void hash256_hex_string(const std::string& src, std::string& hex_str) {
+    hash256_hex_string(src.begin(), src.end(), hex_str);
+}
+
+template <typename InContainer>
+void hash256_hex_string(const InContainer& src, std::string& hex_str) {
+    hash256_hex_string(src.begin(), src.end(), hex_str);
+}
+
+template <typename InContainer>
+std::string hash256_hex_string(const InContainer& src) {
+    return hash256_hex_string(src.begin(), src.end());
+}
+template<typename OutIter>void hash256(std::ifstream& f, OutIter first, OutIter last){
+    hash256(std::istreambuf_iterator<char>(f), std::istreambuf_iterator<char>(), first,last);
+
+}
+}// namespace picosha2
+#endif  // PICOSHA2_H

tests/opatIO/meson.build

@@ -13,7 +13,7 @@ foreach test_file : test_sources
   test_exe = executable(
       exe_name,
       test_file,
-      dependencies: [gtest_dep],
+      dependencies: [gtest_dep, picosha2_dep],
       include_directories: include_directories('../../src/opatIO/public'),
       link_with: libopatIO,  # Link the dobj library
       install_rpath: '@loader_path/../../src'  # Ensure runtime library path resolves correctly

tests/opatIO/opatIOTest.cpp

@@ -5,8 +5,10 @@
 #include <vector>
 #include <set>
 #include <sstream>
+#include <cstring>
+#include "picosha2.h"
 
-std::string EXAMPLE_FILENAME = std::string(getenv("MESON_SOURCE_ROOT")) + "/tests/opatIO/test.opat";
+std::string EXAMPLE_FILENAME = std::string(getenv("MESON_SOURCE_ROOT")) + "/tests/opatIO/synthetic_tables.opat";
 
 /**
  * @file opatIOTest.cpp
@@ -32,6 +34,9 @@ TEST_F(opatIOTest, Constructor) {
     OpatIO opatIO(EXAMPLE_FILENAME);
 }
 
+/**
+ * @test Verify the header is read correctly.
+ */
 TEST_F(opatIOTest, Header) {
     OpatIO opatIO(EXAMPLE_FILENAME);
     Header header = opatIO.getHeader();
@@ -39,42 +44,60 @@ TEST_F(opatIOTest, Header) {
     EXPECT_EQ(header.numTables, 20);
     EXPECT_EQ(header.headerSize, 256);
     EXPECT_EQ(header.indexOffset, 416416);
-    EXPECT_EQ(std::string(header.creationDate), "Feb 16, 2025");
+    EXPECT_EQ(std::string(header.creationDate), "Feb 17, 2025");
-    EXPECT_EQ(std::string(header.sourceInfo), "no source provided by user");
+    EXPECT_EQ(std::string(header.sourceInfo), "utils/opatio/utils/mkTestData.py");
-    EXPECT_EQ(std::string(header.comment), "default header");
+    EXPECT_EQ(std::string(header.comment), "Synthetic Opacity Tables");
+    EXPECT_EQ(header.numIndex, 2);
 }
 
+/**
+ * @test Verify the number of index values is correct. Also check the byte position and index vector
+ */
 TEST_F(opatIOTest, TableIndex) {
     OpatIO opatIO(EXAMPLE_FILENAME);
     std::vector<TableIndex> tableIndex = opatIO.getTableIndex();
     EXPECT_EQ(tableIndex.size(), 20);
-    EXPECT_EQ(tableIndex[0].X, 0.1);
+    EXPECT_EQ(tableIndex[0].index.at(0), 0.1);
-    EXPECT_EQ(tableIndex[0].Z, 0.001);
+    EXPECT_EQ(tableIndex[0].index.at(1), 0.001);
     EXPECT_EQ(tableIndex[0].byteStart, 256);
     EXPECT_EQ(tableIndex[0].byteEnd, 21064);
 }
 
+/**
+ * @test Verify the maxQDepth (for caching)
+ */
 TEST_F(opatIOTest, MaxQDepth) {
     OpatIO opatIO(EXAMPLE_FILENAME);
-    EXPECT_EQ(opatIO.getMaxQDepth(), 10);
+    EXPECT_EQ(opatIO.getMaxQDepth(), 20);
     opatIO.setMaxQDepth(5);
     EXPECT_EQ(opatIO.getMaxQDepth(), 5);
 }
 
+/**
+ * @test Verify the Unload function
+ */
 TEST_F(opatIOTest, Unload){
     OpatIO opatIO(EXAMPLE_FILENAME);
     EXPECT_NO_THROW(opatIO.unload());
     EXPECT_FALSE(opatIO.isLoaded());
 }
 
+/**
+ * @test Verify the lookupTableID function
+ */
 TEST_F(opatIOTest, LookupTableID) {
     OpatIO opatIO(EXAMPLE_FILENAME);
-    EXPECT_EQ(opatIO.lookupTableID(0.321053, 0.0116842), 7);
+    std::vector<double> index = {0.321053, 0.0116842};
+    EXPECT_EQ(opatIO.lookupTableID(index), 7);
 }
 
+/**
+ * @test Verify the GetTable function by checking the first 2x2 square of the table
+ */
 TEST_F(opatIOTest, GetTable) {
     OpatIO opatIO(EXAMPLE_FILENAME);
-    OPATTable tab = opatIO.getTable(0.1, 0.001);
+    std::vector<double> index = {0.1, 0.001};
+    OPATTable tab = opatIO.getTable(index);
     EXPECT_EQ(tab.N_R, 50);
     EXPECT_EQ(tab.N_T, 50);
     EXPECT_DOUBLE_EQ(tab.logR[0], -8.0);
@@ -83,4 +106,28 @@ TEST_F(opatIOTest, GetTable) {
     EXPECT_DOUBLE_EQ(tab.logKappa[0][1], 1.8028572256396647);
     EXPECT_DOUBLE_EQ(tab.logKappa[1][0], 1.8783385110582342);
     EXPECT_DOUBLE_EQ(tab.logKappa[1][1], 1.1005312934444582);
+}
+
+
+
+/**
+ * @test Verify the GetTable function by computing the checksum of the first table and 
+ * comparing it to the stored checksum
+ */
+TEST_F(opatIOTest, Checksum) {
+    OpatIO opatIO(EXAMPLE_FILENAME);
+    std::vector<double> index = {0.1, 0.001};
+    TableIndex tableIndex = opatIO.getTableIndex(index);
+    std::vector<unsigned char> hash = opatIO.computeChecksum(index);
+    std::string hexRepr = picosha2::bytes_to_hex_string(hash);
+
+    std::vector<unsigned char> storedHashVec(tableIndex.sha256, tableIndex.sha256 + 32);
+    std::string storedHexRepr = picosha2::bytes_to_hex_string(storedHashVec);
+    EXPECT_EQ(hexRepr, storedHexRepr);
+}
+
+TEST_F(opatIOTest, ChecksumAll) {
+    OpatIO opatIO(EXAMPLE_FILENAME);
+    opatIO.setMaxQDepth(5);
+    EXPECT_TRUE(opatIO.validateAll());
 }

utils/opatio/src/opatio/opat/opat.py

@@ -24,6 +24,7 @@ class Header:
     creationDate: str    #< Creation date of the file
     sourceInfo: str      #< Source information
     comment: str         #< Comment section
+    numIndex: int        #< Number of values to use when indexing table
     reserved: bytes      #< Reserved for future use
 
 @dataclass
@@ -31,8 +32,7 @@ class TableIndex:
     """
     @brief Structure to hold the index information of a table in an OPAT file.
     """
-    X: float             #< X composition value
+    index: List[float]   #< Index values of the table
-    Z: float             #< Z composition value
     byteStart: int       #< Byte start position of the table
     byteEnd: int         #< Byte end position of the table
     sha256: bytes        #< SHA-256 hash of the table data
@@ -57,7 +57,8 @@ defaultHeader = Header(
     creationDate=datetime.now().strftime("%b %d, %Y"),
     sourceInfo="no source provided by user",
     comment="default header",
-    reserved=b"\x00" * 26
+    numIndex=2,
+    reserved=b"\x00" * 24
 )
 
 class OpatIO:
@@ -166,13 +167,13 @@ class OpatIO:
             raise TypeError(f"{name} must be a non-empty 1D array or iterable")
 
     @staticmethod
-    def compute_checksum(data: bytes) -> bytes:
+    def compute_checksum(data: np.ndarray) -> bytes:
         """
         @brief Compute the SHA-256 checksum of the given data.
         @param data The data to compute the checksum for.
         @return The SHA-256 checksum.
         """
-        return hashlib.sha256(data).digest()
+        return hashlib.sha256(data.tobytes()).digest()
 
     def set_version(self, version: int) -> int:
         """
@@ -206,17 +207,29 @@ class OpatIO:
             raise TypeError(f"comment string ({comment}) is too long ({len(comment)}). Max length is 128")
         self.header.comment = comment
         return self.header.comment
+    
+    def set_numIndex(self, numIndex: int) -> int:
+        """
+        @brief Set the number of values to use when indexing table.
+        @param numIndex The number of values to use when indexing table.
+        @return The set number of values to use when indexing table.
+        """
+        if numIndex < 1:
+            raise ValueError(f"numIndex must be greater than 0! It is currently {numIndex}")
+        self.header.numIndex = numIndex
+        return self.header.numIndex
 
-    def add_table(self, X: float, Z: float, logR: Iterable[float], logT: Iterable[float], logKappa: Iterable[Iterable[float]]):
+    def add_table(self, indicies: Tuple[float], logR: Iterable[float], logT: Iterable[float], logKappa: Iterable[Iterable[float]]):
         """
         @brief Add a table to the OPAT file.
-        @param X The X composition value.
+        @param indicies The index values of the table.
-        @param Z The Z composition value.
         @param logR The logR values.
         @param logT The logT values.
         @param logKappa The logKappa values.
         @throws ValueError if logKappa is not a non-empty 2D array or if logR and logT are not 1D arrays.
         """
+        if len(indicies) != self.header.numIndex:
+            raise ValueError(f"indicies must have length {self.header.numIndex}! Currently it has length {len(indicies)}")
         self.validate_logKappa(logKappa)
         self.validate_1D(logR, "logR")
         self.validate_1D(logT, "logT")
@@ -236,7 +249,7 @@ class OpatIO:
             logKappa = logKappa
         )
     
-        self.tables.append(((X, Z), table))
+        self.tables.append((indicies, table))
         self.header.numTables += 1
 
 
@@ -246,7 +259,7 @@ class OpatIO:
         @return The header as bytes.
         """
         headerBytes = struct.pack(
-            "<4s H I I Q 16s 64s 128s 26s",
+            "<4s H I I Q 16s 64s 128s H 24s",
             self.header.magic.encode('utf-8'),
             self.header.version,
             self.header.numTables,
@@ -255,6 +268,7 @@ class OpatIO:
             self.header.creationDate.encode('utf-8'),
             self.header.sourceInfo.encode('utf-8'),
             self.header.comment.encode('utf-8'),
+            self.header.numIndex,
             self.header.reserved
         )
         return headerBytes
@@ -276,7 +290,7 @@ class OpatIO:
             *logT,
             *logKappa
         )
-        checksum = self.compute_checksum(tableBytes)
+        checksum = self.compute_checksum(logKappa)
         return (checksum, tableBytes)
 
     def _tableIndex_bytes(self, tableIndex: TableIndex) -> bytes:
@@ -286,16 +300,16 @@ class OpatIO:
         @return The table index as bytes.
         @throws RuntimeError if the table index entry does not have 64 bytes.
         """
+        tableIndexFMTString = "<"+"d"*self.header.numIndex+f"QQ"
         tableIndexBytes = struct.pack(
-            '<ddQQ',
+            tableIndexFMTString,
-            tableIndex.X,
+            *tableIndex.index,
-            tableIndex.Z,
             tableIndex.byteStart,
             tableIndex.byteEnd
         )
         tableIndexBytes += tableIndex.sha256
 
-        if len(tableIndexBytes) != 64:
+        if len(tableIndexBytes) != 16+self.header.numIndex*8+32:
             raise RuntimeError(f"Each table index entry must have 64 bytes. Due to an unknown error the table index entry for (X,Z)=({tableIndex.X},{tableIndex.Z}) header has {len(tableIndexBytes)} bytes")
         
         return tableIndexBytes
@@ -313,24 +327,28 @@ class OpatIO:
   creationDate: {self.header.creationDate}
   sourceInfo: {self.header.sourceInfo}
   comment: {self.header.comment}
+  numIndex: {self.header.numIndex}
   reserved: {self.header.reserved}
 )"""
         return reprString
 
-    def _format_table_as_string(self, table: OPATTable, X: float, Z: float) -> str:
+    def _format_table_as_string(self, table: OPATTable, indices: List[float]) -> str:
         """
         @brief Format a table as a string.
         @param table The OPAT table.
-        @param X The X composition value.
+        @indices The index values of the table.
-        @param Z The Z composition value.
         @return The formatted table as a string.
         """
         tableString: List[str] = []
         # fixed width X and Z header per table
-        tableString.append(f"X: {X:<10.4f} Z: {Z:<10.4f}")
+        tableIndexString: List[str] = []
+        for index in indices:
+            tableIndexString.append(f"{index:<10.4f}")
+        tableString.append(" ".join(tableIndexString))
         tableString.append("-" * 80)
         # write logR across the top (reserving one col for where logT will be)
-        logRRow = f"{'':<10}"
+        tableString.append(f"{'':<10}{'logR':<10}")
+        logRRow = f"{'logT':<10}"
         logRRowTrue = "".join(f"{r:<10.4f}" for r in table.logR)
         tableString.append(logRRow + logRRowTrue)
         for i, logT in enumerate(table.logT):
@@ -354,10 +372,19 @@ class OpatIO:
 
         tableRows: List[str] = []
         tableRows.append("\nTable Indexes in OPAT File:\n")
-        tableRows.append(f"{'X':<10} {'Z':<10} {'Byte Start':<15} {'Byte End':<15} {'Checksum (SHA-256)'}")
+        headerString: str = ''
+        for indexID, index in enumerate(table_indexes[0].index):
+            indexKey = f"Index {indexID}"
+            headerString += f"{indexKey:<10}"
+        headerString += f"{'Byte Start':<15} {'Byte End':<15} {'Checksum (SHA-256)'}"
+        tableRows.append(headerString)
         tableRows.append("=" * 80)
         for entry in table_indexes:
-            tableRows.append(f"{entry.X:<10.4f} {entry.Z:<10.4f} {entry.byteStart:<15} {entry.byteEnd:<15} {entry.sha256[:16]}...")
+            tableEntry = ''
+            for index in entry.index:
+                tableEntry += f"{index:<10.4f}"
+            tableEntry += f"{entry.byteStart:<15} {entry.byteEnd:<15} {entry.sha256[:16]}..."
+            tableRows.append(tableEntry)
         return '\n'.join(tableRows)
 
     def save_as_ascii(self, filename: str) -> str:
@@ -370,12 +397,11 @@ class OpatIO:
         currentStartByte: int = 256
         tableIndexs: List[bytes] = []
         tableStrings: List[bytes] = []
-        for (X, Z), table in self.tables:
+        for index, table in self.tables:
             checksum, tableBytes = self._table_bytes(table)
-            tableStrings.append(self._format_table_as_string(table, X, Z) + "\n")
+            tableStrings.append(self._format_table_as_string(table, index) + "\n")
             tableIndex = TableIndex(
-                X = X,
+                index = index,
-                Z = Z,
                 byteStart = currentStartByte,
                 byteEnd = currentStartByte + len(tableBytes),
                 sha256 = checksum
@@ -386,20 +412,27 @@ class OpatIO:
             currentStartByte += len(tableBytes)
         self.header.indexOffset = currentStartByte
         with open(filename, 'w') as f:
-            f.write(f"{self.header.magic}\n")
+            f.write("This is an ASCII representation of an OPAT file, it is not a valid OPAT file in and of itself.\n")
-            f.write(f"Version: {self.header.version}\n")
+            f.write("This file is meant to be human readable and is not meant to be read by a computer.\n")
-            f.write(f"numTables: {self.header.numTables}\n")
+            f.write("The purpose of this file is to provide a human readable representation of the OPAT file which can be used for debugging purposes.\n")
-            f.write(f"headerSize (bytes): {self.header.headerSize}\n")
+            f.write("The full binary specification of the OPAT file can be found in the OPAT file format documentation at:\n")
-            f.write(f"tableIndex Offset (bytes): {self.header.indexOffset}\n")
+            f.write(" https://github.com/4D-STAR/4DSSE/blob/main/specs/OPAT/OPAT.pdf\n")
-            f.write(f"Creation Date: {self.header.creationDate}\n")
+            f.write("="*35 + " HEADER " + "="*36 + "\n")
-            f.write(f"Source Info: {self.header.sourceInfo}\n")
+            f.write(f">> {self.header.magic}\n")
-            f.write(f"Comment: {self.header.comment}\n")
+            f.write(f">> Version: {self.header.version}\n")
-            f.write("="*80 + "\n")
+            f.write(f">> numTables: {self.header.numTables}\n")
+            f.write(f">> headerSize (bytes): {self.header.headerSize}\n")
+            f.write(f">> tableIndex Offset (bytes): {self.header.indexOffset}\n")
+            f.write(f">> Creation Date: {self.header.creationDate}\n")
+            f.write(f">> Source Info: {self.header.sourceInfo}\n")
+            f.write(f">> Comment: {self.header.comment}\n")
+            f.write(f">> numIndex: {self.header.numIndex}\n")
+            f.write("="*37 + " DATA " + "="*37 + "\n")
             f.write("="*80 + "\n")
             for tableString in tableStrings:
                 f.write(tableString)
                 f.write("="*80 + "\n")
-            f.write("="*80 + "\n")
+            f.write("="*36 + " INDEX " + "="*37 + "\n")
             f.write(self.print_table_indexes(tableIndexs))
 
     def save(self, filename: str) -> str:
@@ -417,11 +450,10 @@ class OpatIO:
         currentStartByte: int = 256
         tableIndicesBytes: List[bytes] = []
         tablesBytes: List[bytes] = []
-        for (X, Z), table in self.tables:
+        for index, table in self.tables:
             checksum, tableBytes = self._table_bytes(table)
             tableIndex = TableIndex(
-                X = X,
+                index,
-                Z = Z,
                 byteStart = currentStartByte,
                 byteEnd = currentStartByte + len(tableBytes),
                 sha256 = checksum
@@ -455,7 +487,7 @@ def loadOpat(filename: str) -> OpatIO:
     opat = OpatIO()
     with open(filename, 'rb') as f:
         headerBytes: bytes = f.read(256)
-        unpackedHeader = struct.unpack("<4s H I I Q 16s 64s 128s 26s", headerBytes)
+        unpackedHeader = struct.unpack("<4s H I I Q 16s 64s 128s H 24s", headerBytes)
         loadedHeader = Header(
             magic = unpackedHeader[0].decode().replace("\x00", ""),
             version = unpackedHeader[1],
@@ -465,19 +497,22 @@ def loadOpat(filename: str) -> OpatIO:
             creationDate = unpackedHeader[5].decode().replace("\x00", ""),
             sourceInfo = unpackedHeader[6].decode().replace("\x00", ""),
             comment = unpackedHeader[7].decode().replace("\x00", ""),
-            reserved = unpackedHeader[8]
+            numIndex = unpackedHeader[8],
+            reserved = unpackedHeader[9]
         )
         opat.header = loadedHeader
         f.seek(opat.header.indexOffset)
         tableIndices: List[TableIndex] = []
-        while tableIndexEntryBytes := f.read(32):
+        tableIndexChunkSize = 16 + loadedHeader.numIndex*8
-            unpackedTableIndexEntry = struct.unpack("<ddQQ", tableIndexEntryBytes)
+        tableIndexFMTString = "<"+"d"*loadedHeader.numIndex+"QQ"
+        while tableIndexEntryBytes := f.read(tableIndexChunkSize):
+            unpackedTableIndexEntry = struct.unpack(tableIndexFMTString, tableIndexEntryBytes)
             checksum = f.read(32)
+            index = unpackedTableIndexEntry[:loadedHeader.numIndex]
             tableIndexEntry = TableIndex(
-                X = unpackedTableIndexEntry[0],
+                index = index,
-                Z = unpackedTableIndexEntry[1],
+                byteStart = unpackedTableIndexEntry[loadedHeader.numIndex],
-                byteStart = unpackedTableIndexEntry[2],
+                byteEnd = unpackedTableIndexEntry[loadedHeader.numIndex+1],
-                byteEnd = unpackedTableIndexEntry[3],
                 sha256 = checksum
             )
             tableIndices.append(tableIndexEntry)
@@ -500,5 +535,5 @@ def loadOpat(filename: str) -> OpatIO:
             logT = np.array(unpackedData[N_R: N_R+N_T], dtype=np.float64)
             logKappa = np.array(unpackedData[N_R+N_T:], dtype=np.float64).reshape((N_R, N_T))
 
-            opat.add_table(tableIndex.X, tableIndex.Z, logR, logT, logKappa)
+            opat.add_table(tableIndex.index, logR, logT, logKappa)
     return opat

utils/opatio/utils/mkTestData.py

@@ -0,0 +1,27 @@
+from opatio import OpatIO
+import numpy as np
+np.random.seed(42)
+
+def generate_synthetic_opacity_table(n_r, n_t):
+    logR = np.linspace(-8, 2, n_r, dtype=np.float64)  # log Density grid
+    logT = np.linspace(3, 9, n_t, dtype=np.float64)  # log Temperature grid
+    logK = np.random.uniform(-2, 2, size=(n_r, n_t)).astype(np.float64)  # Synthetic Opacity
+    return logR, logT, logK
+
+if __name__ == "__main__":
+    n_r = 50
+    n_t = 50
+    num_tables = 20
+    XValues = np.linspace(0.1, 0.7, num_tables)
+    ZValues = np.linspace(0.001, 0.03, num_tables)
+    opat = OpatIO()
+    opat.set_comment("Synthetic Opacity Tables")
+    opat.set_source("utils/opatio/utils/mkTestData.py")
+
+    for i in range(num_tables):
+        logR, logT, logK = generate_synthetic_opacity_table(n_r, n_t)
+        opat.add_table((XValues[i], ZValues[i]), logR, logT, logK)
+
+    opat.save("testData/synthetic_tables.opat")
+    opat.save_as_ascii("testData/synthetic_tables_OPAT.ascii")
+