//
// $Source: /cvsroot/gambit/gambit/sources/tools/enumpoly/linrcomb.imp,v $
// $Date: 2006/01/07 05:21:41 $
// $Revision: 1.4 $
//
// DESCRIPTION:
// Implementation of linear combination class
//
// This file is part of Gambit
// Copyright (c) 2002, The Gambit Project
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program 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 General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
//

#include "libgambit/libgambit.h"
#include "linrcomb.h"

//-------------------------------------------------------------------------
// LinearCombination<T>: Constructors, destructors, constructive operators
//-------------------------------------------------------------------------

template <class T> 
LinearCombination<T>::LinearCombination(const Gambit::Matrix<T> &M)
: scrambled(M), weights(M.NumRows()), last_row_is_spanned(true)
{ 
  int K = scrambled.NumRows() - 1;

  int i,j;
  Gambit::Matrix<T> B(K,K);           // Initialize the matrix recording the
  for (i = 1; i <= K; i++) // effect of all row operations.
    for (j = 1; j <= K; j++)
      if (i == j) B(i,j) = 1;
      else        B(i,j) = 0;

  Gambit::Vector<int> PivotRow(K);
  for (i = 1; i <= K; i++) PivotRow[i] = i;
  Gambit::Vector<int> PivotColumn(K);

  int c = 1;
  for (i = 1; i <= K-1; i++) { // Each iteration is a
                                   // Gaussian column clearing

    int i1 = i;        // Find pivot row and column
    while (scrambled(PivotRow[i1],c) == (T)0 && c <= M.NumColumns())
      if (i1 < K) i1++;
      else { c++; i1 = i; }
    if (c > M.NumColumns()) {
      //      gout << "First NumRows() rows not linearly independent in "
      //	<< "LinearCombination(Gambit::Matrix<T> M).\n";
      exit(1);
    }

    int tmp = PivotRow[i];       // Swap
    PivotRow[i] = PivotRow[i1];
    PivotRow[i1] = tmp;

    PivotColumn[i] = c;

    for (int i2 = 1; i2 <= K; i2++) { // Each i2 is a pair of row operations.
      if (i2 != i) {
	T factor = scrambled(PivotRow[i2],PivotColumn[i])
	                             /
		   scrambled(PivotRow[i],PivotColumn[i]);
	AddMultipleOfRowiToRowj(PivotRow[i],PivotRow[i2],-factor,scrambled);
	AddMultipleOfRowiToRowj(PivotRow[i],PivotRow[i2],-factor,B);
      }
    }

    c++;
  }

  if (K > 0) {
    if (K == 1) PivotColumn[K] = 1;
    else        PivotColumn[K] = PivotColumn[K-1]+1;
    while ( PivotColumn[K] <= M.NumColumns() &&
	   scrambled(PivotRow[K],PivotColumn[K]) == (T)0 ) PivotColumn[K]++;
    if (PivotColumn[K] > M.NumColumns()) {
      //      gout << "First NumRows()-1 rows not linearly independent in "
      //	<< "LinearCombination(Gambit::Matrix<T> M).\n";
      // gout << "At the time of failure, scrambled is \n"
      //	<< scrambled << "\n";
      exit(1);
    }

    for (int i2 = 1; i2 < K; i2++) { // Each i2 is a pair of row operations.
      T factor = scrambled(PivotRow[i2],PivotColumn[K])
                                   /
       	         scrambled(PivotRow[K],PivotColumn[K]);
      AddMultipleOfRowiToRowj(PivotRow[K],PivotRow[i2],-factor);
      AddMultipleOfRowiToRowj(PivotRow[K],PivotRow[i2],-factor,B);
    }
  }

  Gambit::Vector<T> xBinverse(K);
  for (int i3 = 1; i3 <= K; i3++) { // Each i3 is a row operation on last row.
    T factor = scrambled(K+1,PivotColumn[i3])
                            /
           scrambled(PivotRow[i3],PivotColumn[i3]);
    xBinverse[PivotRow[i3]] = factor;
    AddMultipleOfRowiToRowj(PivotRow[i3],K+1,-factor);
  }
  Gambit::Vector<T> x(xBinverse * B);

  for (j = 1; j <= M.NumColumns(); j++)
    if (scrambled(K+1,j) != (T)0) last_row_is_spanned = false;

  if (last_row_is_spanned) {
    for (int i4 = 1; i4 <= K; i4++) weights[i4] = -x[i4];
    weights[K+1] = 1;
  }
}

template <class T> 
LinearCombination<T>::LinearCombination(const LinearCombination<T> &L)
: scrambled(L.scrambled), 
  weights(L.weights),
  last_row_is_spanned(L.last_row_is_spanned)
{ }
 
template <class T> LinearCombination<T>::~LinearCombination()
{ }

//-------------------------------------------------------------------------
//                    LinearCombination<T>: Comparison operators
//-------------------------------------------------------------------------

template <class T> 
bool LinearCombination<T>::operator==(const LinearCombination<T> &L) const
{
  if (scrambled == L.scrambled && 
      last_row_is_spanned == L.last_row_is_spanned &&
      weights == L.weights)
    return true;
  else return false;
}

template <class T> 
bool LinearCombination<T>::operator!=(const LinearCombination<T> &L) const
{ return !(*this == L); }

//-------------------------------------------------------------------------
//                    LinearCombination<T>: Information
//-------------------------------------------------------------------------

template <class T> bool LinearCombination<T>::LastRowIsSpanned() const
{
  return last_row_is_spanned;
}

template <class T> Gambit::Vector<T> LinearCombination<T>::LinearDependence() const
{
  return weights;
}

//-------------------------------------------------------------------------
//                  LinearCombination<T>: Private Members
//-------------------------------------------------------------------------

template <class T> 
void LinearCombination<T>::AddMultipleOfRowiToRowj(const int& i, 
						   const int& j,
						   const T& scalar,
						   Gambit::Matrix<T>& B)
{
  for (int k = 1; k <= B.NumColumns(); k++)
    B(j,k) += scalar * B(i,k);
}

template <class T> 
void LinearCombination<T>::AddMultipleOfRowiToRowj(const int& i, 
						   const int& j,
						   const T& scalar)
{
  AddMultipleOfRowiToRowj(i,j,scalar,scrambled);
}