Adjusted initial parameters for maximum prediction power; in particular,

[foosball] / foosrank.cpp

diff --git a/foosrank.cpp b/foosrank.cpp
index 1c4353b8b7749c5413296aa19b83c5bc34ab8d1b..03edcb5b4c1047572a644aa4122af70e288c32f7 100644 (file)
--- a/foosrank.cpp
+++ b/foosrank.cpp
@@ -98,7 +98,7 @@ static double fac(int x)
        return prod;
 }
 
-static void compute_opponent_rating_pdf(int k, int a, double mu2, double sigma2, double winfac, vector<pair<double, double> > &result)
+static void compute_opponent_rating_pdf(int k, int a, double mu2, double sigma2, double winfac, vector<pair<double, double> > *result)
 {
        double binomial_precompute = prodai(k, a) / fac(k-1);
        winfac /= rating_constant;
@@ -106,27 +106,31 @@ static void compute_opponent_rating_pdf(int k, int a, double mu2, double sigma2,
        int sz = (6000.0 - 0.0) / int_step_size;
        double h = (6000.0 - 0.0) / sz;
 
-       fftw_plan f1, f2, b;
-       complex<double> *func1, *func2, *res;
-
-       func1 = reinterpret_cast<complex<double> *>(fftw_malloc(sz*2*sizeof(complex<double>)));
-       func2 = reinterpret_cast<complex<double> *>(fftw_malloc(sz*2*sizeof(complex<double>)));
-       res = reinterpret_cast<complex<double> *>(fftw_malloc(sz*2*sizeof(complex<double>)));
-       f1 = fftw_plan_dft_1d(sz*2,
-               reinterpret_cast<fftw_complex*>(func1),
-               reinterpret_cast<fftw_complex*>(func1),
-               FFTW_FORWARD,
-               FFTW_MEASURE);
-       f2 = fftw_plan_dft_1d(sz*2,
-               reinterpret_cast<fftw_complex*>(func2),
-               reinterpret_cast<fftw_complex*>(func2),
-               FFTW_FORWARD,
-               FFTW_MEASURE);
-       b = fftw_plan_dft_1d(sz*2,
-               reinterpret_cast<fftw_complex*>(res),
-               reinterpret_cast<fftw_complex*>(res),
-               FFTW_BACKWARD,
-               FFTW_MEASURE);
+       static bool inited = false;
+       static fftw_plan f1, f2, b;
+       static complex<double> *func1, *func2, *res;
+
+       if (!inited) {
+               func1 = reinterpret_cast<complex<double> *>(fftw_malloc(sz*2*sizeof(complex<double>)));
+               func2 = reinterpret_cast<complex<double> *>(fftw_malloc(sz*2*sizeof(complex<double>)));
+               res = reinterpret_cast<complex<double> *>(fftw_malloc(sz*2*sizeof(complex<double>)));
+               f1 = fftw_plan_dft_1d(sz*2,
+                       reinterpret_cast<fftw_complex*>(func1),
+                       reinterpret_cast<fftw_complex*>(func1),
+                       FFTW_FORWARD,
+                       FFTW_MEASURE);
+               f2 = fftw_plan_dft_1d(sz*2,
+                       reinterpret_cast<fftw_complex*>(func2),
+                       reinterpret_cast<fftw_complex*>(func2),
+                       FFTW_FORWARD,
+                       FFTW_MEASURE);
+               b = fftw_plan_dft_1d(sz*2,
+                       reinterpret_cast<fftw_complex*>(res),
+                       reinterpret_cast<fftw_complex*>(res),
+                       FFTW_BACKWARD,
+                       FFTW_MEASURE);
+               inited = true;
+       }
        
        // start off by zero
        for (int i = 0; i < sz*2; ++i) {
@@ -154,7 +158,7 @@ static void compute_opponent_rating_pdf(int k, int a, double mu2, double sigma2,
                func2[(i - sz/2 + sz*2)%(sz*2)].real() = prob_score_real(k, a, binomial_precompute, x2*winfac);
        }
 
-       result.reserve(sz*2);
+       result->reserve(sz*2);
 
        // convolve
        fftw_execute(f1);
@@ -164,23 +168,23 @@ static void compute_opponent_rating_pdf(int k, int a, double mu2, double sigma2,
        }
        fftw_execute(b);
 
-       result.reserve(sz);
+       result->reserve(sz);
        for (int i = 0; i < sz; ++i) {
                double r1 = i*h;
-               result.push_back(make_pair(r1, abs(res[i])));
+               result->push_back(make_pair(r1, abs(res[i])));
        }
 }
 
 // normalize the curve so we know that A ~= 1
-static void normalize(vector<pair<double, double> > &curve)
+static void normalize(vector<pair<double, double> > *curve)
 {
        double peak = 0.0;
-       for (vector<pair<double, double> >::const_iterator i = curve.begin(); i != curve.end(); ++i) {
+       for (vector<pair<double, double> >::const_iterator i = curve->begin(); i != curve->end(); ++i) {
                peak = max(peak, i->second);
        }
 
        double invpeak = 1.0 / peak;
-       for (vector<pair<double, double> >::iterator i = curve.begin(); i != curve.end(); ++i) {
+       for (vector<pair<double, double> >::iterator i = curve->begin(); i != curve->end(); ++i) {
                i->second *= invpeak;
        }
 }
@@ -241,7 +245,7 @@ static void solve_matrix(double *A, double *x, double *B)
 
 // Give an OK starting estimate for the least squares, by numerical integration
 // of statistical moments.
-static void estimate_musigma(vector<pair<double, double> > &curve, double &mu_result, double &sigma_result)
+static void estimate_musigma(const vector<pair<double, double> > &curve, double *mu_result, double *sigma_result)
 {
        double h = (curve.back().first - curve.front().first) / (curve.size() - 1);
 
@@ -272,8 +276,8 @@ static void estimate_musigma(vector<pair<double, double> > &curve, double &mu_re
        ex = (h/3.0) * ex / area;
        ex2 = (h/3.0) * ex2 / area;
 
-       mu_result = ex;
-       sigma_result = sqrt(ex2 - ex * ex);
+       *mu_result = ex;
+       *sigma_result = sqrt(ex2 - ex * ex);
 }
        
 // Find best fit of the data in curves to a Gaussian pdf, based on the
@@ -283,7 +287,7 @@ static void estimate_musigma(vector<pair<double, double> > &curve, double &mu_re
 // Note that the algorithm blows up quite hard if the initial estimate is
 // not good enough. Use estimate_musigma to get a reasonable starting
 // estimate.
-static void least_squares(vector<pair<double, double> > &curve, double mu1, double sigma1, double &mu_result, double &sigma_result)
+static void least_squares(const vector<pair<double, double> > &curve, double mu1, double sigma1, double *mu_result, double *sigma_result)
 {
        double A = 1.0;
        double mu = mu1;
@@ -362,18 +366,18 @@ static void least_squares(vector<pair<double, double> > &curve, double mu1, doub
                        break;
        }
 
-       mu_result = mu;
-       sigma_result = sigma;
+       *mu_result = mu;
+       *sigma_result = sigma;
 }
 
-static void compute_new_rating(double mu1, double sigma1, double mu2, double sigma2, int score1, int score2, double &mu, double &sigma, double &probability)
+void compute_new_rating(double mu1, double sigma1, double mu2, double sigma2, int score1, int score2, double *mu, double *sigma, double *probability)
 {
        vector<pair<double, double> > curve;
 
        if (score1 > score2) {
-               compute_opponent_rating_pdf(score1, score2, mu2, sigma2, -1.0, curve);
+               compute_opponent_rating_pdf(score1, score2, mu2, sigma2, -1.0, &curve);
        } else {
-               compute_opponent_rating_pdf(score2, score1, mu2, sigma2, 1.0, curve);
+               compute_opponent_rating_pdf(score2, score1, mu2, sigma2, 1.0, &curve);
        }
 
        // multiply in the gaussian
@@ -418,25 +422,25 @@ static void compute_new_rating(double mu1, double sigma1, double mu2, double sig
                sum /= (sigma2 * sqrt(2.0 * M_PI));
 #endif
 
-               probability = sum;
+               *probability = sum;
        }
 
        double mu_est, sigma_est;
-       normalize(curve);
-       estimate_musigma(curve, mu_est, sigma_est);
+       normalize(&curve);
+       estimate_musigma(curve, &mu_est, &sigma_est);
        least_squares(curve, mu_est, sigma_est, mu, sigma);
 }
 
-static void compute_new_double_rating(double mu1, double sigma1, double mu2, double sigma2, double mu3, double sigma3, double mu4, double sigma4, int score1, int score2, double &mu, double &sigma, double &probability)
+static void compute_new_double_rating(double mu1, double sigma1, double mu2, double sigma2, double mu3, double sigma3, double mu4, double sigma4, int score1, int score2, double *mu, double *sigma, double *probability)
 {
        vector<pair<double, double> > curve, newcurve;
        double mu_t = mu3 + mu4;
        double sigma_t = sqrt(sigma3*sigma3 + sigma4*sigma4);
                        
        if (score1 > score2) {
-               compute_opponent_rating_pdf(score1, score2, mu_t, sigma_t, -1.0, curve);
+               compute_opponent_rating_pdf(score1, score2, mu_t, sigma_t, -1.0, &curve);
        } else {
-               compute_opponent_rating_pdf(score2, score1, mu_t, sigma_t, 1.0, curve);
+               compute_opponent_rating_pdf(score2, score1, mu_t, sigma_t, 1.0, &curve);
        }
 
        newcurve.reserve(curve.size());
@@ -513,12 +517,12 @@ static void compute_new_double_rating(double mu1, double sigma1, double mu2, dou
                sum /= (sigma_t * sqrt(2.0 * M_PI));
 #endif
 
-               probability = sum;
+               *probability = sum;
        }
 
        double mu_est, sigma_est;
-       normalize(newcurve);
-       estimate_musigma(newcurve, mu_est, sigma_est);
+       normalize(&newcurve);
+       estimate_musigma(newcurve, &mu_est, &sigma_est);
        least_squares(newcurve, mu_est, sigma_est, mu, sigma);
 }
 
@@ -543,12 +547,8 @@ int main(int argc, char **argv)
                int score1 = atoi(argv[9]);
                int score2 = atoi(argv[10]);
                double mu, sigma, probability;
-               compute_new_double_rating(mu1, sigma1, mu2, sigma2, mu3, sigma3, mu4, sigma4, score1, score2, mu, sigma, probability);
-               if (score1 > score2) {
-                       printf("%f %f %f\n", mu, sigma, probability);
-               } else {
-                       printf("%f %f %f\n", mu, sigma, probability);
-               }
+               compute_new_double_rating(mu1, sigma1, mu2, sigma2, mu3, sigma3, mu4, sigma4, score1, score2, &mu, &sigma, &probability);
+               printf("%f %f %f\n", mu, sigma, probability);
        } else if (argc > 8) {
                double mu3 = atof(argv[5]);
                double sigma3 = atof(argv[6]);
@@ -561,10 +561,10 @@ int main(int argc, char **argv)
                        double newmu1_1, newmu1_2, newmu2_1, newmu2_2;
                        double newsigma1_1, newsigma1_2, newsigma2_1, newsigma2_2;
                        double probability;
-                       compute_new_double_rating(mu1, sigma1, mu2, sigma2, mu3, sigma3, mu4, sigma4, k, i, newmu1_1, newsigma1_1, probability);
-                       compute_new_double_rating(mu2, sigma2, mu1, sigma1, mu3, sigma3, mu4, sigma4, k, i, newmu1_2, newsigma1_2, probability);
-                       compute_new_double_rating(mu3, sigma3, mu4, sigma4, mu1, sigma1, mu2, sigma2, i, k, newmu2_1, newsigma2_1, probability);
-                       compute_new_double_rating(mu4, sigma4, mu3, sigma3, mu1, sigma1, mu2, sigma2, i, k, newmu2_2, newsigma2_2, probability);
+                       compute_new_double_rating(mu1, sigma1, mu2, sigma2, mu3, sigma3, mu4, sigma4, k, i, &newmu1_1, &newsigma1_1, &probability);
+                       compute_new_double_rating(mu2, sigma2, mu1, sigma1, mu3, sigma3, mu4, sigma4, k, i, &newmu1_2, &newsigma1_2, &probability);
+                       compute_new_double_rating(mu3, sigma3, mu4, sigma4, mu1, sigma1, mu2, sigma2, i, k, &newmu2_1, &newsigma2_1, &probability);
+                       compute_new_double_rating(mu4, sigma4, mu3, sigma3, mu1, sigma1, mu2, sigma2, i, k, &newmu2_2, &newsigma2_2, &probability);
                        printf("%u-%u,%f,%+f,%+f,%+f,%+f\n",
                                k, i, probability, newmu1_1-mu1, newmu1_2-mu2,
                                newmu2_1-mu3, newmu2_2-mu4);
@@ -573,11 +573,10 @@ int main(int argc, char **argv)
                        double newmu1_1, newmu1_2, newmu2_1, newmu2_2;
                        double newsigma1_1, newsigma1_2, newsigma2_1, newsigma2_2;
                        double probability;
-                       compute_new_double_rating(mu1, sigma1, mu2, sigma2, mu3, sigma3, mu4, sigma4, k, i, newmu1_1, newsigma1_1, probability);
-                       compute_new_double_rating(mu1, sigma1, mu2, sigma2, mu3, sigma3, mu4, sigma4, i, k, newmu1_1, newsigma1_1, probability);
-                       compute_new_double_rating(mu2, sigma2, mu1, sigma1, mu3, sigma3, mu4, sigma4, i, k, newmu1_2, newsigma1_2, probability);
-                       compute_new_double_rating(mu3, sigma3, mu4, sigma4, mu1, sigma1, mu2, sigma2, k, i, newmu2_1, newsigma2_1, probability);
-                       compute_new_double_rating(mu4, sigma4, mu3, sigma3, mu1, sigma1, mu2, sigma2, k, i, newmu2_2, newsigma2_2, probability);
+                       compute_new_double_rating(mu1, sigma1, mu2, sigma2, mu3, sigma3, mu4, sigma4, i, k, &newmu1_1, &newsigma1_1, &probability);
+                       compute_new_double_rating(mu2, sigma2, mu1, sigma1, mu3, sigma3, mu4, sigma4, i, k, &newmu1_2, &newsigma1_2, &probability);
+                       compute_new_double_rating(mu3, sigma3, mu4, sigma4, mu1, sigma1, mu2, sigma2, k, i, &newmu2_1, &newsigma2_1, &probability);
+                       compute_new_double_rating(mu4, sigma4, mu3, sigma3, mu1, sigma1, mu2, sigma2, k, i, &newmu2_2, &newsigma2_2, &probability);
                        printf("%u-%u,%f,%+f,%+f,%+f,%+f\n",
                                i, k, probability, newmu1_1-mu1, newmu1_2-mu2,
                                newmu2_1-mu3, newmu2_2-mu4);
@@ -586,28 +585,24 @@ int main(int argc, char **argv)
                int score1 = atoi(argv[5]);
                int score2 = atoi(argv[6]);
                double mu, sigma, probability;
-               compute_new_rating(mu1, sigma1, mu2, sigma2, score1, score2, mu, sigma, probability);
+               compute_new_rating(mu1, sigma1, mu2, sigma2, score1, score2, &mu, &sigma, &probability);
 
-               if (score1 > score2) {
-                       printf("%f %f %f\n", mu, sigma, probability);
-               } else {
-                       printf("%f %f %f\n", mu, sigma, probability);
-               }
+               printf("%f %f %f\n", mu, sigma, probability);
        } else {
                int k = atoi(argv[5]);
 
                // assess all possible scores
                for (int i = 0; i < k; ++i) {
                        double newmu1, newmu2, newsigma1, newsigma2, probability;
-                       compute_new_rating(mu1, sigma1, mu2, sigma2, k, i, newmu1, newsigma1, probability);
-                       compute_new_rating(mu2, sigma2, mu1, sigma1, i, k, newmu2, newsigma2, probability);
+                       compute_new_rating(mu1, sigma1, mu2, sigma2, k, i, &newmu1, &newsigma1, &probability);
+                       compute_new_rating(mu2, sigma2, mu1, sigma1, i, k, &newmu2, &newsigma2, &probability);
                        printf("%u-%u,%f,%+f,%+f\n",
                                k, i, probability, newmu1-mu1, newmu2-mu2);
                }
                for (int i = k; i --> 0; ) {
                        double newmu1, newmu2, newsigma1, newsigma2, probability;
-                       compute_new_rating(mu1, sigma1, mu2, sigma2, i, k, newmu1, newsigma1, probability);
-                       compute_new_rating(mu2, sigma2, mu1, sigma1, k, i, newmu2, newsigma2, probability);
+                       compute_new_rating(mu1, sigma1, mu2, sigma2, i, k, &newmu1, &newsigma1, &probability);
+                       compute_new_rating(mu2, sigma2, mu1, sigma1, k, i, &newmu2, &newsigma2, &probability);
                        printf("%u-%u,%f,%+f,%+f\n",
                                i, k, probability, newmu1-mu1, newmu2-mu2);
                }