discrete_distribution_ii.h

/*
 *  
 *  Copyright (c) 2015, G.A. Kohring
 *  All rights reserved.
 *  
 *  Redistribution and use in source and binary forms, with or without
 *  modification, are permitted provided that the following conditions 
 *  are met:
 *  
 *      * Redistributions of source code must retain the above copyright
 *        notice, this list of conditions and the following disclaimer.
 *  
 *      * Redistributions in binary form must reproduce the above copyright 
 *        notice, this list of conditions and the following disclaimer in the 
 *        documentation and/or other materials provided with the distribution.
 * 
 *  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 
 *  "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
 *  TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR 
 *  PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR 
 *  CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 
 *  EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 
 *  PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
 *  OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
 *  WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
 *  OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
 *  ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 *
 */

#ifndef DISCRETE_DISTRIBUTION_II_H
#define DISCRETE_DISTRIBUTION_II_H

/*
 * A class for producing random integers i in [0,n) with prbability p(i).
 *
 * This implementation uses the square histrogram method in combination with a 
 * small look-up table.
 *
 * (see: "Fast Generation of Discrete Random Variables", by G. Marsaglia, 
 *       W. W. Tsang and J. Wang, Journal of Statistical Software,
 *       July 2004, Vol. 11., Nr. 3. )
 */

#include <algorithm>
#include <cstdio>
#include <cstdlib>
#include <deque>
#include <functional>
#include <map>
#include <cmath>
#include <numeric>
#include <random>
#include <vector>

namespace rng {

template<typename _IntType = int >
class discrete_distribution_30bit {

    static_assert( std::is_integral<_IntType>::value,
                   "template argument not an integral type");


 private:

    std::vector<double> P;
    std::vector<double> V;
    std::vector<int>    K;
    int _J[256];

    int size;

 public:

    typedef _IntType result_type;

    /* 
     * The constructor takes as input iterators over probabilities. The 
     * probabilities should sum to one.
     *
     * Only the first thirty bits in the probabilities are used.
     */

    template <class InputIterator>
        discrete_distribution_30bit( InputIterator first, InputIterator last ): 
                                                          P( first, last ) {
   
        if ( P.size() == 0 ) {
            perror( "no probablities." );
            abort();
        }

        size = P.size();

        K.resize( size, 0 );
        V.resize( size, 0.0 );

        double max_val = static_cast<double>( ( (1 << 30) - 1 ) );
        double norm = 0.0;

        for ( int n = 0; n < size; ++n ) {
            P[n] = std::floor( ( P[n]*max_val ) );
            norm += P[n];
        }

        for ( int n = 0; n < size; ++n ) P[n] /= norm;

        init();
    }

    ~discrete_distribution_30bit(){}

    /* 
     * Generates the next random number in the sequence.
     *
     */
    template<typename _UniformRandomNumberGenerator>
    result_type operator()( _UniformRandomNumberGenerator& _urng ) {

        unsigned uran = _urng();

        int d = _J[ uran & 255 ];

        if ( d >= 0 ) {
            return d;
        } else {
            double U = static_cast<double>(uran)*2.328306437e-10;
            d = static_cast<int>( static_cast<double>( size )*U );
            if ( U < V[d] ) {
                return d;
            } else {
                return K[d];
            }
        }

    }


    void reset() {}

    /* 
     * Returns a vector containing the probabilities used in this
     * instance.
     *
     */
    std::vector<double> probabilities() const {
        return std::vector<double>( P.begin(), P.end() );
    }

    /* 
     * Returns the minimum random integer.
     *
     */
    result_type min() const { 
        return result_type( 0 );
    }

    /* 
     * Returns the maximum random integer.
     *
     */
    result_type max() const { 
        return result_type( size - 1 );
    }

    /* 
     * Check for equality by checking that the probabilities are the same.
     *
     */
    friend bool operator==( const discrete_distribution_30bit& v1, 
                            const discrete_distribution_30bit& v2) {
        return v1.P == v2.P;
    }


private:
    /* 
     * Initialize the look-up table and the histogram.
     *
     */
    void init() {

        std::vector<double> probs( P.begin(), P.end() );

        // ensure the probabilities are normalized
        double sum = 0.0;
        for ( int i = 0; i < size; ++i ) sum += probs[i];
        for ( int i = 0; i < size; ++i ) probs[i] /= sum;


        // fill in the look-up table
        int L = 0;
        double p_sum = 0.0;
        double  a = 1.0/static_cast<double>( size );
        for ( int i = 0; i < size; ++i ) {
            int k = static_cast<int>( 256.0*probs[i] );
            probs[i] = 256.0*probs[i] - static_cast<double>( k );
            p_sum += probs[i];
            for ( int j = 0; j < k; ++j ) _J[L+j] = i;
            L += k;
        }

        // fill empty table slots with -1
        for( int i = L; i < 256; ++i ) _J[i] = -1;

/*
        fprintf(stderr,"Table is %5.2f percent full\n",
                                    100.0*static_cast<double>(L)/256.0 );
*/


        // Initialize K, V
        for ( int j = 0; j < size; ++j ) {
            K[j] = j;
            V[j] = a*static_cast<double>(j + 1);
        }
        
        // normalize new probs

        double p_norm = 1.0/p_sum;
        for ( int i = 0; i < size; ++i ) probs[i] *= p_norm;

        std::multimap<double,int> mp;
        for ( int i = 0; i < size; ++i ) {
            mp.insert( std::pair<double,int>( probs[i], i ) );
        }

        // Apply Robin Hood rule

        for ( int i = 1; i <= size; ++i ) {
            std::multimap<double,int>::iterator imin = mp.begin();
            std::multimap<double,int>::iterator imax = std::prev( mp.end() );

            double min_p = (*imin).first;
            int min_i = (*imin).second;

            double max_p = (*imax).first;
            int max_i = (*imax).second;

            V[min_i] = min_p + static_cast<double>(min_i)*a;
            K[min_i] = max_i;

            mp.erase( std::prev( mp.end() ) );
            mp.erase( mp.begin() );

            max_p = (max_p + min_p) - a;

            mp.insert( std::pair<double,int>( a, min_i ) );
            mp.insert( std::pair<double,int>( max_p, max_i ) );

        }

    }

};

template<typename _IntType>
inline bool operator!=( const rng::discrete_distribution_30bit<_IntType>& v1, 
                        const rng::discrete_distribution_30bit<_IntType>& v2) {
        return !(v1.P == v2.P);
}

};  // namespace rng

#endif // DISCRETE_DISTRIBUTION_II_H