Create a gist now

Instantly share code, notes, and snippets.

What would you like to do?
R/Rpy2 Plotting along a Mathematica 9 Demo
#!/usr/bin/env python
# Copyright: Harald Schilly <harald.schilly@univie.ac.at>
# 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 :-)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment