quiksolv.imp

Gambit 是一个游戏库理论软件

//
// $Source: /home/gambit/CVS/gambit/sources/poly/quiksolv.imp,v $
// $Date: 2002/08/27 17:29:49 $
// $Revision: 1.2 $
//
// DESCRIPTION:
// Implementation of quick-solver classes
//
// 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 "quiksolv.h"
#include "math/gmath.h"

//---------------------------------------------------------------
//                      class: QuikSolv
//---------------------------------------------------------------

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

template <class T> QuikSolv<T>::QuikSolv(const gPolyList<T>& given,
					 gStatus &p_status) 
  : System(given), 
    gDoubleSystem(given.AmbientSpace(),given.TermOrder(),
		  given.NormalizedList()),
    NoEquations( gmin(System.Dmnsn(),System.Length()) ),
    NoInequalities( gmax(System.Length() - System.Dmnsn(),0) ),
    TreesOfPartials(gDoubleSystem),
    HasBeenSolved(false), 
    Roots(),
    isMultiaffine(System.IsMultiaffine()),
    Equation_i_uses_var_j(Eq_i_Uses_j()),
    m_status(p_status)
{ }

template <class T> QuikSolv<T>::QuikSolv(const gPolyList<T>& given,
					 const int&          no_eqs,
					 gStatus &p_status) 
  : System(given), 
    gDoubleSystem(given.AmbientSpace(),given.TermOrder(),
		  given.NormalizedList()),
    NoEquations(no_eqs),
    NoInequalities(System.Length() - no_eqs),
    TreesOfPartials(gDoubleSystem),
    HasBeenSolved(false), 
    Roots(),
    isMultiaffine(System.IsMultiaffine()),
    Equation_i_uses_var_j(Eq_i_Uses_j()),
    m_status(p_status)
{ }

template<class T> QuikSolv<T>::QuikSolv(const QuikSolv& qs)
  : System(qs.System), 
    gDoubleSystem(qs.gDoubleSystem),
    NoEquations(qs.NoEquations),
    NoInequalities(qs.NoEquations),
    TreesOfPartials(qs.TreesOfPartials), 
    HasBeenSolved(qs.HasBeenSolved), 
    Roots(qs.Roots),
    isMultiaffine(qs.isMultiaffine),
    Equation_i_uses_var_j(qs.Equation_i_uses_var_j),
    m_status(qs.m_status)
  
{ }

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

//-------------------------------------------------------
//      Supporting Calculations for the Constructors
//-------------------------------------------------------

template <class T> gRectArray<bool> QuikSolv<T>::Eq_i_Uses_j() const
{
  gRectArray<bool> answer(System.Length(),Dmnsn());
  for (int i = 1; i <= System.Length(); i++)
    for (int j = 1; j <= Dmnsn(); j++)
      if (System[i].DegreeOfVar(j) > 0)
	answer(i,j) = true;
      else
	answer(i,j) = false;
  return answer;
}


//---------------------------
// Find Roots Using Pelican
//---------------------------

template<class T> 
bool QuikSolv<T>::AllRealRootsFromPelican(const gPolyList<gDouble> & sys, 
					  gList<gVector<gDouble> > &ans) const
{
  try {
  PelView pel(sys);

  if (!pel.FoundAllRoots())
    return false;
  ans = pel.RealRoots();
  }
  catch (ErrorInPelican) {
    return false;
  }
  catch (ErrorInQhull) {
    return false;
  }

  return true;
}

template<class T> bool
QuikSolv<T>::PelicanRoots(const gRectangle<T>& r, 
			        gList<gVector<gDouble> > &answer) const
{
  gList<gVector<gDouble> > firstcut;
  bool pelworked = 
    AllRealRootsFromPelican(gDoubleSystem.InteriorSegment(1,Dmnsn()), 
			    firstcut);
  if (!pelworked)
    return false;

  for (int i = 1; i <= firstcut.Length(); i++) {
    firstcut[i] = NewtonPolishOnce(firstcut[i]);
    firstcut[i] = NewtonPolishOnce(firstcut[i]);
    bool isokay(true);
    if (!TogDouble(r).Contains(firstcut[i]))
      isokay = false;
    for (int j = Dmnsn() + 1; isokay && j <= System.Length(); j++)
      if (gDoubleSystem[j].Evaluate(firstcut[i]) < (gDouble)0)
	isokay = false;
    if (isokay)
      answer += firstcut[i];
  }
  
  return true;
}

//---------------------------
// Is a root impossible?
//---------------------------

template<class T> 
bool QuikSolv<T>::SystemHasNoRootsIn(const gRectangle<gDouble>& r,
				     gArray<int>& precedence)        const
{
  for (int i = 1; i <= System.Length(); i++) {

    if ( (precedence[i] <= NoEquations && 
	              TreesOfPartials[precedence[i]].PolyHasNoRootsIn(r)) ||
         (precedence[i] >  NoEquations && 
	   TreesOfPartials[precedence[i]].PolyEverywhereNegativeIn(r)) ) {
      if (i != 1) { // We have found a new "most likely to never vanish"
	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;
}

//--------------------------------------
// Does Newton's method lead to a root?
//--------------------------------------

template <class T>
bool QuikSolv<T>::NewtonRootInRectangle(const gRectangle<gDouble>& r,
					      gVector<gDouble>& point) const
{
  assert (NoEquations == System.Dmnsn());

  gVector<gDouble> zero(Dmnsn());
  for (int i = 1; i <= Dmnsn(); i++) zero[i] = (gDouble)0;

  gVector<gDouble> oldevals = TreesOfPartials.ValuesOfRootPolys(point,
								NoEquations);

  if ( oldevals == zero ) 
    if (r.Contains(point)) return true;
    else                   return false;

  gRectangle<gDouble> bigr = r.SameCenterDoubleSideLengths();

  gVector<gDouble> newpoint(Dmnsn());

  while (true) {
    try {
      newpoint = NewtonPolishOnce(point);
    }
    catch (gSquareMatrix<gDouble>::MatrixSingular) {
      bool nonsingular = false;
      int direction = 1;
      while (direction < Dmnsn() && !nonsingular) {
	gVector<gDouble> perturbed_point(point);
	if (r.UpperBoundOfCoord(direction) > point[direction])
	  perturbed_point[direction] += 
	    (r.UpperBoundOfCoord(direction) - point[direction])/10;
	else
	  perturbed_point[direction] += 
	    (r.LowerBoundOfCoord(direction) - point[direction])/10;
	nonsingular = true;

	try {
	  newpoint = point + 
	    (NewtonPolishOnce(perturbed_point) - perturbed_point);
	}
	catch (gSquareMatrix<gDouble>::MatrixSingular) {
	  nonsingular = false;
	}
	direction++;
      }
      if (!nonsingular) {
	gVector<gDouble> perturbed_point(point);
	if (r.UpperBoundOfCoord(direction) > point[direction])
	  perturbed_point[direction] += 
	    (r.UpperBoundOfCoord(direction) - point[direction])/10;
	else
	  perturbed_point[direction] += 
	    (r.LowerBoundOfCoord(direction) - point[direction])/10;
	newpoint = point + 
	  (NewtonPolishOnce(perturbed_point) - perturbed_point);
      }
    }

    if ( !bigr.Contains(newpoint) ) return false;
    point = newpoint;
    
    gVector<gDouble> newevals = TreesOfPartials.ValuesOfRootPolys(point,
								  NoEquations);

    if (newevals * newevals > oldevals * oldevals) return false;
    if ( newevals == zero )
      if (r.Contains(point)) {
	point = SlowNewtonPolishOnce(point);
	point = SlowNewtonPolishOnce(point);
	return true;
      }
      else                   return false;                

    oldevals = newevals;
  }
}

template <class T>
bool QuikSolv<T>::NewtonRootNearRectangle(const gRectangle<gDouble>& r,
				    	        gVector<gDouble>& point) const
{
  gVector<gDouble> zero(NoEquations);
  for (int i = 1; i <= NoEquations; i++) zero[i] = (gDouble)0;

  gVector<gDouble> oldevals = TreesOfPartials.ValuesOfRootPolys(point,
								NoEquations);

  gRectangle<gDouble> bigr = r.CubeContainingCrcmscrbngSphere();

  if ( oldevals == zero ) 
    if (bigr.Contains(point)) return true;
    else                      return false;

  while (true) {
    gVector<gDouble> newpoint = SlowNewtonPolishOnce(point);
    if ( !bigr.Contains(newpoint) ) 
      return false;
    point = newpoint;
    gVector<gDouble> newevals = TreesOfPartials.ValuesOfRootPolys(point,
								  NoEquations);
    if ( newevals == zero )
      if (r.Contains(point)) return true;
      else                   return false;                

    oldevals = newevals;
  }
}

//------------------------------------
// Is the Newton root the only root?
//------------------------------------

template<class T>
   gDouble QuikSolv<T>::MaxDistanceFromPointToVertexAfterTransformation(
				      const gRectangle<gDouble>& r,
				      const gVector<gDouble>& p,
				      const gSquareMatrix<gDouble>& M)    const
{
  assert( r.Contains(p) );

  gDouble max = (gDouble)0;

  gArray<int> ListOfTwos(Dmnsn());
  for (int i = 1; i <= Dmnsn(); i++) ListOfTwos[i] = 2;
  gIndexOdometer ListOfTopBottoms(ListOfTwos);

  while (ListOfTopBottoms.Turn()) {
    gVector<gDouble> diffs(Dmnsn());
    for (int i = 1; i <= Dmnsn(); i++)
      if (ListOfTopBottoms[i] == 2)
	diffs[i] = r.CartesianFactor(i).UpperBound() - p[i];
      else 
	diffs[i] = p[i] - r.CartesianFactor(i).LowerBound();
    gVector<gDouble> new_diffs = M * diffs;
    gDouble squared_length = new_diffs * new_diffs;
    if (max < squared_length)
       max = squared_length;
  }

  return sqrt((gDouble)max);
}

template<class T>
bool QuikSolv<T>::HasNoOtherRootsIn(const gRectangle<gDouble>& r,
			  const gVector<gDouble>& p,
			  const gSquareMatrix<gDouble>& M)  const
{
  assert (NoEquations == System.Dmnsn());

  gPolyList<gDouble> system1 = gDoubleSystem.TranslateOfSystem(p);

  gPolyList<gDouble> system2 = system1.SystemInNewCoordinates(M);

  gDouble radius = MaxDistanceFromPointToVertexAfterTransformation(r,p,M);

  gDouble max = (gDouble)0;
  for (int i = 1; i <= Dmnsn(); i++) 
    max += system2[i].MaximalValueOfNonlinearPart(radius);

  if (max >= radius) return false;
  else               return true;
}


//--------------------------------------------
// Does Newton's method yield a unique root?
//--------------------------------------------