Merge pull request #8 from tboudreaux/spec/OPAT

OPAT file format

meson.build

@@ -1,4 +1,4 @@
-project('4DSSE', 'cpp', version: '0.0.1a', default_options: ['cpp_std=c++23'])
+project('4DSSE', 'cpp', version: '0.0.1a', default_options: ['cpp_std=c++23'], meson_version: '>=1.6.0')
 
 # Add default visibility for all C++ targets
 add_project_arguments('-fvisibility=default', language: 'cpp')

meson_options.txt

@@ -1 +1 @@
-option('build_tests', type: 'boolean', value: true, description: 'Build tests')
+option('build_tests', type: 'boolean', value: true, description: 'Build tests')

src/meson.build

@@ -3,4 +3,5 @@ subdir('resources')
 
 # Build the main source code
 subdir('dobj')
-subdir('const')
+subdir('const')
+subdir('opatIO')

src/opatIO/meson.build

@@ -0,0 +1,18 @@
+# Define the library
+opatIO_sources = files(
+  'private/opatIO.cpp',
+)
+
+opatIO_headers = files(
+    'public/opatIO.h'
+)
+
+# Define the libopatIO library so it can be linked against by other parts of the build system
+libopatIO = library('opatIO',
+    opatIO_sources, 
+    include_directories: include_directories('public'),
+    cpp_args: ['-fvisibility=default'],
+    install : true)
+
+# Make headers accessible
+install_headers(opatIO_headers, subdir : '4DSSE/opatIO')

src/opatIO/private/opatIO.cpp

@@ -0,0 +1,433 @@
+#include "opatIO.h"
+#include <fstream>
+#include <iostream>
+#include <stdexcept>
+#include <cstring>
+#include <algorithm>
+#include <iomanip>
+#include <map>
+#include <utility>
+#include <cmath>
+#include <limits>
+#include <deque>
+
+// Constructor
+OpatIO::OpatIO() {}
+
+OpatIO::OpatIO(std::string filename) : filename(filename) {
+    load();
+}
+
+// Destructor
+OpatIO::~OpatIO() {
+    unload();
+}
+
+// Load the OPAT file
+void OpatIO::load() {
+    if (loaded) return;
+
+    std::ifstream file(filename, std::ios::binary);
+    if (!file.is_open()) {
+        throw std::runtime_error("Could not open file: " + filename);
+    }
+
+    readHeader(file);
+    readTableIndex(file);
+
+    loaded = true;
+    file.close();
+}
+
+// // Unload the OPAT file
+void OpatIO::unload() {
+    if (!loaded) return;
+
+    tableIndex.clear();
+    while (!tableQueue.empty()) {
+        tableQueue.pop_front();
+    }
+
+    loaded = false;
+}
+
+// Read the header from the file
+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);
+    }
+}
+
+// Read the table index from the 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);
+    }
+    buildTableIDToComposition();
+}
+
+void OpatIO::buildTableIDToComposition(){
+    tableIDToComposition.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);
+        tableID++;
+    }
+    XZLookupEpsilon();
+}
+
+void OpatIO::XZLookupEpsilon() {
+    /* 
+     Get 10% of the minimum spacing between XZ values 
+     in the tableIDToComposition map. This can be used
+     to set the comparison distance when doing a reverse
+     lookup (composition -> tableID)
+     */
+    std::vector<double> Xvalues, Zvalues;
+    double epsilonX, epsilonZ, xgap, zgap;
+
+    // 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;
+        }
+    }
+    // 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 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){
+            return tableID;
+        }
+        tableID++;
+    }
+    return -1;
+}
+
+// Get a table from the queue
+OPATTable OpatIO::getTableFromQueue(int tableID) {
+    for (const auto &table : tableQueue) {
+        if (table.first == tableID) {
+            return table.second;
+        }
+    }
+    throw std::out_of_range("Table not found!");
+}
+
+// Add a table to the queue
+void OpatIO::addTableToQueue(int tableID, OPATTable table) {
+    if (static_cast<int>(tableQueue.size()) >= maxQDepth) {
+        removeTableFromQueue();
+    }
+    std::pair<int, OPATTable> IDTablePair = {tableID, table};
+    tableQueue.push_back(IDTablePair);
+}
+
+// Remove a table from the queue
+void OpatIO::removeTableFromQueue() {
+    if (!tableQueue.empty()) {
+        tableQueue.pop_front();
+    }
+}
+
+// Flush the queue
+void OpatIO::flushQueue() {
+    while (!tableQueue.empty()) {
+        tableQueue.pop_back();
+        tableQueue.pop_front();
+    }
+}
+
+// Get the OPAT version
+uint16_t OpatIO::getOPATVersion() {
+    return header.version;
+}
+
+// Get a table for given X and Z
+OPATTable OpatIO::getTable(double X, double Z) {
+    int tableID = lookupTableID(X, Z);
+    if (tableID == -1) {
+        throw std::out_of_range("X Z Pair Not found!");
+    }
+    try {
+        return getTableFromQueue(tableID);
+    }
+    catch(const std::out_of_range &e) { 
+        return getTable(tableID);
+    }
+}
+
+OPATTable OpatIO::getTable(int tableID) {
+    std::ifstream file(filename, std::ios::binary);
+    if (!file.is_open()) {
+        throw std::runtime_error("Could not open file: " + filename);
+    }
+
+    uint64_t byteStart = tableIndex[tableID].byteStart;
+
+    file.seekg(byteStart, std::ios::beg);
+
+    OPATTable table;
+
+    // Step 1: Read N_R and N_T
+    file.read(reinterpret_cast<char*>(&table.N_R), sizeof(uint32_t));
+    file.read(reinterpret_cast<char*>(&table.N_T), sizeof(uint32_t));
+    
+
+    // Resize vectors to hold the correct number of elements
+    table.logR.resize(table.N_R);
+    table.logT.resize(table.N_T);
+    table.logKappa.resize(table.N_R, std::vector<double>(table.N_T));
+
+    // Step 2: Read logR values
+    file.read(reinterpret_cast<char*>(table.logR.data()), table.N_R * sizeof(double));
+
+    // Step 3: Read logT values
+    file.read(reinterpret_cast<char*>(table.logT.data()), table.N_T * sizeof(double));
+
+    // Step 4: Read logKappa values (flattened row-major order)
+    for (size_t i = 0; i < table.N_R; ++i) {
+
+        file.read(reinterpret_cast<char*>(table.logKappa[i].data()), table.N_T * sizeof(double));
+    }
+
+    if (!file) {
+        throw std::runtime_error("Error reading table from file: " + filename);
+    }
+
+    addTableToQueue(tableID, table);
+    file.close();
+    return table;
+}
+
+// Set the maximum queue depth
+void OpatIO::setMaxQDepth(int depth) {
+    maxQDepth = depth;
+}
+
+int OpatIO::getMaxQDepth() {
+    return maxQDepth;
+}
+
+// Set the filename
+void OpatIO::setFilename(std::string filename) {
+    if (loaded) {
+        throw std::runtime_error("Cannot set filename while file is loaded");
+    }
+    this->filename = filename;
+}
+
+// Check if the file is loaded
+bool OpatIO::isLoaded() {
+    return loaded;
+}
+
+// Print the header
+void OpatIO::printHeader() {
+    std::cout << "Version: " << header.version << std::endl;
+    std::cout << "Number of Tables: " << header.numTables << std::endl;
+    std::cout << "Header Size: " << header.headerSize << std::endl;
+    std::cout << "Index Offset: " << header.indexOffset << std::endl;
+    std::cout << "Creation Date: " << header.creationDate << std::endl;
+    std::cout << "Source Info: " << header.sourceInfo << std::endl;
+    std::cout << "Comment: " << header.comment << std::endl;
+}
+
+// Print the table index
+void OpatIO::printTableIndex() {
+    if (tableIndex.empty()) {
+        std::cout << "No table indexes found." << std::endl;
+        return;
+    }
+
+    // Print table header
+    std::cout << std::left << std::setw(10) << "X"
+        << std::setw(10) << "Z"
+        << std::setw(15) << "Byte Start"
+        << std::setw(15) << "Byte End"
+        << "Checksum (SHA-256)" << std::endl;
+
+    std::cout << std::string(80, '=') << std::endl;  // Separator line
+
+    // 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::setw(15) << index.byteEnd
+            << std::hex; // Switch to hex mode for checksum
+
+        for (int i = 0; i < 8; ++i) {  // Print first 8 bytes of SHA-256 for brevity
+            std::cout << std::setw(2) << std::setfill('0') << (int)index.sha256[i];
+        }
+
+        std::cout << "..." << std::dec << std::setfill(' ') << std::endl;  // Reset formatting
+    }
+}
+
+void OpatIO::printTable(OPATTable table, uint32_t truncateDigits) {
+    int printTo;
+    bool truncate = false;
+    if (table.N_R > truncateDigits) {
+        printTo = truncateDigits;
+        truncate = true;
+    } else {
+        printTo = table.N_R-1;
+    }
+    std::cout << "LogR (size: " << table.logR.size() << "): [";
+    for (int i = 0; i < printTo; ++i) {
+        std::cout << table.logR.at(i) << ", ";
+    }
+    if (truncate) {
+        std::cout << "..., ";
+        for (int i = truncateDigits; i > 1; --i) {
+            std::cout << table.logR.at(table.logR.size() - i) << ", ";
+        }
+    }
+    std::cout << table.logR.back() << "]" << std::endl;
+
+    if (table.N_T > truncateDigits) {
+        printTo = truncateDigits;
+        truncate = true;
+    } else {
+        printTo = table.N_T-1;
+    }
+    std::cout << "LogT (size: " << table.logT.size() << "): [";
+    for (int i = 0; i < printTo; ++i) {
+        std::cout << table.logT.at(i) << ", ";
+    }
+    if (truncate) {
+        std::cout << "..., ";
+        for (int i = truncateDigits; i > 1; --i) {
+            std::cout << table.logT.at(table.logT.size() - i) << ", ";
+        }
+    }
+    std::cout << table.logT.back() << "]" << std::endl;
+
+    bool truncateRow = false;
+    bool truncateCol = false;
+    int printToRow, printToCol;
+    if (table.N_T > truncateDigits) {
+        printToRow = truncateDigits;
+        truncateRow = true;
+    } else {
+        printToRow = table.N_T-1;
+    }
+    if (table.N_R > truncateDigits) {
+        printToCol = truncateDigits;
+        truncateCol = true;
+    } else {
+        printToCol = table.N_R-1;
+    }
+
+    std::cout << "LogKappa (size: " << table.N_R << " x " << table.N_T << "): \n[";
+    for (int rowIndex = 0; rowIndex < printToRow; rowIndex++) {
+        std::cout << "[";
+        for (int colIndex = 0; colIndex < printToCol; colIndex++) {
+            std::cout << table.logKappa.at(rowIndex).at(colIndex) << ", ";
+        }
+        if (truncateRow) {
+            std::cout << "..., ";
+            for (int i = truncateDigits; i > 1; i--) {
+                std::cout << table.logKappa.at(rowIndex).at(table.logKappa.at(rowIndex).size() - i) << ", ";
+            }
+        }
+        std::cout << table.logKappa.at(rowIndex).back() << "],\n";
+    } 
+    if (truncateCol) {
+        std::cout << ".\n.\n.\n";
+        for (int rowIndex = truncateDigits; rowIndex > 1; rowIndex--) {
+            std::cout << "[";
+            for (int colIndex = 0; colIndex < printToCol; colIndex++) {
+                std::cout << table.logKappa.at(rowIndex).at(colIndex) << ", ";
+            }
+            if (truncateRow) {
+                std::cout << "..., ";
+                for (int i = truncateDigits; i > 1; i--) {
+                    std::cout << table.logKappa.at(rowIndex).at(table.logKappa.at(rowIndex).size() - i) << ", ";
+                }
+            }
+            std::cout << table.logKappa.at(rowIndex).back() << "],\n";
+        } 
+        std::cout << "[";
+        for (int colIndex = 0; colIndex < printToCol; colIndex++) {
+            std::cout << table.logKappa.back().at(colIndex) << ", ";
+        }
+        if (truncateRow) {
+            std::cout << "..., ";
+            for (int i = truncateDigits; i > 1; i--) {
+                std::cout << table.logKappa.back().at(table.logKappa.back().size() - i) << ", ";
+            }
+        }
+        std::cout << table.logKappa.back().back() << "]";
+    }
+    std::cout << "]" << std::endl;
+}
+
+void OpatIO::printTable(double X, double Z, uint32_t truncateDigits) {
+    int tableID = lookupTableID(X, Z);
+    OPATTable table = getTable(tableID);
+    printTable(table, truncateDigits);
+}
+
+// Get the table index
+std::vector<TableIndex> OpatIO::getTableIndex() {
+    return tableIndex;
+}
+
+// Get the header
+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 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;
+// }
+
+// // 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;
+// }

src/opatIO/public/opatIO.h

@@ -0,0 +1,99 @@
+#ifndef OPATIO_H
+#define OPATIO_H
+
+#include <iostream>
+#include <fstream>
+#include <string>
+#include <vector>
+#include <deque>
+#include <map>
+#include <utility>
+
+struct Header {
+    char magic[4];
+    uint16_t version;
+    uint32_t numTables;
+    uint32_t headerSize;
+    uint64_t indexOffset;
+    char creationDate[16];
+    char sourceInfo[64];
+    char comment[128];
+    char reserved[26]; 
+};
+
+struct TableIndex {
+    double X;
+    double Z;
+    // double C; // For type 2 OPAL tables. When not set will be 0
+    // double O; // For type 2 OPAL tables. When not set will be 0
+    uint64_t byteStart;
+    uint64_t byteEnd;
+    char sha256[32];
+};
+
+struct OPATTable {
+    uint32_t N_R;
+    uint32_t N_T;
+    std::vector<double> logR;
+    std::vector<double> logT;
+    std::vector<std::vector<double>> logKappa;
+};
+
+class opatIOTest; // Friend for gtest
+
+class OpatIO {
+private:
+    Header header;
+    std::vector<TableIndex> tableIndex;
+    std::deque<std::pair<int, OPATTable>> tableQueue;
+    std::map<int, std::pair<double, double>> tableIDToComposition;
+    std::pair<double, double> XZepsilon;
+    int maxQDepth = 10;
+    std::string filename;
+    bool loaded = false;
+
+    void readHeader(std::ifstream &file);
+    void readTableIndex(std::ifstream &file);
+
+    OPATTable getTableFromQueue(int tableID);
+    void addTableToQueue(int tableID, OPATTable table);
+    void removeTableFromQueue();
+    void flushQueue();
+
+    OPATTable getTable(int tableID);
+    void printTable(OPATTable table, uint32_t truncateDigits=5);
+
+    void XZLookupEpsilon();
+    void buildTableIDToComposition();
+
+public:
+    OpatIO();
+    OpatIO(std::string filename);
+    ~OpatIO();
+
+    OPATTable getTable(double X, double Z);
+    void setMaxQDepth(int depth);
+    int getMaxQDepth();
+    void setFilename(std::string filename);
+    void load();
+    void unload();
+    bool isLoaded();
+
+    void printHeader();
+    void printTableIndex();
+    void printTable(double X, double Z, uint32_t truncateDigits=5);
+
+    std::vector<TableIndex> getTableIndex();
+    Header getHeader();
+
+    std::vector<OPATTable> getClosestXTables(double X, double ZExact, double C=0, double O=0, int numTables=1);
+    std::vector<OPATTable> getClosestZTables(double XExact, double Z, double C=0, double O=0, int numTables=1);
+    std::vector<OPATTable> getClosestTables(double X, double Z, double C=0, double O=0, int numTables=1);  
+
+    int lookupTableID(double X, double Z);
+    int lookupClosestTableID(double X, double Z);
+    uint16_t getOPATVersion();
+};
+
+
+#endif

tests/meson.build

@@ -5,6 +5,7 @@ gtest_nomain_dep = dependency('gtest', main: false, required : true)
 # Subdirectories for unit and integration tests
 subdir('dobj')
 subdir('const')
+subdir('opatIO')
 
 # Subdirectories for sandbox tests
 subdir('dobj_sandbox')

tests/opatIO/meson.build

@@ -0,0 +1,27 @@
+# Test files for opatIO
+test_sources = [
+    'opatIOTest.cpp',
+]
+
+
+
+foreach test_file : test_sources
+  exe_name = test_file.split('.')[0]
+  message('Building test: ' + exe_name)
+
+  # Create an executable target for each test
+  test_exe = executable(
+      exe_name,
+      test_file,
+      dependencies: [gtest_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
+  )
+
+  # Add the executable as a test
+  test(
+    exe_name,
+    test_exe,
+    env: ['MESON_SOURCE_ROOT=' + meson.project_source_root(), 'MESON_BUILD_ROOT=' + meson.project_build_root()])
+endforeach

tests/opatIO/opatIOTest.cpp

@@ -0,0 +1,86 @@
+#include <gtest/gtest.h>
+#include "opatIO.h"
+#include <iostream>
+#include <string>
+#include <vector>
+#include <set>
+#include <sstream>
+
+std::string EXAMPLE_FILENAME = std::string(getenv("MESON_SOURCE_ROOT")) + "/tests/opatIO/example.opat";
+
+/**
+ * @file opatIOTest.cpp
+ * @brief Unit tests for the OpatIO class and associated structs.
+ */
+
+/**
+ * @brief Test suite for the const class.
+ */
+class opatIOTest : public ::testing::Test {};
+
+/**
+ * @test Verify default constructor initializes correctly.
+ */
+TEST_F(opatIOTest, DefaultConstructor) {
+    EXPECT_NO_THROW(OpatIO());
+}
+
+/* 
+ * @test Verify constructor initializes correctly with a file.
+ */
+TEST_F(opatIOTest, Constructor) {
+    OpatIO opatIO(EXAMPLE_FILENAME);
+}
+
+TEST_F(opatIOTest, Header) {
+    OpatIO opatIO(EXAMPLE_FILENAME);
+    Header header = opatIO.getHeader();
+    EXPECT_EQ(header.version, 1);
+    EXPECT_EQ(header.numTables, 20);
+    EXPECT_EQ(header.headerSize, 256);
+    EXPECT_EQ(header.indexOffset, 416416);
+    EXPECT_EQ(std::string(header.creationDate), "Feb 14, 2025");
+    EXPECT_EQ(std::string(header.sourceInfo), "MESA 12700, Synthetic Opacity Data");
+    EXPECT_EQ(std::string(header.comment), "log10 kappa (cm^2/g)");
+}
+
+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].Z, 0.001);
+    EXPECT_EQ(tableIndex[0].byteStart, 256);
+    EXPECT_EQ(tableIndex[0].byteEnd, 21064);
+}
+
+TEST_F(opatIOTest, MaxQDepth) {
+    OpatIO opatIO(EXAMPLE_FILENAME);
+    EXPECT_EQ(opatIO.getMaxQDepth(), 10);
+    opatIO.setMaxQDepth(5);
+    EXPECT_EQ(opatIO.getMaxQDepth(), 5);
+}
+
+TEST_F(opatIOTest, Unload){
+    OpatIO opatIO(EXAMPLE_FILENAME);
+    EXPECT_NO_THROW(opatIO.unload());
+    EXPECT_FALSE(opatIO.isLoaded());
+}
+
+TEST_F(opatIOTest, LookupTableID) {
+    OpatIO opatIO(EXAMPLE_FILENAME);
+    EXPECT_EQ(opatIO.lookupTableID(0.321053, 0.0116842), 7);
+}
+
+TEST_F(opatIOTest, GetTable) {
+    OpatIO opatIO(EXAMPLE_FILENAME);
+    OPATTable tab = opatIO.getTable(0.1, 0.001);
+    EXPECT_EQ(tab.N_R, 50);
+    EXPECT_EQ(tab.N_T, 50);
+    EXPECT_DOUBLE_EQ(tab.logR[0], -8.0);
+    EXPECT_DOUBLE_EQ(tab.logT[0], 3.0);
+    EXPECT_DOUBLE_EQ(tab.logKappa[0][0], -0.50183952461055);
+    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);
+}

utils/opatio/pyproject.toml

@@ -0,0 +1,16 @@
+[build-system]
+requires = ["setuptools", "wheel"]
+build-backend = "setuptools.build_meta"
+
+[project]
+name = "opatio"
+version = "0.1.0a"
+description = "A python module for handling OPAT files"
+readme = "readme.md"
+authors = [{name = "Emily M. Boudreaux", email = "emily.boudreaux@dartmouth.edu"}]
+requires-python = ">=3.8"
+dependencies = ["numpy >= 1.21.1"]
+
+[tool.setuptools]
+packages = ["opatio", "opatio.opat"]
+package-dir = {"" = "src"}

utils/opatio/readme.md

@@ -0,0 +1,46 @@
+# opatIO python module
+This module defines a set of tools to build, write, and read OPAT files. 
+The OPAT fileformat is a custom file format designed to efficiently store
+opacity information for a variety of compositions. 
+
+## Installation
+You can install this module with pip
+```bash
+git clone <repo>
+cd 4DSSE/utils/opat
+pip install .
+```
+
+## General Usage
+The general way that this module is mean to be used is to first build a schema for the opaticy table and then save that to disk. The module will handle all the byte aligment and lookup table construction for you. 
+
+A simple example might look like the following
+
+```python
+from opatio import OpatIO
+
+opacityFile = OpatIO()
+opacityFile.set_comment("This is a sample opacity file")
+opaticyFile.set_source("OPLIB")
+
+# some code to get a logR, logT, and logKappa table
+# where logKappa is of size (n,m) if logR is size n and
+# logT is size m
+
+opacityFile.add_table(X, Z, logR, logT, logKappa)
+opacityFile.save("opacity.opat")
+```
+
+You can also read opat files which have been generated with the loadOpat function
+
+```python
+from opatio import loadOpat
+
+opacityFile = loadOpat("opacity.opat")
+
+print(opacityFile.header)
+print(opaticyFile.tables[0])
+```
+
+## Problems
+If you have problems feel free to either submit an issue to the root github repo (tagged as utils/opatio) or email Emily Boudreaux at emily.boudreaux@dartmouth.edu

utils/opatio/src/opatio/__init__.py

@@ -0,0 +1 @@
+from .opat.opat import OpatIO, loadOpat

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

@@ -0,0 +1,277 @@
+import struct
+import numpy as np
+from datetime import datetime
+
+from dataclasses import dataclass
+
+from typing import Iterable, List, Tuple
+from collections.abc import Iterable as collectionIterable
+
+import hashlib
+
+import os
+
+@dataclass
+class Header:
+    magic: str
+    version: int
+    numTables: int
+    headerSize: int
+    indexOffset: int
+    creationDate: str
+    sourceInfo: str
+    comment: str
+    reserved: bytes
+
+@dataclass
+class TableIndex:
+    X: float
+    Z: float
+    byteStart: int
+    byteEnd: int
+    sha256: bytes
+
+@dataclass
+class OPATTable:
+    N_R: int
+    N_T: int
+    logR: Iterable[float]
+    logT: Iterable[float]
+    logKappa: Iterable[Iterable[float]]
+
+defaultHeader = Header(
+    magic="OPAT",
+    version=1,
+    numTables=0,
+    headerSize=256,
+    indexOffset=0,
+    creationDate=datetime.now().strftime("%b %d, %Y"),
+    sourceInfo="no source provided by user",
+    comment="default header",
+    reserved=b"\x00" * 26
+)
+
+class OpatIO:
+    def __init__(self):
+        self.header: Header = defaultHeader
+        self.tables: List[Tuple[Tuple[float, float], OPATTable]] = []
+
+    @staticmethod
+    def validate_char_array_size(s: str, nmax: int) -> bool:
+        if len(s) > nmax:
+            return False
+        return True
+
+    @staticmethod
+    def validate_logKappa(logKappa):
+        if isinstance(logKappa, np.ndarray):
+            if logKappa.ndim == 2:
+                return
+            else:
+                raise ValueError("logKappa must be a non-empty 2D array")
+
+        if isintance(logKappa, collectionIterable) and all(isinstance(row, collectionIterable) for row in logKappa):
+            try:
+                first_row = next(iter(logKappa))
+                if all(isinstance(x, (int, float)) for x in first_row):
+                    return
+                else:
+                    raise ValueError("logKappa must be fully numeric")
+            except StopIteration:
+                raise ValueError("logKappa must be a non-empty 2D iterable")
+        else:
+            raise TypeError("logKappa must be a non-empty 2D array or iterable")
+
+    @staticmethod
+    def validate_1D(arr, name: str):
+        if isinstance(arr, np.ndarray):
+            if arr.ndim == 1:
+                return
+            else:
+                raise ValueError(f"{name} must be a 1D numpy array")
+        if isinstance(arr, collectionIterable) and not isinstance(arr, (str, bytes)):
+            if all(isinstance(x, (int, float)) for x in arr):
+                return
+            else:
+                raise ValueError(f"{name} must be fully numeric")
+        else:
+            raise TypeError(f"{name} must be a non-empty 2D array or iterable")
+
+    @staticmethod
+    def compute_checksum(data: bytes) -> bytes:
+        return hashlib.sha256(data).digest()
+
+    def set_version(self, version: int) -> int:
+        self.header.version = version
+        return self.header.version
+
+    def set_source(self, source: str) -> str:
+        if not self.validate_char_array_size(source, 64):
+            raise TypeError(f"sourceInfo string ({source}) is too long ({len(source)}). Max length is 64")
+        self.header.sourceInfo = source
+        return self.header.sourceInfo
+
+    def set_comment(self, comment: str) -> str:
+        if not self.validate_char_array_size(comment, 128):
+            raise TypeError(f"comment string ({comment}) is too long ({len(comment)}). Max length is 128")
+        self.header.comment = comment
+        return self.header.comment
+
+    def add_table(self, X: float, Z: float, logR: Iterable[float], logT: Iterable[float], logKappa: Iterable[Iterable[float]]):
+        self.validate_logKappa(logKappa)
+        self.validate_1D(logR, "logR")
+        self.validate_1D(logT, "logT")
+
+        logR = np.array(logR)
+        logT = np.array(logT)
+        logKappa = np.array(logKappa)
+
+        if logKappa.shape != (logR.shape[0], logT.shape[0]):
+            raise ValueError(f"logKappa must be of shape ({len(logR)} x {len(logT)})! Currently logKappa has shape {logKappa.shape}")
+
+        table = OPATTable(
+            N_R = logR.shape[0],
+            N_T = logT.shape[0],
+            logR = logR,
+            logT = logT,
+            logKappa = logKappa
+        )
+    
+        self.tables.append(((X, Z), table))
+        self.header.numTables += 1
+
+
+    def _header_bytes(self) -> bytes:
+        headerBytes = struct.pack(
+            "<4s H I I Q 16s 64s 128s 26s",
+            self.header.magic.encode('utf-8'),
+            self.header.version,
+            self.header.numTables,
+            self.header.headerSize,
+            self.header.indexOffset,
+            self.header.creationDate.encode('utf-8'),
+            self.header.sourceInfo.encode('utf-8'),
+            self.header.comment.encode('utf-8'),
+            self.header.reserved
+        )
+        return headerBytes
+
+    def _table_bytes(self, table: OPATTable) -> Tuple[bytes, bytes]:
+        logR = table.logR.flatten()
+        logT = table.logT.flatten()
+        logKappa = table.logKappa.flatten()
+        tableBytes = struct.pack(
+            f"<II{table.N_R}d{table.N_T}d{table.N_R*table.N_T}d",
+            table.N_R,
+            table.N_T,
+            *logR,
+            *logT,
+            *logKappa
+        )
+        checksum = self.compute_checksum(tableBytes)
+        return (checksum, tableBytes)
+
+    def _tableIndex_bytes(self, tableIndex: TableIndex) -> bytes:
+        tableIndexBytes = struct.pack(
+            '<ddQQ',
+            tableIndex.X,
+            tableIndex.Z,
+            tableIndex.byteStart,
+            tableIndex.byteEnd
+        )
+        tableIndexBytes += tableIndex.sha256
+
+        if len(tableIndexBytes) != 64:
+            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
+
+    def save(self, filename: str) -> str:
+        tempHeaderBytes = self._header_bytes()
+
+        if len(tempHeaderBytes) != 256:
+            raise RuntimeError(f"Header must have 256 bytes. Due to an unknown error the header has {len(tempHeaderBytes)} bytes")
+
+        currentStartByte: int = 256
+        tableIndicesBytes: List[bytes] = []
+        tablesBytes: List[bytes] = []
+        for (X, Z), table in self.tables:
+            checksum, tableBytes = self._table_bytes(table)
+            tableIndex = TableIndex(
+                X = X,
+                Z = Z,
+                byteStart = currentStartByte,
+                byteEnd = currentStartByte + len(tableBytes),
+                sha256 = checksum
+            )
+            tableIndexBytes = self._tableIndex_bytes(tableIndex)
+            tablesBytes.append(tableBytes)
+            tableIndicesBytes.append(tableIndexBytes)
+
+            currentStartByte += len(tableBytes)
+        self.header.indexOffset = currentStartByte
+        headerBytes = self._header_bytes()
+
+        with open(filename, 'wb') as f:
+            f.write(headerBytes)
+            for tableBytes in tablesBytes:
+                f.write(tableBytes)
+            for tableIndexBytes in tableIndicesBytes:
+                f.write(tableIndexBytes)
+
+        if os.path.exists(filename):
+            return filename
+
+
+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)
+        loadedHeader = Header(
+            magic = unpackedHeader[0].decode(),
+            version = unpackedHeader[1],
+            numTables = unpackedHeader[2],
+            headerSize = unpackedHeader[3],
+            indexOffset = unpackedHeader[4],
+            creationDate = unpackedHeader[5].decode(),
+            sourceInfo = unpackedHeader[6].decode(),
+            comment = unpackedHeader[7].decode(),
+            reserved = unpackedHeader[8]
+        )
+        opat.header = loadedHeader
+        f.seek(opat.header.indexOffset)
+        tableIndices: List[TableIndex] = []
+        while tableIndexEntryBytes := f.read(32):
+            unpackedTableIndexEntry = struct.unpack("<ddQQ", tableIndexEntryBytes)
+            checksum = f.read(32)
+            tableIndexEntry = TableIndex(
+                X = unpackedTableIndexEntry[0],
+                Z = unpackedTableIndexEntry[1],
+                byteStart = unpackedTableIndexEntry[2],
+                byteEnd = unpackedTableIndexEntry[3],
+                sha256 = checksum
+            )
+            tableIndices.append(tableIndexEntry)
+        
+        currentStartByte = opat.header.headerSize
+        f.seek(currentStartByte)
+        for tableIndex in tableIndices:
+            f.seek(tableIndex.byteStart)
+            byteLength = tableIndex.byteEnd - tableIndex.byteStart
+            tableBytes = f.read(byteLength)
+
+            nr_nt_fmt = "<II"
+            nr_nt_size = struct.calcsize(nr_nt_fmt)
+            N_R, N_T = struct.unpack(nr_nt_fmt, tableBytes[:nr_nt_size])
+
+            dataFormat = f"<{N_R}d{N_T}d{N_R*N_T}d"
+            unpackedData = struct.unpack(dataFormat, tableBytes[nr_nt_size:])
+
+            logR = np.array(unpackedData[:N_R], dtype=np.float64)
+            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)
+    return opat
+