/*

This is the C computer program "drawprobs.c", by Jeffrey S. Rosenthal, 2022.

See the explanation and simulations at: http://probability.ca/fdraw/

This program compares probabilities and expected values of functionals of
the FIFA World Cup group draws, using Rejection-Sampler Uniform draws,
versus using FIFA's Sequential Draws.  It was written in conjunction with
the research paper "Football Group Draw Probabilities and Corrections" by
Gareth O. Roberts and Jeffrey S. Rosenthal, 2022.

The parameters NUMREJ, NUMSEQ, and TT at the top of the file, and also the
FUNCTIONALS h() at the bottom, maybe changed by the user.

To run this program, type "cc drawprobs.c -lm" to compile it, and then
type "a.out" to run it.

Comments and questions welcome: http://probability.ca/jeff/contact.html

*/


// USER-DEFINED PARAMETERS -- MAY BE CHANGED:
#define NUMREJ 10000    // Number of Rejection Samplers; must be positive.
#define NUMSEQ NUMREJ  // Number of Sequential Samplers; must be positive.
#define TT 1    // 1 for 2022 World Cup, or 0 for 2018 World Cup.
// MAY ALSO CHANGE THE FUNCTIONALS h() AT THE *END* OF THIS FILE.



#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include <sys/time.h>
#include <sys/types.h>
#include <unistd.h>

#define VERBOSE 0
#define EUROCHECK 1
#define MAXNUMH 99
int NUMH = 0;

#define NUMPOTS 4
#define POTSIZE 8
#define MAXNAMELENGTH 20
#define NUMGROUPS 8
#define GROUPSIZE 4

#if TT == 1

char *teamname[NUMPOTS][POTSIZE] = {
  { "Qatar", "Belgium", "Brazil", "France",
                "Argentina", "England", "Portugal", "Spain" },
  { "Denmark", "Netherlands", "Germany", "Switzerland",
                "Croatia", "Mexico", "USA", "Uruguay" },
  { "Iran", "Serbia", "Japan", "Senegal",
                "Tunisia", "Poland", "SKorea", "Morocco" },
  { "EurB", "AsSA", "OcNCA", "SaudiArabia",
                "Cameroon", "Ecuador", "Canada", "Ghana" },
};

// 0=Europe, 1=SouthAmerica, 2=NCamer, 3=AFC, 4=CAF, 5=OFC (Oceanian)
// 0=EU, 1=SA, 2=NA, 3=AS, 4=AF, 5=OC
int region[NUMPOTS][POTSIZE] = {
  { 3,0,1,0,1,0,0,0 },
  { 0,0,0,0,0,2,2,1 },
  { 3,0,3,4,4,0,3,4 },
  { 0,1,2,3,4,1,2,4 },
};

int extraregion[NUMPOTS][POTSIZE] = {
  { -1,-1,-1,-1,-1,-1,-1,-1 },
  { -1,-1,-1,-1,-1,-1,-1,-1 },
  { -1,-1,-1,-1,-1,-1,-1,-1 },
  { -1,3,5,-1,-1,-1,-1,-1 },
};

#else

char *teamname[NUMPOTS][POTSIZE] = {
  { "Russia", "Germany", "Brazil", "Portugal",
  		"Argentina", "Belgium", "Poland", "France" },
  { "Spain", "Peru", "Switzerland", "England",
  		"Colombia", "Mexico", "Uruguay", "Croatia" },
  { "Denmark", "Iceland", "CostaRica", "Sweden",
  		"Tunisia", "Egypt", "Senegal", "Iran" },
  { "Serbia", "Nigeria", "Australia", "Japan",
  		"Morocco", "Panama", "KoreaRep", "SaudiArabia" },
};

int region[NUMPOTS][POTSIZE] = {
  { 0,0,1,0,1,0,0,0 },
  { 0,1,0,0,1,2,1,0 },
  { 0,0,2,0,4,4,4,3 },
  { 0,4,3,3,4,2,3,3 },
};

int extraregion[NUMPOTS][POTSIZE] = {
  { -1,-1,-1,-1,-1,-1,-1,-1 },
  { -1,-1,-1,-1,-1,-1,-1,-1 },
  { -1,-1,-1,-1,-1,-1,-1,-1 },
  { -1,-1,-1,-1,-1,-1,-1,-1 },
};

#endif

#define NUMREGIONS 6
int regmax[NUMREGIONS] = { 2,1,1,1,1,1 };
int group[NUMPOTS][POTSIZE];

int newperm[99], taken[99], numacc, attemptsum, numundo;

/* BEGIN MAIN PROGRAM. */

int main(int argc, char **argv)
{

    int i,j, itnum;
    int ifloor(), isvalid();
    void randperm(), rejectionsamp(), printgroups(), printgrouplist();
    void seedrnd(), randassign(), seqassign(), displayteam();
    double theh[MAXNUMH], hsum[MAXNUMH], hav1[MAXNUMH], hav2[MAXNUMH];
    // void h(double *);
    void h();
    struct timeval starttime, endtime;

    // Initialisations
    seedrnd();
    gettimeofday(&starttime, NULL);

    printf("\nTeam Names (Regions):\n");
    for (i=0; i<NUMPOTS; i++) {
      printf("Pot %d: ", i);
      for (j=0; j<POTSIZE; j++)
        // printf("%s(%d) ", teamname[i][j], region[i][j]);
        displayteam(i,j);
      printf("\n");
    }
    printf("\n");

    printf("Running %d rejection samplers ...\n\n", NUMREJ);
    attemptsum = 0;
    for (i=0; i<NUMH; i++)
      hsum[i] = 0.0;
    for (itnum=1; itnum<=NUMREJ; itnum++) {
      if (VERBOSE)
        printf("Rej #%d/%d: ", itnum, NUMREJ);
      rejectionsamp();
      h(theh);  // Compute all the functional values.
      for (i=0; i<NUMH; i++)
        hsum[i] = hsum[i] + theh[i];
    }
    for (i=0; i<NUMH; i++)
      hav1[i] = hsum[i] / NUMREJ;
    printf("\nAchieved %d samples in %d attempts (avr=%.2f).\n",
    		NUMREJ, attemptsum, 1.0 * attemptsum / NUMREJ);
    printf("\nRejection sampler (uniform) final groups:\n");
    printgroups();
    printgrouplist();

    printf("Running %d sequential samplers ...\n", NUMSEQ);
    for (i=0; i<NUMH; i++)
      hsum[i] = 0.0;
    for (itnum=1; itnum<=NUMSEQ; itnum++) {
      if (VERBOSE)
        printf("Seq #%d/%d: ", itnum, NUMSEQ);
      seqassign();
      h(theh);  // Compute all the functional values.
      for (i=0; i<NUMH; i++)
        hsum[i] = hsum[i] + theh[i];
    }
    printf("\nSequential sampler (biased) final groups:\n");
    printgroups();
    for (i=0; i<NUMH; i++)
      hav2[i] = hsum[i] / NUMSEQ;

    printf("Average functional values, based on %d and %d samples:\n",
						NUMREJ, NUMSEQ);
    for (i=0; i<NUMH; i++)
      printf("h%d:  %f (uniform)  versus  %f (sequential)\n",
      		i, hav1[i], hav2[i]);

    /* Finish up. */
    gettimeofday(&endtime, NULL);
    printf("\nTotal elapsed time = %f seconds.\n\n",
    		endtime.tv_sec - starttime.tv_sec
        		+ (endtime.tv_usec - starttime.tv_usec) / 1e6);
    return(0);
}

void displayteam(int iii, int jjj) {
    int ex = extraregion[iii][jjj];
    printf("%s[%d-%d](%d", teamname[iii][jjj], iii, jjj, region[iii][jjj]);
    if (ex >=0 )
      printf(",%d", ex);
    printf(") ");
}

void rejectionsamp() {
    int i, isvalid();
    void randassign();

    i = 0;
    while (1) {
      randassign();
      i++;
      if (isvalid())
	break;
    }
    if (VERBOSE)
      printf("Success after %d attempts.\n", i+1);
    attemptsum = attemptsum + i+1;
}

void printgroups() {
    int i,j,g;
    void displayteam();

    for (g=0; g<NUMGROUPS; g++) {
      printf("Group %d: ", g);
      for (i=0; i<NUMPOTS; i++) {
        for (j=0; j<POTSIZE; j++) {
          if (group[i][j]==g)
            // printf("%s[%d-%d](%d) ", teamname[i][j], i, j, region[i][j]);
            displayteam(i,j);
        }
      }
      printf("\n");
    }
    printf("\n");
}

void printgrouplist() {
    int i,j;
    for (i=0; i<NUMPOTS; i++) {
        for (j=0; j<POTSIZE; j++)
          printf("%d, ", group[i][j]);
        printf("\n");
    }
    printf("\n");
}

// https://en.wikipedia.org/wiki/Fisher%E2%80%93Yates_shuffle
void randperm(int len, int fixfirst) {
  int ii, newval, tmp, ifloor();
  for (ii=0; ii<len; ii++)
    newperm[ii] = ii;
  for (ii=len-1; ii>=1; ii--) {
        newval = fixfirst + ifloor((ii-fixfirst+1)*drand48());
        tmp = newperm[ii];
        newperm[ii] = newperm[newval];
        newperm[newval] = tmp;
  }
}

void randassign() {
  int ii, jj;
  void randperm();

  for (ii=0; ii<GROUPSIZE; ii++) {
    randperm(NUMGROUPS, (ii==0));
    for (jj=0; jj<NUMGROUPS; jj++) {
      group[ii][jj] = newperm[jj];
// printf("Team %d-%d assigned to group %d.\n", ii, jj, group[ii][jj]);
    }
  }
}

void seqassign() {
  int ii, jj;
  // , nn, okay, tryspot;
  // int isvalid(), kk, lastplace;
  // int numskips = 0;
  int placerest();

  // printf("Trying sequential assignment ... ");

  // Initialise all group choices to -1 so they don't count, etc:
  for (ii=0; ii<NUMPOTS; ii++)
    for (jj=0; jj<POTSIZE; jj++)
      group[ii][jj] = -1;
  numundo = 0;

  // Assign the host team [0,0].
  group[0][0] = 0;
  if (placerest(0)) {
      if (VERBOSE)
        printf("Success; total number of undo was %d.\n", numundo);
      ;
  } else {
      printf("Failure!\n");
      exit(1);
  }

}

int placerest (int row) {
  int gg, jj;
  int isvalid(), ifloor(), valcheck;
  int numalreadyplaced = 0;
  int placedat[NUMGROUPS];

  for (gg=0; gg<NUMGROUPS; gg++)
      placedat[gg] = 0;

  for (jj=0; jj<POTSIZE; jj++) {
    if ( (gg=group[row][jj]) >= 0) {
      numalreadyplaced++;
      placedat[gg] = 1;
    }
  }

  if (numalreadyplaced >= NUMGROUPS) {  // Finished this row.

    if (row >= GROUPSIZE-1)  // Finished the last row.
      return(1);
    else
      return( placerest(row+1) );

  } else {

    // Choose a team jj from the current row who haven't yet been placed.
    while (group[row][jj=ifloor(drand48()*POTSIZE)] >= 0)
      // printf("Try another jj. numalreadyplaced=%d.\n", numalreadyplaced);
      ;
if (VERBOSE)
  printf("On row %d, chose team jj=%d, where numal=%d, group[%d][%d]=%d.\n",
		    row, jj, row, jj, numalreadyplaced, group[row][jj]);

    // Place team [row][jj] somewhere:
    for (gg=0; gg<NUMGROUPS; gg++) {
      if (!placedat[gg]) {
	  // Try placing team [row][jj] at group gg.
//if ((row==0)&&(
if (VERBOSE)
//if (numundo>=100000)
  printf("AA. numundo=%d. Trying to place team %d-%d at group %d.\n",
  numundo, row, jj, gg);
  //sleep(1);
          group[row][jj] = gg;
	  //placedat[row][gg] = 1;
	  valcheck = isvalid();
	  if ( valcheck && (placerest(row)) ) {
//if (row==0)
if (VERBOSE)
//if (numundo>=100000)
  printf("Placed.\n");
	    return(1);
          } else {
	    // Undo placement at gg.
//if (row==0)
	    if (valcheck) {
	      numundo++;
if (VERBOSE) {
//if (numundo>100000) {
  printf("numundo=%d. Trying %s[%d-%d] at group %d, numal=%d, ", numundo, teamname[row][jj], row,
jj, gg, numalreadyplaced);
  printf("Not placed.\n");
  printgroups();
}
	    }
            group[row][jj] = -1;
	    //placedat[row][gg] = 0;
	  }
      }
    }
    // The for loop has ended, so we have failed to place the rest.
    return(0);

  }
}

int isvalid () {
  int ii, jj, regcount[NUMGROUPS][NUMREGIONS], g;
  int totassigned, complete, reg, ex;
  int europefreegroups = 0;
  int europetocome = 0;

  // Initialisations:
  totassigned = 0;
  for (ii=0; ii<NUMGROUPS; ii++)
    for (jj=0; jj<NUMREGIONS; jj++)
      regcount[ii][jj] = 0;

  // Count the group region counts so far:
  for (ii=0; ii<NUMPOTS; ii++) {
    for (jj=0; jj<POTSIZE; jj++) {
      g = group[ii][jj];
      if ( (g >= 0) && (g < NUMGROUPS) ) {
        reg = region[ii][jj];
        ex = extraregion[ii][jj];
        regcount[g][reg]++;
      	totassigned++;
	if (ex >= 0)
          regcount[g][ex]++;
      }
      if ( (g<0) && (region[ii][jj]==0) )
        europetocome++;
    }
  }
  complete = (totassigned == NUMPOTS * POTSIZE);

if ( EUROCHECK && (!complete) ) {
  // Special check for not-enough-European-teams-left:
  for (ii=0; ii<NUMGROUPS; ii++)
    if (regcount[ii][0]==0)
      europefreegroups++;
  if (europetocome < europefreegroups) {
    if (VERBOSE) {
      if ( (numundo>10000) && (totassigned<30) ) {
        printf("totassigned=%d, numundo=%d. Insufficient Europes (%d) for void groups (%d).\n", totassigned, numundo, europetocome, europefreegroups);
        printgroups();
      }
    }
    return(0);
  }
}

  // Check for invalid group/region totals.
  for (ii=0; ii<NUMGROUPS; ii++) {
    if ( complete && (regcount[ii][0]==0) ) {
// printf("Group %d has zero teams from Europe.\n", ii);
      return(0);
    }
    for (jj=0; jj<NUMREGIONS; jj++) {
      if (regcount[ii][jj] > regmax[jj]) {
// printf("isvalid: Group %d Region %d count=%d is too high.\n",
// ii, jj, regcount[ii][jj]);
        return(0);
      }
    }
  }

  // The assignment passed all the checks:
  return(1);
}


/* IFLOOR */
int ifloor(double x)  /* returns floor(x) as an integer */
{
    return((int)floor(x));
}

/* SEEDRND. */
void seedrnd()
{
    int seed;
    struct timeval tmptv;
    gettimeofday (&tmptv, (struct timezone *)NULL);
    seed = (int) (tmptv.tv_usec - 1000000 * ifloor ( ( (double) tmptv.tv_usec )
                                  / ( (double) 1000000.0 ) ) );
    srand48(seed);
}


/* FUNCTIONALS h(). */
void h(double *hval) {
  int ii, jj, regcount[NUMGROUPS][NUMREGIONS], g;
  for (ii=0; ii<NUMGROUPS; ii++)
    for (jj=0; jj<NUMREGIONS; jj++)
      regcount[ii][jj] = 0;
  for (ii=0; ii<NUMPOTS; ii++)
    for (jj=0; jj<POTSIZE; jj++)
      if ( ((g=group[ii][jj]) >= 0) && (g < NUMGROUPS) )
        regcount[ group[ii][jj] ][ region[ii][jj] ]++;

  if (TT) {  // FUNCTIONALS FOR THE 2022 WORLD CUP -- MAY BE CHANGED:
    NUMH = 10;
    hval[0] = 1.0 * (regcount[0][0] >= 1);     // European with Qatar.
    hval[1] = 1.0 * (regcount[group[0][5]][0]-1); // European with England.
    hval[2] = 1.0 * (group[0][5]==group[1][2]); // England with Germany.
    hval[3] = 1.0 * (group[1][6]==0);       // USA in Group A with Qatar.
    hval[4] = 1.0 * (group[3][6]==0);       // Canada in Group A with Qatar.
    hval[5] = 1.0 * (group[1][2]==0);       // Germany in Group A with Qatar.
    hval[6] = 1.0 * (group[3][0]==0);       // EurB in Group A with Qatar.
    hval[7] = 1.0 * (group[3][6]==group[0][1]);   // Canada with Belgium.
    hval[8] = 1.0 * (group[3][6]==group[1][4]);   // Canada with Croatia.
    hval[9] = 1.0 * (group[3][6]==group[2][7]);   // Canada with Morocco.

  } else {  // FUNCTIONALS FOR THE 2018 WORLD CUP -- MAY BE CHANGED:
    NUMH = 3;
    hval[0] = 1.0 * (regcount[0][0]-1);     // European with Russia.
    hval[1] = 1.0 * (regcount[group[1][3]][0]-1); // European with England.
    hval[2] = 1.0 * (group[1][3]==group[0][1]); // England with Germany.
  }

}