gpartltr.imp

Gambit 是一个游戏库理论软件

//
// $Source: /home/gambit/CVS/gambit/sources/poly/gpartltr.imp,v $
// $Date: 2002/08/27 17:29:46 $
// $Revision: 1.2 $
//
// DESCRIPTION:
// Implementation of TreeOfPartials
//
// 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 "gpartltr.h"

//---------------------------------------------------------------
//                   class: TreeOfPartials
//---------------------------------------------------------------

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

template <class T> TreeOfPartials<T>::TreeOfPartials(const gPoly<T>& given) 
: PartialTree(given)
{
  TreeOfPartialsRECURSIVE(PartialTree, PartialTree.RootNode());  
}

//-------------------------------------------------------------------------
// Recursive generation of all partial derivatives of the root polynomial
//-------------------------------------------------------------------------

template <class T>
void TreeOfPartials<T>::TreeOfPartialsRECURSIVE(gTree<gPoly<T> >& t,
					  gTreeNode<gPoly<T> >* n) const
{
  if (n->GetData().Degree() >= 1) {
    for (int i = 1; i <= n->GetData().Dmnsn(); i++) {
      t.InsertAt(n->GetData().PartialDerivative(i),n);
      TreeOfPartialsRECURSIVE(t,n->GetYoungest());
    }
  }
}

template<class T> TreeOfPartials<T>::TreeOfPartials(const TreeOfPartials& qs)
: PartialTree(qs.PartialTree)
{
}

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

template <class T> 
T TreeOfPartials<T>::ValueOfPartialOfRootPoly(const int& coord, 
					      const gVector<T>& p) const
{
  if (RootPoly().Degree() <= 0) return (T)0;
  else {
    int i = 1; 
    gTreeNode<gPoly<T> >* node = RootNode()->GetEldest();
    while (i < coord) {
      i++;
      node = node->GetNext();
    }
    T answer = node->GetData().Evaluate(p);
    return answer;
  }
}

template <class T>
gVector<T> TreeOfPartials<T>::VectorOfPartials(const gVector<T>& point) const
{
  gVector<T> answer(Dmnsn());
  for (int i = 1; i <= Dmnsn(); i++)
    answer[i] = ValueOfPartialOfRootPoly(i,point);
  return answer;
}

template <class T> T
TreeOfPartials<T>::MaximalNonconstantContributionRECURSIVE(
					     const gTreeNode<gPoly<T> >* n,
					     const gVector<T>& p,
					     const gVector<T>& halvesoflengths,
					     gVector<int>& wrtos) const
{
  T answer = (T)0;

  if (n->GetEldest() != NULL) {
    gList<gTreeNode<gPoly<T> >*> children = PartialTree.Children(n);
    for (int i = 1; i <= children.Length(); i++) {
      wrtos[i]++;

      T increment = children[i]->GetData().Evaluate(p);
      if (increment < (T)0) increment = -increment;
     	  
      for (int j = 1; j <= p.Length(); j++) 
	for (int k = 1; k <= wrtos[j]; k++) {
	  increment *= halvesoflengths[j];
	  increment /= (T)k;
	}

      answer += increment;

      answer += MaximalNonconstantContributionRECURSIVE(children[i],
							p,
							halvesoflengths,
							wrtos);
      
      wrtos[i]--;
    }
  }

  return answer;
}

template <class T> T
TreeOfPartials<T>::MaximalNonconstantDifferenceRECURSIVE(
					     const gTreeNode<gPoly<T> >* n1,
					     const gTreeNode<gPoly<T> >* n2,
					     const gVector<T>& p,
					     const gVector<T>& halvesoflengths,
					     gVector<int>& wrtos) const
{
  T answer = (T)0;

  if (n1->GetEldest() != NULL && n2->GetEldest() != NULL) {
    gList<gTreeNode<gPoly<T> >*> children1 = PartialTree.Children(n1);
    gList<gTreeNode<gPoly<T> >*> children2 = PartialTree.Children(n2);
    for (int i = 1; i <= children1.Length(); i++) {
      wrtos[i]++;

      T increment = children1[i]->GetData().Evaluate(p) - 
	            children2[i]->GetData().Evaluate(p);
      if (increment < (T)0) increment = -increment;
     	  
      for (int j = 1; j <= p.Length(); j++) 
	for (int k = 1; k <= wrtos[j]; k++) {
	  increment *= halvesoflengths[j];
	  increment /= (T)k;
	}

      answer += increment;

      answer += MaximalNonconstantDifferenceRECURSIVE(children1[i],
						      children2[i],
						      p,
						      halvesoflengths,
						      wrtos);
      
      wrtos[i]--;
    }
  }

  else if (n1->GetEldest() != NULL && n2->GetEldest() == NULL) {
    gList<gTreeNode<gPoly<T> >*> children1 = PartialTree.Children(n1);
    for (int i = 1; i <= children1.Length(); i++) {
      wrtos[i]++;

      T increment = children1[i]->GetData().Evaluate(p);
      if (increment < (T)0) increment = -increment;
     	  
      for (int j = 1; j <= p.Length(); j++) 
	for (int k = 1; k <= wrtos[j]; k++) {
	  increment *= halvesoflengths[j];
	  increment /= (T)k;
	}

      answer += increment;

      answer += MaximalNonconstantContributionRECURSIVE(children1[i],
							p,
							halvesoflengths,
							wrtos);
      
      wrtos[i]--;
    }
  }

  else if (n1->GetEldest() == NULL && n2->GetEldest() != NULL) {
    gList<gTreeNode<gPoly<T> >*> children2 = PartialTree.Children(n2);
    for (int i = 1; i <= children2.Length(); i++) {
      wrtos[i]++;

      T increment = children2[i]->GetData().Evaluate(p);
      if (increment < (T)0) increment = -increment;
     	  
      for (int j = 1; j <= p.Length(); j++) 
	for (int k = 1; k <= wrtos[j]; k++) {
	  increment *= halvesoflengths[j];
	  increment /= (T)k;
	}

      answer += increment;

      answer += MaximalNonconstantContributionRECURSIVE(children2[i],
							p,
							halvesoflengths,
							wrtos);
      
      wrtos[i]--;
    }
  }

  return answer;
}


template <class T> T
TreeOfPartials<T>::MaximalNonconstantContribution(const gVector<T>& p,
						  const gVector<T>& 
						     halvesoflengths) const
{
  gVector<int> WithRespectTos(p.Length());
  for (int i = 1; i <= p.Length(); i++) WithRespectTos[i] = 0;

  return MaximalNonconstantContributionRECURSIVE(RootNode(),
						 p,
						 halvesoflengths,
						 WithRespectTos); 
}

template <class T> T
TreeOfPartials<T>::MaximalNonconstantDifference(const TreeOfPartials<T>& other,
						const gVector<T>& p,
						const gVector<T>& 
						   halvesoflengths) const
{
  gVector<int> WithRespectTos(p.Length());
  for (int i = 1; i <= p.Length(); i++) WithRespectTos[i] = 0;

  return MaximalNonconstantDifferenceRECURSIVE(other.RootNode(),
					       RootNode(),
					       p,
					       halvesoflengths,
					       WithRespectTos); 
}


template <class T>
T TreeOfPartials<T>::EvaluateRootPoly(const gVector<T>& point) const 
{
  return RootNode()->GetData().Evaluate(point); 
}

template <class T>
T TreeOfPartials<T>::ValueOfRootPoly(const gVector<T>& point) const 
{ 
  return RootPoly().Evaluate(point); 
}


template <class T> gOutput& operator << (gOutput& output, 
					 const TreeOfPartials<T>& /* x */)
{
  // Changed to just return silently.
  return output;
}

template <class T> bool
TreeOfPartials<T>::PolyHasNoRootsIn(const gRectangle<T>& r) const
{
  if (this->RootNode()->GetData().IsMultiaffine())
    return MultiaffinePolyHasNoRootsIn(r);
  else {
    gVector<T> center = r.Center();
    
    T constant = this->RootNode()->GetData().Evaluate(center);
    if (constant < (T)0) constant = - constant;
    
    gVector<T> HalvesOfSideLengths = r.SideLengths();
    for (int k = 1; k <= Dmnsn(); k++) HalvesOfSideLengths[k] /= (T)2;
    
    T max = 
      this->MaximalNonconstantContribution(center, HalvesOfSideLengths);
  
    if (max >= constant) 
      return false;
    else 
      return true;
  }
}

template <class T> bool 
TreeOfPartials<T>::MultiaffinePolyHasNoRootsIn(const gRectangle<T>& r) const
{
  int sign;
  if (this->RootNode()->GetData().Evaluate(r.Center()) > (T)0)
    sign = 1;
  else
    sign = -1;
  
  gArray<int> zeros(Dmnsn());
  gArray<int> ones(Dmnsn());
  for (int j = 1; j <= Dmnsn(); j++) {
    zeros[j] = 0; 
    if (this->RootNode()->GetData().DegreeOfVar(j) > 0) ones[j] = 1;
    //Equation_i_uses_var_j(index,j)) ones[j] = 1; 
    else                                ones[j] = 0;
  }
  gIndexOdometer topbottoms(zeros,ones);
  
  while (topbottoms.Turn()) {
    gVector<T> point(Dmnsn());
    for (int i = 1; i <= Dmnsn(); i++)
      if (topbottoms[i] == 0)
	point[i] = r.LowerBoundOfCoord(i);
      else
	point[i] = r.UpperBoundOfCoord(i);
    if ((T)sign * this->RootNode()->GetData().Evaluate(point) <=  (T)0)
      return false;
  }
  
  return true;
}



template<class T> bool
TreeOfPartials<T>::MultiaffinePolyEverywhereNegativeIn(
					        const gRectangle<T>& r) const
{ 
  if (Dmnsn() == 0) {
    gVector<T> point(Dmnsn());
    if (this->RootNode()->GetData().Evaluate(point) >= (T)0)
      return false;
    else
      return true;
  }

  gArray<int> zeros(Dmnsn());
  gArray<int> ones(Dmnsn());
  for (int j = 1; j <= Dmnsn(); j++) {
    zeros[j] = 0; 
    if (this->RootNode()->GetData().DegreeOfVar(j) > 0) ones[j] = 1; 
    else                                                ones[j] = 0;
  }
  gIndexOdometer topbottoms(zeros,ones);
  
  while (topbottoms.Turn()) {
    gVector<T> point(Dmnsn());
    for (int i = 1; i <= Dmnsn(); i++)
      if (topbottoms[i] == 0)
	point[i] = r.LowerBoundOfCoord(i);
      else
	point[i] = r.UpperBoundOfCoord(i);
    if (this->RootNode()->GetData().Evaluate(point) >= (T)0)
      return false;
  }
  
  return true;
}

template<class T> bool
TreeOfPartials <T>::PolyEverywhereNegativeIn(const gRectangle<T>& r) const
{
  if (this->RootNode()->GetData().IsMultiaffine())
    return MultiaffinePolyEverywhereNegativeIn(r);
  else {
    gVector<T> center = r.Center();

    T constant = 
      this->RootNode()->GetData().Evaluate(center);
    if (constant >= (T)0) return false;
    
    gVector<T> HalvesOfSideLengths = r.SideLengths();
    for (int k = 1; k <= Dmnsn(); k++) HalvesOfSideLengths[k] /= (T)2;
    
    T max = this->MaximalNonconstantContribution(center, HalvesOfSideLengths);

    if (max + constant >= (T)0) return false;
    else                     return true;
  }
}




//---------------------------------------------------------------
//                    class: ListOfPartialTrees
//---------------------------------------------------------------

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

template <class T> 
ListOfPartialTrees<T>::ListOfPartialTrees(const gList<gPoly<T> >& given) 
: PartialTreeList()
{
  for (int i = 1; i <= given.Length(); i++) 
    PartialTreeList += TreeOfPartials<T>(given[i]);
}

template <class T> 
ListOfPartialTrees<T>::ListOfPartialTrees(const gPolyList<T>& given) 
: PartialTreeList()
{
  for (int i = 1; i <= given.Length(); i++) 
    PartialTreeList += TreeOfPartials<T>(given[i]);
}

template<class T> 
ListOfPartialTrees<T>::ListOfPartialTrees(const ListOfPartialTrees& qs)
: PartialTreeList(qs.PartialTreeList)
{
}

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

template <class T> bool 
ListOfPartialTrees<T>::operator ==(const ListOfPartialTrees<T>& rhs) const
{
  if (Length() != rhs.Length()) return false;
  for (int i = 1; i <= Length(); i++)
    if ((*this)[i] != rhs[i])
      return false;
  return true;
}

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

template <class T> 
gMatrix<T> ListOfPartialTrees<T>::DerivativeMatrix(const gVector<T>& p) const
{
  gMatrix<T> answer(Length(),Dmnsn());
  int i;
  for (i = 1; i <= Length(); i++)
    for (int j = 1; j <= Dmnsn(); j++)
      answer(i,j) = (*this)[i].ValueOfPartialOfRootPoly(j,p);
  
  return answer;
}

template <class T> gMatrix<T> 
ListOfPartialTrees<T>::DerivativeMatrix(const gVector<T>& p,
					const int& NoEquations) const
{
  gMatrix<T> answer(NoEquations,Dmnsn());
  int i;
  for (i = 1; i <= NoEquations; i++)
    for (int j = 1; j <= Dmnsn(); j++)
      answer(i,j) = (*this)[i].ValueOfPartialOfRootPoly(j,p);
  
  return answer;
}

template <class T> gSquareMatrix<T> 
ListOfPartialTrees<T>::SquareDerivativeMatrix(const gVector<T>& p) const
{
  assert (Length() >= Dmnsn());
  gSquareMatrix<T> answer(Dmnsn());
  int i;
  for (i = 1; i <= Dmnsn(); i++)
    for (int j = 1; j <= Dmnsn(); j++)
      answer(i,j) = (*this)[i].ValueOfPartialOfRootPoly(j,p);
 
  return answer;
}

template <class T> gVector<T> 
ListOfPartialTrees<T>::ValuesOfRootPolys(const gVector<T>& point,
					 const int& NoEquations) const
{
  gVector<T> answer(NoEquations);
  for (int i = 1; i <= NoEquations; i++)
    answer[i] = PartialTreeList[i].EvaluateRootPoly(point);

  return answer;
}

template <class T> T  
ListOfPartialTrees<T>::MaximalNonconstantDifference(const int& i,
						    const int& j,
						    const gVector<T>& point, 
				  const gVector<T>& halvesoflengths)    const
{
  return PartialTreeList[i].MaximalNonconstantDifference(PartialTreeList[j],
							 point,
							 halvesoflengths);
}

template <class T> gOutput& operator << (gOutput& output, 
					 const ListOfPartialTrees<T>& /* x */)
{
  // Changed to just return silently.
  return output;
}