//
// $Source: /cvsroot/gambit/gambit/sources/tools/enumpoly/ineqsolv.imp,v $
// $Date: 2006/01/07 05:41:26 $
// $Revision: 1.4 $
//
// DESCRIPTION:
// Implementation of IneqSolv
//
// 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 "ineqsolv.h"

//---------------------------------------------------------------
//                      class: IneqSolv
//---------------------------------------------------------------

//---------------------------
// Constructors / Destructors
//---------------------------

template <class T> IneqSolv<T>::IneqSolv(const gPolyList<T>& given)
  : System(given), 
    TreesOfPartials(given),
    Epsilon((T)1/(T)1000000)
    //    Epsilon((T)0),
    //    HasBeenSolved(false), 
    //    HasASolution(triUNKNOWN), 
    //    Sample(given.Dmnsn()),

{ }

template<class T> IneqSolv<T>::IneqSolv(const IneqSolv& qs)
  : System(qs.System), 
    TreesOfPartials(qs.TreesOfPartials), 
    Epsilon(qs.Epsilon)
    //    HasBeenSolved(qs.HasBeenSolved), 
    //    HasASolution(qs.HasASolution), 
    //    Sample(qs.Sample),
{ }

template<class T> IneqSolv<T>::~IneqSolv()
{ }

//----------------------------------
//        Operators
//----------------------------------

template<class T> 
IneqSolv<T>& IneqSolv<T>::operator=(const IneqSolv<T> & rhs)
{
  assert (System == rhs.System);

  if (*this != rhs) {
    Epsilon = rhs.Epsilon;
  }
  return *this;
}

template<class T>  
bool IneqSolv<T>::operator==(const IneqSolv<T> & rhs) const
{
    if (System  != rhs.System || Epsilon != rhs.Epsilon)
         return false;
    else return true;
}

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


//---------------------------
//      Calculations
//---------------------------

template<class T> 
bool IneqSolv<T>::IsASolution(const Gambit::Vector<T>& v) const
{
  bool answer(true);
  for (int i = 1; i <= System.Length() && answer; i++)
    if ( System[i].Evaluate(v) < -Epsilon )
      answer = false;
  return answer;
}

template<class T> 
bool IneqSolv<T>::SystemHasNoSolutionIn(const gRectangle<T>& r,
					      Gambit::Array<int>& precedence) const
{
  for (int i = 1; i <= System.Length(); i++) {

    if ( TreesOfPartials[precedence[i]].PolyEverywhereNegativeIn(r) ) {

      /*
      //DEBUG
      gout << "The polynomial " << System[precedence[i]] 
	   << " has no roots in ";
      gRectangle<T> newrect(r);
      gout << newrect;
      gout << ".\n";
      */

      if (i != 1) { // We have found a new "most likely to never be positive"
	int tmp = precedence[i];
	for (int j = 1; j <= i-1; j++) 
	  precedence[i - j + 1] = precedence[i - j];
	precedence[1] = tmp;
      }
      return true; 
    }
  }
  return false;
}

template <class T> const bool
IneqSolv<T>::ASolutionExistsRecursion(const gRectangle<T>& r, 
				            Gambit::Vector<T>& sample,
				            Gambit::Array<int>& precedence) const
{
  /*
  //DEBUG
  gout << "The rectangle is\n";
  for (int i = 1; i <= r.Dmnsn() - 1; i++)
    gout << r.CartesianFactor(i) << "x";
  if (r.Dmnsn() > 0) 
    gout << r.CartesianFactor(r.Dmnsn()) << "\n";
  */

  // Check for user interrupt
  if ( IsASolution(r.Center()) )
    return true;

  if ( SystemHasNoSolutionIn(r, precedence) ) 
    return false;

  int N = r.NumberOfCellsInSubdivision();
  for (int i = 1; i <= N; i++) 
    if (ASolutionExistsRecursion(r.SubdivisionCell(i),
				 sample,
				 precedence))
      return true;
  return false; 
}

template <class T> const bool
IneqSolv<T>::ASolutionExists (const gRectangle<T>& r, Gambit::Vector<T>& sample)
{
     // precedence orders search for everywhere negative poly
  Gambit::Array<int> precedence(System.Length());  
  for (int i = 1; i <= System.Length(); i++) precedence[i] = i;

  /*
  //DEBUG
  gout << "The system is\n" << System << "\n";
  */

  bool answer = ASolutionExistsRecursion(r, sample, precedence);

  return answer;
}