Aug 122015
 

R的诡异错误 (3)
A pitfall in R (3)

下面的代码有两行,第一行计算Genotype的Minor Allele Frequency,第二行计算Weight matrix (用SKAT paper建议的Beta分布的weight)。
下面的代码有问题吗?

# G is a sample by marker matrix
maf <- colMeans(G) / 2
w <- diag(dbeta(maf, 1, 25))

粗粗一看,上面的代码的思路是正确的,不过这里有两个坑。
第一个坑是缺失数据(Missing data),如果Genotype矩阵G有缺失数据,那么maf和w都会有缺失的数据;
第二个坑比较难看出来:当G只有一列的时候,maf也会是一个标量(Scalar),同样dbeta函数的返回值是标量。
你可能想diag函数作用在一个标量上,其返回值是一个1×1的矩阵。
但是,diag()的返回的是一个….不确定的结果,看个例子吧:

diag(3)            => 3x3 单位矩阵
diag(3.9)          => 还是3x3 单位矩阵
diag(0)            => 是0x0 矩阵
diag(-1)           => 报错
diag(c(1, 2, 3)    => 3x3 矩阵, 对角线是1,2,3

总之,第二个坑告诉我们,如果dbeta返回的是一个长度为1的数值,别指望w一定是1×1的矩阵。

继续吐槽:
diag()这个函数设计的不好用。在R内部,标量就是被实现为长度为1的向量。
在现在的设计下,diag要求使用者同时考虑变量的类型(矩阵,向量,缺失)以及变量的长度和值(向量长度为1且大于等于0,向量长度为1且小于0,向量长度大于1)。
对于一个适应了使用有类型系统的编程语言的程序员,这种要求会增加使用者的脑力负担,同时阅读者在语义上也很容易出现混淆。
我建议的解决方案是:diag(向量)永远输出一个对角线是该向量的矩阵,同时定义一个新函数(可以叫diag_len,类比seq_len;或者叫eye,类比Matlab里的eye函数)专门用来生成单位矩阵。

最后总结一下,我觉得这样写比较

# G is a sample by marker matrix
maf <- colMeans(G, na.rm = TRUE) / 2
m <- length(maf)
w <- matrix(0, nr = m, nc = m) 
diag(w) <- dbeta(maf, 1, 25) 
# 或者:w <- diag(dbeta(maf, 1, 25), nrow = m, ncol = m)
Jun 162015
 

不要用C++11特性开发R包
Don’t use C++11 features to develop R packages

这是个教训。在我开发seqminer的时候浪费太多时间了。希望其他人不要重复我的错误。

C++11有不少有意思的特性,但为什么不要在开发R包的时候使用C++11的特性呢(这里指R 3.2)?下面有一些原因。

1. 没有清晰的文档
R对于C++11的支持写的很模糊。你需要看两个文档,Wring R ExtensionsR Installation and Administration
第一个文档提到你需要CXX_STD,第二个文档提到CXX1XFLAGS。这时候你知道应该设置一些变量,但怎么设置?没有文档告诉你。

这里我找到一个方案,可以凑合用。你需要改configure.ac文件。不过你仍会碰到下面提到的第3个问题。

dnl https://github.com/cran/dfcomb/blob/master/configure.ac
dnl copy the following lines to configure.ac

AC_DEFUN([AX_GET_R_CONF_VAR],[$1=`"${R_HOME}/bin/R" CMD config $1`])
AX_GET_R_CONF_VAR([CC])
AX_GET_R_CONF_VAR([CFLAGS])
AX_GET_R_CONF_VAR([CPPFLAGS])
AX_GET_R_CONF_VAR([CXXCPP])
AX_GET_R_CONF_VAR([CXX1X])
AX_GET_R_CONF_VAR([CXX1XSTD])
AX_GET_R_CONF_VAR([CXX1XFLAGS])
AX_GET_R_CONF_VAR([CXX1XXPICFLAGS])

CXX="${CXX1X} ${CXX1XSTD}"
CXXFLAGS="${CXX1XFLAGS} ${CXX1XPICFLAGS}"
CXXPICFLAGS="${CXX1XXPICFLAGS}"

2. 跨平台支持很艰难

当你把包提交到CRAN时,很有可能在Solaris系统上出现编译或运行的错误。
如果你使用了C++11的功能,你可能想在自己安装Solaris试试。
从我的经验来看,Solaris11 + Solaris Studio 12.4环境比较好搭建。
但CRAN上用的是很老的Solaris10 + Solaris Studio 12.3。
如果你想搭建这个系统,你需要从源代码编译R,然后你会发现,编译不了!
我发现blas库无法被链接。因此你在提交R包到CRAN时,最好明确说,我不支持Solaris。

3. 支持Linux系统也会出问题

在Linux系统上,也许gcc的版本很新,可以支持C++11,但如果安装R的时候没有设置C++11需要的一些变量,你还是用不了。
比如Ubuntu系统里, R是用apt安装的,这个版本的R是默认不知道你有没有支持C++11的编译器的。
这就是说,你在开发时使用的R,CRAN上的R都可以安装你的包,但不保证你的用户能安装。

怎么知道你的R是否支持C++11呢?你可以运行“R CMD config”:

  CXX1X         C++ compiler command for C++11 code
  CXX1XSTD      flag used to enable C++11 support
  CXX1XFLAGS    C++11 compiler flags
  CXX1XXPICFLAGS
                special flags for compiling C++11 code to be turned into
                a shared library

然后看一下CXX1X的值, “R CMD config CXX1X”,如果这个值是空的。默认情况下没法编译安装你写的R包(install.packages())。
顺便说一下: CXX1XXPICFLGS是一个笔误,正确的拼写是CXX1XPICFLAGS (不是两个XX)

解决方法:
先在shell下:

export CXX1X=”g++ -std=c++11″
export SHLIB_CXX1XLD=”g++ -std=c++11 -shared”

然后R里面就能用install.packages()安装了。
顺便说一下,SHLIB_CXX1XLD 没有写在手册里……

现在总结一下,如果你用了C++11的特性写R包,你会花很多时间读手册让R去选用合适的编译器,很难去支持CRAN上的老旧的Solaris,不能保证最终用户能够编译安装你的R包。
所以,我的建议的是,如果想让更多的人用你的R包,别折腾C++11的特性去适应R,还是写最portable的C++代码吧。

Apr 012015
 

R的诡异错误(2)
A pitfall in R

R里面的陷阱一个接一个,继上次说的诡异错误,最近我又跌进了另一个坑。

先创造一个叫df的data.frame, 如下所示:

options(stringsAsFactors = FALSE)
df <- data.frame(
  group = c("group1", "group2"),
  val = c(100, 9)
  )
df
# df content:
#
#   group val
#1 group1 100
#2 group2   9

我们统计一下val列在每组(group)的均值:

## get mean value per group
library(plyr)
res <- ddply(df, .(group), function(x) {
  c(group = unique(x$group), 
    val = mean(x$val))
  })
res

res的结果如下:

> res
   group val
1 group1  100
2 group2  9
> res[1,2] > res[2,2]
[1] FALSE

困惑的地方来了,res[1,2]显示的是“100”, res[2,2]现实的是“9”, 为什么res[1,2] > res[2,2]的结果是FALSE?
R肯定是认为100比9大。这显然是不对的。

那上面哪个地方有陷阱呢?再仔细看一下高亮显示的两行,看出问题了么?

## get mean value per group
res <- ddply(df, .(group), function(x) {
  c(group = unique(x$group), 
    val = mean(x$val))
  })

陷阱在于第3-4行,当用c()来构造一个向量时,val被自动从数值型转成字符型。也就是说,存储的是“100”和“9”,不是100和9。
在字符型数值的比较中,“9”是会比“100”大的。因此order输出的是2 1.

从这个坑里爬出来,不由得吐槽一下,R的自动类型转换能减少程序的代码量,但是这种类型转换很容易引入错误。
也就是说R编程时不需要程序员去检查类型,但在底层R会自动的(有时也是自作聪明的)转换类型。
这就要求程序员额外费心考虑类型之间的自动转换,实际上加重了编程的负担。
如果不然,程序员放任R做自动的类型匹配和转换,就可能得到错误的结果,这在科学研究中实在是很危险。

既然R的坑这么多,我想过去试试新出的Julia或者Go语言。可惜现在这些语言都没有很好的支持统计计算,比如Julia目前缺少对缺失数据(missing data)的支持,
比较有希望的NumPy写起来语法繁杂,近期吸引不到大多数的统计学家(统计的用户太少,机器学习的用户多)。
要是能有一种计算工具,可以让统计学家很容易的使用,同时又有较强的类型系统,还能把大多数R包迁移过去,那这种工具必有远大前程。

Mar 052015
 

设置代理服务器
Setting up proxy

出于对患者信息的保护,学校内网接入Internet时必须使用代理(Proxy)。
这令很多软件的安装使用过程变得复杂。

这个帖子会列出给常用软件设置代理的方法。

R

最直接的设置代理的方法是设置环境变量: http_proxy, https_proxy, 例如:

Sys.setenv(http_proxy="http://proxy.swmed.edu:3128")
Sys.setenv(https_proxy="http://proxy.swmed.edu:3128")

如果上述方法不起作用,也有可能需要:

# In R, a partial solution is to use:
options(download.file.method="wget")

这是在curl不能正常使用代理,但wget可以用的时候的解决方案。

还有一种情况是用devtools安装R的扩展包,可以用:

# For devtools
library(httr)
set_config(
  use_proxy(url="proxy.swmed.edu", port=3128)
)

如果Mac电脑上没有安装X11,装package时需要指定repository(来源),比如:

install.packages('RMySQL', repos='http://cran.us.r-project.org')

Docker

很多软件提供了Dockerfile。用户需要使用docker build来制作自己的容器(docker)。
如果在防火墙后面, 简单的方法是用(link):

docker build 
  --build-arg https_proxy=$HTTP_PROXY --build-arg http_proxy=$HTTP_PROXY 
  --build-arg HTTP_PROXY=$HTTP_PROXY --build-arg HTTPS_PROXY=$HTTP_PROXY 
  --build-arg NO_PROXY=$NO_PROXY  --build-arg no_proxy=$NO_PROXY -t java .

另外一个方法是手动修改Dockerfile,加入ENV语句,例如:

ENV HTTP_PROXY http://proxy.swmed.edu
ENV HTTPS_PROXY http://proxy.swmed.edu
ENV http_proxy http://proxy.swmed.edu
ENV https_proxy http://proxy.swmed.edu
....
Feb 032015
 

R函数式编程
Functional Programming in R

R语言的循环很慢,同时语法有很复杂,陷阱很多。举例来说,逐行遍历一个矩阵有两种写法:

(1)

for (i in seq_len(nrow(m))){print(i)}

(2)

for (i in 1:nrow(m)){print(i)}

这两个写法那个对?或者两个都对?
熟悉R的用户自然会用(1),因为当m矩阵的行数为零时,只有(1)的结果是正确的:

> m <- matrix(, nr = 0, nc = 0)
> for (i in seq_len(nrow(m))){print(i)}
> for (i in 1:nrow(m)){print(i)}
[1] 1
[1] 0

为了避免手动写循环,一种折中的解决方案是用Filter, Reduce, Map等函数式编程的概念。
这些函数可以比较自然的处理列表(list)或者向量(vector)类型。举个例子:

l <- as.list(seq(9, 1, -1))
d <- data.frame(a = c(1, 2), b = c(3, 4))
m <- matrix(seq(9, 1, -1), 3, 3)
v <- seq(9, 1, -1)

isEven <- function(x) { x %% 2 == 0 }
isEven(1:3)

Filter(isEven, l)
## Filter(isEven, d) ## this does not work
Filter(isEven, m)
Filter(isEven, v)

Position(isEven, l)
# Position(isEven, d) ## this does not work
Position(isEven, m)
Position(isEven, v)

Find(isEven, l)
# Find(isEven, d) ## this does not work
Find(isEven, m)
Find(isEven, v)

可以看到,对于列表(list)或者向量(vector)类型,R内建的函数式函数(更多的见[1])还是很方便的。
如果想把函数式编程玩出更多花样,可以看看Hadley的新书Adv-R(参见[2])。

这里还有个问题,这种函数式的写法速度如何?下面的代码告诉我们函数式编程节省了写程序的时间,在某些情况下也许也能节省计算时间。

library(microbenchmark)
v <- seq(10000)
m1 <- function(v) {
  v[v%%2 == 0]
}
m2 <- function(v) {
  res <- vector()
  for(i in v) {
    if (i %% 2 == 0) {
      res <- c(res, i)
    }
  }
  res
}
m3 <- function(v) {
  Filter(isEven, v)
}

microbenchmark(res1 <- m1(v),
               res2 <- m2(v),
               res3 <- m3(v), times = 1000)

下面是结果:

Unit: microseconds
          expr       min        lq       mean    median         uq       max neval
 res1 <- m1(v)   341.038   348.526   529.9475   356.225   380.6965  58359.46  1000
 res2 <- m2(v) 33687.987 38433.148 48014.3656 39791.000 41751.3955 303604.96  1000
 res3 <- m3(v)  6178.547  6648.369  8493.9365  6999.976  7445.0725 124781.90  1000

在上面的例子中,
第一种是向量化的赋值,速度最快。
第二种是用循环,速度最慢,
第三种是函数式,速度居中。

不过据我的经验,当v的维数很小时(比如只有100个数),函数式相对于循坏的优势越来越小,甚至比循环要慢。

题外话:据说R的核心是类似于Scheme,也是一种函数式语言。但是在速度上对函数式编程的支持看来还有很大的提升空间。

[1]https://stat.ethz.ch/R-manual/R-devel/library/base/html/funprog.html
[2]Advanced R http://adv-r.had.co.nz/Functional-programming.html

Apr 132014
 

R package的UBSAN测试
How to perform UBSAN tests on R packages

最近需要提交一个R package (seqminer)到CRAN,我们需要保证这个包代码通过UBSAN 测试(Undefined bahavior sanitizer tests,详见这个链接).
这个测试可以有效的检测结果不确定的指令(例如打印INT_MAX+1)。

目前只有Clang(版本3.3以后)才能支持UBSAN。为了对R package进行UBSAN测试,我们首先需要一个特殊的R版本,即用打开UBSAN支持(-fsanitize=undefined)的Clang编译的R。参考了网上一篇Blog,写一下简单步骤:

1. 下载编译工具

sudo apt-get install valgrind subversion r-base-dev clang-3.4 texlive-fonts-extra texlive-latex-extra
sudo apt-get build-dep cran-base

2. 下载R源代码(devel版本)

svn co https://svn.r-project.org/R/trunk ~/R-devel

3. 替换R-devel/config.site

CC="clang -std=gnu99 -fsanitize=undefined"
CFLAGS="-fno-omit-frame-pointer -Wall -pedantic -mtune=native"
F77="gfortran"
LIBnn="lib64"
LDFLAGS="-L/usr/local/lib64 -L/usr/local/lib"
CXX="clang++ -std=c++11 -fsanitize=undefined"
CXXFLAGS="-fno-omit-frame-pointer -Wall -pedantic -mtune=native"
FC=${F77}

4. 编译

cd R-devel
./configure --with-x=no --without-recommended-packages
make

这里–without-recommended-packages会省略一些R的packages(例如boot,nlme,Matrix) 以提升编译速度。如果不用这个选项,我发现Matrix package总是没法编译

5. 建立~/.R/Makevars

CC = clang -std=gnu99 -fsanitize=undefined -fno-omit-frame-pointer
CXX = clang++ -fsanitize=undefined -fno-omit-frame-pointer
PKG_LIBS = /usr/lib/llvm-3.4/lib/clang/3.4/lib/linux/libclang_rt.ubsan_cxx-x86_64.a

这里CXX 一定要用clang++,如果用clang的话会出现链接错误,比如找不到这个符号: _ZTVN10__cxxabiv117__class_type_infoE

这里还需要指定PKG_LIBS包括LLVM的库文件,不然就找不到UBSAN的符号。

6. 测试package

R CMD build seqminer/
R CMD check seqminer_2.7.tar.gz

UBSAN测试的结果会显示在:seqminer.Rcheck/tests/tests.Rout
如果有错误,你会看很长的一行, runtime error….

其他技巧

检测内存还可以用valgrind:

R CMD check --use-valgrind seqminer_2.7.tar.gz
R -d valgrind --vanilla < mypkg-Ex.R
R -d "valgrind --tool=memcheck --leak-check=full" --vanilla < mypkg-Ex.R
Mar 092014
 

A pitfall in R

先看下面的代码:

freq <- 0.29
altAllele <- freq * 200 * 2
print(altAllele) ## this will output 116
print(length(rep(1, altAllele))) ## what is the output of this??

最后一行输出的是115, 不是116。为什么?
答案是:

altAllele是一个numeric类型的数,尽管打印出来是116,但在电脑里存储的是一个比整数116稍小的浮点数。因此rep(1, altAllele)是一个长度是115的vector。

R的程序有很多陷阱。这个陷阱是因为R不需要定义数据类型,但它的底层还是有类型的。因此写R程序一定要注意这点。解决方法如下:

freq <- 0.29
altAllele <- as.integer(freq * 200 * 2)
print(altAllele) ## this will output 116
print(length(rep(1, altAllele))) ## what is the output of this??

Python里有没有类似的陷阱?结论是一样有 (Python 2.7.5 NumPy 1.8.0):

freq = .29
altAllele = freq * 200 * 2
len(np.repeat(1, altAllele)) ## output 115 here.

但是Python里打印altAllele的值不是116,而是115.99999999999999。这总比R打印出16要好一点。

Oct 182013
 

R语言如何用索引
How to use index in R

R语言有三种索引方式: [[, [, $
本文基于R语言的定义,用例子里说明这三种索引方式的使用条件和它们之间的细微差别。

首先,R里面有下面的索引方式:

x[i]
x[i, j]
x[[i]]
x[[i, j]]
x$a
x$"a"

总结并翻译一下:

1. 如果x是一个vector (包括matrix和array),用[i]或者[[i]]得到的是第i个元素。而且一般用[而不是[[,用[[的话会失去变量的名字(name和dimnames)。比如:

aVec <- c(1,2,3)
names(aVec) <- c("x1", "x2", "x3")
> aVec[1]
x1 
 1 
> aVec[[1]]
[1] 1

2. 如果x是一个list,那么[i]得到的是一个list,用[[i]]得到的是第i个元素。aList[1]得到的是一个list,aList[[1]]得到的一个vector。

aList <- list(1, 2)
> aList[1]
$x1
[1] 1

> aList[[1]]
[1] 1

3. 如果索引是一个vector(比如一次选取多个元素),只能用[,不能用[[。选出来的元素类型,索引一个或者多个是一样效果。

4. 用$加一个变量来索引list,这个变量会按字面选取,不会先求值再选取。比如: x<-1; aList[[x]] 和 aList[[1]]不一样。 5. 用$来索引list之外的元素,得到的是NULL 关于各种index的测试,还可以参考这个链接

Jul 082013
 

让R更快
Speed up R

用R 3.0.1来做基准,用Rscript运行R-benchmark-25.R计时。

1) 基准时间:

197.16u 1.74s 199.71r 1703184kB 0 Rscript R-benchmark-25.R

之后下载OpenBLAS,编译,然后把编译好的libopenblas.so 替换R自带的libRblas.so.
2) 测试单线程的速度:

Single thread OpenBLAS

105.83u 108.53s 58.04r 3072256kB 0 Rscript R-benchmark-25.R

 

3) 下面试试四个线程:

Four threaded OpenBLAS

export OPENBLAS_NUM_THREADS=4

79.39u 23.67s 58.54r 2122688kB 0 Rscript R-benchmark-25.R
可以看出使用四个线程和一个线程耗时差不多一样。但如果细看矩阵计算时间(这里忽略),四个线程还是能提速不少的。
但是考虑到使用方便,一个线程在大多数情况下应该就够用了。