Cdf.h

#ifndef _GK_CDF_H
#define _GK_CDF_H

#include "Sampler.h"


namespace gk {
    
//! construction d'une fonction de repartition et choix aleatoire d'un element en fonction de la fonction de repartition. 
class Cdf
{
    std::vector<float> m_cdf;
    float m_total;
    
public:
    //! constructeur, indique le nombre d'elements (optionnel).
    Cdf( const int reserve= 0 ) 
        :
        m_cdf(),
        m_total(0.f)
    {
        if(reserve > 0)
            m_cdf.reserve(reserve);
    }
    
    //! insere un element dans la fonction de repartition.
    void insert( const float v )
    {
        assert( v!= 0.f );      // la pdf doit etre strictement positive

        // accumule les valeurs pour 'construire' la fonction de repartition.
        m_total= m_total + v;
        m_cdf.push_back(m_total);        
    }
    
#if 0
    //! choisit un element, renvoie son indice et la probabilite de l'avoir choisi.
    //! utilise une recherche lineaire.
    float sampleSlow( int& id ) const
    {
        assert(m_cdf.empty() == false);

        // choisit une valeur aleatoire entre 0 et 1
        const float r= drand48() * m_total;
        // recherche l'element i tel que cdf[i-1] < r*total < cdf[i]
        int  i;
        const int count= (int) m_cdf.size();
        for(i= 0; i < count; i++)
            if(r < m_cdf[i])
                break;
        
        if(i < count)
            id= i;
        else
            id= i -1;
        
        // retrouve la valeur associee a l'element et renvoie la probabilite de l'avoir selectionne
        if(id > 0)
            return (m_cdf[id] - m_cdf[id -1]) / m_total;
        else
            return m_cdf[0] / m_total;
    }
#endif

    //! choisit un element, renvoie son indice et la probabilite de l'avoir choisi.
    //! utilise une recherche dichotomique dans l'ensemble de valeurs strictement croissant.
    float sample( Sampler& sampler, int& id ) const
    {
        assert(m_cdf.empty() == false);

        // choisit une valeur aleatoire entre 0 et 1
        const float r= sampler.uniformFloat() * m_total;
        // recherche l'element i tel que cdf[i-1] < r*total < cdf[i]
        int p= 0;
        int q= (int) m_cdf.size();
        while(p < q)
        {
            const int m= (p+q) / 2;
            if(m_cdf[m] < r)
                p= m+1;
            else
                q= m;
        }
        assert(p >= 0 && p < (int) m_cdf.size());
        
        // retrouve la valeur associee a l'element et renvoie la probabilite de l'avoir selectionne
        id= p;
        if(p > 0) 
            return (m_cdf[p] - m_cdf[p -1]) / m_total;
        else
            return m_cdf[0] / m_total;
    }
    
    //! renvoie la probabilite d'avoir selectionne un element.
    float pdf( const int id ) const
    {
        assert(m_cdf.empty() == false);
        
        if(id == 0)
            return m_cdf[0] / m_total;
        else
            return (m_cdf[id] - m_cdf[id -1]) / m_total;
    }
    
    //! choisit un element uniformement, renvoie son indice et la probabilite de l'avoir choisi.
    float sampleUniform( Sampler& sampler, int &id ) const
    {
        assert(m_cdf.empty() == false);
        // choisit une valeur aleatoire entre 0 et n-1
        id= drand48() * ((int) m_cdf.size() -1);
        assert(id < (int) m_cdf.size());
        // renvoie la probabilite d'avoir selectionne la valeur : 1 / n
        return 1.f / (float) m_cdf.size();
    }
    
    //! renvoie la probabilite d'avoir selection un element.
    float pdfUniform( const int id ) const
    {
        assert(m_cdf.empty() == false);
        return 1.f / (float) m_cdf.size();
    }
};

}

#endif

---