ExampleSparseMatrix.cpp

/*
This file is part of the BIAS library (Basic ImageAlgorithmS).

Copyright (C) 2003-2009    (see file CONTACT for details)
Multimediale Systeme der Informationsverarbeitung
Institut fuer Informatik
Christian-Albrechts-Universitaet Kiel


BIAS is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version.

BIAS is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License
along with BIAS; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/



/** @example ExampleSparseMatrix.cpp
    @relates SparseMatrix
    @brief Test implementation of sparse matrix operations.
    @see BIAS::SparseMatrix
    @ingroup g_examples
    @author esquivel
    @date 11/2008
*/

#include <iostream>
#include <Base/Math/Random.hh>
#include <Base/Math/SparseMatrix.hh>
#include <MathAlgo/SVD.hh>

using namespace std;
using namespace BIAS;

/** brief Helper function to get rid of small matrix entries. */
BIAS::Matrix<double> CropMatrix(const BIAS::Matrix<double> &M,
                                double epsilon = 1e-15)
{
  BIAS::Matrix<double> A(M);
  const unsigned int n = A.GetRows();
  const unsigned int m = A.GetCols();
  for (unsigned int i = 0; i < n; i++) {
    for (unsigned int j = 0; j < m; j++) {
      if (fabs(A[i][j]) < epsilon) A[i][j] = 0.0;
    }
  }
  return A;
}

int main()
{

  cout << "-- Testing base operations -- " << endl;

  BIAS::Matrix<double> A(6, 6, MatrixIdentity);
  A[2][1] = 17;
  A[3][5] = -0.1;
  A[4][4] = 0.0;
  BIAS::SparseMatrix S(A);

  cout << "- Dense is A = " << A << endl;
  cout << "- Sparse is S = " << endl;
  S.Print("S", cout);

  cout << "- Modifying S (consistency check)... " << flush;
  S.AddToDiagonal(2.0);
  cout << (S.CheckConsistency() ? "OK, " : "Not OK!, ") << flush;
  S.MultiplyDiagonalBy(0.5);
  cout << (S.CheckConsistency() ? "OK, " : "Not OK!, ") << flush;
  S.SetElement(1, 1, 0.0);
  cout << (S.CheckConsistency() ? "OK, " : "Not OK!, ") << flush;
  S.SetElement(1, 4, 1.0);
  cout << (S.CheckConsistency() ? "OK" : "Not OK!")
       << endl << flush;
//
//  S.SetElement(6,2,1.0);
//  S.SetElement(6,2,2.0);
//  cout << (S.CheckConsistency() ? "OK" : "Not OK!")
//  << endl << flush;
  
  cout << "sparse S after setting single elements is " << S << endl << flush;
  
  SparseMatrix SAddTest(2,2);
  SAddTest.SetElement(0,0,1.0);
  SAddTest.SetElement(1,1,1.0);
  cout << "s before adding " << SAddTest << endl;
  Matrix<double> addBottom(2,3,MatrixIdentity);
  addBottom[0][2] = 2;
  Matrix<double> addRight(3,2,MatrixIdentity);
  addRight[2][0] = 3;
  SparseMatrix addBottomSparse(addBottom);
  SparseMatrix addRightSparse(addRight);
  
  cout << "adding to bottom " << addBottomSparse  << addBottomSparse.GetRows() << " " << addBottomSparse.GetCols() << endl;
  cout << "adding to right " << addRightSparse << addRightSparse.GetRows() << " " << addRightSparse.GetCols() <<endl;
  SAddTest.AppendMatrixBottom(addBottomSparse);
  cout << "s after adding to bottom " << SAddTest << " size " << SAddTest.GetRows() << " " << SAddTest.GetCols() << endl;
  SAddTest.AppendMatrixRight(addRightSparse);
  cout << "s after adding to right " << SAddTest << " size " << SAddTest.GetRows() << " " << SAddTest.GetCols() << endl;
  
  cout << (SAddTest.CheckConsistency() ? "OK" : "Not OK!")
  << endl << flush;
  
  BIAS::Matrix<double> B;
  S.GetAsDense(B);

  cout << "- Modified sparse is S = " << endl;
  S.Print("S", cout);
  cout << "- Conversion to dense is B = " << B << endl;

  BIAS::Matrix<double> BT, BTB;
  BT = B.Transpose();
  BTB = BT * B;
  BIAS::SparseMatrix ST(S.GetColsNum(), S.GetRowsNum()), STS(ST.GetRowsNum(), S.GetColsNum());
  S.Transpose(ST);
  ST.Multiply(S, STS);

  cout << "- Sparse multiplication S'S = " << endl;
  STS.Print("S'S", cout);
  cout << "- Dense multiplication B'B = " << BTB << endl;

  double m;
  unsigned int i, j;
  m = STS.GetMaxDiagonalElement(i);
  cout << "- Maximal diagonal element of S'S is [" << i << "][" << i
       << "] = " << m << endl;
  j = 4;
  m = S.GetMaxRowElement(j, i);
  cout << "- Maximal element of row " << j << " in S is [" << j << "][" << i
       << "] = " << m << endl;
  m = ST.GetMaxColumnElement(j, i);
  cout << "- Maximal element of column " << j << " of S' is [" << i
       << "][" << j << "] = " << m << endl;

  cout << "- Checking consistency of S... ";
  cout << (S.CheckConsistency() ? "OK!" : "failed!") << endl;
  cout << "- Checking consistency of S'... ";
  cout << (ST.CheckConsistency() ? "OK!" : "failed!") << endl;
  cout << "- Checking consistency of S'S... ";
  cout << (STS.CheckConsistency() ? "OK!" : "failed!") << endl;


  cout << endl << "-- Testing matrix inversion -- " << endl << endl;

  BIAS::Matrix<double> M(6, 6, MatrixZero);
  M[0][0] =  2.0;
  M[1][1] =  4.0;
  M[2][2] = -1.0;
  M[3][3] = -3.0;
  M[4][4] =  0.5;
  M[5][5] =  1.0;
  M[3][4] =  1.5;
  M[2][4] =  0.1;
  M[3][1] =  2.5;
  BIAS::SparseMatrix N(M);
  BIAS::SparseMatrix Nold(M);
  N.SetElement(5, 5, 1.0);
  cout << "- Dense is M = " << M << endl;
  cout << "- Sparse is N = " << endl;
  N.Print("N", cout);
  cout << "- Checking consistency of N... ";
  cout << (N.CheckConsistency() ? "OK!" : "failed!") << endl;

  BIAS::SVD svd;
  BIAS::Matrix<double> Minv = svd.Invert(M);
  cout << "SVD result of inverse is M^-1 = " << CropMatrix(Minv) << endl;

  BIAS::SparseMatrix Ninv(N.GetRowsNum(), N.GetColsNum());
  N.Invert(Ninv);

  cout << "- Sparse after inversion by Gauss-Jordan is N = " << endl;
  N.Print("N", cout);
  cout << "- Inverse is N^-1 = " << endl;
  Ninv.Print("N^-1", cout);
  cout << "- Checking consistency of N... ";
  cout << (N.CheckConsistency() ? "OK!" : "failed!") << endl;
  cout << "- Checking consistency of N^-1... ";
  cout << (Ninv.CheckConsistency() ? "OK!" : "failed!") << endl;

  cout << "- M * M^-1 = " << CropMatrix(M*Minv) << endl;

  BIAS::SparseMatrix NNinv(Nold.GetRowsNum(), Ninv.GetColsNum());
  Nold.Multiply(Ninv, NNinv);
  cout << "- N * N^-1 = " << endl;
  NNinv.Print("N*N^-1", cout);
  cout << "- Checking consistency of N*Ninv... ";
  cout << (NNinv.CheckConsistency() ? "OK!" : "failed!") << endl;


  cout << endl << "-- Testing linear equation systems -- " << endl << endl;

  BIAS::SparseMatrix L(M);
  BIAS::SparseMatrix Lold(M);
  BIAS::Vector<double> b(6), x(6), b2(6), x2(6);
  b[0] =  1.0;
  b[1] =  2.0;
  b[2] =  3.0;
  b[3] = -4.0;
  b[4] =  5.0;
  b[5] = -6.0;

  cout << "- Solving linear equation system L*x = b via dense SVD and sparse "
       << "matrix Gauss-Jordan approach..." << endl;
  svd.Solve(M, b, x);
  L.Solve(b, x2);
  cout << "- Sparse LES matrix after Gauss-Jordan is L = " << endl;
  L.Print("L", cout);
  cout << "- Checking consistency of L... ";
  cout << (L.CheckConsistency() ? "OK!" : "failed!") << endl;

  cout << "- Result of SVD is x = " << x << endl;
  cout << "- Result of sparse Gauss-Jordan is x = " << x2 << endl;
  Lold.Multiply(x2, b2);
  cout << "- Evaluation of |L*x - b| yields "
       << BIAS::Vector<double>(b-b2).NormL2()
       << " (compare with dense SVD: "
       << BIAS::Vector<double>(b-(M*x)).NormL2()
       << ")" << endl << endl;

  return 0;
}