TestGaussJordan.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 */

/** @file TestGaussJordan.cpp
 *  @ingroup g_tests 
    @brief Test GaussJordanAlgorithm computing the reduced row echelon form
    @relates Matrix
    @author woelk 05/2008 (c) www. vision-n.de */

#include <Base/Math/Matrix.hh>
#include <Base/Math/Random.hh>
#include <Base/Common/CompareFloatingPoint.hh>

using namespace BIAS;
using namespace std;

Random myrand;

void GenerateRandomMatrix(Matrix<double>& A);

bool TestForReducedRowEchelonForm(Matrix<double>& A);


int main(int argc, char *argv[])
{
  const bool verbose = false;
  // check num_sizes different matric sizes
  const int num_sizes = 100;
  for (int t=0; t<num_sizes; t++){
    const int numr=myrand.GetUniformDistributedInt(1,10), 
      numc=myrand.GetUniformDistributedInt(1,10);
    Matrix<double> A(numr, numc), Aorig;


    // check differently filled matrices
    const int num_tries = 1000;
    for (int i=0; i<num_tries; i++){
      GenerateRandomMatrix(A);
      Aorig = A;
      //if (verbose) cout << "------------------------------------\nA: "<<A;
      A.GaussJordan();
      if (verbose) cout <<"\nrref(A): "<<A;
      if (!TestForReducedRowEchelonForm(A)){
        if (!verbose) 
          BIASERR("error: matrix is not in reduced row echelon form: "<<A);
        //if (verbose) 
          cout << "\nA: "<<Aorig;
          //if (verbose) 
          cout << "  !!FAILED!!\n";
        return -1;
      }
      if (verbose) cout << "passed!\n";
    }
  }

  return 0;
}

void GenerateRandomMatrix(Matrix<double>& A)
{
  static int num_calls = 0;
  const double range=1.;
  const int numr=A.num_rows(), numc=A.num_cols();
  for (int r=0; r<numr; r++){
    for (int c=0; c<numc; c++){
      A[r][c] = myrand.GetUniformDistributed(-range, range);
      // insert some zeroes in every second generation
      //if (abs(A[r][c])<0.5) A[r][c] = 0.;
      if ((num_calls%2)==0 && abs(A[r][c])<0.5) A[r][c] = 0.;
    }
  }
  num_calls++;
}

bool TestForReducedRowEchelonForm(Matrix<double>& A)
{
  const int numr=A.num_rows(), numc=A.num_cols();
  int lead=-1, r, c;
  for (r=0; r<numr; r++){
    // look for the first non-zero element in the row
    for (c=0; c<numc && Equal(A[r][c], 0.0); c++){ }
    // check if it is right of the leading elemnt from the previous row 
    if (c<=lead && lead<numc){ cout << "lead pos\n"; return false; }
    // remember the leading element
    lead = c;
    if (lead==numc) continue;
    // check if the leading element is one
    if (!Equal(A[r][lead], 1.0)) { cout << "lead == 1.\n"; return false; }
    // check if the leading element id the only nonzero element in its column
    for (int rr=0; rr<numr; rr++){
      if (rr==r) continue;
      if (!Equal(A[rr][lead], 0.0)) { cout << "top == 0.\n"; return false; }
    }
  }
  return true;
}