Last active
March 4, 2025 18:42
-
-
Save unrealwill/d1bc68d1f5c7ee6fe72b76dc578985e9 to your computer and use it in GitHub Desktop.
Revisions
-
unrealwill revised this gist
Mar 4, 2025 . 1 changed file with 88 additions and 81 deletions.There are no files selected for viewing
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode charactersOriginal file line number Diff line number Diff line change @@ -84,6 +84,13 @@ double probaInverseSelection( int32_t result, double* probas, int32_t* flip, int return p; } bool initCombination( int32_t* values, int k) { for( int32_t i = 0 ; i < k ; i++) values[i] = i+1; return true; } //next_combination Copy pasted from https://gist.github.com/shaunlebron/2843980 int32_t next_combination(int32_t *values, int32_t k, int32_t n) { @@ -176,7 +183,6 @@ double hammingNeighborhoodComputeExpectedValue( int32_t bet, double* payout,doub //we iterate over all bets which only differ on the selected matches int32_t possibleresult = bet; //We use the fact that we can factor the sum and we use the fact that the proba for an event sum to 1 double probaresult = probaInverseSelection<NbMatch,NbResByMatch >(possibleresult, probas,indDiffMatch,hd); val += payout[hd]*probaresult; } @@ -189,50 +195,48 @@ double hammingNeighborhoodComputeExpectedValue( int32_t bet, double* payout,doub template<int32_t NbMatch, int32_t NbResByMatch > void getDigits( int32_t input, int32_t* digits ) { for( int32_t i = 0 ; i < NbMatch ; i++) { int32_t nextinput = input / NbResByMatch; digits[i] = input - nextinput*NbResByMatch ; input = nextinput; } } template<int32_t NbMatch, int32_t NbResByMatch > int32_t fromDigits( int32_t* digits ) { int32_t out = 0; for( int32_t i = 0 ; i < NbMatch ;i++) { out *= NbResByMatch ; out += digits[NbMatch-1-i]; } return out; } template<int32_t NbMatch, int32_t NbResByMatch > int32_t computeNeighbor( int32_t possibleresult, int32_t j, int32_t * indDiffMatch, int32_t hd) { int32_t digits[NbMatch]; int jj = j; getDigits<NbMatch,NbResByMatch>( possibleresult, digits); for( int kk = 0 ; kk < hd ;kk++) { int32_t nextjj = jj / (NbResByMatch-1); int32_t offsetDigit = jj - nextjj * (NbResByMatch-1); jj = nextjj; int dig = digits[indDiffMatch[kk]-1] + 1 + offsetDigit; digits[indDiffMatch[kk]-1] = dig < NbResByMatch ? dig: dig-NbResByMatch; } int32_t neighbor = fromDigits<NbMatch,NbResByMatch>( digits ); return neighbor; } int32_t main() { @@ -300,8 +304,45 @@ int32_t main() std::chrono::steady_clock::time_point begin = std::chrono::steady_clock::now(); //For each possible result Compute average number of winners by category aka hamming distance == 0, hamming distance == 1, hamming distance ==2, .., hamming_distance=maxpayouthd //int32_t allmatch = pow( nbresultbymatch , nbmatch); /* double* avgWinnerByCat = new double[ allmatch * maxpayouthd]; for( int i = 0 ; i < allmatch*maxpayouthd ;i++) { avgWinnerByCat[i] = 0.0; } double val = 0.0; for( int32_t possibleresult= 0; possibleresult < allmatch ; possibleresult++) { double probaresult = proba<nbmatch,nbresultbymatch >(possibleresult, probas); for( int32_t possibleBet = 0 ; possibleBet < allmatch ; possibleBet ++) { if( hammingdistance<nbmatch,nbresultbymatch >( possibleresult,possibleBet)< maxpayouthd ) { //cout << "here" << endl; int hd = hammingdistance<nbmatch,nbresultbymatch >( possibleresult,possibleBet); avgWinnerByCat[ possibleresult*maxpayouthd + hd] += probaresult; } } } cout << "avgWinnerByCat" << endl; for( int i = 0 ; i < allmatch*maxpayouthd ; i++) cout << "avgWinnerByCat["<<i<<"] = " << avgWinnerByCat[i] << endl; */ bool useNaive = false; bool neighborInBetSpace = false; //Example for computing the expected value of a single bet { srand(42); @@ -331,6 +372,8 @@ bool useNaive = false; begin = std::chrono::steady_clock::now(); if( neighborInBetSpace == false) if( useNaive ) for( int32_t j = 0 ; j < nbdistinctbets; j++) { @@ -344,12 +387,12 @@ bool useNaive = false; betExpectedValue[j] = hammingNeighborhoodComputeExpectedValue<nbmatch,nbresultbymatch,maxpayouthd>(j,payout,probas); } //We can also permute the order of the neighborhood but some tricks don't apply which make them slower //It is slower but allows you to have a payout which is a function of the possible result for example //when the payout is dependent on the number of winner of the bet if( neighborInBetSpace == true) //We now compute the expected value for all possible distinct bets so that we can pick the best one by doing a bruteforce max if( useNaive ) { @@ -368,63 +411,27 @@ bool useNaive = false; { //Hamming Neighbor iteration int32_t indDiffMatch[maxpayouthd+1]; for( int32_t possibleresult= 0; possibleresult < allmatch ; possibleresult++) for( int32_t hd = 0 , co = pow(nbresultbymatch-1,hd) ; hd < maxpayouthd ; hd++, co = pow(nbresultbymatch-1,hd) ) for( bool firstComb = initCombination(indDiffMatch,hd);firstComb || next_combination(indDiffMatch,hd,nbmatch) ; firstComb = false ) for (int32_t j = 0 ; j < co ;j++ ) { double probaresult = proba<nbmatch,nbresultbymatch>(possibleresult, probas); double scatteredValue = payout[ hd ] * probaresult ; int32_t neighbor = computeNeighbor<nbmatch,nbresultbymatch>( possibleresult, j, indDiffMatch, hd); betExpectedValue[neighbor] += scatteredValue ; } } //TODO: dynamic programming solution //Split betExpectedValue by number of matches, and hamming distances //Then partition the combination to expose Pascal Triangle formula int32_t bestbet = 0; double bestvalue = betExpectedValue[0]; -
Mar 1, 2025 .There are no files selected for viewing
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode charactersOriginal file line number Diff line number Diff line change @@ -0,0 +1,445 @@ #include <iostream> #include <math.h> #include <random> #include <chrono> using namespace std; //Find the best quiniela bet //Compile with g++ -O3 quiniela.cpp -o quiniela //On my machine with g++ version 9.4.0 the code find the best bet in 155 s //When compiled with g++ version 10.5.0 the code is much slower and need 390 s //Single threaded, and no simd optimization template<int32_t NbMatch, int32_t NbResByMatch > int32_t hammingdistance( int32_t result, int32_t pred) { int32_t sum = 0; //We compare in base NbResByMatch for( int32_t i = 0 ; i < NbMatch ; i++) { int32_t nextresult = result / NbResByMatch ; int32_t nextpred = pred / NbResByMatch ; int32_t remainder1 = result - NbResByMatch *nextresult; int32_t remainder2 = pred - NbResByMatch *nextpred; sum+= (remainder1 != remainder2); result = nextresult; pred = nextpred; } return sum; } template<int32_t NbMatch, int32_t NbResByMatch > double proba( int32_t result, double* probas) { double p = 1.0; for( int32_t i = 0 ; i < NbMatch ; i++) { int32_t nextresult = result/ NbResByMatch ; int32_t remainder1 = result - nextresult*NbResByMatch ; p *= probas[NbResByMatch *i+remainder1]; result = nextresult; } return p; } template<int32_t NbMatch, int32_t NbResByMatch > double logproba( int32_t result, double* logprobas) { double logp = 0.0; for( int32_t i = 0 ; i < NbMatch ; i++) { int32_t nextresult = result/ NbResByMatch ; int32_t remainder1 = result - nextresult*NbResByMatch ; logp += logprobas[NbResByMatch *i+remainder1]; result = nextresult; } return logp; } //flip is given in ascending order template<int32_t NbMatch, int32_t NbResByMatch > double probaInverseSelection( int32_t result, double* probas, int32_t* flip, int32_t nbflip) { double p = 1.0; int32_t indflip = 0; for( int32_t i = 0 ; i < NbMatch ; i++) { int32_t nextresult = result/ NbResByMatch ; int32_t remainder1 = result - nextresult*NbResByMatch ; if( indflip< nbflip && flip[indflip] == (i+1) ) { p*=(1- probas[NbResByMatch *i+remainder1]); indflip ++; } else { p *= probas[NbResByMatch *i+remainder1]; } result = nextresult; } return p; } //next_combination Copy pasted from https://gist.github.com/shaunlebron/2843980 int32_t next_combination(int32_t *values, int32_t k, int32_t n) { // identify the rightmost index that can be shifted to the right // e.g. 11000111 // ^ int32_t i = k-1; if (values[i] == n) { i--; while (i >= 0 && values[i] == values[i+1]-1) i--; } // exit if no more combinations can be made // e.g. 00011111 if (i==-1) return 0; // shift chosen index to the right // e.g. // (before: 11000111) // ^ // (after: 10100111) // ^ values[i]++; // left-collapse any following indexes // (before: 10100111) // ^ *** // (after: 10111100) // ^*** i++; while (i < k) { values[i] = values[i-1]+1; i++; } return 1; } template<int32_t NbMatch, int32_t NbResByMatch > double naiveComputeExpectedValue( int32_t bet, double* payout,double* probas) { int32_t allmatch = pow( NbResByMatch , NbMatch); double val = 0.0; for( int32_t possibleresult= 0; possibleresult < allmatch ; possibleresult++) { double probaresult = proba<NbMatch,NbResByMatch >(possibleresult, probas); val += payout[ hammingdistance<NbMatch,NbResByMatch >( possibleresult,bet) ] * probaresult ; } return val; } //Slower than doing the multiplication when NbMatch is small template<int32_t NbMatch, int32_t NbResByMatch > double naiveComputeExpectedValueLogProba( int32_t bet, double* payout,double* logprobas) { int32_t allmatch = pow( NbResByMatch , NbMatch); double val = 0.0; for( int32_t possibleresult= 0; possibleresult < allmatch ; possibleresult++) { double logprobaresult = logproba<NbMatch,NbResByMatch >(possibleresult, logprobas); val += payout[ hammingdistance<NbMatch,NbResByMatch >( possibleresult,bet) ] * exp(logprobaresult) ; } return val; } template<int32_t NbMatch, int32_t NbResByMatch , int32_t maxpayouthd> double hammingNeighborhoodComputeExpectedValue( int32_t bet, double* payout,double* probas) { int32_t allmatch = pow( NbResByMatch , NbMatch); double val = 0.0; //Partition on hamming distance for( int32_t hd = 0 ; hd < maxpayouthd ; hd++) { //We pick the index of the match which will be different //We iterate over n choose k only //can be done in NbMatch ! / ( (NbMatch-hd)!* hd!) iterations int32_t indDiffMatch[maxpayouthd+1]; //The convention from the copy pasted code is that indices start at 1 for( int32_t i = 0 ; i < hd ; i++) indDiffMatch[i] = i+1; do { //we iterate over all bets which only differ on the selected matches int32_t possibleresult = bet; //We use the fact that we can factor the sum and we use the fact that the proba for an event sum to 1 //so we can replace pow(NbResByMatch ,hd) by a single sum double probaresult = probaInverseSelection<NbMatch,NbResByMatch >(possibleresult, probas,indDiffMatch,hd); val += payout[hd]*probaresult; } while( next_combination(indDiffMatch,hd,NbMatch) ); } return val; } /* template<int32_t NbMatch, int32_t NbResByMatch > int32_t zeroDigits( int32_t input, int32_t* digitstozero, int32_t size ) { int32_t digits[ NbMatch ]; for( int32_t i = 0 ; i < NbMatch ; i++) { digits[i] = input % NbResByMatch ; input /= NbResByMatch ; } for( int32_t i = 0 ; i < size ; i++) digits[digitstozero[i]] = 0; int32_t out = 0; for( int32_t i = 0 ; i < NbMatch ;i++) { out += digits[NbMatch-1-i]; out *= NbResByMatch ; } return out; } template<int32_t NbMatch, int32_t NbResByMatch > int32_t setDigit( int32_t input, int32_t digit, int32_t val ) { int32_t digits[ NbMatch ]; for( int32_t i = 0 ; i < NbMatch ; i++) { digits[i] = input % NbResByMatch ; input /= NbResByMatch ; } digits[digit] = val; int32_t out = 0; for( int32_t i = 0 ; i < NbMatch ;i++) { out *= NbResByMatch ; out += digits[NbMatch-1-i]; } return out; } */ int32_t main() { //Only valid when pow(nbresultbymatch,nbmatch ) can be stored in a int32_t const int32_t nbmatch = 14; const int32_t nbresultbymatch = 3; const int32_t maxpayouthd = 5; //the max hamming distance is nbmatch double* payout = new double[nbmatch]; for( int32_t i = 0 ; i < nbmatch ; i++) payout[i] = 0.0; //payout structure : we win if hammind distance is less or equal than 4 payout[0] = 1.0; payout[1] = 1.0; payout[2] = 1.0; payout[3] = 1.0; payout[4] = 1.0; double* probas = new double[nbmatch*nbresultbymatch]; double* logprobas = new double[nbmatch*nbresultbymatch]; std::mt19937 generator(42); std::uniform_real_distribution<double> distribution(0.0,1.0); for( int32_t i = 0 ; i < nbmatch*nbresultbymatch ; i++) { probas[i] = distribution(generator); } //we normalize so that they sum to 1 for( int32_t i = 0 ; i < nbmatch ; i++) { double sum = 0.0; for( int32_t j = 0 ; j < nbresultbymatch ; j++) { sum += probas[nbresultbymatch*i+j]; } for( int32_t j = 0 ; j < nbresultbymatch ; j++) { probas[nbresultbymatch*i+j] /= sum; } } for( int32_t i = 0 ; i < nbmatch ;i++ ) for( int32_t j = 0 ; j < nbresultbymatch ; j++) cout << "probas["<<i << "," << j << "] = " << probas[i*nbresultbymatch+j] << endl; for( int32_t i = 0 ; i < nbmatch * nbresultbymatch ; i++) { logprobas[i] = log( probas[i]); } int32_t allmatch = pow( nbresultbymatch, nbmatch); const int32_t nbdistinctbets = pow(nbresultbymatch,nbmatch); double * betExpectedValue = new double[nbdistinctbets]; for( int32_t i = 0 ; i < nbdistinctbets; i++) betExpectedValue[i] = 0.0; std::chrono::steady_clock::time_point begin = std::chrono::steady_clock::now(); bool useNaive = false; //Example for computing the expected value of a single bet { srand(42); int32_t singlebet = rand()%nbdistinctbets; cout<< "naive for single bet " << singlebet << endl; double val = naiveComputeExpectedValue<nbmatch,nbresultbymatch>(singlebet,payout,probas); std::chrono::steady_clock::time_point end = std::chrono::steady_clock::now(); cout << "naive val " << val << endl; std::cout << "Time difference = " << std::chrono::duration_cast<std::chrono::microseconds>(end - begin).count() / 1.0e6 << "[s]" << std::endl; begin = std::chrono::steady_clock::now(); double valhamming = hammingNeighborhoodComputeExpectedValue<nbmatch,nbresultbymatch,maxpayouthd>(singlebet,payout,probas); end = std::chrono::steady_clock::now(); cout << "hamming val " << valhamming << endl; std::cout << "Time difference = " << std::chrono::duration_cast<std::chrono::microseconds>(end - begin).count() / 1.0e6 << "[s]" << std::endl; begin = std::chrono::steady_clock::now(); double valnaivelp = naiveComputeExpectedValueLogProba<nbmatch,nbresultbymatch>(singlebet,payout,logprobas); end = std::chrono::steady_clock::now(); cout << "valnaivelp " << valnaivelp << endl; std::cout << "Time difference = " << std::chrono::duration_cast<std::chrono::microseconds>(end - begin).count() / 1.0e6 << "[s]" << std::endl; } //return 0; begin = std::chrono::steady_clock::now(); if( useNaive ) for( int32_t j = 0 ; j < nbdistinctbets; j++) { if( j%10000 == 0) cout << j << " / " << nbdistinctbets << endl; betExpectedValue[j] = naiveComputeExpectedValue<nbmatch,nbresultbymatch>(j,payout,probas); } else for( int32_t j = 0 ; j < nbdistinctbets; j++) { if( j%10000 == 0) cout << j << " / " << nbdistinctbets << endl; betExpectedValue[j] = hammingNeighborhoodComputeExpectedValue<nbmatch,nbresultbymatch,maxpayouthd>(j,payout,probas); } //We can also permute the order of the two loops but some tricks don't apply which make them slower /* //We now compute the expected value for all possible distinct bets so that we can pick the best one by doing a bruteforce max if( useNaive ) { //Naive iteration for( int32_t possibleresult= 0; possibleresult < allmatch ; possibleresult++) { if( possibleresult % 1000 == 0) cout<< possibleresult << " / " << allmatch << endl; double probaresult = proba<nbmatch,nbresultbymatch>(possibleresult, probas); for( int32_t j = 0 ; j < nbdistinctbets; j++) betExpectedValue[j] += payout[ hammingdistance<nbmatch,nbresultbymatch>( possibleresult,j) ] * probaresult ; } } else { //Hamming Neighbor iteration for( int32_t possibleresult= 0; possibleresult < allmatch ; possibleresult++) { if( possibleresult % 1000 == 0) cout<< possibleresult << " / " << allmatch << endl; double probaresult = proba<nbmatch,nbresultbymatch>(possibleresult, probas); //Partition on hamming distance for( int32_t hd = 0 ; hd < maxpayouthd ; hd++) { //We pick the index of the match which will be different //We iterate over choose istinct hd in nbmatch //We pick with replacement and we will filter on hamming distance later int32_t indDiffMatch[maxpayouthd+1]; for( int32_t i = 0 ; i < hd ; i++) indDiffMatch[i] = i+1; do { //we iterate over all bets which only differ on the selected matches int32_t neighbor = possibleresult; //we zero the selected match digit in base nbresultbymatch //zeroDigits( neighbor, indDiffMatch,hd ); //We can't apply the same trick using probaInverseSelection because the values are scattered int32_t co = pow(nbresultbymatch,hd); for (int32_t j = 0 ; j < co ;j++ ) { int32_t jj = j; //TODO optimization : we should skip if the remainder == digit( possibleResult, indDiffMatch[kk] ) //This would allow skip the check on hammingDistance //But this replacement should be faster than a hammingdistance check which is already fast for( int32_t kk = 0 ; kk < hd ; kk++) { int32_t remainder = jj % nbresultbymatch; jj /= nbresultbymatch; //cout << " indDiffMatch["<<kk <<"] " << indDiffMatch[kk] << endl; //the minus 1 is the next combination permutation neighbor = setDigit<nbmatch,nbresultbymatch>( neighbor, indDiffMatch[kk] -1,remainder); } if( hammingdistance<nbmatch,nbresultbymatch>( possibleresult,neighbor) == hd) { //cout << "neighbor " << neighbor << " nbdistinctbets " << nbdistinctbets << endl; betExpectedValue[neighbor] += payout[ hd ] * probaresult ; } } }while( next_combination(indDiffMatch,hd,nbmatch)); } } } */ int32_t bestbet = 0; double bestvalue = betExpectedValue[0]; for( int32_t j = 1 ; j < nbdistinctbets; j++) { //cout << j << " : " << betExpectedValue[j] << "curbestval " << bestvalue << endl; if( betExpectedValue[j] > bestvalue) { bestbet = j; bestvalue = betExpectedValue[j]; } } std::chrono::steady_clock::time_point end = std::chrono::steady_clock::now(); cout << "bestbet " << bestbet << " best value " << bestvalue << endl; std::cout << "Time difference = " << std::chrono::duration_cast<std::chrono::microseconds>(end - begin).count() / 1.0e6 << "[s]" << std::endl; }