3
\$\begingroup\$

Am simply wondering if I made any egregious mistakes while implementing the Alias Method as an IEnumerator<TElement>; would also like to know if there are any sorts of general improvements that can be made to the design of the class.

Usage:

var seed = SecureRandom.Default.NextUInt32();
var rng = Pcg32XshRr.New(seed, 1U);
var generator = ProbabilisticEnumerator
    .New(
        elementWeights: new Dictionary<string, int> {
            { "A", 1 },
            { "B", 2 },
            { "C", 4 }
        },
        randomNumberGenerator: rng
    );
var summary = generator
    .Take(100000)
    .GroupBy(item => item)
    .Select(item => new {
        Element = item.Key,
        Count = item.Count(),
    })
    .OrderBy(item => item.Element);

foreach (var item in summary) {
    Console.WriteLine($"{item.Element} | {item.Count}");
}

Code:

Nuget Package

Source Repository

/// <summary>
/// Represents an enumerator that yields elements in accordance with the rules descibed by a probability table; relies on Michael D Vose's implementation of <a href="https://en.wikipedia.org/wiki/Alias_method">Walker's Alias Method.</a>
/// </summary>
/// <typeparam name="TElement">The type of elements encapsulated by the enumerator.</typeparam>
/// <remarks>
/// Derived from https://github.com/BlueRaja/Weighted-Item-Randomizer-for-C-Sharp.
/// </remarks>
public class ProbabilisticEnumerator<TElement> : IEnumerable<TElement>, IEnumerator<TElement>
{
    private readonly struct ElementMetadata
    {
        public int ActualIndex { get; }
        public int AliasedIndex { get; }
        public int Threshold { get; }

        public ElementMetadata(int actualIndex, int aliasedIndex, int biasAreaSize) {
            ActualIndex = actualIndex;
            AliasedIndex = aliasedIndex;
            Threshold = biasAreaSize;
        }
    }

    private readonly ElementMetadata[] m_elementMetadata;
    private readonly TElement[] m_elements;
    private readonly int m_heightPerRectangle;
    private readonly IUniformlyDistributedRandomNumberGenerator m_randomNumberGenerator;

    /// <summary>
    /// Gets the next random element.
    /// </summary>
    public TElement Current {
        get {
            var elementMetadata = m_elementMetadata;
            var elements = m_elements;
            var heightPerRectangle = m_heightPerRectangle;
            var randomNumberGenerator = m_randomNumberGenerator;
            var randomHeight = randomNumberGenerator.NextInt32(0, heightPerRectangle);
            var randomMetadata = elementMetadata[randomNumberGenerator.NextInt32(0, (elementMetadata.Length - 1))];

            return ((randomHeight <= randomMetadata.Threshold) ? elements[randomMetadata.ActualIndex] : elements[randomMetadata.AliasedIndex]);
        }
    }
    /// <summary>
    /// Gets the next random element.
    /// </summary>
    object IEnumerator.Current => Current;

    /// <summary>
    /// Initializes a new instance of the <see cref="ProbabilisticEnumerator{TElement}"/> class.
    /// </summary>
    /// <param name="elementWeights">The collection of element-to-weight pairs that defines the rules for the table.</param>
    /// <param name="randomNumberGenerator">The source of random numbers that will be used to perform extract elements from the table.</param>
    private ProbabilisticEnumerator(IReadOnlyDictionary<TElement, int> elementWeights, IUniformlyDistributedRandomNumberGenerator randomNumberGenerator) {
        if (elementWeights.IsNull()) {
            throw new ArgumentNullException(paramName: nameof(elementWeights));
        }

        var count = unchecked((ulong)elementWeights.Count);
        var elements = new TElement[count];
        var index = 0;
        var totalWeight = 0UL;

        foreach (var kvp in elementWeights) {
            var element = kvp.Key;
            var weight = kvp.Value;

            if (0 > weight) {
                throw new ArgumentOutOfRangeException(actualValue: weight, message: "weight must be a positive integer", paramName: nameof(weight));
            }

            elements[index++] = element;
            totalWeight += unchecked((ulong)weight);
        }

        var gcd = BitwiseHelpers.GreatestCommonDivisor(count, totalWeight);
        var heightPerRectangle = checked((int)(totalWeight / gcd));
        var weightMultiplier = checked((int)(count / gcd));

        m_elementMetadata = InitializeMetadata(elementWeights, weightMultiplier, heightPerRectangle);
        m_elements = elements;
        m_heightPerRectangle = heightPerRectangle;
        m_randomNumberGenerator = randomNumberGenerator;
    }

    /// <summary>
    /// Releases all resources used by this <see cref="ProbabilisticEnumerator{TElement}"/> instance.
    /// </summary>
    public void Dispose() { }
    /// <summary>
    /// Returns an enumerator that yields a random element from the table.
    /// </summary>
    public IEnumerator<TElement> GetEnumerator() => this;
    /// <summary>
    /// Returns an enumerator that yields a random element from the table.
    /// </summary>
    IEnumerator IEnumerable.GetEnumerator() => GetEnumerator();
    /// <summary>
    /// Returns true.
    /// </summary>
    public bool MoveNext() => true;
    /// <summary>
    /// Throws <see cref="NotSupportedException"/>.
    /// </summary>
    public void Reset() => new NotSupportedException();

    private static ElementMetadata[] InitializeMetadata(IReadOnlyDictionary<TElement, int> elementWeights, int weightMultiplier, int heightPerRectangle) {
        var count = elementWeights.Count;
        var elementMetadata = new ElementMetadata[count];
        var index = 0;
        var stackLarge = new Stack<KeyValuePair<int, int>>();
        var stackSmall = new Stack<KeyValuePair<int, int>>();

        foreach (var kvp in elementWeights) {
            var newWeight = (kvp.Value * weightMultiplier);

            if (newWeight > heightPerRectangle) {
                stackLarge.Push(new KeyValuePair<int, int>(index++, newWeight));
            }
            else {
                stackSmall.Push(new KeyValuePair<int, int>(index++, newWeight));
            }
        }

        while (0 < stackLarge.Count) {
            var largeItem = stackLarge.Pop();
            var smallItem = stackSmall.Pop();

            largeItem = new KeyValuePair<int, int>(largeItem.Key, (largeItem.Value - (heightPerRectangle - smallItem.Value)));

            if (largeItem.Value > heightPerRectangle) {
                stackLarge.Push(largeItem);
            }
            else {
                stackSmall.Push(largeItem);
            }

            elementMetadata[--count] = new ElementMetadata(smallItem.Key, largeItem.Key, smallItem.Value);
        }

        while (0 < stackSmall.Count) {
            var smallItem = stackSmall.Pop();

            elementMetadata[--count] = new ElementMetadata(smallItem.Key, smallItem.Key, heightPerRectangle);
        }

        return elementMetadata;
    }

    /// <summary>
    /// Initializes a new instance of the <see cref="ProbabilisticEnumerator{TElement}"/> class.
    /// </summary>
    /// <param name="elementWeights">The collection of element-to-weight pairs that defines the rules for the table.</param>
    /// <param name="randomNumberGenerator">The source of random numbers that will be used to perform extract elements from the table.</param>
    public static ProbabilisticEnumerator<TElement> New(IReadOnlyDictionary<TElement, int> elementWeights, IUniformlyDistributedRandomNumberGenerator randomNumberGenerator) => new ProbabilisticEnumerator<TElement>(elementWeights, randomNumberGenerator);
}

/// <summary>
/// A collection of methods that directly or indirectly augment the <see cref="ProbabilisticEnumerator{TElement}"/> class.
/// </summary>
public static class ProbabilisticEnumerator
{
    /// <summary>
    /// Initializes a new instance of the <see cref="ProbabilisticEnumerator{TElement}"/> class.
    /// </summary>
    /// <param name="elementWeights">The collection of element-to-weight pairs that defines the rules for the table.</param>
    /// <param name="randomNumberGenerator">The source of random numbers that will be used to perform extract elements from the table.</param>
    public static ProbabilisticEnumerator<TElement> New<TElement>(IReadOnlyDictionary<TElement, int> elementWeights, IUniformlyDistributedRandomNumberGenerator randomNumberGenerator) => ProbabilisticEnumerator<TElement>.New(elementWeights, randomNumberGenerator);
}
\$\endgroup\$
  • \$\begingroup\$ Where can I find some of the missing types like BitwiseHelpers or the IUniformlyDistributedRandomNumberGenerator? Are these azure-related are they some utilities of yours? I didn't find them in your repository... \$\endgroup\$ – t3chb0t Apr 22 at 13:28
  • 1
    \$\begingroup\$ @t3chb0t Sorry, I didn't include them since they didn't seem necessary for the review. They're in the same Azure DevOps project (just different repositories), here ya go: BitwiseHelpers, UniformDistribution. \$\endgroup\$ – Kittoes0124 Apr 22 at 16:09
  • \$\begingroup\$ Thanks; for the review probably not but I'd like to try to run it and maybe find other interesting stuff - perfect, it's working ;-) \$\endgroup\$ – t3chb0t Apr 22 at 17:39
  • \$\begingroup\$ @t3chb0t Ah, NP. Just FYI, the NuGet package will include those dependencies for you but one totally understands if you don't trust it and want to cobble things together yourself. \$\endgroup\$ – Kittoes0124 Apr 22 at 17:45
2
\$\begingroup\$
  public static ulong GreatestCommonDivisor(ulong x, ulong y)
    {
      if (x == 0) { return y; }
      if (y == 0) { return x; }

      var g = ((int)CountTrailingZeros(x | y));

      x >>= ((int)CountTrailingZeros(x));

      do
      {
        y >>= ((int)CountTrailingZeros(y));

        if (x > y)
        {
          var z = x;

          x = y;
          y = z;
        }

        y = (y - x);
      } while (0 != y);

      return (x << g);
    }

This can be done a lot easier and with fewer iterations by using Euclid's algorithm:

static ulong gcd(ulong x, ulong y)
{
  if (x == 0) { return y; }
  if (y == 0) { return x; }

  ulong d = 0;

  while (x > 0)
  {
    d = x % y;
    if (d == 0)
      return y;
    x = y;
    y = d;
  }

  return 1;
}

I haven't studied your implementation of the initialization in details, but at first sight it looks a lot more complicated than the implementation provided in this paper. It seems that you're trying to avoid floating point numbers?


int count = 10000;
  Dictionary<string, long> stats = new Dictionary<string, long>
  {
    {"A", 0},
    {"B", 0},
    {"C", 0},
  };

  for (int i = 0; i < count; i++)
  {

     var rng = Pcg32XshRr.New(0, 1);


    var generator = ProbabilisticEnumerator
        .New(
            elementWeights: new Dictionary<string, int> {
        { "A", 1 },
        { "B", 2 },
        { "C", 4 }
            },
            randomNumberGenerator: rng
        )
        .Take(500);

    var summary = generator
        .GroupBy(item => item)
        .Select(item => new
        {
          Element = item.Key,
          Count = item.Count(),
        })
        .OrderBy(item => item.Element);

    foreach (var item in summary)
    {
      stats[item.Element] += item.Count;
      //Console.WriteLine($"{item.Element} | {item.Count}");
    }
  }

  Console.WriteLine();
  foreach (var entry in stats)
  {
    Console.WriteLine($"{entry.Key} : {entry.Value / count}");
  }

When I run this distribution 100000, I get an average distribution as:

A : 70
B : 169
C : 261

I initialize like this: Pcg32XshRr.New(0, 1); which caused it to start the same sequence each time, but trying with this initialization: var rng = Pcg32XshRr.New(DateTime.Now.Ticks, DateTime.Now.Ticks / 10000); it gets worse:

A : 83
B : 145
C : 270

I would expect it to be more like:

A: 71 (1 / 7) * 500
B: 142
C: 285  

Or maybe I misunderstand the concept?.


  var count = unchecked((ulong)elementWeights.Count);

This seems strange. By default a C# assemblies are compiled unchecked, so it should not be necessary, unless you compile with the check flag set? (But if I try, running it in a checked environment, Pcg32XshRr.Sample() throws an OverflowException in this line:

uint threshold = ((((uint)(-exclusiveHigh)) % exclusiveHigh));


IEnumerable<TElement>, IEnumerator<TElement>

There is rarely reasons for implementing both these interfaces, and I don't see the need in this class either. IEnumerable<T> should be sufficient and can cover most needs.


            elementWeights: new Dictionary<string, int> {
        { "A", 1 },
        { "B", 2 },
        { "C", 4 }

Requiring input data in this way may be a little cumbersome in real life, because you most often will have the data and the probabilities in separate sets, so I would take two arguments - a data set and its corresponding probabilities. You should then of course check the length etc...


I think, the one class do too much, and I would split it all into more classes in order to make it all a little more SOLID. A design could be as the below, but there may be others as well:

public interface IUniformRandomGenerator
{
  double Next(int max);
}

public interface IBiasedRandomGenerator
{
  int Next { get; }
}

public class IBiasedRandomGenerator : IRandomGenerator
{
  public BiasedRandomGenerator(IList<double> probabilities, IUniformRandomGenerator uniformGenerator)
  {
    // TODO Initialize
  }

  public int Next
  {
    get
    {
      return default;
    }
  }
}

public class BiasedRandomEnumerator<TElement> : IEnumerable<TElement>
{
  private readonly IList<TElement> m_elements;
  private readonly IBiasedRandomGenerator m_random;

  public BiasedRandomEnumerator(IList<TElement> elements, IList<double> probabilities, int seed)
  {
    ValidateInput(elements, probabilities);

    m_elements = elements;
    m_random = new BiasedRandomGenerator(probabilities, new UniformRandomGenerator(seed));
  }

  public BiasedRandomEnumerator(IList<TElement> elements, IBiasedRandomGenerator random)
  {
    m_elements = elements;
    m_random = random;
  }

  private void ValidateInput(IList<TElement> elements, IList<double> probabilities)
  {
    // TODO
  }

  public IEnumerator<TElement> GetEnumerator()
  {
    while (true)
    {
      yield return m_elements[m_random.Next];
    }
  }

  IEnumerator IEnumerable.GetEnumerator()
  {
    return GetEnumerator();
  }
}

In this way each class has only one responsibility, and the interfaces secures loose couplings between them. I have experimented a little with the names, but I won't defend them to the end of times.

\$\endgroup\$
  • \$\begingroup\$ Can you share the code you used to test the distribution; in addition to the possible bugs, I can imagine a few legitimate ways for you to have ended up with a "fixed random distribution" instead of a "random random distribution." \$\endgroup\$ – Kittoes0124 Apr 22 at 19:30
  • \$\begingroup\$ Good call on the GCD function, it is a remnant of some research that I was doing; the current implementation is far from ideal since I don't have access to a single instruction CTZ function. Just did some performance testing that I have been meaning to do and confirmed that my binary version is 5x slower than the one that simply relies on mod. \$\endgroup\$ – Kittoes0124 Apr 22 at 19:53
  • \$\begingroup\$ @Kittoes0124: See my update of the section. You're right I was only running a "random distribution". Please tell me how to do it right. \$\endgroup\$ – Henrik Hansen Apr 22 at 19:54
  • 1
    \$\begingroup\$ I made an edit that generates a random seed. Also, it looks like there is a bug in my translation since the result gets even more skewed with randomized seed; I suspect it is some sort of "off by 1" error that probably happens in the initialization. Have confirmed that the bug isn't in the Pcg32XshRr generator by using SecureRandom.New() (which is much slower, but guaranteed random). \$\endgroup\$ – Kittoes0124 Apr 22 at 20:08
  • 1
    \$\begingroup\$ Found it, randomHeight = randomNumberGenerator.NextInt32(0, heightPerRectangle) should actually be randomHeight = randomNumberGenerator.NextInt32(1, heightPerRectangle). Most references depend on a generator with an exclusive upper bound; my implementation is inclusive and I failed to properly adjust. \$\endgroup\$ – Kittoes0124 Apr 22 at 21:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.