LinearAlgebra.h

Go to the documentation of this file.

#pragma once
 
#include <limits>
#include <iostream>
#include <vector>
#include <cmath>
 
#include <TMatrixD.h>
#include <TVectorD.h>
#include <TDecompSVD.h>
 
 
const double EPSILON = 1e-9;
 
 
namespace LA
{
    //___________________________________________________________________________
    std::vector<double> ROOTEquationSolver(const double& m1, const double& b1,  
                                           const double& m2, const double& b2,  
                                           const double& m3, const double& b3,  
                                           bool debugMode)
    {
        std::vector<std::vector<double>> mat = {
            {-m1,   0,  1},
            {  0, -m2,  1},
            {-m3,   1,  0}
        };
        
        // Construct a TMatrixD equivalent of the 2D std::vector
        const int rows = mat.size();
        const int cols = mat.empty() ? 0 : mat[0].size();
        TMatrixD rootMatrix(rows, cols);
        for (int i = 0; i < rows; ++i) {
            for (int j = 0; j < cols; ++j) {
                rootMatrix(i, j) = mat[i][j];
            }
        }
        
        TDecompSVD svd(rootMatrix);
        TVectorD b(3); b(0) = b1; b(1) = b2; b(2) = b3;
        Bool_t is_good_solution;
        
        TVectorD solution(b.GetNrows());
 
        solution = svd.Solve(b, is_good_solution);
 
        if (!is_good_solution)  
        {
            if (debugMode)
                std::cerr<<"\n\nWARNING: The system cannot be solved or is ill-conditioned.\n"; 
            return {}; 
        }
 
        std::vector<double> solutionVec = VectorUtilities::TVectorDToStdVector(solution);
        return solutionVec;
    }
    
    
    
    //___________________________________________________________________________
    void GaussianElimination(std::vector<std::vector<double>>& mat, bool debugMode = false)
    {
        if (debugMode)
        {
            std::cout<<"\nBefore Gaussian elimination\n";
            VectorUtilities::printVectorOfVectors(mat);
            std::cout<<"\n";
        }
        
        int n = mat.size();
 
        for (int i = 0; i < n; ++i)
        {
            int pivotRow = i;
            for (int row = i + 1; row < n; ++row)
            {
                if (std::fabs(mat[row][i]) > std::fabs(mat[pivotRow][i]))
                    pivotRow = row;
            }
            std::swap(mat[pivotRow], mat[i]);
 
            if (std::fabs(mat[i][i]) < EPSILON) continue;
 
            for (int col = i + 1; col <= n; ++col)
            {
                mat[i][col] /= mat[i][i];
            }
            mat[i][i] = 1.0;
 
            for (int row = i + 1; row < n; ++row)
            {
                double factor = mat[row][i];
                for (int col = i; col <= n; ++col) {
                    mat[row][col] -= factor * mat[i][col];
                }
            }
        }
        
        if (debugMode)
        {
            std::cout<<"\nAfter Gaussian elimination\n";
            VectorUtilities::printVectorOfVectors(mat);
            std::cout<<"\n";
        }
    }
 
 
    //___________________________________________________________________________
    int GetMatrixRank(const std::vector<std::vector<double>>& mat)
    {
        int rank = 0;
        for (const auto& row : mat) {
            for (double val : row) {
                if (std::fabs(val) > EPSILON) {
                    ++rank;
                    break;
                }
            }
        }
        return rank;
    }
 
 
    //___________________________________________________________________________
    int CheckSolutions(const std::vector<std::vector<double>>& mat, std::vector<double>& solution)
    {
        
        int n = mat.size();
        int rankA = GetMatrixRank(mat);
 
        std::vector<std::vector<double>> augmentedMat(mat);
        for (int i = 0; i < n; ++i) {
            augmentedMat[i].pop_back();  
        }
        
        int rankAugmented = GetMatrixRank(augmentedMat);
 
        if (rankA < rankAugmented)
            return -1;               
        else if (rankA < n)
            return 0;                
 
        solution.assign(n, 0);
        for (int i = n - 1; i >= 0; --i) {
            solution[i] = mat[i][n];
            
            for (int j = i + 1; j < n; ++j) {
                solution[i] -= mat[i][j] * solution[j];
            }
            
            if (std::fabs(mat[i][i]) > EPSILON) {
                solution[i] /= mat[i][i];
            }
            
            else if (std::fabs(solution[i]) > EPSILON) {
                return -1; 
            }
        }
        return 1; 
    }
 
 
    //___________________________________________________________________________
    bool ValidateSolution(const std::vector<std::vector<double>>& origMat,
                          const std::vector<double>& solution,
                          bool debugMode)
    {
        
        int n = origMat.size();
        for (int i = 0; i < n; ++i)
        {
            double sum = 0.0;
            for (int j = 0; j < n; ++j) {
                sum += origMat[i][j] * solution[j];
            }
            
            debugMode = false;
            
            if (debugMode)
                std::cout<<"For equation "<<i<<", sum(0, n-1) = "<<sum<<", origMat[i][n] = "<<origMat[i][n]<<std::endl;
            if (std::fabs(sum - origMat[i][n]) > EPSILON)
                return false; 
        }
        return true; 
    }
 
 
    //___________________________________________________________________________
    bool IsSlopeUndefined(double slope)
    {
        return std::abs(slope) > std::numeric_limits<double>::max();
    }
 
 
    //___________________________________________________________________________
    std::vector<double> SolveSystemOfThreeEquations(const double& m1, const double& b1,  
                                                    const double& m2, const double& b2,  
                                                    const double& m3, const double& b3,  
                                                    bool debugMode)
    {
        std::vector<std::vector<double>> mat = {
            {-m1,   0,  1,  b1},
            {  0, -m2,  1,  b2},
            {-m3,   1,  0,  b3}
        };
 
        // Store the original matrix for validation
        std::vector<std::vector<double>> origMat = mat;
 
        // Perform Gaussian elimination
        GaussianElimination(mat);
 
        std::vector<double> solution;
        int result = CheckSolutions(mat, solution);
        
        if (result == 1)
        {
            if (debugMode)
            {
                std::cout<<"\nThe system has a unique solution:"<<std::endl;
                for (double s : solution) {
                    std::cout << s << " ";
                }
                std::cout<<std::endl;
            }
 
            // Validate the solution
            if (ValidateSolution(origMat, solution, debugMode))
            {
                if (debugMode)
                    std::cout<<"The solution is valid."<<std::endl;
                return solution;
            }
 
            else
            {
                if (debugMode)
                    std::cout<<"\n\nERROR: The solution is invalid.\nThis should never happen! Please check logic!\n\n"<<std::endl;
                return {-9.0};
            }
        }
        
        else if (result == 0)
        {
            if (debugMode)
                std::cout<<"\nThe system has infinite solutions."<<std::endl;
            return {1.0};
        }
        
        else
        {
            if (debugMode)
                std::cout<<"\nThe system has no solution."<<std::endl;
            return {0.0};
        }
    }
}