//
// $Source: /cvsroot/gambit/gambit/sources/tools/enummixed/lrsnash.cc,v $
// $Date: 2006/11/13 06:16:24 $
// $Revision: 1.3 $
//
// DESCRIPTION:
// Compute Nash equilibria via enumerating extreme points, lrslib version
//
// This file is part of Gambit
// Copyright (c) 2006, The Gambit Project
// Based on the implementation in lrslib 4.2b, which is
// Copyright (c) 1995-2005, David Avis
//
// 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.
//
//
// This file is based heavily on nash.c distributed with lrslib 4.2b.
// I have tried where possible to maintain as much of the original, and
// make a minimal set of changes to interface with Gambit's game
// representation library. This way, hopefully, if future versions of lrslib
// fix bugs or create enhancements that are relevant, it will be fairly
// easy to port them to this program. -- TLT, 23.viii.2006
//
#include <iostream>
#include <stdio.h>
#include <string.h>
// The order of these next includes is important, because of macro definitions
#include "libgambit/libgambit.h"
using namespace Gambit;
extern "C" {
#include "lrslib.h"
}
//==========================================================================
// Building the problem representations
//==========================================================================
//
// These two functions are based upon the program setupnash.c from the
// lrslib distribution, and the user's guide documentation.
// There are two separate functions, one for each player's problem.
// According to the user's guide, the ordering of the constraint rows
// is significant, and differs between the players; for player 1's problem
// the nonnegativity constraints come first, whereas for player 2's problem
// they appear later. Experiments suggest this is in fact true, and
// reversing them breaks something.
//
void FillNonnegativityRows(lrs_dic *P, lrs_dat *Q,
int firstRow, int lastRow, int n)
{
const int MAXCOL = 1000; /* maximum number of columns */
long num[MAXCOL], den[MAXCOL];
for (long row = firstRow; row <= lastRow; row++) {
num[0] = 0;
den[0] = 1;
for (long col = 1; col < n; col++) {
num[col] = (row-firstRow+1 == col) ? 1 : 0;
den[col] = 1;
}
lrs_set_row(P, Q, row, num, den, GE);
}
}
void FillConstraintRows(lrs_dic *P, lrs_dat *Q,
const StrategySupport &p_support,
int p1, int p2, int firstRow)
{
const int MAXCOL = 1000; /* maximum number of columns */
long num[MAXCOL], den[MAXCOL];
Game game = p_support.GetGame();
Rational min = game->GetMinPayoff() - Rational(1);
PureStrategyProfile cont(game);
for (long row = firstRow; row < firstRow + p_support.NumStrategies(p1);
row++) {
num[0] = 0;
den[0] = 1;
cont.SetStrategy(p_support.GetStrategy(p1, row - firstRow + 1));
for (long st = 1; st <= p_support.NumStrategies(p2); st++) {
cont.SetStrategy(p_support.GetStrategy(p2, st));
Rational x = cont.GetPayoff<Rational>(p1) - min;
num[st] = -x.numerator().as_long();
den[st] = x.denominator().as_long();
}
num[p_support.NumStrategies(p2)+1] = 1;
den[p_support.NumStrategies(p2)+1] = 1;
lrs_set_row(P, Q, row, num, den, GE);
}
}
void FillLinearityRow(lrs_dic *P, lrs_dat *Q, int m, int n)
{
const int MAXCOL = 1000; /* maximum number of columns */
long num[MAXCOL], den[MAXCOL];
num[0] = -1;
den[0] = 1;
for (int i = 1; i < n-1; i++) {
num[i] = 1;
den[i] = 1;
}
num[n-1] = 0;
den[n-1] = 1;
lrs_set_row(P, Q, m, num, den, EQ);
}
//
// Build the H-representation for player p1
//
void BuildRep(lrs_dic *P, lrs_dat *Q, const StrategySupport &p_support,
int p1, int p2)
{
long m=Q->m; /* number of inequalities */
long n=Q->n;
if (p1 == 1) {
FillConstraintRows(P, Q, p_support, p1, p2, 1);
FillNonnegativityRows(P, Q, p_support.NumStrategies(p1) + 1,
p_support.MixedProfileLength(), n);
}
else {
FillNonnegativityRows(P, Q, 1, p_support.NumStrategies(p2), n);
FillConstraintRows(P, Q, p_support, p1, p2,
p_support.NumStrategies(p2) + 1);
}
FillLinearityRow(P, Q, m, n);
}
//========================================================================
// Forward declarations of useful functions
//========================================================================
/* lrs driver, argv[2]= 2nd input file for nash equilibria */
long nash2_main (lrs_dic *P1, lrs_dat *Q1, lrs_dic *P2orig,
lrs_dat *Q2, long *numequilib,
lrs_mp_vector output1, lrs_mp_vector output2,
const StrategySupport &p_support);
long lrs_getfirstbasis2 (lrs_dic ** D_p, lrs_dat * Q, lrs_dic *P2orig,
lrs_mp_matrix * Lin, long no_output);
long getabasis2 (lrs_dic * P, lrs_dat * Q, lrs_dic * P2orig, long order[]);
// This is a modified version of lrs_output from the original, in which
// we output the equilibria found in Gambit format.
void nashoutput(lrs_dat *Q1, lrs_mp_vector output1,
lrs_dat *Q2, lrs_mp_vector output2,
const StrategySupport &p_support);
//
// This is the main function, based on main() from lrslib's 'nash' driver.
//
void LrsSolve(const StrategySupport &p_support)
{
lrs_dic *P1,*P2; /* structure for holding current dictionary and indices */
lrs_dat *Q1,*Q2; /* structure for holding static problem data */
lrs_mp_vector output1; /* holds one line of output; ray,vertex,facet,linearity */
lrs_mp_vector output2; /* holds one line of output; ray,vertex,facet,linearity */
lrs_mp_matrix Lin; /* holds input linearities if any are found */
lrs_dic *P2orig; /* we will save player 2's dictionary in getabasis */
long col; /* output column index for dictionary */
long startcol = 0;
long prune = FALSE; /* if TRUE, getnextbasis will prune tree and backtrack */
long numequilib=0; /* number of nash equilibria found */
long oldnum=0;
/* global variables lrs_ifp and lrs_ofp are file pointers for input and output */
/* they default to stdin and stdout, but may be overidden by command line parms. */
/***************************************************
Step 0:
Do some global initialization that should only be done once,
no matter how many lrs_dat records are allocated. db
***************************************************/
if (!lrs_init("")) {
return;
}
/*********************************************************************************/
/* Step 1: Allocate lrs_dat, lrs_dic and set up the problem */
/*********************************************************************************/
Q1 = lrs_alloc_dat ("LRS globals"); /* allocate and init structure for static problem data */
if (Q1 == NULL) {
return;
}
Q1->nash=TRUE;
Q1->n = p_support.NumStrategies(1) + 2;
Q1->m = p_support.MixedProfileLength() + 1;
P1 = lrs_alloc_dic (Q1); /* allocate and initialize lrs_dic */
if (P1 == NULL) {
return;
}
BuildRep(P1, Q1, p_support, 2, 1);
output1 = lrs_alloc_mp_vector (Q1->n + Q1->m); /* output holds one line of output from dictionary */
/* allocate and init structure for player 2's problem data */
Q2 = lrs_alloc_dat ("LRS globals");
if (Q2 == NULL) {
return;
}
Q2->nash=TRUE;
Q2->n = p_support.NumStrategies(2) + 2;
Q2->m = p_support.MixedProfileLength() + 1;
P2 = lrs_alloc_dic (Q2); /* allocate and initialize lrs_dic */
if (P2 == NULL) {
return;
}
BuildRep(P2, Q2, p_support, 1, 2);
output2 = lrs_alloc_mp_vector (Q2->n + Q2->m); /* output holds one line of output from dictionary */
P2orig = lrs_getdic(Q2); /* allocate and initialize lrs_dic */
if (P2orig == NULL)
return;
copy_dict(Q2,P2orig,P2);
/*********************************************************************************/
/* Step 2: Find a starting cobasis from default of specified order */
/* P1 is created to hold active dictionary data and may be cached */
/* Lin is created if necessary to hold linearity space */
/* Print linearity space if any, and retrieve output from first dict. */
/*********************************************************************************/
if (!lrs_getfirstbasis (&P1, Q1, &Lin, TRUE))
return;
if (Q1->dualdeg)
{
printf("\n*Warning! Dual degenerate, ouput may be incomplete");
printf("\n*Recommendation: Add dualperturb option before maximize in first input file\n");
}
if (Q1->unbounded)
{
printf("\n*Warning! Unbounded starting dictionary for p1, output may be incomplete");
printf("\n*Recommendation: Change/remove maximize option, or include bounds \n");
}
/* Pivot to a starting dictionary */
/* There may have been column redundancy */
/* If so the linearity space is obtained and redundant */
/* columns are removed. User can access linearity space */
/* from lrs_mp_matrix Lin dimensions nredundcol x d+1 */
if (Q1->homogeneous && Q1->hull)
startcol++; /* col zero not treated as redundant */
for (col = startcol; col < Q1->nredundcol; col++) /* print linearity space */
lrs_printoutput (Q1, Lin[col]); /* Array Lin[][] holds the coeffs. */
/*********************************************************************************/
/* Step 3: Terminate if lponly option set, otherwise initiate a reverse */
/* search from the starting dictionary. Get output for each new dict. */
/*********************************************************************************/
/* We initiate reverse search from this dictionary */
/* getting new dictionaries until the search is complete */
/* User can access each output line from output which is */
/* vertex/ray/facet from the lrs_mp_vector output */
/* prune is TRUE if tree should be pruned at current node */
do
{
prune=lrs_checkbound(P1,Q1);
if (!prune && lrs_getsolution (P1, Q1, output1, col))
{
oldnum=numequilib;
nash2_main(P1,Q1,P2orig,Q2,&numequilib,output1,output2,p_support);
if (numequilib > oldnum || Q1->verbose)
{
if(Q1->verbose)
prat(" \np2's obj value: ",P1->objnum,P1->objden);
//fprintf (lrs_ofp, "\n");
}
}
}
while (lrs_getnextbasis (&P1, Q1, prune));
lrs_clear_mp_vector(output1, Q1->m + Q1->n);
lrs_clear_mp_vector(output2, Q2->m + Q2->n);
Q2->Qhead = P2; /* reset this or you crash free_dic */
lrs_free_dic (P2,Q2); /* deallocate lrs_dic */
lrs_free_dat (Q2); /* deallocate lrs_dat */
lrs_free_dic (P1,Q1); /* deallocate lrs_dic */
lrs_free_dat (Q1); /* deallocate lrs_dat */
lrs_close ("");
fclose(lrs_ofp);
}
/*********************************************/
/* end of nash driver */
/*********************************************/
/**********************************************************/
/* nash2_main is a second driver used in computing nash */
/* equilibria on a second polytope interleaved with first */
/**********************************************************/
long nash2_main (lrs_dic *P1, lrs_dat *Q1, lrs_dic *P2orig,
lrs_dat *Q2, long *numequilib,
lrs_mp_vector output1, lrs_mp_vector output2,
const StrategySupport &p_support)
{
lrs_dic *P2; /* This can get resized, cached etc. Loaded from P2orig */
lrs_mp_matrix Lin; /* holds input linearities if any are found */
long col; /* output column index for dictionary */
long startcol = 0;
long prune = FALSE; /* if TRUE, getnextbasis will prune tree and backtrack */
long nlinearity;
long *linearity;
static long firstwarning=TRUE; /* FALSE if dual deg warning for Q2 already given */
static long firstunbounded=TRUE; /* FALSE if dual deg warning for Q2 already given */
long i,j;
/* global variables lrs_ifp and lrs_ofp are file pointers for input and output */
/* they default to stdin and stdout, but may be overidden by command line parms. */
/*********************************************************************************/
/* Step 1: Allocate lrs_dat, lrs_dic and set up the problem */
/*********************************************************************************/
P2=lrs_getdic(Q2);
copy_dict(Q2,P2,P2orig);
/* Here we take the linearities generated by the current vertex of player 1*/
/* and append them to the linearity in player 2's input matrix */
/* next is the key magic linking player 1 and 2 */
/* be careful if you mess with this! */
linearity=Q2->linearity;
nlinearity=0;
for(i=Q1->lastdv+1;i <= P1->m; i++)
{
if (!zero(P1->A[P1->Row[i]][0]))
{
j = Q1->inequality[P1->B[i]-Q1->lastdv];
if (j < Q1->linearity[0])
linearity[nlinearity++]= j;
}
}
/* add back in the linearity for probs summing to one */
linearity[nlinearity++]= Q1->linearity[0];
/*sort linearities */
for (i = 1; i < nlinearity; i++)
reorder (linearity, nlinearity);
if(Q2->verbose)
{
fprintf(lrs_ofp,"\np2: linearities %ld",nlinearity);
for (i=0;i < nlinearity; i++)
fprintf(lrs_ofp," %ld",linearity[i]);
}
Q2->nlinearity = nlinearity;
Q2->polytope = FALSE;
/*********************************************************************************/
/* Step 2: Find a starting cobasis from default of specified order */
/* P2 is created to hold active dictionary data and may be cached */
/* Lin is created if necessary to hold linearity space */
/* Print linearity space if any, and retrieve output from first dict. */
/*********************************************************************************/
if (!lrs_getfirstbasis2 (&P2, Q2, P2orig, &Lin, TRUE))
goto sayonara;
if (firstwarning && Q2->dualdeg)
{
firstwarning=FALSE;
printf("\n*Warning! Dual degenerate, ouput may be incomplete");
printf("\n*Recommendation: Add dualperturb option before maximize in second input file\n");
}
if (firstunbounded && Q2->unbounded)
{
firstunbounded=FALSE;
printf("\n*Warning! Unbounded starting dictionary for p2, output may be incomplete");
printf("\n*Recommendation: Change/remove maximize option, or include bounds \n");
}
/* Pivot to a starting dictionary */
/* There may have been column redundancy */
/* If so the linearity space is obtained and redundant */
/* columns are removed. User can access linearity space */
/* from lrs_mp_matrix Lin dimensions nredundcol x d+1 */
if (Q2->homogeneous && Q2->hull)
startcol++; /* col zero not treated as redundant */
/* for (col = startcol; col < Q2->nredundcol; col++)*/ /* print linearity space */
/*lrs_printoutput (Q2, Lin[col]);*/ /* Array Lin[][] holds the coeffs. */
/*********************************************************************************/
/* Step 3: Terminate if lponly option set, otherwise initiate a reverse */
/* search from the starting dictionary. Get output for each new dict. */
/*********************************************************************************/
/* We initiate reverse search from this dictionary */
/* getting new dictionaries until the search is complete */
/* User can access each output line from output which is */
/* vertex/ray/facet from the lrs_mp_vector output */
/* prune is TRUE if tree should be pruned at current node */
do
{
prune=lrs_checkbound(P2,Q2);
col=0;
if (!prune && lrs_getsolution (P2, Q2, output2, col))
{
(*numequilib)++;
if (Q2->verbose)
prat(" \np1's obj value: ",P2->objnum,P2->objden);
nashoutput(Q1, output1, Q2, output2, p_support);
}
}
while (lrs_getnextbasis (&P2, Q2, prune));
sayonara:
lrs_free_dic(P2,Q2);
return 0;
}
/*********************************************/
/* end of nash2_main */
/*********************************************/
/* In lrs_getfirstbasis and lrs_getnextbasis we use D instead of P */
/* since the dictionary P may change, ie. &P in calling routine */
#define D (*D_p)
long
lrs_getfirstbasis2 (lrs_dic ** D_p, lrs_dat * Q, lrs_dic * P2orig, lrs_mp_matrix * Lin, long no_output)
/* gets first basis, FALSE if none */
/* P may get changed if lin. space Lin found */
/* no_output is TRUE supresses output headers */
{
long i, j, k;
/* assign local variables to structures */
lrs_mp_matrix A;
long *B, *C, *Row, *Col;
long *inequality;
long *linearity;
long hull = Q->hull;
long m, d, lastdv, nlinearity, nredundcol;
static long ocount=0;
m = D->m;
d = D->d;
lastdv = Q->lastdv;
nredundcol = 0L; /* will be set after getabasis */
nlinearity = Q->nlinearity; /* may be reset if new linearity read */
linearity = Q->linearity;
A = D->A;
B = D->B;
C = D->C;
Row = D->Row;
Col = D->Col;
inequality = Q->inequality;
/* default is to look for starting cobasis using linearies first, then */
/* filling in from last rows of input as necessary */
/* linearity array is assumed sorted here */
/* note if restart/given start inequality indices already in place */
/* from nlinearity..d-1 */
for (i = 0; i < nlinearity; i++) /* put linearities first in the order */
inequality[i] = linearity[i];
k = 0; /* index for linearity array */
if (Q->givenstart)
k = d;
else
k = nlinearity;
for (i = m; i >= 1; i--)
{
j = 0;
while (j < k && inequality[j] != i)
j++; /* see if i is in inequality */
if (j == k)
inequality[k++] = i;
}
if (Q->debug)
{
fprintf (lrs_ofp, "\n*Starting cobasis uses input row order");
for (i = 0; i < m; i++)
fprintf (lrs_ofp, " %ld", inequality[i]);
}
if (!Q->maximize && !Q->minimize)
for (j = 0; j <= d; j++)
itomp (ZERO, A[0][j]);
/* Now we pivot to standard form, and then find a primal feasible basis */
/* Note these steps MUST be done, even if restarting, in order to get */
/* the same index/inequality correspondance we had for the original prob. */
/* The inequality array is used to give the insertion order */
/* and is defaulted to the last d rows when givenstart=FALSE */
if (!getabasis2 (D, Q,P2orig, inequality))
return FALSE;
if(Q->debug)
{
fprintf(lrs_ofp,"\nafter getabasis2");
printA(D, Q);
}
nredundcol = Q->nredundcol;
lastdv = Q->lastdv;
d = D->d;
/********************************************************************/
/* now we start printing the output file unless no output requested */
/********************************************************************/
if (!no_output || Q->debug)
{
fprintf (lrs_ofp, "\nV-representation");
/* Print linearity space */
/* Don't print linearity if first column zero in hull computation */
k = 0;
if (nredundcol > k)
{
fprintf (lrs_ofp, "\nlinearity %ld ", nredundcol - k); /*adjust nredundcol for homog. */
for (i = 1; i <= nredundcol - k; i++)
fprintf (lrs_ofp, " %ld", i);
} /* end print of linearity space */
fprintf (lrs_ofp, "\nbegin");
fprintf (lrs_ofp, "\n***** %ld rational", Q->n);
} /* end of if !no_output ....... */
/* Reset up the inequality array to remember which index is which input inequality */
/* inequality[B[i]-lastdv] is row number of the inequality with index B[i] */
/* inequality[C[i]-lastdv] is row number of the inequality with index C[i] */
for (i = 1; i <= m; i++)
inequality[i] = i;
if (nlinearity > 0) /* some cobasic indices will be removed */
{
for (i = 0; i < nlinearity; i++) /* remove input linearity indices */
inequality[linearity[i]] = 0;
k = 1; /* counter for linearities */
for (i = 1; i <= m - nlinearity; i++)
{
while (k <= m && inequality[k] == 0)
k++; /* skip zeroes in corr. to linearity */
inequality[i] = inequality[k++];
}
} /* end if linearity */
if (Q->debug)
{
fprintf (lrs_ofp, "\ninequality array initialization:");
for (i = 1; i <= m - nlinearity; i++)
fprintf (lrs_ofp, " %ld", inequality[i]);
}
if (nredundcol > 0)
{
*Lin = lrs_alloc_mp_matrix (nredundcol, Q->n);
for (i = 0; i < nredundcol; i++)
{
if (!(Q->homogeneous && Q->hull && i == 0)) /* skip redund col 1 for homog. hull */
{
lrs_getray (D, Q, Col[0], D->C[0] + i - hull, (*Lin)[i]); /* adjust index for deletions */
}
if (!removecobasicindex (D, Q, 0L))
return FALSE;
}
} /* end if nredundcol > 0 */
if (Q->verbose)
{
fprintf (lrs_ofp, "\nNumber of pivots for starting dictionary: %ld",Q->count[3]);
ocount=Q->count[3];
}
/* Do dual pivots to get primal feasibility */
if (!primalfeasible (D, Q))
{
if ( Q->verbose )
{
fprintf (lrs_ofp, "\nNumber of pivots for feasible solution: %ld",Q->count[3]);
fprintf (lrs_ofp, " - No feasible solution");
ocount=Q->count[3];
}
return FALSE;
}
if (Q->verbose)
{
fprintf (lrs_ofp, "\nNumber of pivots for feasible solution: %ld",Q->count[3]);
ocount=Q->count[3];
}
/* Now solve LP if objective function was given */
if (Q->maximize || Q->minimize)
{
Q->unbounded = !lrs_solvelp (D, Q, Q->maximize);
/* check to see if objective is dual degenerate */
j = 1;
while (j <= d && !zero (A[0][j]))
j++;
if (j <= d)
Q->dualdeg = TRUE;
}
else
/* re-initialize cost row to -det */
{
for (j = 1; j <= d; j++)
{
copy (A[0][j], D->det);
storesign (A[0][j], NEG);
}
itomp (ZERO, A[0][0]); /* zero optimum objective value */
}
/* reindex basis to 0..m if necessary */
/* we use the fact that cobases are sorted by index value */
if (Q->debug)
printA (D, Q);
while (C[0] <= m)
{
i = C[0];
j = inequality[B[i] - lastdv];
inequality[B[i] - lastdv] = inequality[C[0] - lastdv];
inequality[C[0] - lastdv] = j;
C[0] = B[i];
B[i] = i;
reorder1 (C, Col, ZERO, d);
}
if (Q->debug)
{
fprintf (lrs_ofp, "\n*Inequality numbers for indices %ld .. %ld : ", lastdv + 1, m + d);
for (i = 1; i <= m - nlinearity; i++)
fprintf (lrs_ofp, " %ld ", inequality[i]);
printA (D, Q);
}
if (Q->restart)
{
if (Q->debug)
fprintf (lrs_ofp, "\nPivoting to restart co-basis");
if (!restartpivots (D, Q))
return FALSE;
D->lexflag = lexmin (D, Q, ZERO); /* see if lexmin basis */
if (Q->debug)
printA (D, Q);
}
/* Check to see if necessary to resize */
if (Q->inputd > D->d)
*D_p = resize (D, Q);
return TRUE;
}
/********* end of lrs_getfirstbasis ***************/
long
getabasis2 (lrs_dic * P, lrs_dat * Q, lrs_dic * P2orig, long order[])
/* Pivot Ax<=b to standard form */
/*Try to find a starting basis by pivoting in the variables x[1]..x[d] */
/*If there are any input linearities, these appear first in order[] */
/* Steps: (a) Try to pivot out basic variables using order */
/* Stop if some linearity cannot be made to leave basis */
/* (b) Permanently remove the cobasic indices of linearities */
/* (c) If some decision variable cobasic, it is a linearity, */
/* and will be removed. */
{
long i, j, k;
/* assign local variables to structures */
lrs_mp_matrix A = P->A;
long *B = P->B;
long *C = P->C;
long *Row = P->Row;
long *Col = P->Col;
long *linearity = Q->linearity;
long *redundcol = Q->redundcol;
long m, d, nlinearity;
long nredundcol = 0L; /* will be calculated here */
static long firsttime=TRUE;
static long *linindex;
m = P->m;
d = P->d;
nlinearity = Q->nlinearity;
if(firsttime)
{
firsttime = FALSE;
linindex = (long int *) calloc ((m + d + 2), sizeof (long));
}
else /* after first time we update the change in linearities from the last time, saving many pivots */
{
for(i=1;i<=m+d;i++)
linindex[i]=FALSE;
if(Q->debug)
fprintf(lrs_ofp,"\nlindex =");
for(i=0;i<nlinearity;i++)
{
linindex[d+linearity[i]]=TRUE;
if(Q->debug)
fprintf(lrs_ofp," %ld",d+linearity[i]);
}
for(i=1;i<=m;i++)
{
if(linindex[B[i]]) /* pivot out unwanted linearities */
{
k=0;
while(k<d && (linindex[C[k]] || zero (A[Row[i]][Col[k]])))
k++;
if (k < d)
{
j=i; /* note this index changes in update, cannot use i!)*/
if(C[k] > B[j]) /* decrease i or we may skip a linearity */
i--;
pivot (P, Q, j, k);
update (P, Q, &j, &k);
}
else
if(Q->debug || Q->verbose)
fprintf(lrs_ofp,"\n*Couldn't remove linearity i=%ld B[i]=%ld",i,B[i]);
/* this is not necessarily an error, eg. two identical rows/cols in payoff matrix */
}
}
goto hotstart;
}
/* standard lrs processing is done on only the first call to getabasis2 */
if (Q->debug)
{
fprintf (lrs_ofp, "\ngetabasis from inequalities given in order");
for (i = 0; i < m; i++)
fprintf (lrs_ofp, " %ld", order[i]);
}
for (j = 0; j < m; j++)
{
i = 0;
while (i <= m && B[i] != d + order[j])
i++; /* find leaving basis index i */
if (j < nlinearity && i > m) /* cannot pivot linearity to cobasis */
{
if (Q->debug)
printA (P, Q);
#ifndef LRS_QUIET
fprintf (lrs_ofp, "\nCannot find linearity in the basis");
#endif
return FALSE;
}
if (i <= m)
{ /* try to do a pivot */
k = 0;
while (C[k] <= d && zero (A[Row[i]][Col[k]]))
k++;
if (C[k] <= d)
{
pivot (P, Q, i, k);
update (P, Q, &i, &k);
}
else if (j < nlinearity)
{ /* cannot pivot linearity to cobasis */
if (zero (A[Row[i]][0]))
{
#ifndef LRS_QUIET
fprintf (lrs_ofp, "\n*Input linearity in row %ld is redundant--skipped", order[j]);
#endif
linearity[j] = 0;
}
else
{
if (Q->debug)
printA (P, Q);
if (Q->verbose)
fprintf (lrs_ofp, "\nInconsistent linearities");
return FALSE;
}
} /* end if j < nlinearity */
} /* end of if i <= m .... */
} /* end of for */
/* update linearity array to get rid of redundancies */
i = 0;
k = 0; /* counters for linearities */
while (k < nlinearity)
{
while (k < nlinearity && linearity[k] == 0)
k++;
if (k < nlinearity)
linearity[i++] = linearity[k++];
}
nlinearity = i;
/* column dependencies now can be recorded */
/* redundcol contains input column number 0..n-1 where redundancy is */
k = 0;
while (k < d && C[k] <= d)
{
if (C[k] <= d) /* decision variable still in cobasis */
redundcol[nredundcol++] = C[k] - Q->hull; /* adjust for hull indices */
k++;
}
/* now we know how many decision variables remain in problem */
Q->nredundcol = nredundcol;
Q->lastdv = d - nredundcol;
/* if not first time we continue from here after loading dictionary */
hotstart:
if (Q->debug)
{
fprintf (lrs_ofp, "\nend of first phase of getabasis2: ");
fprintf (lrs_ofp, "lastdv=%ld nredundcol=%ld", Q->lastdv, Q->nredundcol);
fprintf (lrs_ofp, "\nredundant cobases:");
for (i = 0; i < nredundcol; i++)
fprintf (lrs_ofp, " %ld", redundcol[i]);
printA (P, Q);
}
/* here we save dictionary for use next time, *before* we resize */
copy_dict(Q,P2orig,P);
/* Remove linearities from cobasis for rest of computation */
/* This is done in order so indexing is not screwed up */
for (i = 0; i < nlinearity; i++)
{ /* find cobasic index */
k = 0;
while (k < d && C[k] != linearity[i] + d)
k++;
if (k >= d)
{
if(Q->debug || Q->verbose)
fprintf (lrs_ofp, "\nCould not remove cobasic index");
/* not neccesarily an error as eg., could be repeated row/col in payoff */
}
else
{
removecobasicindex (P, Q, k);
d = P->d;
}
}
if (Q->debug && nlinearity > 0)
printA (P, Q);
/* set index value for first slack variable */
/* Check feasability */
if (Q->givenstart)
{
i = Q->lastdv + 1;
while (i <= m && !negative (A[Row[i]][0]))
i++;
if (i <= m)
fprintf (lrs_ofp, "\n*Infeasible startingcobasis - will be modified");
}
return TRUE;
} /* end of getabasis2 */
//==========================================================================
// Outputting equilibria
//==========================================================================
// This is just a modified 'prat' from lrsmp.c
void
printrat (char name[], lrs_mp Nin, lrs_mp Din) /*reduce and print Nin/Din */
{
lrs_mp Nt, Dt;
long i;
fprintf (stdout, "%s", name);
/* reduce fraction */
copy (Nt, Nin);
copy (Dt, Din);
reduce (Nt, Dt);
/* print out */
if (sign (Nin) * sign (Din) == NEG)
fprintf (stdout, "-");
// else
// fprintf (stdout, "");
fprintf (stdout, "%lu", Nt[length (Nt) - 1]);
for (i = length (Nt) - 2; i >= 1; i--)
fprintf (stdout, FORMAT, Nt[i]);
if (!(Dt[0] == 2 && Dt[1] == 1)) /* rational */
{
fprintf (stdout, "/");
fprintf (stdout, "%lu", Dt[length (Dt) - 1]);
for (i = length (Dt) - 2; i >= 1; i--)
fprintf (stdout, FORMAT, Dt[i]);
}
}
void
nashoutput(lrs_dat *Q1, lrs_mp_vector output1,
lrs_dat *Q2, lrs_mp_vector output2,
const StrategySupport &p_support)
{
std::cout << "NE";
// The -1 is because the last entry in the vector is the payoff
// of the other player
long i = 1;
GamePlayer player1 = p_support.GetGame()->GetPlayer(1);
for (int j = 1; j <= player1->NumStrategies(); j++) {
if (p_support.Contains(player1->GetStrategy(j))) {
std::cout << ",";
printrat("", output1[i++], output1[0]);
}
else {
std::cout << ",0";
}
}
i = 1;
GamePlayer player2 = p_support.GetGame()->GetPlayer(2);
for (int j = 1; j <= player2->NumStrategies(); j++) {
if (p_support.Contains(player2->GetStrategy(j))) {
std::cout << ",";
printrat("", output2[i++], output2[0]);
}
else {
std::cout << ",0";
}
}
std::cout << std::endl;
fflush(stdout);
}
syntax highlighted by Code2HTML, v. 0.9.1