title: Python Custom Zipf Number Generator Performing Poorly
body: <p>I needed a custom Zipf-like number generator because <code>numpy.random.zipf</code> function doesn't achieve what I need. Firstly, its <code>alpha</code> must be greater than <code>1.0</code> and I need an alpha of <code>0.5</code>. Secondly, its cardinality is directly related to the sample size and I need to make more samples than the cardinality, e.g. make a list of 1000 elements from a Zipfian distribution of only 6 unique values.</p>
  
  <p>@stanga posted <a href="https://stackoverflow.com/a/8788662/3310838">a great solution</a> to this. </p>
  
  <pre><code>import random 
  import bisect 
  import math 
  
  class ZipfGenerator: 
  
      def __init__(self, n, alpha): 
          # Calculate Zeta values from 1 to n: 
          tmp = [1. / (math.pow(float(i), alpha)) for i in range(1, n+1)] 
          zeta = reduce(lambda sums, x: sums + [sums[-1] + x], tmp, [0]) 
  
          # Store the translation map: 
          self.distMap = [x / zeta[-1] for x in zeta] 
  
      def next(self): 
          # Take a uniform 0-1 pseudo-random value: 
          u = random.random()  
  
          # Translate the Zipf variable: 
          return bisect.bisect(self.distMap, u) - 1
  </code></pre>
  
  <p>The <code>alpha</code> can be less than <code>1.0</code> and the sampling can be infinite for a fixed cardinality <code>n</code>. The problem is that it runs too slow.</p>
  
  <pre><code># Calculate Zeta values from 1 to n: 
  tmp = [1. / (math.pow(float(i), alpha)) for i in range(1, n+1)] 
  zeta = reduce(lambda sums, x: sums + [sums[-1] + x], tmp, [0])
  </code></pre>
  
  <p>These two lines are the culprits. When I choose <code>n=50000</code> I can generate my list in ~10 seconds. I need to execute this when <code>n=5000000</code> but it's not feasible. I don't fully understand why this is performing so slow because (I think) it has linear complexity and the floating point operations seem simple. I am using Python 2.6.6 on a good server.</p>
  
  <p>Is there an optimization I can make or a different solution altogether that meet my requirements?</p>
  
  <p><strong>EDIT</strong>: I'm updating my question with a possible solution using modifications recommended by @ev-br . I've simplified it as a subroutine that returns the entire list. @ev-br was correct to suggest changing <code>bisect</code> for <code>searchssorted</code> as the former proved to be a bottleneck as well.</p>
  
  <pre><code>def randZipf(n, alpha, numSamples): 
      # Calculate Zeta values from 1 to n: 
      tmp = numpy.power( numpy.arange(1, n+1), -alpha )
      zeta = numpy.r_[0.0, numpy.cumsum(tmp)]
      # Store the translation map: 
      distMap = [x / zeta[-1] for x in zeta]
      # Generate an array of uniform 0-1 pseudo-random values: 
      u = numpy.random.random(numSamples)
      # bisect them with distMap
      v = numpy.searchsorted(distMap, u)
      samples = [t-1 for t in v]
      return samples
  </code></pre>