R/Rpy2 Plotting along a Mathematica 9 Demo

mma9rpy2-2.png

mma9rpy2.png

mma9rpy2.py

#!/usr/bin/env python
# Copyright: Harald Schilly <[email protected]>
# License: Apache 2.0
# This script demos how rpy2 is capable to interface with R
# Inspired by: http://www.wolfram.com/mathematica/new-in-9/built-in-integration-with-r/hierarchical-clustering.html
# Dependencies (in Linux Ubuntu): python-numpy r-recommended r-cran-ggplot2 python-rpy2

# get the R environment
import rpy2
from rpy2 import robjects as robj
from rpy2.robjects import r
from rpy2.robjects.packages import importr

# R's base package, base.<TAB> has completion
base = importr('base')

# lattice plot package
lattice = importr('lattice')

# package containing lot of datasets
datasets = importr('datasets')

# activate seamless conversions of numpy ndarrays to R objects
from rpy2.robjects import numpy2ri
numpy2ri.activate()

# despite Python's "print", we can also expose R's "print"
rprint = rpy2.robjects.globalenv.get("print")

### the actual fun starts here

# Data straight from numpy
import numpy as np
mdata = np.random.randn(10, 5)

# Python's print vs. R's print
print mdata
rprint(mdata) # different
print r.summary(mdata) # rprint should be the same

# The following are two different ways, how the dimnames can be set
if False:
  # Case 1: let's do what Mathematica does
  # y in the global variable scope
  robj.globalenv['y'] = mdata

  # evaluation without ugly escaping
  r("""
  dimnames(y) <- list(
                  paste("g", 1:10, sep=""),
                  paste("t", 1:5,  sep=""))
  y
  """)

  # show what's there:
  mdata = y = robj.globalenv['y']
  rprint(y)

else:
  # Case 2: robject.do_slot("...") and rinterface.StrSexpVector
  # see: http://rpy.sourceforge.net/rpy2/doc-2.1/html/rinterface.html#pass-by-value-paradigm
  # mdata converted to an R object
  mdata = numpy2ri.numpy2ri(mdata)
  import rpy2.rinterface as ri
  # get dimensions automatically (as an extra)
  dims = list(mdata.do_slot("dim"))
  descr = ri.baseenv["list"](
           ri.StrSexpVector(['g%s'%_ for _ in range(dims[0])]),
           ri.StrSexpVector(['t%s'%_ for _ in range(dims[1])]))
  mdata.do_slot_assign("dimnames", descr)
  # Note: Laurent Gautier pointed out, that one can use ro.ListVector
  # not sure how, though ...
  # http://blog.harald.schil.ly/2012/11/mathematica-9s-r-integration-vs-rpy2.html?showComment=1354737231460#c7980673096108868874

# correlation and distance matrices
corrm = r.cor(r.t(mdata), method="spearman")
stats = importr("stats")

# no idea how to do the 1-matrix automagically
robj.globalenv['corrm'] = corrm
distm = stats.as_dist(r("1-corrm"))

rprint(corrm)
rprint(distm)

# do the clustering
hr = stats.hclust(distm, method = "complete", members = r("NULL"))

# plotting can be done, too
# http://rpy.sourceforge.net/rpy2/doc-2.2/html/graphics.html

# fire up ggplot2
ggplot2 = importr('ggplot2')
# we won't use it, see rpy2 documentation for working examples

# plotting to PNG, instead of the usual X11
grdevices = importr('grDevices')
grdevices.png(file="mma9rpy2.png", width=512, height=300)
try:
  r.par(mfrow = r.c(1,2))
  r.plot(hr, hang = 0.1)
  r.plot(hr, hang = -0.1)
finally:
  # end of plotting, closing png device
  grdevices.dev_off()

# feel free to replace png with pdf (function and filename)
grdevices.png(file="mma9rpy2-2.png", width = 512, height = 512)
try:
  r.heatmap(mdata)
finally:
  grdevices.dev_off()

# Note: the artefacts in the plot are known, workaraound in the pipeline :-)