lptab.imp

Gambit 是一个游戏库理论软件

//
// $Source: /home/gambit/CVS/gambit/sources/numerical/lptab.imp,v $
// $Date: 2002/09/26 17:50:55 $
// $Revision: 1.1.2.1 $
//
// DESCRIPTION:
// Implementation of LP tableaus
//
// 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 "lptab.h"
#include "base/base.h"

// ---------------------------------------------------------------------------
//                   LPTableau member definitions 
// ---------------------------------------------------------------------------


template <class T>
LPTableau<T>::LPTableau(const gMatrix<T> &A, const gVector<T> &b)
  : Tableau<T>(A,b), dual(A.MinRow(),A.MaxRow()),
    unitcost(A.MinRow(),A.MaxRow()), cost(A.MinCol(),A.MaxCol())
{ };

template <class T>
LPTableau<T>::LPTableau(const gMatrix<T> &A, const gBlock<int> &art, 
			const gVector<T> &b)
  : Tableau<T>(A,art,b), dual(A.MinRow(),A.MaxRow()),
    unitcost(A.MinRow(),A.MaxRow()), cost(A.MinCol(),A.MaxCol()+art.Length())
{ };

template <class T>
LPTableau<T>::LPTableau(const LPTableau<T> &orig)
  : Tableau<T>(orig), dual(orig.dual),  unitcost(orig.unitcost),
    cost(orig.cost)
{ }

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


template <class T>
LPTableau<T>& LPTableau<T>::operator=(const LPTableau<T> &orig)
{
  Tableau<T>::operator=(orig);
  if(this!= &orig) {
    dual=orig.dual;
    unitcost= orig.unitcost;
    cost= orig.cost;
  }
  return *this;
}

// cost-based functions

TEMPLATE_SPECIALIZATION()
void LPTableau<gRational>::SetCost(const gVector<gRational>& c)
{
  int i;
  if(cost.First()==c.First() && cost.Last()==c.Last()) {
    for(i=cost.First();i<=cost.Last();i++) cost[i] = c[i]*(gRational)Tableau<gRational>::TotDenom();
    for(i=unitcost.First();i<=unitcost.Last();i++) unitcost[i] = (gRational)0;
    Refactor();
    SolveDual();
    return;
  }
  if(c.First()!=cost.First()) throw BadDim();
  if(c.Last()!=(cost.Last()+unitcost.Length())) throw BadDim();
  for(i=c.First();i<=cost.Last();i++)
    cost[i]=c[i]*(gRational)Tableau<gRational>::TotDenom();
  for(i=unitcost.First();i<=unitcost.Last();i++)
    unitcost[i]=c[cost.Length()+i-unitcost.First()+1];
  // gout << "\nc: " << c.First() << " " << c.Last() << " " << c;
  // gout << "\ncost: " << cost.First() << " " << cost.Last() << " " << cost;
  // gout << "\nunit: " << unitcost.First() << " " << unitcost.Last() << " " << unitcost;
  //** added for gRational
  Refactor();
  SolveDual();
}

TEMPLATE_SPECIALIZATION()
void LPTableau<double>::SetCost(const gVector<double>& c)
{
  int i;
  if(cost.First()==c.First() && cost.Last()==c.Last()) {
    for(i=cost.First();i<=cost.Last();i++) cost[i] = c[i];
    for(i=unitcost.First();i<=unitcost.Last();i++) unitcost[i] = (double)0;
    Refactor();
    SolveDual();
    return;
  }
  if(c.First()!=cost.First()) throw BadDim();
  if(c.Last()!=(cost.Last()+unitcost.Length())) throw BadDim();
  for(i=c.First();i<=cost.Last();i++)
    cost[i]=c[i];
  for(i=unitcost.First();i<=unitcost.Last();i++)
    unitcost[i]=c[cost.Length()+i-unitcost.First()+1];
  // gout << "\nc: " << c.First() << " " << c.Last() << " " << c;
  // gout << "\ncost: " << cost.First() << " " << cost.Last() << " " << cost;
  // gout << "\nunit: " << unitcost.First() << " " << unitcost.Last() << " " << unitcost;
  //** added for gRational
  Refactor();
  SolveDual();
}

template <class T>
void LPTableau<T>::SetCost(const gVector<T> &uc, const gVector<T> &c)
{
  if(cost.First()!=c.First() || cost.Last()!=c.Last()) 
    throw BadDim();
  if(unitcost.First()!=uc.First() || unitcost.Last()!=uc.Last()) 
    throw BadDim();
  int i;
  for(i=cost.First();i<=cost.Last();i++) cost[i]=c[i];
  for(i=unitcost.First();i<=unitcost.Last();i++) unitcost[i]=uc[i];
  SolveDual();
}


template <class T>
gVector<T> LPTableau<T>::GetCost(void) const
{
  gVector<T> x(cost.First(),cost.Last());
  for(int i=x.First();i<=x.Last();i++)x[i]=cost[i];
  return x;
  //  return cost; 
}


template <class T>
gVector<T> LPTableau<T>::GetUnitCost(void) const
{
  gVector<T> x(unitcost.First(),unitcost.Last());
  for(int i=x.First();i<=x.Last();i++)x[i]=unitcost[i];
  return x;
}

TEMPLATE_SPECIALIZATION() double LPTableau<double>::TotalCost(void)
{
  gVector<double> tmpcol(MinRow(),MaxRow());
  BasisSelect(unitcost,cost,tmpcol);
  return tmpcol*solution;
}

TEMPLATE_SPECIALIZATION() gRational LPTableau<gRational>::TotalCost(void)
{
  gVector<gRational> tmpcol(MinRow(),MaxRow());
  gVector<gRational> sol(MinRow(),MaxRow());
  BasisSelect(unitcost,cost,tmpcol);
  BasisVector(sol);
  for(int i=tmpcol.First();i<=tmpcol.Last();i++)
    if(Label(i)>0)tmpcol[i]/=(gRational)Tableau<gRational>::TotDenom();

  return tmpcol*sol;
}

TEMPLATE_SPECIALIZATION()
void LPTableau<double>::DualVector(gVector<double> &L) const
{
  L= dual;
}

TEMPLATE_SPECIALIZATION()
void LPTableau<gRational>::DualVector(gVector<gRational> &out) const
{
  out= dual;
  //  for(int i=out.First();i<=out.Last();i++) 
  //  if(Label(i)>=0) out[i]*=TotDenom();
}

template <class T>
T LPTableau<T>::RelativeCost(int col) const
{
  gVector<T> tmpcol(MinRow(),MaxRow());
  if( col<0 ) {
    return unitcost[-col] - dual[-col];
  }
  else {
    GetColumn(col, (gVector<T> &)tmpcol);
    return cost[col] - dual*tmpcol;
  }
}

/*

template <class T>
void LPTableau<T>::RelativeCostVector(gVector<T> &relunitcost,
				      gVector<T> &relcost) const
{

  if(!A->CheckColumn(relunitcost)) throw BadDim();
  if(!A->CheckRow(relcost)) throw BadDim();
  
  relunitcost= unitcost - dual;
  relcost= cost - dual*A;  // pre multiplication not defined?  
}
*/


template <class T>
void LPTableau<T>::SolveDual()
{
  gVector<T> tmpcol1(MinRow(),MaxRow());
  BasisSelect(unitcost,cost,tmpcol1);
  SolveT(tmpcol1,dual);
}

// Redefined functions

template <class T>
void LPTableau<T>::Refactor()
{
  // gout << "\nIn LPTableau<T>::Refactor()";  
  Tableau<T>::Refactor();
  SolveDual();
}

template <class T>
void LPTableau<T>::Pivot(int outrow,int col)
{
  // gout << "\nIn LPTableau<T>::Pivot() ";
  // gout << "outrow: " << outrow << " col: " << col;
  // BigDump(gout);
  Tableau<T>::Pivot(outrow,col);
  SolveDual();
}

template <class T>
void LPTableau<T>::ReversePivots(gList<gArray<int> > &PivotList)
{
  gVector<T> tmpcol(MinRow(),MaxRow());
  // gout << "\nIn LPTableau<T>::ReversePivots";
  bool flag;
  int i,j,k,enter;
  T ratio,a_ij,a_ik,b_i,b_k,c_j,c_k,c_jo,x;
  gList<int> BestSet;
  gArray<int> pivot(2);
  gVector<T> tmpdual(MinRow(),MaxRow());

  gVector<T> solution(tmpcol);  //$$
  BasisVector(solution);        //$$

  // BigDump(gout);
  // gout << "\ncost: " << GetCost();
  // gout << "\nunitcost: " << GetUnitCost() << "\n";
  // for(i=MinCol();i<=MaxCol();i++) gout << " " << RelativeCost(i);
  // for(i=MinRow();i<=MaxRow();i++) gout << " " << RelativeCost(-i);

  for(j=-MaxRow();j<=MaxCol();j++) if(j && !Member(j)  && !IsBlocked(j)) {
    SolveColumn(j,tmpcol);
    // gout << "\nColumn " << j;
    // gout << "\nPivCol = " << tmpcol;
    // gout << "\ncurrentSolCol = " << solution;
    
    // find all i where prior tableau is primal feasible
    
    BestSet.Flush();
    for(i=MinRow();i<=MaxRow();i++)
      if(GtZero(tmpcol[i])) BestSet.Append(i);
    if(BestSet.Length()>0) {
      ratio = solution[BestSet[1]]/tmpcol[BestSet[1]];
      // find max ratio
      for(i=2;i<=BestSet.Length();i++) {
	x = solution[BestSet[i]]/tmpcol[BestSet[i]];
	if(GtZero(x-ratio)) ratio = x;
      }
      // eliminate nonmaximizers
      for(i=BestSet.Length();i>=1;i--) {
	x = solution[BestSet[i]]/tmpcol[BestSet[i]];
	if(LtZero(x-ratio)) BestSet.Remove(i);
      }	

      // check that j would be the row to exit in prior tableau

      // first check that prior pivot entry > 0 
      for(i=BestSet.Length();i>=1;i--) {
	a_ij = (T)1/tmpcol[BestSet[i]];
	if(LeZero(a_ij)) {
	  // gout << "\nj not row to exit in prior tableau: a_ij <= 0";
	  BestSet.Remove(i);
	}
	else {
	  // next check that prior pivot entry attains max ratio
	  b_i = solution[BestSet[i]]/tmpcol[BestSet[i]];
	  ratio = b_i/a_ij;
  
	  flag = 0;
	  for(k=tmpcol.First();k<=tmpcol.Last() && !flag;k++) 
	    if(k!=BestSet[i]) {
	      a_ik = - a_ij * tmpcol[k];
	      b_k = solution[k] - b_i*tmpcol[k];
	      if(GtZero(a_ik) && GtZero(b_k/a_ik -ratio)) {
		// gout << "\nj not row to exit in prior tableau: ";
		// gout << "higher ratio at row= " << k;
		BestSet.Remove(i);
		flag = 1;
	      }
	      else if(GtZero(a_ik) && EqZero(b_k/a_ik-ratio) && Label(k)<j) {
		// gout << "\nj not row to exit in prior tableau: ";
		// gout << "same ratio,lower lex at k= " << k;
		BestSet.Remove(i);
		flag = 1;
	      }
	    }
	}
      }
    }
    // gout << "\nafter checking rows, BestSet = ";
    // BestSet.Dump(gout);

    // check that i would be the column to enter in prior tableau

    for(i=BestSet.Length();i>=1;i--) {
      enter = Label(BestSet[i]);
      // gout << "\nenter = " << enter;
      
      tmpcol = (T)0;
      tmpcol[BestSet[i]]=(T)1;
      // gout << "\ntmpcol, loc 1: " << tmpcol;
      SolveT(tmpcol,tmpdual);
      // gout << "\ntmpcol, loc 2: " << tmpcol;
      // gout << "\ntmpdual, loc 1: " << tmpdual;
      
/*      if( j<0 )
	{ tmpcol=(T)0; tmpcol[-j]=(T)1; }
      else
	A->GetColumn(j,tmpcol);
*/
      GetColumn(j,tmpcol);      
      // gout << "\ncol " << j << ": " << tmpcol;
      a_ij = tmpdual*tmpcol;
      c_j = RelativeCost(j);
      if(EqZero(a_ij)) {
	// gout << "\ni not col to enter in prior tableau: ";
	// gout << "a_ij=0";
	BestSet.Remove(i);
      }
      else {
	ratio = c_j/a_ij;
	// gout << " ratio: " << ratio;
	if(enter<0) 
	  a_ik = tmpdual[-enter];
	else {
	  GetColumn(enter,tmpcol);
//	  A->GetColumn(enter,tmpcol);
	  a_ik = tmpdual*tmpcol;
	}
	c_k = RelativeCost(enter);
	c_jo = c_k - a_ik * ratio; 
	// gout << "\ntmpdual = " << tmpdual << "\n";
	// gout << " c_j:" << c_j; 
	// gout << " c_k:" << c_k; 
	// gout << " c_jo:" << c_jo; 
	// gout << " a_ij:" << a_ij; 
	// gout << " a_ik:" << a_ik; 
	if(GeZero(c_jo)) {
	  // gout << "\ni not col to enter in prior tableau: ";
	  // gout << "c_jo<0";
	  BestSet.Remove(i);
	}
	else {
	  flag=0;
	  for(k=-b->Last();k<enter && !flag;k++) if(k!=0) {
	    if(k<0) 
	      a_ik=tmpdual[-k];
	    else {
//	      A->GetColumn(k,tmpcol);
	      GetColumn(k,tmpcol);
	      a_ik = tmpdual*tmpcol;
	    }
	    c_k = RelativeCost(k);