[Bug c++/57925] New: discrete_distribution can be improved to O(1) per sampling

[Bug c++/57925] New: discrete_distribution can be improved to O(1) per sampling

[yangzhe1990 at gmail dot com]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=57925

            Bug ID: 57925
           Summary: discrete_distribution can be improved to O(1) per
                    sampling
           Product: gcc
           Version: 4.8.1
            Status: UNCONFIRMED
          Severity: enhancement
          Priority: P3
         Component: c++
          Assignee: unassigned at gcc dot gnu.org
          Reporter: yangzhe1990 at gmail dot com

Current implementation of discrete_distribution employs a strait-forward
algorithm by partial_sum and std::lower_bound, which is an O(log n) per sample
algorithm. It can be improved to an O(1) per sampling, by using the pair and
alias method, just like what GSL did.

The (log n) factor of std::lower_bound is large even for n = 3. 

The algorithm is relatively simple, and you will definitely enjoy implementing
it. Here's my sample code. Although it's quite informal, but I hope you can be
convinced that the effort of this improvement is small but its so great to have
it in the STL.

class my_discrete_distribution {
private:
    vector<double> paired;
    vector<double> weight;
    uniform_real_distribution<double> u;
public:
    template<class InputIterator>
    my_discrete_distribution(InputIterator begin, InputIterator end)
    {
        for (; begin != end; ++begin)
            weight.push_back(*begin);
        int size = weight.size();
        vector<int> small(size);
        int small_cnt = 0;
        for (int i = 0; i < size; ++i) {
            weight[i] *= size;
            if (weight[i] <= 1)
                small[small_cnt++] = i;
        }
        paired.resize(size);
        for (int i = 0; i < size; ++i)
            paired[i] = i;
        for (int i = 0; i < size; ++i) {
            double w = weight[i];
            if (w > 1) {
                int j; int p;
                for (j = small_cnt - 1; j >= 0 && w > 1; --j) {
                    p = small[j];
                    paired[p] = i;
                    w -= (1 - weight[p]);
                }
                small[j + 1] = i;
                small_cnt = j + 2;
                weight[i] = w;
            }
        }
        u.param(uniform_real_distribution<double>(0.0, size + 0.0).param());
    }

    template<class RNG>
    int operator() (RNG &rng) {
        const double a = u(rng);
        int int_part = (int)a;
        if (a - int_part >= weight[int_part])
            return paired[int_part];
        else
            return int_part;
    }
};

[Bug libstdc++/57925] discrete_distribution can be improved to O(1) per sampling

[paolo.carlini at oracle dot com]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=57925

--- Comment #1 from Paolo Carlini <paolo.carlini at oracle dot com> ---
Ok, but since you don't seem to have a Copyright assignment on file, we can't
really use your code and we would have to, eg, adapt from scratch the GSL code
or something entirely different. Are you willing to assign the Copyright to
FSF? I can send you the questionnaire to request the official forms.

[Bug libstdc++/57925] discrete_distribution can be improved to O(1) per sampling

[yangzhe1990 at gmail dot com]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=57925

--- Comment #2 from yangzhe1990 at gmail dot com ---
(In reply to Paolo Carlini from comment #1)
> Ok, but since you don't seem to have a Copyright assignment on file, we
> can't really use your code and we would have to, eg, adapt from scratch the
> GSL code or something entirely different. Are you willing to assign the
> Copyright to FSF? I can send you the questionnaire to request the official
> forms.

Cool. Feel free to use or modify my code, or do anything on this code. It's
public domain.

BTW, boost also implemented this algorithm.

[Bug libstdc++/57925] discrete_distribution can be improved to O(1) per sampling

[paolo.carlini at oracle dot com]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=57925

Paolo Carlini <paolo.carlini at oracle dot com> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
             Status|UNCONFIRMED                 |ASSIGNED
   Last reconfirmed|                            |2013-07-22
           Assignee|unassigned at gcc dot gnu.org      |paolo.carlini at oracle dot com
   Target Milestone|---                         |4.9.0
     Ever confirmed|0                           |1

--- Comment #3 from Paolo Carlini <paolo.carlini at oracle dot com> ---
Handling this.

[Bug libstdc++/57925] discrete_distribution can be improved to O(1) per sampling

[jakub at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=57925

Jakub Jelinek <jakub at gcc dot gnu.org> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
   Target Milestone|4.9.0                       |4.9.1

--- Comment #4 from Jakub Jelinek <jakub at gcc dot gnu.org> ---
GCC 4.9.0 has been released

[Bug libstdc++/57925] discrete_distribution can be improved to O(1) per sampling

[jakub at gcc dot gnu.org]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=57925

Jakub Jelinek <jakub at gcc dot gnu.org> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
   Target Milestone|4.9.1                       |4.9.2

--- Comment #5 from Jakub Jelinek <jakub at gcc dot gnu.org> ---
GCC 4.9.1 has been released.

[Bug libstdc++/57925] discrete_distribution can be improved to O(1) per sampling

[jakub at gcc dot gnu.org]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=57925

Jakub Jelinek <jakub at gcc dot gnu.org> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
   Target Milestone|4.9.2                       |4.9.3

--- Comment #6 from Jakub Jelinek <jakub at gcc dot gnu.org> ---
GCC 4.9.2 has been released.

[Bug libstdc++/57925] discrete_distribution can be improved to O(1) per sampling

[jakub at gcc dot gnu.org]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=57925

--- Comment #7 from Jakub Jelinek <jakub at gcc dot gnu.org> ---
GCC 4.9.3 has been released.

[Bug libstdc++/57925] discrete_distribution can be improved to O(1) per sampling

[jakub at gcc dot gnu.org]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=57925

Jakub Jelinek <jakub at gcc dot gnu.org> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
   Target Milestone|4.9.3                       |4.9.4

[Bug libstdc++/57925] discrete_distribution can be improved to O(1) per sampling

[pinskia at gcc dot gnu.org]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=57925

Andrew Pinski <pinskia at gcc dot gnu.org> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
   Target Milestone|4.9.4                       |---