Find the best quiniela bet for specific probabilities

quiniela.cpp

#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;
}

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)
{
    // 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
            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 >
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()
{
    //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();


//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);
    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( neighborInBetSpace == false)
    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 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 )
    {
    //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
    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];
    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;

}