/* ALLEQ.C - A program that computes all the Nash equilibria.in 
   2-player GSCs. Copyright (C) 2002 Federico Echenique

   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.
 */

/*
   Date: Nov. 22, 2002
   Author: Federico Echenique
   Email: fede@hss.caltech.edu
   Website:
   http://www.hss.caltech.edu/~fede/
*/

/* This program is an implementation of the Algorithm in the
   paper "Finding All Equilibria," by Federico Echenique (Caltech
   Working Social Sciences Working Paper #1147

   This is implemented for particular funcional forms, defined     
   in the C functions u1 and u2. The functional forms have          
   parameters A1, A2, B1, B2. You can change the functional forms, 
   if you need more parameters you will have to add them.           */

#include <math.h>
#include <stdio.h>
#include <time.h>

#define RAND_INT(l,h) (((int)(random() * ((double)(h)-(l)+1))) + (l))

/* Function prototypes */

double random (void);  
void rand_seed (unsigned int);
double u1(int s1,int s2);
double u2(int s1,int s2);
int BR1(int s2, int myS1);
int uBR1(int s2, int myS1);
int BR2(int s1, int myS2);
int uBR2(int s1, int myS2);
void topkis(int myS1, int myS2, int *eq);
void topkis_largest(int *largest_eq);
int CHECK(int myS1,int myS2, int small_s1, int small_s2, int * proposed_eq);
int norepetition(int h, int x, int y, int * ptr1, int * ptr2);
int norepetitionglq(int h, int x, int y, int * ptr1, int * ptr2);
int no_smaller(int h, int g, int * ptr1, int * ptr2) ;

void *malloc(size_t size);
void *realloc(void *ptr,size_t size);
void free(void *ptr);

static unsigned int SEED = 93186752;

/* The following are the parameters used in the payoff functions. The
values are set in the main function below.*/

double A1=0.5;
double A2=0.5;
double B1=0.5;
double B2=0.5;

/* The number of strategies available for each player is N. This is also set */
/* in the main function. The variables largest_eq are used as an upper bound */
/* on best-response calculations. They should start out as =N. */

int N = 100;
int largest_eq1 = 100;
int largest_eq2 = 100;

/* ---- MAIN ---- */

int main(void)
{
  int *MS1, *MS2, *MS1P, *MS2P;
  int eq[2], smallest_eq[2], proposed_eq[2];
  int S1, S2, i, numbeq, k, kp, f, g, count, change_smallest;
  double showme1, showme2;
  time_t tinic, tfinal;

  /*  Size of the game (number of strategies of each p.).*/

  printf("Enter the number of strategies: ");
  scanf("%i", &N);

  /*  Get parameter values for payoffs. */

/*          A1=(double)random(); A2= (double)random(); */
/*          B1=(double)random(); B2= (double)random();   */

     A1=1.0;
     A2=0.0;
     B1=1.0; 
     B2=0.0;

  largest_eq1 = N;
  largest_eq2 = N;


  /*  Initialize time. */

  tinic =time(NULL);

/*   Initialize the number of equilibria to 1. The smallest equilibrium */
/*   is not counted, so need to start the count at 1. */

  numbeq=1;

  /*  First get the largest equilibrium of the game. */

  topkis_largest(eq);

  largest_eq1=*eq;
  largest_eq2=*(eq+1);

/* Second get the smallest equilibrium of the game. Initialize proposed_eq  */
/* to the smallest eq. */

  topkis(0,0,eq);

  *smallest_eq = *eq; 
  *(smallest_eq + 1) = *(eq + 1);
  *proposed_eq = *eq; 
  *(proposed_eq + 1) = *(eq + 1);

/* Print out the smallest and largest equilibria. */

  showme1= (double) 100*smallest_eq[0]/N;
  showme2= (double) 100*smallest_eq[1]/N;
  printf("The smallest equilibrium is %f,%f \n",showme1,showme2);

  showme1= (double) 100*largest_eq1/N;
  showme2= (double) 100*largest_eq2/N;
  printf("The largest equilibrium is %f,%f \n",showme1,showme2);


 /*   --------  Initialize state of the algorithm ---- */

/*   (MS1, MS2) will be the "current" state of the algorith, */
/*   (MS1P, MS2P) will be the "next" state of the algorithm */

  k=1; /* k is the length of the MS1 and MS2 arrays */

 kp=0; /* k is the length of the MS1 and MS2 arrays */

/*   Allocating memory to MS1 and MS2: */

  if((MS1=malloc(sizeof(int)*k))==NULL){
    printf("Memory allocation error. \n");
    exit(0);
  }

  if((MS2=malloc(sizeof(int)*k))==NULL){
    printf("Memory allocation error. \n");
    exit(0);
  }

/*    Initial state of the algorithm is the smallest eq.: */

  *MS1 = *smallest_eq;
  *MS2 = *(smallest_eq+1);

/*    Allocating memory to MS1P and MS2P: */

  if((MS1P=malloc(sizeof(int)*1))==NULL){
    printf("Memory allocation error. \n");
    exit(0);
  }

  if((MS2P=malloc(sizeof(int)*1))==NULL){
    printf("Memory allocation error. \n");
    exit(0);
  }

/*   Now the recusion that takes a state (MS1,MS2) and gets  */
/*   a new state (MS1P,MS2) */

  while (0 < k){

    for(f=0; f < k; f++){

      S1 = *(MS1 + f) + 1;
      S2 = *(MS2 + f);

    if((S1 <= largest_eq1)&&(S2 <= largest_eq2)){

        topkis(S1, S2, eq); 

      kp++;
     
      if((MS1P=realloc(MS1P,sizeof(int)*kp))==NULL){
        printf("Memory allocation error, damn malloc! \n");
        exit(0);
      }

      if((MS2P=realloc(MS2P,sizeof(int)*kp))==NULL){
        printf("Memory allocation error, damn malloc! \n");
        exit(0);
      }

      *(MS1P + (kp-1)) = *eq;
      *(MS2P + (kp-1)) = *(eq + 1);

    }

      S1=*(MS1 + f);
      S2=*(MS2 + f) + 1;

    if((S1 <= largest_eq1)&&(S2 <= largest_eq2)){

        topkis(S1, S2, eq); 

      kp++;
     
      if((MS1P=realloc(MS1P,sizeof(int)*kp))==NULL){
        printf("Memory allocation error, damn malloc! \n");
        exit(0);
      }

      if((MS2P=realloc(MS2P,sizeof(int)*kp))==NULL){
        printf("Memory allocation error, damn malloc! \n");
        exit(0);
      }

      *(MS1P + (kp-1)) = *eq;
      *(MS2P + (kp-1)) = *(eq + 1);

    }
    }
   

/* Now we have the new state. We shall copy it back to (MS1,MS2) */
/* avoiding repetitions and leaving out non-minimal elements. */

/* The size of (MS1,MS2) is now kp, need to realloc memory. */

    k=kp;

    if (k > 0) {

      if((MS1=realloc(MS1,sizeof(int)*k))==NULL){
        printf("Memory allocation error, damn malloc! \n");
        exit(0);
      }

      if((MS2=realloc(MS2,sizeof(int)*k))==NULL){
        printf("Memory allocation error, damn malloc! \n");
        exit(0);
      }

/* Copy the new state to (MS1,MS2) avoiding repetitions. */

      count = 0;

      for (g=0; g < k; g++) {
	if (norepetition(g, *(MS1P+g), *(MS2P+g), MS1, MS2)) {

	  *(MS1 + count) = *(MS1P + g);
	  *(MS2 + count) = *(MS2P + g);
	  count++;
	}
      }

/* Copy (MS1,MS2) back to (MS1P,MS2P). */

    k = count;

      for (g=0; g < k; g++) {

	  *(MS1P + g) = *(MS1 + g);
	  *(MS2P + g) = *(MS2 + g);
	
      }

/* Finally, copy (MS1P,MS2P) to (MS1,MS2) retaining only the */
/* minimal elemements. */

      count = 0;

      for (g=0; g < k; g++) {
	if (no_smaller(k, g, MS1P, MS2P)) {
	  *(MS1 + count) = *(MS1P + g);
	  *(MS2 + count) = *(MS2P + g);
	  count++;
	}
      }

/* Now we shall check wich of these are equilibria, for that we */
/* shall first see if we can use the proposed_eq as a lower bound on */
/* the "check" calculations. If we can, then proposed_eq will be
   copied to smallest_eq.*/

      k=count;

      change_smallest = 0;
    for (g=0; g < k; g++) {

      if (( *proposed_eq < *(MS1 + g) ) && ( *(proposed_eq+1) < *(MS2 + g) ))

	change_smallest++;
    }

    if (change_smallest == k)
      {*smallest_eq = * proposed_eq;
      *(smallest_eq + 1) = *(proposed_eq + 1);
      }


    /* Now we go through all the elments in (MS1,MS2) CHECKing if
       the are equilibria */

    for (g=0; g < k; g++) {

      numbeq += CHECK(*(MS1 + g) , *(MS2 + g), *smallest_eq, *(smallest_eq+1),
		      proposed_eq);
    }

    }

    kp=0;

  }

  /* Now we are done with the algorithm. Compute how long it took and
     report the number of equilibria. */

  tfinal =time(NULL);

  printf("numb. of eq. %d \t time: %f seconds \n", numbeq, 
	 difftime(tfinal,tinic));

  return 0;
}

/*  This function checks if a proposed equilibrium (S1,S2) is an eq. */
int CHECK(int myS1, int myS2, int small_s1, int small_s2, int * proposed_eq)
{
  int br1, br2, iseq, leq1, leq2;
  iseq=0;
  leq1 = largest_eq1;
  leq2 = largest_eq2;
  largest_eq1 = myS1;
  largest_eq2 = myS2;

  br1=BR1(myS2, small_s1); br2=BR2(myS1, small_s2);

  if((br1 == myS1) && (br2 == myS2)) {
    iseq=1;
    *proposed_eq = myS1;
    *(proposed_eq + 1) = myS2;
  }

  largest_eq1 = leq1;
  largest_eq2 = leq2;

  return iseq;

}

/* This function computes the smallest equilibrium when players */
/* are restricted to strategies larger than (S1,S2) */
void topkis(int S1, int S2, int *eq)
{
  int s1_init, s2_init, br1, br2;
  s1_init = S1; 
  s2_init = S2;
  br1=BR1(s2_init, S1); 
  br2=BR2(s1_init, S2);

  while( (s1_init != br1) || (s2_init != br2) )
  {
    s1_init=br1; 
    s2_init=br2;
    br1=BR1(s2_init, s1_init); 
    br2=BR2(s1_init, s2_init);
  }

  *eq = s1_init; 
  *(eq + 1) = s2_init;
}

/*  This function computes the largest equilibrium of the game*/
void topkis_largest(int *eq)
{
  int s1, s2, br1, br2;
  s1=N; 
  s2=N;
  br1=uBR1(s2,0); 
  br2=uBR2(s1,0);

  while( (s1 != br1) || (s2 != br2) )
  {
    s1=br1; 
    s2=br2;
    br1=uBR1(s2,0); 
    br2=uBR2(s1,0);
  }
  *eq = s1; 
  *(eq+1) = s1; 
}

/* This function finds the best-response by 1 to a strategy s2 by 2 */
/*  when 1 is restricted to strategies greater than or equal to S1 */
int BR1(int mys2, int myS1)
{
  int x, curmax;
  double payoff, curmaxv;

  curmax=myS1; curmaxv=u1(myS1,mys2);

  for(x=myS1; x <= largest_eq1; x++)
  {
    payoff = u1(x,mys2);
    if(payoff > curmaxv) 
      {
      curmaxv=payoff;
      curmax=x;
      }
  }
  return curmax;
}

/* This function finds the best-response by 2 to a strategy s1 by 1 */
/*  when 2 is restricted to strategies greater than or equal to S2 */
int BR2(int mys1, int myS2)
{
  int x, curmax;
  double payoff, curmaxv;

  curmax=myS2; curmaxv=u2(mys1,myS2);

  for(x=myS2; x <= largest_eq2; x++)
  {
    payoff = u2(mys1,x);
    if(payoff > curmaxv) 
      {
      curmaxv=payoff;
      curmax=x;
      }
  }
  return curmax;
}

/*  This function finds the upper-selection best-response by 1 to a  */
/*  strategy s2 by 2 when 1 is restricted to strategies greater than  */
/*  or equal to S1 */

int uBR1(int mys2, int myS1)
{
  int x, curmax;
  double payoff, curmaxv;

  curmax = largest_eq1; 
  curmaxv = u1(largest_eq1,mys2);

  for(x=largest_eq1; myS1 <= x; x--)
  {
    payoff = u1(x,mys2);
    if(payoff > curmaxv) 
      {
      curmaxv=payoff;
      curmax=x;
      }
  }
  return curmax;
}


/*  This function finds the upper-selection best-response by 2 to a  */
/*  strategy s1 by 1 when 2 is restricted to strategies greater than  */
/*  or equal to S2 */

int uBR2(int mys1, int myS2)
{
  int x, curmax;
  double payoff, curmaxv;

  curmax = largest_eq1; 
  curmaxv = u2(mys1,largest_eq1);

  for(x=largest_eq1; myS2 <= x; x--)
  {
    payoff = u2(mys1,x);
    if(payoff > curmaxv) 
      {
      curmaxv=payoff;
      curmax=x;
      }
  }
  return curmax;
}

/* This function is the payoff function of player 1*/
double u1(int s1,int s2)
{
  double value, s1here, s2here, ls1;
  s1here= (double) 100*s1/N;
  s2here= (double) 100*s2/N;
  value= -0.01*A1*(s1here-s2here)*(s1here-s2here)+2*A2*sin(100*s1here)
+ 0.0001*( (1-A1)*s1here*(1+s2here) - 0.01*(0.5-A2)*s1here*s1here );
  return value;
}

/* This function is the payoff function of player 2*/
double u2(int s1,int s2)
{
  double value, s1here, s2here;
  s1here= (double) 100*s1/N;
  s2here= (double) 100*s2/N;
  value= -0.01*B1*(s1here-s2here)*(s1here-s2here)+2*B2*sin(100*s2here)
+ 0.0001*( (1-B1)*s1here*(1+s2here) - 0.01*(0.5-B2)*s2here*s2here );
  return value;
}


int norepetition(int h, int x, int y, int * ptr1, int * ptr2) {
  int norep = 1;
  int i;
  if (h > 0) {
    for (i=0; i < h; i++) {
      if( (*(ptr1 + i) == x) && (*(ptr2 + i) == y)) norep = 0;
    }
  }

  return norep;

}


int norepetitionglq(int h, int x, int y, int * ptr1, int * ptr2) {
  int norep = 1;
  int i;
  if (h > 0) {
    for (i=0; i < h; i++) {
      if( (*(ptr1 + i) <= x) && (*(ptr2 + i) <= y)) norep = 0;
    }
  }

  return norep;

}

int no_smaller(int h, int g, int * ptr1, int * ptr2) {
  int nosmall = 1;
  int i;
    for (i=0; i < h; i++) {
      if(  
	 ( (*(ptr1 + i) < *(ptr1 + g)) && (*(ptr2 + i) <= *(ptr2 + g)) ) ||
	 ( (*(ptr1 + i) <= *(ptr1 + g)) && (*(ptr2 + i) < *(ptr2 + g)) )
	  ) nosmall = 0;
    }

  return nosmall;

}

double random ()  
{
   static unsigned int a = 1588635695, m = 4294967291U, q = 2, r = 1117695901;
   SEED = a*(SEED % q) - r*(SEED / q);
   return ((double)SEED / (double)m);
 }

void rand_seed (unsigned int init)   {if (init != 0) SEED = init;}