Matrix.cpp

 
/* SPDX-License-Identifier: CC0-1.0 */
#include <vector>
#include <iostream>
#include <string>
#include "jama/jama_svd.h"
#include "jama/jama_eig.h"
#include "jama/jama_lu.h"
 
#include "Constants.h"
#include "IException.h"
#include "IString.h"
#include "Matrix.h"
 
namespace Isis {
  Matrix::Matrix(const int n, const int m, const double value) {
    if(n < 1 || m < 1) {
      std::string m = "Invalid matrix dimensions";
      throw IException(IException::Programmer, m, _FILEINFO_);
    }
    p_matrix = TNT::Array2D<double>(n, m, value);
  }
 
  Matrix::Matrix(TNT::Array2D<double> matrix) {
    p_matrix = matrix.copy();
  }
 
  Matrix::~Matrix() {}
 
  Matrix Matrix::Identity(const int n) {
    if(n < 1) {
      std::string m = "Invalid matrix dimensions";
      throw IException(IException::Programmer, m, _FILEINFO_);
    }
    Matrix identity(n, n);
    for(int i = 0; i < identity.Rows(); i++) {
      identity[i][i] = 1.0;
    }
    return identity;
  }
 
  double Matrix::Determinant() {
    if(Rows() != Columns()) {
      std::string m = "Unable to calculate the determinant, the matrix is not square.";
      throw IException(IException::Programmer, m, _FILEINFO_);
    }
    JAMA::LU<double> lu(p_matrix);
    return lu.det();
  }
 
  double Matrix::Trace() {
    if(Rows() != Columns()) {
      std::string m = "Unable to calculate the trace, the matrix is not square.";
      throw IException(IException::Programmer, m, _FILEINFO_);
    }
    double trace = 0.0;
    for(int i = 0; i < Rows(); i++) {
      trace += p_matrix[i][i];
    }
    return trace;
  }
 
  Matrix Matrix::Multiply(Matrix &matrix) {
    if(Columns() != matrix.Rows()) {
      std::string m = "Incompatible matrix dimensions, cannot multiply the matrices.";
      throw IException(IException::Programmer, m, _FILEINFO_);
    }
    TNT::Array2D<double> m(matrix.Rows(), matrix.Columns());
    for(int i = 0; i < m.dim1(); i++) {
      for(int j = 0; j < m.dim2(); j++) {
        m[i][j] = matrix[i][j];
      }
    }
    return Matrix(TNT::matmult(p_matrix, m));
  }
 
  Matrix Matrix::Add(Matrix &matrix) {
    if(Rows() != matrix.Rows() || Columns() != matrix.Columns()) {
      std::string m = "Incompatible matrix dimensions, cannot add the matrices.";
      throw IException(IException::Programmer, m, _FILEINFO_);
    }
    TNT::Array2D<double> m(matrix.Rows(), matrix.Columns());
    for(int i = 0; i < m.dim1(); i++) {
      for(int j = 0; j < m.dim2(); j++) {
        m[i][j] = matrix[i][j];
      }
    }
    return Matrix(p_matrix + m);
  }
 
  Matrix Matrix::Subtract(Matrix &matrix) {
    if(Rows() != matrix.Rows() || Columns() != matrix.Columns()) {
      std::string m = "Incompatible matrix dimensions, cannot subtract the matrices.";
      throw IException(IException::Programmer, m, _FILEINFO_);
    }
    TNT::Array2D<double> m(matrix.Rows(), matrix.Columns());
    for(int i = 0; i < m.dim1(); i++) {
      for(int j = 0; j < m.dim2(); j++) {
        m[i][j] = matrix[i][j];
      }
    }
 
    return Matrix(p_matrix - m);
  }
 
  Matrix Matrix::MultiplyElementWise(Matrix &matrix) {
    if(Rows() != matrix.Rows() || Columns() != matrix.Columns()) {
      std::string m = "Incompatible matrix dimensions, cannot multiply the matrices.";
      throw IException(IException::Programmer, m, _FILEINFO_);
    }
    TNT::Array2D<double> m(matrix.Rows(), matrix.Columns());
    for(int i = 0; i < m.dim1(); i++) {
      for(int j = 0; j < m.dim2(); j++) {
        m[i][j] = matrix[i][j];
      }
    }
    return Matrix(p_matrix * m);
  }
 
  Matrix Matrix::Multiply(double scalar) {
    Matrix product(Rows(), Columns());
    for(int i = 0; i < Rows(); i++) {
      for(int j = 0; j < Columns(); j++) {
        product[i][j] = p_matrix[i][j] * scalar;
      }
    }
    return product;
  }
 
  Matrix Matrix::Transpose() {
    TNT::Array2D<double> transpose(p_matrix.dim2(), p_matrix.dim1());
    for(int i = 0; i < transpose.dim1(); i++) {
      for(int j = 0; j < transpose.dim2(); j++) {
        transpose[i][j] = p_matrix[j][i];
      }
    }
    return Matrix(transpose);
  }
 
  Matrix Matrix::Inverse() {
    if(Rows() != Columns()) {
      std::string m = "Unable to calculate the inverse, the matrix is not square.";
      throw IException(IException::Programmer, m, _FILEINFO_);
    }
    TNT::Array2D<double> id(p_matrix.dim1(), p_matrix.dim2(), 0.0);
    for(int i = 0; i < p_matrix.dim1(); i++) id[i][i] = 1;
 
    JAMA::LU<double> lu(p_matrix);
    if(lu.det() == 0.0) {
      std::string m = "Cannot take the inverse of the matrix";
      throw IException(IException::Programmer, m, _FILEINFO_);
    }
 
    return Matrix(lu.solve(id));
  }
 
  std::vector<double> Matrix::Eigenvalues() {
    if(Rows() != Columns()) {
      std::string m = "Unable to calculate eigenvalues, the matrix is not square.";
      throw IException(IException::Programmer, m, _FILEINFO_);
    }
    JAMA::Eigenvalue<double> E(p_matrix);
    TNT::Array2D<double> D;
    E.getD(D);
 
    std::vector<double> eigenvalues(D.dim1());
    for(int i = 0; i < D.dim1(); i++) {
      eigenvalues[i] = D[i][i];
    }
 
    return eigenvalues;
  }
 
  Matrix Matrix::Eigenvectors() {
    if(Rows() != Columns()) {
      std::string m = "Unable to calculate eigenvectors, the matrix is not square.";
      throw IException(IException::Programmer, m, _FILEINFO_);
    }
    JAMA::Eigenvalue<double> E(p_matrix);
    TNT::Array2D<double> V;
    E.getV(V);
 
    return Matrix(V);
  }
 
  ostream &operator<<(ostream &os, Matrix &matrix) {
    for(int i = 0; i < matrix.Rows(); i++) {
      for(int j = 0; j < matrix.Columns(); j++) {
        os << matrix[i][j];
        if(j < matrix.Columns() - 1) os << " ";
      }
      if(i < matrix.Rows() - 1) os << endl;
    }
    return os;
  }
} // end namespace isis