feat(NetworkJacobian): Rank caching, regularization, and CSV export

Comptuting the rank of a large matrix with QR factorization can be
expensive and the rank is often never needed. We have implimented a
caching system so that the rank is only computed when asked for, and
then only once. Further, the regularization method which was previously
in an anonomous namespace inside of a single translation unit has been
moved to the jacobian header and implimentation file. This allows all
parts of GridFire to access the same regularization method. Finally a
small CSV output method has been added which is useful for debugging

src/include/gridfire/engine/types/jacobian.h

@@ -1,14 +1,18 @@
 #pragma once
 
 #include "fourdst/atomic/atomicSpecies.h"
+#include "fourdst/composition/composition_abstract.h"
+#include "quill/Logger.h"
 #include <Eigen/SparseCore>
 #include <Eigen/SparseQR>
 
 #include <tuple>
 #include <functional>
 #include <unordered_map>
+#include <optional>
 
 namespace gridfire {
+    constexpr double MIN_ABUNDANCE_TO_CONTRIBUTE_TO_JACOBIAN = 1e-100;
     using JacobianEntry = std::pair<std::pair<fourdst::atomic::Species, fourdst::atomic::Species>, double>;
 
     class NetworkJacobian {
@@ -60,9 +64,18 @@ namespace gridfire {
         [[nodiscard]] Eigen::SparseMatrix<double> data() const;
         [[nodiscard]] const std::unordered_map<fourdst::atomic::Species, size_t>& mapping() const;
 
+        void to_csv(const std::string& filename) const;
+
     private:
         Eigen::SparseMatrix<double> m_jacobianMatrix;
         std::unordered_map<fourdst::atomic::Species, size_t> m_speciesToIndexMap;
-        size_t m_rank;
+
+        mutable std::optional<size_t> m_rank = std::nullopt;
     };
+
+    NetworkJacobian regularize_jacobian(
+        const NetworkJacobian& jacobian,
+        const fourdst::composition::CompositionAbstract& comp,
+        std::optional<quill::Logger*> logger = std::nullopt
+    );
 }

src/lib/engine/types/jacobian.cpp

@@ -1,7 +1,15 @@
 #include "gridfire/engine/types/jacobian.h"
+
+#include <format>
+#include <fstream>
+#include <print>
+#include <chrono>
+#include <ranges>
 #include <Eigen/SparseCore>
 #include <Eigen/SparseQR>
 
+#include "quill/LogMacros.h"
+
 namespace gridfire {
     NetworkJacobian::NetworkJacobian(
         const Eigen::SparseMatrix<double>& jacobianMatrix,
@@ -12,13 +20,6 @@ namespace gridfire {
             m_speciesToIndexMap[species] = i;
         }
 
-        if (m_jacobianMatrix.rows() == 0 || m_jacobianMatrix.cols() == 0) {
-            m_rank = 0;
-        } else {
-            Eigen::SparseQR<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>> solver;
-            solver.compute(m_jacobianMatrix);
-            m_rank = solver.rank();
-        }
     }
 
     NetworkJacobian::NetworkJacobian(
@@ -57,7 +58,7 @@ namespace gridfire {
 
     double NetworkJacobian::operator()(const size_t i, const size_t j) const {
         if (i >= m_jacobianMatrix.rows() || j >= m_jacobianMatrix.cols()) {
-            throw std::out_of_range("Index out of bounds in NetworkJacobian operator().");
+            throw std::out_of_range(std::format("Index ({}, {}) out of bounds in NetworkJacobian operator() for jacobian of shape ({}, {}).", i, j, m_jacobianMatrix.rows(), m_jacobianMatrix.cols()));
         }
         return m_jacobianMatrix.coeff(i, j);
     }
@@ -73,7 +74,7 @@ namespace gridfire {
 
     void NetworkJacobian::set(const size_t i, const size_t j, const double value) {
         if (i >= m_jacobianMatrix.rows() || j >= m_jacobianMatrix.cols()) {
-            throw std::out_of_range("Index out of bounds in NetworkJacobian set().");
+            throw std::out_of_range(std::format("Index ({}, {}) out of bounds in NetworkJacobian set() for jacobian of shape ({}, {}).", i, j, m_jacobianMatrix.rows(), m_jacobianMatrix.cols()));
         }
         m_jacobianMatrix.coeffRef(i, j) = value;
     }
@@ -91,7 +92,17 @@ namespace gridfire {
     }
 
     size_t NetworkJacobian::rank() const {
-        return m_rank;
+        if (m_rank) {
+            return m_rank.value();
+        }
+        if (m_jacobianMatrix.rows() == 0 || m_jacobianMatrix.cols() == 0) {
+            m_rank = 0;
+        } else {
+            Eigen::SparseQR<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>> solver;
+            solver.compute(m_jacobianMatrix);
+            m_rank = solver.rank();
+        }
+        return m_rank.value();
     }
 
     bool NetworkJacobian::singular() const {
@@ -157,4 +168,75 @@ namespace gridfire {
     const std::unordered_map<fourdst::atomic::Species, size_t> & NetworkJacobian::mapping() const {
         return m_speciesToIndexMap;
     }
+
+    void NetworkJacobian::to_csv(const std::string &filename) const {
+        std::ofstream file(filename);
+        if (file.is_open()) {
+            const size_t numSpecies = m_jacobianMatrix.rows();
+            size_t i = 0;
+            file << "# Species: ";
+            for (size_t k = 0; k < numSpecies; ++k) {
+                fourdst::atomic::Species species = std::ranges::find_if(
+                    m_speciesToIndexMap,
+                    [k](const auto& pair) {
+                        return pair.second == k;
+                    })->first;
+                file << species.name();
+                if (i < numSpecies - 1) {
+                    file << ",";
+                }
+                i++;
+            }
+            const std::chrono::system_clock::time_point now = std::chrono::system_clock::now();
+            const std::time_t now_c = std::chrono::system_clock::to_time_t(now);
+
+            file << "\n";
+            file << "# Rows: " << numSpecies << ", Columns: " << numSpecies << "\n";
+            file << "# Generated on " << std::ctime(&now_c) << "\n";
+            file << "# Generated by GridFire NetworkJacobian::to_csv\n";
+
+            for (size_t row = 0; row < numSpecies; ++row) {
+                for (size_t col = 0; col < numSpecies; ++col) {
+                    file << this->operator()(row, col);
+                    if (col < numSpecies - 1) {
+                        file << ",";
+                    }
+                }
+                file << "\n";
+            }
+            file.close();
+        }
+    }
+
+    NetworkJacobian regularize_jacobian(
+        const NetworkJacobian& jacobian,
+        const fourdst::composition::CompositionAbstract& comp,
+        const std::optional<quill::Logger*> logger
+    ) {
+        const std::vector<JacobianEntry> infs = jacobian.infs();
+        const std::vector<JacobianEntry> nans = jacobian.nans();
+        if (infs.size() == 0 && nans.size() == 0) {
+            return jacobian;
+        }
+
+        NetworkJacobian newJacobian = jacobian;
+        for (const auto& [iSp, dSp] : infs | std::views::keys) {
+            if (comp.getMolarAbundance(iSp) < MIN_ABUNDANCE_TO_CONTRIBUTE_TO_JACOBIAN || comp.getMolarAbundance(dSp) < MIN_ABUNDANCE_TO_CONTRIBUTE_TO_JACOBIAN) {
+                newJacobian.set(iSp, dSp, 0.0);
+                if (logger) {
+                    LOG_TRACE_L1(logger.value(), "Regularized Jacobian entry ({}, {}) from inf to 0.0 due to low abundance.", iSp.name(), dSp.name());
+                }
+            }
+        }
+        for (const auto& [iSp, dSp] : nans | std::views::keys) {
+            if (comp.getMolarAbundance(iSp) < MIN_ABUNDANCE_TO_CONTRIBUTE_TO_JACOBIAN || comp.getMolarAbundance(dSp) < MIN_ABUNDANCE_TO_CONTRIBUTE_TO_JACOBIAN) {
+                newJacobian.set(iSp, dSp, 0.0);
+                if (logger) {
+                    LOG_TRACE_L1(logger.value(), "Regularized Jacobian entry ({}, {}) from inf to 0.0 due to low abundance.", iSp.name(), dSp.name());
+                }
+            }
+        }
+
+        return newJacobian;
+    }
 }