]> git.sesse.net Git - wloh/commitdiff
Add a Monte Carlo simulator.
authorSteinar H. Gunderson <Steinar H. Gunderson sesse@debian.org>
Sat, 17 Mar 2012 13:38:25 +0000 (14:38 +0100)
committerSteinar H. Gunderson <Steinar H. Gunderson sesse@debian.org>
Sat, 17 Mar 2012 13:38:25 +0000 (14:38 +0100)
mcwordfeud.cpp [new file with mode: 0644]
ziggurat.cpp [new file with mode: 0644]
ziggurat.hpp [new file with mode: 0644]

diff --git a/mcwordfeud.cpp b/mcwordfeud.cpp
new file mode 100644 (file)
index 0000000..7470af3
--- /dev/null
@@ -0,0 +1,188 @@
+#include <stdio.h>
+#include <math.h>
+#include <string.h>
+#include <stdlib.h>
+
+#include <map>
+#include <vector>
+#include <string>
+#include <algorithm>
+
+#include "ziggurat.hpp"
+
+using namespace std;
+
+#define MAX_PLAYERS 16
+
+struct player {
+       int player_index;
+       int points, margin;
+};
+struct highest_ranking {
+public:
+       bool operator() (const player &a, const player &b) const {
+               if (a.points != b.points) {
+                       return a.points > b.points;
+               }
+               return a.margin > b.margin;
+       }
+};
+
+#if 0
+float draw_unit()
+{
+       return rand() * (1.0f / RAND_MAX);
+}
+
+// FIXME: replace rejection method with a better method! this one can only ever go to [-5z, +5z]
+// and is very slow.
+float draw_gaussian(float stddev)
+{
+       for ( ;; ) {
+               float z = draw_unit() * 10.0f - 5.0f;
+               float y = draw_unit();
+               if (y < exp(-z*z)) {
+                       return z * stddev;
+               }
+       }       
+}
+#else
+float draw_gaussian(float mu, float stddev)
+{
+       static bool inited = false;
+       static long unsigned seed = 123456789;
+       int kn[128];
+       float fn[128], wn[128];
+       if (!inited) {
+               r4_nor_setup(kn, fn, wn);
+               inited = true;
+       }
+
+       return r4_nor(&seed, kn, fn, wn) * stddev + mu;
+}
+#endif
+
+int main(int argc, char **argv)
+{
+       int trials = atoi(argv[1]);
+
+       int num_players;
+       if (scanf("%d", &num_players) != 1) {
+               fprintf(stderr, "Could't read number of players\n");
+               exit(1);
+       }
+
+       if (num_players > MAX_PLAYERS) {
+               fprintf(stderr, "Max %d players supported\n");
+               exit(1);
+       }
+
+       vector<string> players;
+       map<string, int> player_map;
+       float ratings[MAX_PLAYERS];
+       bool has_scores[MAX_PLAYERS][MAX_PLAYERS];
+       for (int pl1 = 0; pl1 < num_players; ++pl1) {
+               for (int pl2 = 0; pl2 < num_players; ++pl2) {
+                       has_scores[pl1][pl2] = false;
+               }
+       }
+
+       int scores[MAX_PLAYERS][MAX_PLAYERS];
+       for (int i = 0; i < num_players; ++i) {
+               char buf[256];
+               float rating;
+               int ret = scanf("%s %f", buf, &rating);
+               if (ret < 1) {
+                       fprintf(stderr, "Couldn't read player %d\n", i);
+                       exit(1);
+               }
+               if (ret != 2) {
+                       rating = 1500.0f;
+               }
+
+               players.push_back(buf);
+               player_map[buf] = i;
+               ratings[i] = rating;
+       }
+
+       for ( ;; ) {
+               char pl1[256], pl2[256];
+               int score1, score2;
+
+               if (scanf("%s %s %d %d", pl1, pl2, &score1, &score2) != 4) {
+                       break;
+               }
+
+               if (player_map.count(pl1) == 0) {
+                       fprintf(stderr, "Unknown player '%s'\n", pl1);
+                       exit(1);
+               }
+               if (player_map.count(pl2) == 0) {
+                       fprintf(stderr, "Unknown player '%s'\n", pl2);
+                       exit(1);
+               }
+
+               scores[player_map[pl1]][player_map[pl2]] = score1;
+               scores[player_map[pl2]][player_map[pl1]] = score2;
+               has_scores[player_map[pl1]][player_map[pl2]] = true;
+               has_scores[player_map[pl2]][player_map[pl1]] = true;
+       }
+
+       int placements[MAX_PLAYERS][MAX_PLAYERS];
+       for (int i = 0; i < num_players; ++i) {
+               for (int j = 0; j < num_players; ++j) {
+                       placements[i][j] = 0;
+               }
+       }
+
+       for (int i = 0; i < trials; ++i) {
+               // draw the missing matches
+               for (int pl1 = 0; pl1 < num_players; ++pl1) {
+                       for (int pl2 = pl1 + 1; pl2 < num_players; ++pl2) {
+                               if (has_scores[pl1][pl2]) {
+                                       continue;
+                               }
+
+                               float mu = ratings[pl1] - ratings[pl2];
+                               
+                               int score = lrintf(draw_gaussian(mu, 82.9f));
+                               scores[pl1][pl2] = score;
+                               scores[pl2][pl1] = -score;
+                       }
+               }
+
+               // sum up the points and margin for each player
+               player stats[MAX_PLAYERS];
+               for (int pl1 = 0; pl1 < num_players; ++pl1) {
+                       stats[pl1].player_index = pl1;
+                       stats[pl1].points = 0;
+                       stats[pl1].margin = 0;
+                       for (int pl2 = 0; pl2 < num_players; ++pl2) {
+                               if (pl1 == pl2) {
+                                       continue;
+                               }
+                               stats[pl1].margin += scores[pl1][pl2];
+                               stats[pl1].margin -= scores[pl2][pl1];
+                               if (scores[pl1][pl2] > scores[pl2][pl1]) {
+                                       stats[pl1].points += 2;
+                               } else if (scores[pl1][pl2] == scores[pl2][pl1]) {
+                                       stats[pl1].points++;
+                               }
+                       }
+               }
+
+               // order by points and then margin
+               sort(stats, stats + num_players, highest_ranking());
+               for (int j = 0; j < num_players; ++j) {
+                       ++placements[stats[j].player_index][j];
+               }
+       }
+
+       for (int i = 0; i < num_players; ++i) {
+               printf("%-15s", players[i].c_str());
+               for (int j = 0; j < num_players; ++j) {
+                       printf(" %5d", placements[i][j]);
+               }
+               printf("\n");
+       }
+}
diff --git a/ziggurat.cpp b/ziggurat.cpp
new file mode 100644 (file)
index 0000000..dbf21aa
--- /dev/null
@@ -0,0 +1,488 @@
+# include <cstdlib>
+# include <iostream>
+# include <iomanip>
+# include <cmath>
+# include <ctime>
+
+using namespace std;
+
+# include "ziggurat.hpp"
+
+//****************************************************************************80
+
+float r4_exp ( unsigned long int *jsr, int ke[256], float fe[256], 
+  float we[256] )
+
+//****************************************************************************80
+//
+//  Purpose:
+//
+//    R4_EXP returns an exponentially distributed single precision real value.
+//
+//  Discussion:
+//
+//    The underlying algorithm is the ziggurat method.
+//
+//    Before the first call to this function, the user must call R4_EXP_SETUP
+//    to determine the values of KE, FE and WE.
+//
+//  Licensing:
+//
+//    This code is distributed under the GNU LGPL license. 
+//
+//  Modified:
+//
+//    08 December 20080
+//
+//  Author:
+//
+//    John Burkardt
+//
+//  Reference:
+//
+//    George Marsaglia, Wai Wan Tsang,
+//    The Ziggurat Method for Generating Random Variables,
+//    Journal of Statistical Software,
+//    Volume 5, Number 8, October 2000, seven pages.
+//
+//  Parameters:
+//
+//    Input/output, unsigned long int *JSR, the seed.
+//
+//    Input, int KE[256], data computed by R4_EXP_SETUP.
+//
+//    Input, float FE[256], WE[256], data computed by R4_EXP_SETUP.
+//
+//    Output, float R4_EXP, an exponentially distributed random value.
+//
+{
+  int iz;
+  int jz;
+  float value;
+  float x;
+
+  jz = shr3 ( jsr );
+  iz = ( jz & 255 );
+
+  if ( abs ( jz  ) < ke[iz] )
+  {
+    value = ( float ) ( abs ( jz ) ) * we[iz];
+  }
+  else
+  {
+    for ( ; ; )
+    {
+      if ( iz == 0 )
+      {
+        value = 7.69711 - log ( r4_uni ( jsr ) );
+        break;
+      }
+
+      x = ( float ) ( abs ( jz ) ) * we[iz];
+
+      if ( fe[iz] + r4_uni ( jsr ) * ( fe[iz-1] - fe[iz] ) < exp ( - x ) )
+      {
+        value = x;
+        break;
+      }
+
+      jz = shr3 ( jsr );
+      iz = ( jz & 255 );
+
+      if ( abs ( jz ) < ke[iz] )
+      {
+        value = ( float ) ( abs ( jz ) ) * we[iz];
+        break;
+      }
+    }
+  }
+  return value;
+}
+//****************************************************************************80
+
+void r4_exp_setup ( int ke[256], float fe[256], float we[256] )
+
+//****************************************************************************80
+//
+//  Purpose:
+//
+//    R4_EXP_SETUP sets data needed by R4_EXP.
+//
+//  Licensing:
+//
+//    This code is distributed under the GNU LGPL license. 
+//
+//  Modified:
+//
+//    08 December 2008
+//
+//  Author:
+//
+//    John Burkardt
+//
+//  Reference:
+//
+//    George Marsaglia, Wai Wan Tsang,
+//    The Ziggurat Method for Generating Random Variables,
+//    Journal of Statistical Software,
+//    Volume 5, Number 8, October 2000, seven pages.
+//
+//  Parameters:
+//
+//    Output, int KE[256], data needed by R4_EXP.
+//
+//    Output, float FE[256], WE[256], data needed by R4_EXP.
+//
+{
+  double de = 7.697117470131487;
+  int i;
+  const double m2 = 2147483648.0;
+  double q;
+  double te = 7.697117470131487;
+  const double ve = 3.949659822581572E-03;
+
+  q = ve / exp ( - de );
+
+  ke[0] = ( int ) ( ( de / q ) * m2 );
+  ke[1] = 0;
+
+  we[0] = ( float ) ( q / m2 );
+  we[255] = ( float ) ( de / m2 );
+
+  fe[0] = 1.0;
+  fe[255] = ( float ) ( exp ( - de ) );
+
+  for ( i = 254; 1 <= i; i-- )
+  {
+    de = - log ( ve / de + exp ( - de ) );
+    ke[i+1] = ( int ) ( ( de / te ) * m2 );
+    te = de;
+    fe[i] = ( float ) ( exp ( - de ) );
+    we[i] = ( float ) ( de / m2 );
+  }
+  return;
+}
+//****************************************************************************80
+
+float r4_nor ( unsigned long int *jsr, int kn[128], float fn[128], 
+  float wn[128] )
+
+//****************************************************************************80
+//
+//  Purpose:
+//
+//    R4_NOR returns a normally distributed single precision real value.
+//
+//  Discussion:
+//
+//    The value returned is generated from a distribution with mean 0 and 
+//    variance 1.
+//
+//    The underlying algorithm is the ziggurat method.
+//
+//    Before the first call to this function, the user must call R4_NOR_SETUP
+//    to determine the values of KN, FN and WN.
+//
+//  Licensing:
+//
+//    This code is distributed under the GNU LGPL license. 
+//
+//  Modified:
+//
+//    08 December 2008
+//
+//  Author:
+//
+//    John Burkardt
+//
+//  Reference:
+//
+//    George Marsaglia, Wai Wan Tsang,
+//    The Ziggurat Method for Generating Random Variables,
+//    Journal of Statistical Software,
+//    Volume 5, Number 8, October 2000, seven pages.
+//
+//  Parameters:
+//
+//    Input/output, unsigned long int *JSR, the seed.
+//
+//    Input, int KN[128], data computed by R4_NOR_SETUP.
+//
+//    Input, float FN[128], WN[128], data computed by R4_NOR_SETUP.
+//
+//    Output, float R4_NOR, a normally distributed random value.
+//
+{
+  int hz;
+  int iz;
+  const float r = 3.442620;
+  float value;
+  float x;
+  float y;
+
+  hz = shr3 ( jsr );
+  iz = ( hz & 127 );
+
+  if ( abs ( hz ) < kn[iz] )
+  {
+    value = ( float ) ( hz ) * wn[iz];
+  }
+  else
+  {
+    for ( ; ; )
+    {
+      if ( iz == 0 )
+      {
+        for ( ; ; )
+        {
+          x = - 0.2904764 * log ( r4_uni ( jsr ) );
+          y = - log ( r4_uni ( jsr ) );
+          if ( x * x <= y + y );
+          {
+            break;
+          }
+        }
+
+        if ( hz <= 0 )
+        {
+          value = - r - x;
+        }
+        else
+        {
+          value = + r + x;
+        }
+        break;
+      }
+
+      x = ( float ) ( hz ) * wn[iz];
+
+      if ( fn[iz] + r4_uni ( jsr ) * ( fn[iz-1] - fn[iz] ) < exp ( - 0.5 * x * x ) )
+      {
+        value = x;
+        break;
+      }
+
+      hz = shr3 ( jsr );
+      iz = ( hz & 127 );
+
+      if ( abs ( hz ) < kn[iz] )
+      {
+        value = ( float ) ( hz ) * wn[iz];
+        break;
+      }
+    }
+  }
+
+  return value;
+}
+//****************************************************************************80
+
+void r4_nor_setup ( int kn[128], float fn[128], float wn[128] )
+
+//****************************************************************************80
+//
+//  Purpose:
+//
+//    R4_NOR_SETUP sets data needed by R4_NOR.
+//
+//  Licensing:
+//
+//    This code is distributed under the GNU LGPL license. 
+//
+//  Modified:
+//
+//    04 May 2008
+//
+//  Author:
+//
+//    John Burkardt
+//
+//  Reference:
+//
+//    George Marsaglia, Wai Wan Tsang,
+//    The Ziggurat Method for Generating Random Variables,
+//    Journal of Statistical Software,
+//    Volume 5, Number 8, October 2000, seven pages.
+//
+//  Parameters:
+//
+//    Output, int KN[128], data needed by R4_NOR.
+//
+//    Output, float FN[128], WN[128], data needed by R4_NOR.
+//
+{
+  double dn = 3.442619855899;
+  int i;
+  const double m1 = 2147483648.0;
+  double q;
+  double tn = 3.442619855899;
+  const double vn = 9.91256303526217E-03;
+
+  q = vn / exp ( - 0.5 * dn * dn );
+
+  kn[0] = ( int ) ( ( dn / q ) * m1 );
+  kn[1] = 0;
+
+  wn[0] = ( float ) ( q / m1 );
+  wn[127] = ( float ) ( dn / m1 );
+
+  fn[0] = 1.0;
+  fn[127] = ( float ) ( exp ( - 0.5 * dn * dn ) );
+
+  for ( i = 126; 1 <= i; i-- )
+  {
+    dn = sqrt ( - 2.0 * log ( vn / dn + exp ( - 0.5 * dn * dn ) ) );
+    kn[i+1] = ( int ) ( ( dn / tn ) * m1 );
+    tn = dn;
+    fn[i] = ( float ) ( exp ( - 0.5 * dn * dn ) );
+    wn[i] = ( float ) ( dn / m1 );
+  }
+  return;
+}
+//****************************************************************************80
+
+float r4_uni ( unsigned long int *jsr )
+
+//****************************************************************************80
+//
+//  Purpose:
+//
+//    R4_UNI returns a uniformly distributed real value.
+//
+//  Licensing:
+//
+//    This code is distributed under the GNU LGPL license. 
+//
+//  Modified:
+//
+//    20 May 2008
+//
+//  Author:
+//
+//    John Burkardt
+//
+//  Reference:
+//
+//    George Marsaglia, Wai Wan Tsang,
+//    The Ziggurat Method for Generating Random Variables,
+//    Journal of Statistical Software,
+//    Volume 5, Number 8, October 2000, seven pages.
+//
+//  Parameters:
+//
+//    Input/output, unsigned long int *JSR, the seed.
+//
+//    Output, float R4_UNI, a uniformly distributed random value in
+//    the range [0,1].
+//
+{
+  unsigned long int jsr_input;
+  float value;
+
+  jsr_input = *jsr;
+
+  *jsr = ( *jsr ^ ( *jsr <<   13 ) );
+  *jsr = ( *jsr ^ ( *jsr >>   17 ) );
+  *jsr = ( *jsr ^ ( *jsr <<    5 ) );
+
+  value = fmod ( 0.5 + ( float ) ( jsr_input + *jsr ) / 65536.0 / 65536.0, 1.0 );
+
+  return value;
+}
+//****************************************************************************80
+
+unsigned long int shr3 ( unsigned long int *jsr )
+
+//****************************************************************************80
+//
+//  Purpose:
+//
+//    SHR3 evaluates the SHR3 generator for integers.
+//
+//  Licensing:
+//
+//    This code is distributed under the GNU LGPL license. 
+//
+//  Modified:
+//
+//    08 December 2008
+//
+//  Author:
+//
+//    John Burkardt
+//
+//  Reference:
+//
+//    George Marsaglia, Wai Wan Tsang,
+//    The Ziggurat Method for Generating Random Variables,
+//    Journal of Statistical Software,
+//    Volume 5, Number 8, October 2000, seven pages.
+//
+//  Parameters:
+//
+//    Input/output, unsigned long int *JSR, the seed, which is updated 
+//    on each call.
+//
+//    Output, unsigned long int SHR3, the new value.
+//
+{
+  unsigned long int value;
+
+  value = *jsr;
+
+  *jsr = ( *jsr ^ ( *jsr <<   13 ) );
+  *jsr = ( *jsr ^ ( *jsr >>   17 ) );
+  *jsr = ( *jsr ^ ( *jsr <<    5 ) );
+
+  value = value + *jsr;
+
+  return value;
+}
+//****************************************************************************80
+
+void timestamp ( )
+
+//****************************************************************************80
+//
+//  Purpose:
+//
+//    TIMESTAMP prints the current YMDHMS date as a time stamp.
+//
+//  Example:
+//
+//    31 May 2001 09:45:54 AM
+//
+//  Licensing:
+//
+//    This code is distributed under the GNU LGPL license. 
+//
+//  Modified:
+//
+//    24 September 2003
+//
+//  Author:
+//
+//    John Burkardt
+//
+//  Parameters:
+//
+//    None
+//
+{
+# define TIME_SIZE 40
+
+  static char time_buffer[TIME_SIZE];
+  const struct tm *tm;
+  size_t len;
+  time_t now;
+
+  now = time ( NULL );
+  tm = localtime ( &now );
+
+  len = strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );
+
+  cout << time_buffer << "\n";
+
+  return;
+# undef TIME_SIZE
+}
diff --git a/ziggurat.hpp b/ziggurat.hpp
new file mode 100644 (file)
index 0000000..af8272c
--- /dev/null
@@ -0,0 +1,9 @@
+float r4_exp ( unsigned long int *jsr, int ke[256], float fe[256], 
+  float we[256] );
+void r4_exp_setup ( int ke[256], float fe[256], float we[256] );
+float r4_nor ( unsigned long int *jsr, int kn[128], float fn[128], 
+  float wn[128] );
+void r4_nor_setup ( int kn[128], float fn[128], float wn[128] );
+float r4_uni ( unsigned long int *jsr );
+unsigned long int shr3 ( unsigned long int *jsr );
+void timestamp ( );