Add construct and destroy methods to custom_allocator

nhance custom allocator for Boost compatibility and improve numerical stability in polyfit function

omnn/extrapolator/Extrapolator.cpp

@@ -7,33 +7,43 @@
 #include "math/Integer.h"
 #include "math/Sum.h"
 
+#include <boost/numeric/ublas/storage.hpp>
+
+
 using namespace omnn;
 using namespace math;
+using namespace boost::numeric::ublas;
 
 
-auto det_fast(extrapolator_base_matrix matrix)
-{
-    using T = extrapolator_base_matrix::value_type;
-    ublas::permutation_matrix<std::size_t> pivots(matrix.size1());
+// Ensure that matrix_t, vector_t, and permutation_matrix_t are using unbounded_array_wrapper with custom_allocator
+using matrix_t = boost::numeric::ublas::matrix<Valuable, boost::numeric::ublas::basic_row_major<>, unbounded_array_wrapper<Valuable, custom_allocator<Valuable>>>;
+using vector_t = boost::numeric::ublas::vector<Valuable, unbounded_array_wrapper<Valuable, custom_allocator<Valuable>>>;
+using permutation_matrix_t = boost::numeric::ublas::permutation_matrix<std::size_t, unbounded_array_wrapper<std::size_t, custom_allocator<std::size_t>>>;
+
+// This function calculates the determinant of a matrix using LU factorization
+auto det_fast(matrix_t matrix) {
+    // Use the matrix's size1 directly to avoid calling size on the allocator
+    auto matrix_size = matrix.size1();
+    permutation_matrix_t pivots(matrix_size);
 
-    auto isSingular = ublas::lu_factorize(matrix, pivots);
+    // Perform LU factorization on a copy of the matrix to preserve the original matrix
+    matrix_t matrix_copy(matrix);
+    auto isSingular = ublas::lu_factorize(matrix_copy, pivots);
     if (isSingular)
-        return T(0);
+        return Valuable(0);
 
-    T det = 1;
-    for (std::size_t i = 0; i < pivots.size(); ++i)
-    {
+    Valuable det(1);
+    for (std::size_t i = 0; i < matrix_size; ++i) {
         if (pivots(i) != i)
-            det *= static_cast<double>(-1);
+            det *= Valuable(-1);
 
-        det *= matrix(i, i);
+        det *= matrix_copy(i, i);
     }
 
     return det;
 }
 
-bool Extrapolator::Consistent(const extrapolator_base_matrix& augment)
-{
+bool Extrapolator::Consistent(const matrix_t& augment) {
     return Determinant() == det_fast(augment);
 }
 
@@ -97,45 +107,44 @@ Extrapolator::Extrapolator(std::initializer_list<std::vector<T>> dependancy_matr
         ++r;
     }
 }
-Extrapolator::solution_t Extrapolator::Solve(const ublas::vector<T>& augment) const
+Extrapolator::solution_t Extrapolator::Solve(const vector_t& augment) const
 {
-    auto e = *this;
+    auto e = static_cast<const matrix_t&>(*this);
     auto sz1 = size1();
     auto sz2 = size2();
 
-    ublas::vector<T> a(sz2);
-    const ublas::vector<T>* au = &augment;
+    // Initialize vector with size and default value
+    vector_t a(sz2, Valuable(0));
+    const vector_t* au = &augment;
     if (sz1 > sz2 + 1 /*augment*/) {
         // make square matrix to make it solvable by boost ublas
         e = Extrapolator(sz2, sz2);
         // sum first equations
-        a[0] = 0;
         for (auto i = sz2; i--;) {
-            e(0, i) = 0;
+            e.insert_element(0, i, Valuable(0));
         }
         auto d = sz1 - sz2;
         for (auto i = d; i--;) {
             for (auto j = sz2; j--;) {
                 e(0, j) += operator()(i, j);
             }
-            a[0] += augment[i];
+            a(i - d) = (*au)(i); // Use vector's operator() provided by Boost Ublas for element access
         }
 
         for (auto i = d; i < sz1; ++i) {
             for (auto j = sz2; j--;) {
                 e(i - d, j) = operator()(i, j);
             }
-            a[i - d] = augment[i];
+            a(i - d) = (*au)(i); // Use vector's operator() provided by Boost Ublas for element access
         }
         au = &a;
     }
     solution = ublas::solve(e, *au, ublas::upper_tag());
     return solution;
 }
 
-Extrapolator::T Extrapolator::Determinant() const
-{
-    ublas::permutation_matrix<std::size_t> pivots(this->size1());
+Extrapolator::T Extrapolator::Determinant() const {
+    permutation_matrix_t pivots(this->size1());
     auto mLu = *this;
     auto isSingular = ublas::lu_factorize(mLu, pivots);
     if (isSingular)
@@ -152,12 +161,11 @@ Extrapolator::T Extrapolator::Determinant() const
     return det;
 }
 
-bool Extrapolator::Consistent(const ublas::vector<T>& augment)
-{
+bool Extrapolator::Consistent(const vector_t& augment) {
     auto sz = augment.size();
-    extrapolator_base_matrix augmentedMatrix(sz, 1);
+    matrix_t augmentedMatrix(sz, 1);
     for (auto i = sz; i-- > 0;) {
-        augmentedMatrix(i, 0) = augment[i];
+        augmentedMatrix(i, 0) = augment(i);
     }
     return Consistent(augmentedMatrix);
 }

omnn/extrapolator/Extrapolator.h

@@ -23,7 +23,7 @@ namespace omnn::math {
     namespace ublas = boost::numeric::ublas;
     //using extrapolator_value_type = a_rational;
     using extrapolator_value_type = Valuable;
-    using extrapolator_base_matrix = boost::numeric::ublas::matrix<extrapolator_value_type>;
+    using extrapolator_base_matrix = ublas::matrix<extrapolator_value_type>;
 
 
 class Extrapolator
@@ -37,21 +37,30 @@ class Extrapolator
     mutable solution_t solution;
 
 public:
+    using vector_t = ublas::vector<T, custom_allocator<T>>;
+    using matrix_t = ublas::matrix<T, ublas::row_major, custom_allocator<T>>;
+
     using base::base;
 
 	Extrapolator() = default;
 
     Extrapolator(std::initializer_list<std::vector<T>> dependancy_matrix);
 
-    solution_t Solve(const ublas::vector<T>& augment) const;
+    solution_t Solve(const vector_t& augment) const;
 
     T Determinant() const;
 
     /**
      * If possible, make the matrix consistent
      * @return true if it is
      * **/
-    bool Consistent(const ublas::vector<T>& augment);
+    bool Consistent(const vector_t& augment);
+
+    /**
+     * If possible, make the matrix consistent
+     * @return true if it is
+     * **/
+    bool Consistent(const matrix_t& augment);
 
     /**
      * If possible, make the matrix consistent

omnn/math/Polyfit.h

@@ -40,16 +40,27 @@ auto polyfit( const std::vector<T>& oX, const std::vector<T>& oY, size_t nDegree
 {
     using namespace boost::numeric::ublas;
 
+    // Define a custom allocator type for the Boost Numeric Ublass types
+    using custom_alloc_t = custom_allocator<T>;
+
+    // Use the custom allocator for the matrix and vector types
+    using matrix_t = matrix<T, basic_row_major<>, unbounded_array<T, custom_alloc_t>>;
+    using vector_t = vector<T, unbounded_array<T, custom_alloc_t>>;
+
+    // Define a custom allocator type for int
+    using custom_int_alloc_t = custom_allocator<int>;
+
+    // Use the custom allocator for the permutation_matrix type
+    using permutation_matrix_t = permutation_matrix<int, unbounded_array<int, custom_int_alloc_t>>;
     if ( oX.size() != oY.size() )
         throw std::invalid_argument( "X and Y vector sizes do not match" );
 
     // more intuitive this way
     nDegree++;
 
     auto nCount =  oX.size();
-    matrix<T> oXMatrix( nCount, nDegree );
-    matrix<T> oYMatrix( nCount, 1 );
-
+    matrix_t oXMatrix( nCount, nDegree );
+    matrix_t oYMatrix( nCount, 1 );
     // copy y matrix
     for ( size_t i = 0; i < nCount; i++ )
     {
@@ -67,22 +78,125 @@ auto polyfit( const std::vector<T>& oX, const std::vector<T>& oY, size_t nDegree
         }
     }
 
+    // Scale input data
+    T xMean = std::accumulate(oX.begin(), oX.end(), T(0)) / oX.size();
+    T xStd = std::sqrt(std::inner_product(oX.begin(), oX.end(), oX.begin(), T(0)) / oX.size() - xMean * xMean);
+    T yMean = std::accumulate(oY.begin(), oY.end(), T(0)) / oY.size();
+    T yStd = std::sqrt(std::inner_product(oY.begin(), oY.end(), oY.begin(), T(0)) / oY.size() - yMean * yMean);
+
+    std::cout << "Debug: xMean = " << xMean << ", xStd = " << xStd << ", yMean = " << yMean << ", yStd = " << yStd << std::endl;
+
+    for (size_t i = 0; i < nCount; ++i) {
+        for (size_t j = 0; j < nDegree; ++j) {
+            oXMatrix(i, j) = std::pow(oX[i] - xMean, j) / std::pow(xStd, j);
+        }
+        oYMatrix(i, 0) = (oY[i] - yMean) / yStd;
+    }
     // transpose X matrix
-    matrix<T> oXtMatrix( trans(oXMatrix) );
+    matrix_t oXtMatrix( trans(oXMatrix) );
     // multiply transposed X matrix with X matrix
-    matrix<T> oXtXMatrix( prec_prod(oXtMatrix, oXMatrix) );
+    matrix_t oXtXMatrix( prec_prod(oXtMatrix, oXMatrix) );
     // multiply transposed X matrix with Y matrix
-    matrix<T> oXtYMatrix( prec_prod(oXtMatrix, oYMatrix) );
+    matrix_t oXtYMatrix( prec_prod(oXtMatrix, oYMatrix) );
+
+    // Add regularization to handle near-singular matrices
+    const T lambda = 1e-8; // Adjusted regularization parameter
+    for (size_t i = 0; i < oXtXMatrix.size1(); ++i) {
+        oXtXMatrix(i, i) += lambda;
+    }
+
+    std::cout << "Debug: Regularization parameter lambda = " << lambda << std::endl;
+
+    // Solve the system using LU decomposition
+    permutation_matrix<std::size_t> pm(oXtXMatrix.size1());
+    lu_factorize(oXtXMatrix, pm);
+    lu_substitute(oXtXMatrix, pm, oXtYMatrix);
+
+    std::cout << "Debug: Coefficients:" << std::endl;
+    for (size_t i = 0; i < oXtYMatrix.size1(); ++i) {
+        std::cout << oXtYMatrix(i, 0) << " ";
+    }
+    std::cout << std::endl;
+
+    // Check for near-zero determinant
+    T det = 1;
+    for (size_t i = 0; i < oXtXMatrix.size1(); ++i) {
+        det *= oXtXMatrix(i, i);
+    }
+    std::cout << "Debug: Determinant before LU decomposition: " << det << std::endl;
 
+    // Check for singularity after regularization
+    det = 1;
+    for (size_t i = 0; i < oXtXMatrix.size1(); ++i) {
+        det *= oXtXMatrix(i, i);
+    }
+    std::cout << "Debug: Determinant after regularization: " << det << std::endl;
+
+    // Check input data
+    std::cout << "Debug: Input data (X, Y):" << std::endl;
+    for (size_t i = 0; i < nCount; ++i) {
+        std::cout << "(" << oX[i] << ", " << oY[i] << ") ";
+    }
+    std::cout << std::endl;
+
+    // Scale input data
+    T x_mean = std::accumulate(oX.begin(), oX.end(), T(0)) / nCount;
+    T x_std = std::sqrt(std::inner_product(oX.begin(), oX.end(), oX.begin(), T(0)) / nCount - x_mean * x_mean);
+    T y_mean = std::accumulate(oY.begin(), oY.end(), T(0)) / nCount;
+    T y_std = std::sqrt(std::inner_product(oY.begin(), oY.end(), oY.begin(), T(0)) / nCount - y_mean * y_mean);
+
+    std::vector<T> scaled_oX(oX), scaled_oY(oY);
+    for (size_t i = 0; i < nCount; ++i) {
+        scaled_oX[i] = (oX[i] - x_mean) / x_std;
+        scaled_oY[i] = (oY[i] - y_mean) / y_std;
+    }
+
+    // Tikhonov regularization parameter
+    // lambda is already declared earlier, so we remove the redeclaration here
+
+    // Add regularization term to oXtXMatrix
+    for (size_t i = 0; i < oXtXMatrix.size1(); ++i) {
+        oXtXMatrix(i, i) += lambda;
+    }
     // lu decomposition
-    permutation_matrix<int> pert(oXtXMatrix.size1());
+    permutation_matrix_t pert(oXtXMatrix.size1());
     auto singular = lu_factorize(oXtXMatrix, pert);
-    // must be singular
+    std::cout << "Debug: LU factorization result: " << singular << std::endl;
+
+    // Check matrix after LU decomposition
+    std::cout << "Debug: oXtXMatrix after LU decomposition:" << std::endl;
+    for (size_t i = 0; i < oXtXMatrix.size1(); ++i) {
+        for (size_t j = 0; j < oXtXMatrix.size2(); ++j) {
+            std::cout << oXtXMatrix(i, j) << " ";
+        }
+        std::cout << std::endl;
+    }
+
+    // must not be singular
     BOOST_ASSERT( singular == 0 );
 
+    std::cout << "Debug: oXtYMatrix before substitution:" << std::endl;
+    for (size_t i = 0; i < oXtYMatrix.size1(); ++i) {
+        std::cout << oXtYMatrix(i, 0) << " ";
+    }
+    std::cout << std::endl;
     // backsubstitution
     lu_substitute(oXtXMatrix, pert, oXtYMatrix);
 
+    std::cout << "Debug: Resulting coefficients:" << std::endl;
+    for (size_t i = 0; i < oXtYMatrix.size1(); ++i) {
+        std::cout << oXtYMatrix(i, 0) << " ";
+    }
+    std::cout << std::endl;
+
+    // Unscale the coefficients
+    for (size_t i = 0; i < oXtYMatrix.size1(); ++i) {
+        T x_std_power = T(1);
+        for (size_t j = 0; j < i; ++j) {
+            x_std_power *= x_std;
+        }
+        oXtYMatrix(i, 0) = oXtYMatrix(i, 0) * y_std / x_std_power;
+    }
     // copy the result to coeff
     return std::vector<T>( oXtYMatrix.data().begin(), oXtYMatrix.data().end() );
 }