RNG benchmarks in Python vs R

Active Analytics Ltd: posted 23 July 2014 by Chibisi Chima-Okereke

Introduction

Today we provide a helpful script inspired by R's microbenchmark package for benchmarking in Python. We then use this script to compare random number generation (RNG) in Python vs R. You may want to view a previous blog on R vs Python.

At the moment if you want to benchmark functions in Python, there are limited tools that allow you to do this apart from basic clock time. The timeit module is available but it doesn't have the most convenient usage syntax. So after using the microbenchmark package, we were inspired to outline a similar simple script in Python.

The code in this blog entry is run/compiled on Ubuntu 14.04, core i7 processor having 16 GB RAM, using the g++ version 4.8.2 compiler, R version 3.1.0, Python version 2.7.6. The script is now located in github.

Preparations

We load a few python modules that will be used in the script:

import collections as col
import numpy as np
import time as tt
import tabulate as tab

Next we create a SimpleTable object, which is literally a simple table type object that inherits from an OrderedDict in the collections library. It is rudimentary and designed to hold the data from the benchmarks but does not have any checking or advanced features. The keys of the dict object represent the column headers and column data is held in the values:

class SimpleTable(col.OrderedDict):
    """
        SimpleTable is just a OrderedDict with some convenience functions
        It is not designed to be full featured nor does it do any checking
        of the table structure.

        The keys represent the column headers and the values represent the
        columns.
    """
    def nrow(self):
        """
        Returns the number of rows in the table
        """
        return self[list(self.viewkeys())[0]].__len__()
    def ncol(self):
        """
        Returns the number of columns in the table
        """
        return self.__len__()
    def rbind(self, row):
        """
        Binds any suitable dict, OrderedDict, or SimpleTable
        To the current table

        :row: suitable dict or SimpleTable of rows to be bound
        """
        for k in self.viewkeys():
            for i in np.arange(row[k].__len__()):
                self[k].append(row[k][i])
    def __repr__(self):
        return tab.tabulate(self, headers = "keys")

Below we have a function that formats and rounds units of time from seconds to either "ms" miliseconds, "us" microsceonds, "ns" nanoseconds, "s" seconds (rounding), or simply return the "raw" timing output from `tt.time` function.

# Function to transform and round time units
def timeUnits(times, unit = "ms", ndigits = 5):
    """
    Function formats times for the bench function
    :param times: list of times to be formatted
    :param unit: string the units of time: "ms": (milliseconds), "us": (microseconds),
                    "ns": nanoseconds, "s": time in seconds, "raw": raw time in seconds
    :param ndigits: number of decimal places to round down the times
    :return: list of formatted times
    """
    if unit == "ms":
        return [round(i * 1E3, ndigits) for i in times]
    elif unit == "us":
        return [round(i * 1E6, ndigits) for i in times]
    elif unit == "ns":
        return [round(i * 1E9, ndigits) for i in times]
    elif unit == "s":
        return [round(i, ndigits) for i in times]
    elif unit == "raw":
        return times

Benchmarking functions

Next we provide three functions for benchmarking, bench is the workhorse function that carries out the benchmark of a single case defined in a string, the lbenchmark function carries out benchmarks of a list of of strings to be evaluated, and benchmark allows multiple strings to be entered one after the other as different arguments in the function. So for common use, the lbenchmark and benchmark functions are the available - take your pick.

The bench function is given below:

def bench(sExpr, neval=100, units = "ms", ndigits = 5):
    """
    :param expr: string expression to be evaluated
    :param neval: number of evaluations that the statistics will be calculated from
    :param units: string the units of time: "ms": (milliseconds), "us": (microseconds),
                    "ns": nanoseconds, "s": time in seconds, "raw": raw time in seconds
    :param ndigits: number of decimal places to round down the times
    :return: SimpleTable of times min, lower, mid, upper quartiles and max time, and the expression run
    """
    times = np.ndarray(shape=(neval,), dtype = np.float64)
    expr = compile(sExpr, "<string>", "eval")
    for i in np.arange(neval):
        start = tt.time()
        out = eval(expr)
        end = tt.time()
        times[i] = end - start
    times = np.percentile(times, [0, 25, 50, 75, 100])
    times = timeUnits(times, units, ndigits)
    summ = SimpleTable([('expr', [sExpr]), ('min', [times[0]]), ('lq', [times[1]]), ('median', [times[2]]),
                            ('uq', [times[3]]), ('max', [times[4]])], ('neval', [neval]))
    return summ

The lbenchmark function is given below:

def lbenchmark(lExpr, **kwargs):
    """
    List version of benchmark, takes in a list of string expressions as lExpr

    :param lExpr: list of strings to be evaluated
    :param neval: number of evaluations that the statistics will be calculated from
    :param units: string the units of time: "ms": (milliseconds), "us": (microseconds),
                    "ns": nanoseconds, "s": time in seconds, "raw": raw time in seconds
    :param ndigits: number of decimal places to round down the times
    :return: SimpleTable of times min, lower, mid, upper quartiles and max time, and the expression run
    """
    nExpr = lExpr.__len__()
    out = bench(lExpr[0], **kwargs)
    if nExpr > 1:
        for i in np.arange(1, nExpr):
            out.rbind(bench(lExpr[i], **kwargs))
    return out

The benchmark function is given below:

def benchmark(*args, **kwargs):
    """
    Function running bench underneath for multiple strings

    :param args: the strings to be evaluated
    :param kwargs: arguments passed to bench
    :return: SimpleTable of times min, lower, mid, upper quartiles and max time, and the expression run
    """
    nExpr = args.__len__()
    out = bench(args[0], **kwargs)
    if nExpr > 1:
        for i in np.arange(1, nExpr):
            out.rbind(bench(args[i], **kwargs))
    return out

RNG R vs. Python

Here we give examples that benchmark a few functions for sampling from statistical distributions and compare them to their R counterparts. Both R and Python use the Mersenne Twister as the RNG engine:

Examples

In Python

# Units are microseconds ...

>>> evalStrings = ["np.random.uniform(size=1000)", "np.random.normal(size=1000)",
               "np.random.binomial(size=1000, n=10, p = .5)", "np.random.poisson(size=1000, lam=1)"]

>>> bench(evalStrings[0], units = "us")
expr                              min       lq    median      uq      max    neval
----------------------------  -------  -------  --------  ------  -------  -------
np.random.uniform(size=1000)  23.8419  24.0803    25.034  25.034  46.9685      100

>>> print benchmark("np.random.uniform(size=1000)", "np.random.normal(size=1000)",
               "np.random.binomial(size=1000, n=10, p = .5)", "np.random.poisson(size=1000, lam=1)",
               units="us")
expr                                             min       lq    median        uq       max    neval
-------------------------------------------  -------  -------  --------  --------  --------  -------
np.random.uniform(size=1000)                 15.974   16.9277   16.9277   17.345    36.9549      100
np.random.normal(size=1000)                  56.982   57.9357   57.9357   58.8894  148.058       100
np.random.binomial(size=1000, n=10, p = .5)  57.9357  58.8894   64.1346  147.104   190.02        100
np.random.poisson(size=1000, lam=1)          61.9888  63.1809   66.0419   71.0487  105.143       100

>>> lbenchmark(evalStrings, units = "us")
expr                                             min       lq    median       uq       max    neval
-------------------------------------------  -------  -------  --------  -------  --------  -------
np.random.uniform(size=1000)                 24.7955  25.9876   25.9876  28.1334   56.0284      100
np.random.normal(size=1000)                  76.0555  78.9166   85.1154  86.0691  168.085       100
np.random.binomial(size=1000, n=10, p = .5)  59.1278  64.8499   66.5188  96.0827  120.163       100
np.random.poisson(size=1000, lam=1)          64.8499  69.6778   70.0951  71.0487   97.99        100

In R


> require(microbenchmark)
Loading required package: microbenchmark

> n = 1000

> microbenchmark(runif(n), rnorm(n), rbinom(n, p = .5, size = 10), rpois(n, 1))
Unit: microseconds
                          expr    min      lq median      uq     max neval
                      runif(n) 37.317 39.6890 40.914 42.2355  70.507   100
                      rnorm(n) 82.480 85.4910 86.836 88.2970 154.761   100
 rbinom(n, p = 0.5, size = 10) 92.181 94.5645 95.781 97.2435 153.588   100
                   rpois(n, 1) 39.534 41.3360 42.541 43.7010  85.183   100

Conclusion

We have created a few easy to use functions here for quickly carrying out benchmarks. These tools allow programmers to have an idea of what the performance of different functions will be and we provided some examples of usage. The scripts brings some R-like sugar to Python.

The timed RNG tests reveals a little better RNG performance from Python overall (particularly for the uniform distribution at more than x2 faster). We will probably come back to this topic at some later blog.

Data Science Consulting & Software Training

Active Analytics Ltd. is a data science consultancy, and Open Source Statistical Software Training company. Please contact us for more details or to comment on the blog.

Dr. Chibisi Chima-Okereke, R Training, Statistics and Data Analysis.