lptab.imp

Gambit 是一个游戏库理论软件

	    c_jo = c_k - a_ik * ratio; 
	    
	    if(LtZero(c_jo)) { 
	      // gout << "\ni not col to enter in prior tableau: ";
	      // gout << "c_jo < 0 for k = " << k;
	      BestSet.Remove(i);
	      flag=1;
	    }
	  }
	}
      }
    }
    // gout << "\nafter checking cols, BestSet = ";
    // BestSet.Dump(gout);

    if(BestSet.Length()>0) 
      for(i=1;i<=BestSet.Length();i++) {
	pivot[1] = BestSet[i];
	pivot[2] = j;
	PivotList.Append(pivot);
      }
  }
}

template <class T>
bool LPTableau<T>::IsReversePivot(int i, int j)
{
  gVector<T> tmpcol(MinRow(),MaxRow());

  // first check that pivot preserves primal feasibility
  
  // gout << "\nin IsReversePivot, i= " << i << " j = "<< j;
  SolveColumn(j,tmpcol);
  gVector<T> solution(tmpcol);  //$$
  BasisVector(solution);        //$$
  // gout << "\ncurrentPivCol = " << tmpcol;
  // gout << "\ncurrentSolCol = " << solution;
  if(LeZero(tmpcol[i])) { 
    // gout << "\nPrior tableau not primal feasible: currentPivCol[i] <= 0";
    return 0;
  }
  int k;
  T ratio = solution[i]/tmpcol[i];
  // gout << "\nratio = " << ratio;
  
  for(k=tmpcol.First();k<=tmpcol.Last();k++)
    if(GtZero(tmpcol[k]) && GtZero(solution[k]/tmpcol[k]-ratio)) {
      // gout << "\nPrior tableau not primal feasible: i not min ratio";
      return 0;
    }
  // check that j would be the row to exit in prior tableau
  
  T a_ij,a_ik,b_i,b_k,c_j,c_k,c_jo;

  a_ij = (T)1/tmpcol[i];
  if(LeZero(a_ij)) {
    // gout << "\nj not row to exit in prior tableau: a_ij <= 0";
    return 0;
  }
  b_i = solution[i]/tmpcol[i];
  ratio = b_i/a_ij;
  
  for(k=tmpcol.First();k<=tmpcol.Last();k++) 
    if(k!=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;
	return 0;
      }
      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;
	return 0;
      }
    }

  // check that i would be the column to enter in prior tableau
  
  int enter = Label(i);
  // gout << "\nenter = " << enter;
  
  gVector<T> tmpdual(MinRow(),MaxRow());
  tmpcol = (T)0;
  tmpcol[i]=(T)1;
  SolveT(tmpcol,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";
    return 0;
  }
  ratio = c_j/a_ij;
  
  if(enter<0) 
    a_ik = tmpdual[-enter];
  else {
//    A->GetColumn(enter,tmpcol);
    GetColumn(enter,tmpcol);
    a_ik = tmpdual*tmpcol;
  }
  c_k = RelativeCost(enter);
  c_jo = c_k - a_ik * ratio; 
  if(GeZero(c_jo)) {
    // gout << "\ni not col to enter in prior tableau: ";
    // gout << "c_jo<0";
    return 0;
  }

  for(k=-b->Last();k<enter;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);
    c_jo = c_k - a_ik * ratio; 
    
    if(LtZero(c_jo)) { 
      // gout << "\ni not col to enter in prior tableau: ";
      // gout << "c_jo < 0 for k = " << k;
      return 0;
    }
  }
  // gout << "\nValid Reverse pivot at i = " << i << " j =  " << j;
  return 1;
}

template <class T>
void LPTableau<T>::DualReversePivots(gList<gArray<int> > &/*list*/)
{
}

template <class T>
bool LPTableau<T>::IsDualReversePivot(int i, int j)
{
  // first check that pivot preserves dual feasibility
  
  // gout << "\nin IsDualReversePivot, i= " << i << " j = "<< j;

  int k;

  gVector<T> tmpcol(MinRow(),MaxRow());
  gVector<T> tmpdual(MinRow(),MaxRow());
  tmpcol = (T)0;
  tmpcol[i]=(T)1;
  SolveT(tmpcol,tmpdual);

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

  // gout << "\ncurrentPivCol = " << tmpcol;
  // gout << "\ncurrentSolCol = " << solution;

  T a_ij,a_ik,c_j,c_k,ratio;

/*  if( j<0 )
    { tmpcol=(T)0; tmpcol[-j]=(T)1; }
  else
    A->GetColumn(j,tmpcol);
  */

    GetColumn(j,tmpcol);
  
  a_ij = tmpdual*tmpcol;
  c_j = RelativeCost(j);
  if(GeZero(a_ij)) {
    // gout << "\nPrior tableau not dual feasible: ";
    // gout << "a_ij>=0";
    return 0;
  }
  ratio = c_j/a_ij;
  
  for(k=-b->Last();k<=cost.Last();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);
    
    if(LtZero(a_ik) && GtZero(c_k/a_ik-ratio)) { 
      // gout << "\nPrior tableau not dual feasible: ";
      // gout << "\nhigher ratio for k = " << k;
      return 0;
    }
  }

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

  int enter = Label(i);
  // gout << "\nenter = " << enter;

  if(enter<0) 
    a_ik = tmpdual[-enter];
  else {
//    A->GetColumn(enter,tmpcol);
    GetColumn(enter,tmpcol);
    a_ik = tmpdual*tmpcol;
  }
  a_ik = a_ik/a_ij;
  c_k = RelativeCost(enter);
  c_k -= a_ik * c_j; 

  if(GeZero(a_ik)) {
    // gout << "\ni not col to enter in prior tableau: ";
    // gout << "a_ik>=0";
    return 0;
  }
  ratio = c_k/a_ik;

  for(k=-b->Last();k<=cost.Last();k++) if(k!=0) {
    if(k<0) 
      a_ik=tmpdual[-k];
    else {
//    A->GetColumn(k,tmpcol);
    GetColumn(k,tmpcol);
    a_ik = tmpdual*tmpcol;
    }
    a_ik = a_ik/a_ij;
    c_k = RelativeCost(k);
    c_k -= a_ik * c_j; 
    
    if(LtZero(a_ik) && GtZero(c_k/a_ik- ratio)) { 
      // gout << "\ni not col to enter in prior tableau: ";
      // gout << "\nhigher ratio for k = " << k;
      return 0;
    }
    if(k<enter && LtZero(a_ik) && EqZero(c_k/a_ik - ratio)) { 
      // gout << "\ni not col to enter in prior tableau: ";
      // gout << "\nsame ratio and lower lex for k = " << k;
      return 0;
    }
  }

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

  SolveColumn(j,tmpcol);
  // gout << "\ncurrentPivCol = " << tmpcol;
  // gout << "\ncurrentSolCol = " << solution;

  T b_k,b_i;
  b_i= solution[i]/tmpcol[i];
  if(LeZero(b_i)) {
    // gout << "\nj not row to exit in prior tableau: ";
    // gout << "b_i<=0";
    return 0;
    }

  
  for(k=b->First();k<=b->Last();k++) 
    if(k!=i) {
      b_k = solution[k] -  b_i * tmpcol[k];
      if(GtZero(b_k) && Label(k)<j) {
	// gout << "\nj not row to exit in prior tableau: ";
	// gout << "same ratio,lower lex at k= " << k;
	return 0;
      }
    }
  // gout << "\nValid Reverse pivot at i = " << i << " j =  " << j;
  return 1;
}

/*
template <class T>
BFS<T> LPTableau<T>::DualBFS() const
{
  BFS<T> cbfs((T) 0);
  for(int i=MinRow();i<=MaxRow();i++) {
    if(!Member(-i)) 
      cbfs.Define(-i,dual[i]);
  }
  return cbfs;
}
*/

template <class T>
BFS<T> LPTableau<T>::DualBFS() const
{
  BFS<T> cbfs((T) 0);
  gVector<T> d(MinRow(),MaxRow());
  DualVector(d);

  for(int i=MinRow();i<=MaxRow();i++) {
    if(!Member(-i)) 
      cbfs.Define(-i,d[i]);
  }
  // gout << "\ndual: " << d;
  return cbfs;
}

template <class T>
int LPTableau<T>::LastLabel( void )
{
  return artificial.Last();
}

template <class T>
void LPTableau<T>::BigDump(gOutput &to)
{
  //   Tableau<T>::BigDump(to);

  to << "\nBasis:";
  basis.Dump(to);
  to << "\nMinCol(): " << MinCol(); 
  to << "\nMaxCol(): " << MaxCol(); 
  to << "\nLastlabel(): " << LastLabel(); 
  to << "\nNumRows(): " << (*A).NumRows(); 
  gMatrix<T> AA(MinRow(),MaxRow(),MinCol(),LastLabel()+(*A).NumRows());
  gVector<T> bb(MinRow(), MaxRow());

  BasisVector(bb);
  to << "\nBasisVector:\n" << bb;

  for(int j=MinCol();j<=LastLabel();j++) {
    SolveColumn(j, bb);
    for(int i=AA.MinRow();i<=AA.MaxRow();i++) 
      AA(i,j) = bb[i];
  }
  for(int j=MinRow();j<=MaxRow();j++) {
    SolveColumn(-j, bb);
    for(int i=AA.MinRow();i<=AA.MaxRow();i++) 
      AA(i,LastLabel()+j) = bb[i];
  }
  to << "\nTableau:\n" << AA;

  to << "\nCost vector: " << GetCost();
  to << "\nUnit Cost  : " << GetUnitCost();
  to << "\nTotal Cost: " << TotalCost();
}

template <class T>
void LPTableau<T>::BasisSelect(const gBlock<T> &rowv, gVector<T> &colv) const
{
  for(int i=basis.First(); i<=basis.Last(); i++) {
    if(basis.Label(i)<0)
      colv[i]= 0;
    else
      colv[i]= rowv[basis.Label(i)];
  }
}

template <class T>
void LPTableau<T>::BasisSelect(const gBlock<T> &unitv,
			   const gBlock<T> &rowv,
			   gVector<T> &colv ) const
{
  for(int i=basis.First(); i<=basis.Last(); i++) {
    if(basis.Label(i)<0)
      colv[i]= unitv[-basis.Label(i)];
    else
      colv[i]= rowv[basis.Label(i)];
  }
}

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

template <class T> gText LPTableau<T>::BadPivot::Description(void) const
{
  return "Bad Pivot in LPTableau";
}

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

template <class T> gText LPTableau<T>::BadDim::Description(void) const
{
  return "Bad Dimension in LPTableau";
}