/*****************************************************************************
** banzhaf.c = code module for Banzhaf power index computations
**
** Original version: Mark A. Livingston,  24 Dec 95
** Dept. of Computer Science, Univ. of North Carolina
**
** Modification: 08 May 2000 (More efficient use of pseudo-random bits)
** Hewlett-Packard Laboratories, Palo Alto, California
**
** Modification: 05 May 2000 (Added qualified majority switch)
** Hewlett-Packard Laboratories, Palo Alto, California
**
** Modification: June 1996 (Corrected Banzhaf formulas)
** Dept. of Computer Science, Univ. of North Carolina
**
** Copyright (C) 1995--2000, Mark A. Livingston, The University of North
** Carolina at Chapel Hill, Hewlett-Packard Company
**
** Routines for reading data from file and writing output to screen, and for
** generating coalitions, determining critical members, and computing power.
*****************************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <malloc.h>
#include <limits.h>
#include <sys/types.h>
#include <string.h>
#include <time.h>
#include "banzhaf.h"

/*
** These two lines define the behavior of the program for large numbers of
** blocks.  The first line literally defines how "large" is "large," and the
** second says how many coalitions we generate (at random) for large numbers
** of blocks.  That is, how many Monte Carlo trials we use.
*/
#define MAX_BLOCKS_FOR_BRUTE_FORCE 32
#define MAX_RANDOM_COALITIONS 4290000000U

/* Internal variables to help with directing computational effort */
static int total,    /* total number of votes */
           majority, /* number of votes to win */
           max,      /* number of votes of largest block */
           yea;      /* number of votes currently in favor */
static unsigned long coalitionsGenerated;
static Block *LastBlock;

/*
** Opens file, reads data, and closes file.  Allocates and fills in an array
** of structures describing the various voting blocks.  Returns negative value
** on error, or number of blocks read.
*/
int
banzhafReadBlocks( char *fname, Block **blp )
{
    int ct, i;
    char name[ NAMELEN ], buf[ BUFLEN ];
    FILE *fp;

    /* Open file and go to the beginning */
    if( fname == NULL )
      fp = stdin;
    else {
      if( ( fp = fopen( fname, "r" ) ) == NULL ) {
	fprintf( stderr, "Unable to open data file '%s'\n", fname );
	return -1;
      }
      rewind( fp );
    }

    /* Read the number of blocks and allocate space */
    fgets( buf, BUFLEN - 1, fp );
    sscanf( buf, "%d", &ct );
    if( ct <= 0 ) {
      fprintf( stderr, "Invalid number of blocks (%d)\n", ct );
      if( fp != stdin )
	fclose( fp );
      return -1;
    }
    *blp = (Block *) malloc( ( ct + 1 ) * sizeof( Block ) );
    if( *blp == NULL ) {
      fprintf( stderr, "Unable to allocate space for %d blocks\n", ct );
      if( fp != stdin )
	fclose( fp );
      return -1;
    }

    /* read block data line-by-line */
    total = max = 0;
    for( i = 0; i < ct; i++ ) {
      /* read the line */
      if( fgets( buf, BUFLEN - 1, fp ) == NULL ) {
	fprintf( stderr, "Did not find description of block %d (EOF)\n", i );
	if( fp != stdin )
	  fclose( fp );
	return -1;
      }

      /* break up the static fields and store them */
      sscanf( buf, "%s %d %d", name, &((*blp)[i].pop), &((*blp)[i].weight) );
      (*blp)[i].name = (char *) malloc(( strlen( name )+1 ) * sizeof( char ));
      if( (*blp)[i].name == NULL ) {
	fprintf( stderr, "Unable to allocate space for name, block %d\n", i );
	if( fp != stdin )
	  fclose( fp );
	return -1;
      }
      strcpy( (*blp)[i].name, name );

      /* update agrgegate stats */
      total += (*blp)[i].weight;
      if( max < (*blp)[i].weight )
	max = (*blp)[i].weight;

      /* initialize the dynamic and computed fields */
      (*blp)[i].member = False;
      (*blp)[i].power = 0;
      (*blp)[i].ppc = 0.0;
    }

    /* more aggregate stats */
    yea = 0;
    majority = total / 2 + 1; /* integer division */
    printf( "Total = %d, majority = %d, max = %d\n", total, majority, max );
    LastBlock = &((*blp)[ct]);

    if( fp != stdin )
      fclose( fp );

    srand48( (time_t) time( NULL ));
    coalitionsGenerated = 0L;
    return ct;
}

void
banzhafSetMajority( int value )
{
    majority = value;
    printf( "Total = %d, majority = %d, max = %d\n", total, majority, max );
    return;
}

static Boolean
NextBit( void )
{
    static int first = True, pos;
    static long int bitstream;
    static Boolean ret;

    if( first ) {
      first = False;
      bitstream = mrand48( );
      pos = 47;
    }

    if(( bitstream >> ( pos )) & 0x01 )
      ret = True;
    else
      ret = False;

    if( --pos < 0 ) {
      bitstream = mrand48( );
      pos = 47;
    }

    return ret;
}

/*
** If the number of blocks is too many to generate all coalitions, then we
** generate a random coalition for the Monte Carlo approach.  The number of
** blocks that is too many is based on a guess, and is defined above.  I
** suggest that you NOT raise that number.  The algorithm is exponential in
** asymptotic running time.  The number of random coalitions to generate is
** also a guess and defined above.  That one you might feel like changing;
** the running time is linear in that variable.  It returns true if we reach
** the maximum number of random coalitions.
**
** If the number of blocks is not too many, then this routine generates the
** next coalition, similarly to counting in binary (though the number has the
** low-order bits first, in this analogy).  It returns true if all possible
** coalitions have been generated.  */
Boolean
banzhafNextCoalition( Block *bl, int ct, Boolean verbose )
{
    int i, j;
    Boolean ret;

    if( ct > MAX_BLOCKS_FOR_BRUTE_FORCE ) {
      /* too many blocks to generate all coalitions, generate a random one */
      if( verbose )
	printf( "coalition = " );
      for( i = 0; i < ct; i++ ) {
	/*
	** What we want is a string of ct bits to define a coalition.  There
	** really isn't a function to generate that, so I get a pseudo-random
	** double-precision number and test whether is is less then 0.5 or not.
	** Hopefully this isn't introducing any bias, although this is not the
	** greatest pseudo-random number generator anyway.
        **
        ** That method is rather inefficient, so I wrote a different version
        ** that would use every bit from the random number generator.  This
        ** produced rather different output, and thus must have violated the
        ** assumptions I made in my variance analysis, so I am leaving the
        ** default method as the slower but more accurate version. 01 Dec 2000
	*/
#ifdef FASTER_BUT_LESS_ACCURATE
	bl[i].member = NextBit( );
#else
	bl[i].member = ( drand48( ) < 0.5 ) ? 1 : 0;
#endif
	if( verbose )
	  printf( "%c", ( bl[i].member ) ? '1' : '0' );
      }
      if( verbose )
	printf( "\n" );
      ret = ( coalitionsGenerated < MAX_RANDOM_COALITIONS - 1 ) ? False : True;
    } else { /* generate the next coalition in sequence */

      for( i = 0; i < ct && bl[i].member; i++ ) {
	bl[i].member = False;
	yea -= bl[i].weight;
      }

      if( i == ct )
	ret = True;
      else {
	bl[i].member = True;
	yea += bl[i].weight;
	if( i == ct - 4 && !verbose ) {
	  printf( "." );
	  fflush( stdout );
	}
	ret = False;
      }

      if( verbose ) {
	for( j = 0; j < ct; j++ )
	  printf( "%d", ( bl[j].member ) ? 1 : 0 );
	printf( "\n" );
      }
    }

    coalitionsGenerated++;
    return ret;
}

/* Determine if any block is critical - that is, if the current coalition is
   a majority and has no more than TOTAL / 2 + MAX votes, where MAX is the
   weight of the largest block.  If there are more votes than that, then even
   if the largest block of all leaves the coalition, the coalition still
   carries.  Thus no blocks can be critical. */
Boolean
banzhafHasCritical( Block *bl, int ct )
{
    int yea, i;

    /* count 'yea' votes */
    for( yea = 0, i = 0; i < ct; i++ )
      if( bl[i].member )
	yea += bl[i].weight;

    /* determine if this is a majority */
    if( yea >= majority ) {
      /* it is a majority - now does it have too many? */
      if( yea <= majority + max - 1 )
	/* not too many */
	return True;
      else
	return False;
    } else
      return False;
}

void
banzhafCountCritical( Block *bl, int ct )
{
    int yea, i;

    /* count 'yea' votes */
    yea = 0;
    for( i = 0; i < ct; i++ )
      if( bl[i].member )
	yea += bl[i].weight;

    /* a block is critical if the coalition falls without it */
    for( ; bl != LastBlock; bl++ )
      if( bl->member && ( yea - bl->weight < majority ) )
	bl->power++;

    return;
}

#define M_1_SQRT_2_PI 0.39894228040143267794

/*
** Until this call is made, all we've been doing is counting the power of each
** block.  What we want to do now is convert the power points for each block
** into a probability that the block can change the outcome by changing its
** "vote."  The assumption, then, is that each person within a given block has
** the same probability of changing the block's vote with his/her individual
** vote, which is in theory true in a secret ballot, assuming that no one's
** vote (or decision to cast a vote) is affected by anyone else's vote.  That
** assumption has been the subject of some debate over the years, of course.
** 
** At any rate, what this routine does is divide the number of times the block
** was deemed to be "critical" in Banzhaf's terminology (see below) by the
** total number of coalitions.  This is the probability that the block can
** change the election.  Then we assume that the probability of each person
** within the block is the same.  So the probability of a single person
** changing the election is
**
** P( changing his/her block's vote ) * P( block changing outcome of election )
**
** which is equal to
**
**                                      # of times being critical in coalition
** P( changing his/her block's vote ) * --------------------------------------
**                                              total # of coalitions
**
** The first factor must be computed with Stirling's approximation for blocks
** with populations in the millions, such as U.S. states.  The exact formula
** is the probability that all other votes are exactly split in half.  Then the
** last voter decides the election with his/her vote.  The formula requires
** some simple combinatorial formulas, given on my Banzhaf Power Index web page
** <http://www.cs.unc.edu/~livingst/Banzhaf/>.
**
** But what we are really interested in is the ratio of this value for citizens
** of different states.  So we don't really need to divide by the total number
** of coalitions.  You might also think that we don't need to know the first
** factor, only that the value is equal for each voter within a single block.
** That isn't right, though, since the probability of changing a block's
** decision is not linearly related to the population.  Banzhaf and others
** gave brief proofs of this.
**
** Being Critical
** --------------
** Banzhaf defines power as the probability of changing the election with your
** vote.  Thus the power of a block is the probability of turning a majority
** into a minority by leaving the majority.  Thus if the block has enough votes
** to defeat the majority by leaving it, and thereby leaving the majority with
** too few votes to be a majority, it was given a power point in the above
** banzhafCountCritical() routine.
*/
void
banzhafComputeRatios( Block *bl, int ct )
{
    int i, j, tmpI;
    double minPPC = 1e+300, tmpD;
    char *tmpC;

    printf( "coalitions generated = %lu\n", coalitionsGenerated );
    for( i = 0; i < ct; i++ ) {
      /* first compute citizen's probability of changing the block's vote */
      bl[i].Pp_chng = 1.0 / sqrt( (double) bl[i].pop ) * M_1_SQRT_2_PI;
      /* now probability of block changing election's result */
      bl[i].Pb_chng = (double) bl[i].power / (double) coalitionsGenerated;
      /* now each person's power in this block */
      bl[i].ppc = bl[i].Pp_chng * bl[i].Pb_chng;

      /* Keep track of lowest power-per-capita for conversion to ratios */
      if( i == 0 || bl[i].ppc < minPPC )
	minPPC = bl[i].ppc;
    }

    /* convert power-per-capita for each block to normalized power ratio */
    for( i = 0; i < ct; i++ ) {
      bl[i].pr = bl[i].ppc / minPPC;
    }

    /*
    ** Now we want to sort the blocks into descending order of power.
    ** Given that the idea here is minimal programming effort, and given
    ** that the blocks are probably already nearly sorted for long lists
    ** like the U.S. Electoral College (which are likely to be sorted by
    ** either population or Electoral votes, and we happen to know that
    ** such a sorting will put the list in nearly the final order),
    ** bubble sort is quite appropriate here.
    */
    for( i = 0; i < ct - 1; i++ )
      for( j = ct - 1; j > i; j-- )
	if( bl[j].ppc > bl[j-1].ppc ) {
	  Swap( bl[j].name,   bl[j-1].name,   tmpC );
	  Swap( bl[j].pop,    bl[j-1].pop,    tmpI );
	  Swap( bl[j].weight, bl[j-1].weight, tmpI );
	  Swap( bl[j].ppc,    bl[j-1].ppc,    tmpD );
	  Swap( bl[j].power,  bl[j-1].power,  tmpI );
	  Swap( bl[j].pr,     bl[j-1].pr,     tmpD );
	}

    return;
}

void
banzhafPrintTable( Block *bl, int ct )
{
    int i;

    if( coalitionsGenerated > 0L ) {
      printf( "Block   \tPopulation    Weight    Power      BPR  \n" );
      printf( "--------\t------------  ------  ---------  -------\n" );
      for( i = 0; i < ct; i++ )
	printf( "%-16s%12d  %6d  %8ld   %6.3f\n",
		bl[i].name, bl[i].pop, bl[i].weight, bl[i].power, bl[i].pr );
    }

    return;
}