/******************************************************************************
* MODULE     : matrix_test.cc
* DESCRIPTION: Test matrices
* COPYRIGHT  : (C) 2004  Joris van der Hoeven and Gregoire Lecerf
*******************************************************************************
* This software falls under the GNU general public license and comes WITHOUT
* ANY WARRANTY WHATSOEVER. See the file $TEXMACS_PATH/LICENSE for more details.
* If you don't have this file, write to the Free Software Foundation, Inc.,
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
******************************************************************************/

#include <algebramix/algebramix-config.hpp>
#include <numerix/floating.hpp>
#include <algebramix/matrix.hpp>
#include <algebramix/matrix_int.hpp>
#include <algebramix/matrix_integer.hpp>
#include <algebramix/matrix_modular.hpp>
#include <algebramix/matrix_modular_int.hpp>
#include <algebramix/matrix_modular_integer.hpp>

using namespace mmx;
/*
#define I signed char
#define Modulus modulus<I, modulus_int_naive<3> >
#define Scalar modular<Modulus, modular_fixed<I, 7> >
*/
#define Scalar rational
//#define Scalar integer
#define V matrix_naive
#define Matrix matrix<Scalar,V>
#define max_size 7
  
Matrix
sample1 (nat r, nat c) {
  Matrix M (1, r, c);
  for (nat i = 0; i < r; i++)
    for (nat j = 0; j < c; j++)
      M(i,j)= (-i+j)/2;
  return M;
}

Matrix
sample2 (nat r, nat c) {
  Matrix M (1, r, c);
  for (nat i = 0; i < r; i++)
    for (nat j = 0; j < c; j++)
      M(i,j)= rand ();
  return M;
}

// floatings
matrix<floating<> >
sample_floating (nat r, nat c) {
  matrix<floating<> > M (1, r, c);
  for (nat i = 0; i < r; i++)
    for (nat j = 0; j < c; j++)
      M(i,j)= (i+j)/2;
  return M;
}

// integers
integer random (const integer& m) {
  static gmp_randstate_t state;
  static bool initialized = false;
  if (!initialized) {
    gmp_randinit_default (state);
    initialized = true;
  }
  integer x;
  mpz_urandomm (*x, state, *m);
  if (rand () & 1) neg (x);
  return x;
}

matrix<integer>
sample_integer (nat r, nat c, nat s) {
  matrix<integer> M (1, r, c);
  integer m= pow (integer (2), s);
  for (nat i=0; i<r; i++)
    for (nat j=0; j<c; j++) {
      M (i, j)= random (m);
      if (rand () & 1) neg (M (i, j));
    }
  return M;
}

void
quick_test () {
  Matrix M (1, 10, 10);
  mmout << "M =\n" << M << "\n";
  Matrix N (3, 10, 10);
  mmout << "N =\n" << N << "\n";
  mmout << "N + M =\n" << M+N << "\n";
  mmout << "3 N=\n" << Scalar(3) * N << "\n";
  mmout << "quo (N, 2)=\n" << quo (N, Scalar(2)) << "\n";
  N -= M;
  mmout << "N -= M -> \n" << N << "\n";
  N *= 3;
  mmout << "N *= 3 -> \n" << N << "\n";
  N += N;
  mmout << "N += N -> \n" << N << "\n";
  N /= 6;
  mmout << "N /= 6 -> \n" << N << "\n";
  N -= M;
  mmout << "(N-=M) == M -> " << (N==M) << "\n \n";
  //  mmout << "det (N)= " << det (N) << "\n";
  Scalar* o= mmx_new<Scalar> (200);
  Scalar* t= o;
  for (nat i = 0; i < 200; i++) {
    *t = Scalar (i);
    t++;
  }
  Matrix O (o, 20, 10, format<Scalar > ());
  mmout << "O =\n" << O << "\n";
  mmout << "range(O,5,2,8,10) =\n" << range (O, 5, 2, 8, 10) << "\n";
  N *= 10;
  mmout << "N =\n" << N << "\n";
  mmout << "transpose(O) =\n" << transpose (O) << "\n";
  mmout << "N*transpose(O) =\n" << N * transpose (O) << "\n";
  mmout << "O= \n" << O << "\n";
  mmout << "O...=\n" << iterate (O) << "\n";
  swap_row (M, 0, 6);
  mmout << "row(0,M) <-> row(6,M) -> \n" << M << "\n";
}

template<typename W> void
test_mul () {
  Matrix M, A, B;
  for (nat i = 1; i < max_size; i++) {
    for (nat j = 1; j < max_size; j++) {
      for (nat k = 1; k < max_size; k++) {
        mmout << "Testing product in size "
              << i << " x " << j << " x " << k << "  \r" << flush_now;
        A= sample2 (i, j);
        B= sample2 (j, k);
        M= A * B;
        if (as<matrix<Scalar, W> > (M) !=
            as<matrix<Scalar, W> > (A) * as<matrix<Scalar, W> > (B)) {
          mmout << "A = " << A << "\n";
          mmout << "B = " << B << "\n";
          mmout << "\nA * B = "
                << as<matrix<Scalar, W> > (A) * as<matrix<Scalar, W> > (B)
                << "\n";
          mmout << "instead of\n";
          mmout << A * B << "\n";
          exit (-1);
        }
      }
    }
  }
  mmout << "\n";
}
  
template<typename W> void
test_mul_integer (nat s) {
  Matrix M, A, B;
  for (nat l = 1; l < s; l++) 
    for (nat i = 1; i < max_size; i++) {
      for (nat j = 1; j < max_size; j++) {
        for (nat k = 1; k < max_size; k++) {
          mmout << "Testing product in size "
                << i << " x " << j << " x " << l << " x " << l 
                << "  \r" << flush_now;
          A= sample_integer (1 << i, 1 << j, 2 << l);
          B= sample_integer (1 << j, 1 << k, 1 << l);
          M= A * B;
        if (as<matrix<Scalar, W> > (M) !=
            as<matrix<Scalar, W> > (A) * as<matrix<Scalar, W> > (B)) {
          mmout << "A = " << A << "\n";
          mmout << "B = " << B << "\n";
          mmout << "\nA * B = "
                << as<matrix<Scalar, W> > (A) * as<matrix<Scalar, W> > (B)
                << "\n";
          mmout << "instead of\n";
          mmout << A * B << "\n";
          exit (-1);
        }
      }
    }
  }
  mmout << "\n";
}

void
test_row_echelon () {
  Matrix M, R, K;
  for (nat i = 1; i < max_size; i++) {
    for (nat j = 1; j < max_size; j++) {
      mmout << "Testing row echelon in size "
            << i << " x " << j << "  \r" << flush_now;
      M= sample1 (i, j);
      R= row_echelon (M, K);
      if (K * M != R) {
        mmout << "\nM = " << M << "\n";
        mmout << "K = " << K << "\n";
        mmout << "R = " << R << "\n";
        exit (-1);
      }
    }
  }
  mmout << "\n";
}

void
test_col_echelon () {
  Matrix M, C, K;
  for (nat i = 1; i < max_size; i++) {
    for (nat j = 1; j < max_size; j++) {
      mmout << "Testing column echelon in size "
            << i << " x " << j << "  \r" << flush_now;
      M= sample1 (i, j);
      C= column_echelon (M, K);
      if (M * K != C) {
        mmout << "\nM = " << M << "\n";
        mmout << "K = " << K << "\n";
        mmout << "C = " << C << "\n";
        exit (-1); }
    }
  }
  mmout << "\n";
}

void
test_det () {
  Matrix M;
  for (nat i = 1; i < max_size; i++) {
      mmout << "Testing matrix determinant in size "
            << i << " x " << i << "  \r" << flush_now;
    M= sample2 (i, i);
    (void) det (M);
  }
  mmout << "\n";
}

void
test_kernel () {
  Matrix M, K;
  for (nat i = 1; i < max_size; i++) {
    for (nat j = 1; j < max_size; j++) {
      mmout << "Testing matrix kernel in size "
            << i << " x " << j << "  \r" << flush_now;
      M= sample1 (i, j);
      K= kernel (M);
      if (M * K != 0) {
        mmout << "\nM = " << M << "\n";
        mmout << "kernel (M) = " << K << "\n";
        exit (-1); }
    }
  }
  mmout << "\n";
}

void
test_image () {
  Matrix M, K;
  for (nat i = 1; i < max_size; i++) {
    for (nat j = 1; j < max_size; j++) {
      mmout << "Testing matrix image in size "
            << i << " x " << j << "  \r" << flush_now;
      M= sample1 (i, j);
      K= image (M);
    }
  }
  mmout << "\n";
}

void
test_rank () {
  Matrix M, K;
  for (nat i = 1; i < max_size; i++) {
    for (nat j = 1; j < max_size; j++) {
      mmout << "Testing matrix rank in size "
            << i << " x " << j << "  \r" << flush_now;
      M= sample1 (i, j);
      (void) rank (M);
    }
  }
  mmout << "\n";
}

void
test_invert () {
  Matrix M, K;
  for (nat i = 1; i < max_size; i++) {
      M= sample2 (i, i);
      if (det (M) != 0) {
        mmout << "Testing matrix inverse in size "
              << i << " x " << i << "  \r" << flush_now;
        K= invert (M);
        if (M * K != Matrix (Scalar(1), i, i)) {
          mmout << "\nM = " << M << "\n";
          mmout << "M^(-1) = " << K << "\n";
          mmout << "M * M^(-1) = " << M * K << "\n";
          exit (-1); }
      }
  }
  mmout << "\n";
}

void
test_col_orthogonalization () {
  Matrix M, Q, L;
  for (nat i = 1; i < max_size; i++) {
    for (nat j = 1; j < max_size; j++) {
      mmout << "Testing column orthogonalization "
            << i << " x " << j << "  \r" << flush_now;
      M= sample1 (i, j);
      Q= column_orthogonalization (M, L);
      if (Q * L != M) {
        mmout << "\nM = " << M << "\n";
        mmout << "Q = " << Q << "\n";
        mmout << "L = " << L << "\n";
        exit (-1); }
    }
  }
  mmout << "\n";
}

void
test_col_orthonormalization () {
  matrix<floating<> > M, Q, L;
  for (nat i = 1; i < max_size; i++) {
    for (nat j = 1; j < max_size; j++) {
      mmout << "Testing column orthonormalization "
            << i << " x " << j << "  \r" << flush_now;
      M= sample_floating (i, j);
      Q= column_orthonormalization (M, L);
      if (big_max (abs (Q * L - M)) > floating<> (0.000001)) {
        mmout << "\nM = " << M << "\n";
        mmout << "Q = " << Q << "\n";
        mmout << "L = " << L << "\n";
        exit (-1); }
    }
  }
  mmout << "\n";
}

int
main () {
  try {
    //quick_test ();
    //test_mul<matrix_modular_int> ();
    //test_mul<matrix_integer> ();
    //test_mul_integer<matrix_integer> (10);
    test_row_echelon ();
    test_col_echelon ();
    test_det ();
    test_kernel ();
    test_image ();
    test_rank ();
    test_invert ();
    test_col_orthogonalization ();
    test_col_orthonormalization ();
  }
  catch (exception& e) {
    mmout << e << "\n";
    return -1;
  }
  return 0;
}