Apr 202016
 

Some Notes about Makefile

记录一下最近用Makefile的经验和教训
内容比较凌乱。放在这里供以后参考。

1. Target-specific variable的作用域

Target-specific variable指的是变量在某些特定目标下的值。用法是:

target … : variable-assignment

注意这里的变量作用域只是在Recipe部分,而不能在目标(Target)或者依赖(Prerequisites)部分。

target … : prerequisites …
recipe

就是说这个Target-specific变量的值仅在recipe中有效。

2. 除错的方法

用warning 函数,比如$(warning msg)。它可以出现在各个部分: 目标, 依赖,recipe,顶级变量, 都行 (参见[2]):

$(warning A top-level warning)

FOO := $(warning Right-hand side of a simple variable)bar
BAZ = $(warning Right-hand side of a recursive variable)boo

$(warning A target)target: $(warning In a prerequisite list)makefile $(BAZ)
$(warning In a command script)
ls
$(BAZ):

输出:
$ make
makefile:1: A top-level warning
makefile:2: Right-hand side of a simple variable
makefile:5: A target
makefile:5: In a prerequisite list
makefile:5: Right-hand side of a recursive variable
makefile:8: Right-hand side of a recursive variable
makefile:6: In a command script
ls
makefile

还可以配合这些参数–just-print (-n), –print-data-base (-p), and –warn-undefined-variables

3. 管理大型项目

大型项目的代码有层级结构。xieMakefile要注意代码精简,善用函数,尽量自动化,避免重复代码,检查变量的有效性。

例如:
create-output-directories := \
$(shell for f in $(modules); \
do \
$(TEST) -d $$f || $(MKDIR) $$f; \
done)

例如用call, eval, foreach等函数:

# use foreach and test
.PHONY: validate_build
validate_build:
test $(foreach f,$(RELEASE_FILES),-s $f -a) -e .

local_src := $(addprefix $(subdirectory)/,playlist.y scanner.l)

# create library
define make-library
libraries += $(BINARY_DIR)/$1
sources += $2

$(BINARY_DIR)/$1: $(call source-dir-to-binary-dir, \
$(subst .c,.o,$(filter %.c,$2)) \
$(subst .y,.o,$(filter %.y,$2)) \
$(subst .l,.o,$(filter %.l,$2)))
$(AR) $(ARFLAGS) $$@ $$^
endef
$(eval $(call make-library, $(subdirectory)/libdb.a, $(local_src)))


参考
[1]Target-specific variables. GNU make. (link)
[2]Debugging Makefiles, Managing Projects with GNU Make, 3rd Edition, Chapter 12. (link)
[3]C and C++, Managing Projects with GNU Make, 3rd Edition, Chapter 8. (Link)

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包迁移过去,那这种工具必有远大前程。

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

Dec 112014
 

加快C++编译速度
Speed up C++ compiling speed

在开发rvtests时,常常需要重新编译整个项目,但是这是个耗费时间的事。

比如,编译最耗费时间的regression库,用“make”,需要50秒。
在多CPU的机器上,打开并行编译选项后,“make -j”,需要11秒(系统时间是11秒,用户时间是64秒)。

这就是说重新编译regression库,最快也要11秒。

今天我发现有一个好用的工具,ccache
可以把这个时间缩小到5秒 (加速编译速度一倍)。
方法是用:

make CXX="ccache g++"

ccache首次使用会稍微慢一点。因为它会缓存(cache)源程序。但之后编译的速度就变得飞快。
比如上面的方法用5秒。如果还想更快,打开“-j”选项,只需要不到1秒。

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的测试,还可以参考这个链接

Sep 082013
 

矩阵最大特征值的求法

Largest Eigenvalue of Matrix

 

矩阵的特征值分解是在科学计算常见的线性代数方法。在一些情况下,我们不需要全部特征值,只需要一部分特征值。或者因为矩阵的一些性质,我们只对一部分特征值感兴趣。比如对于矩阵 B = A’ * A, 如果A是m by n的实矩阵,m很大,n很小。根据线性代数理论,A, B矩阵的Rank最多是n,并且B矩阵是半正定矩阵。对矩阵B做特征值分解时,其特征值都是非负实数,并且最多有n个正数。因此,如果求解B矩阵的特征值,只需要知道其最大的n个特征值,因为其余特征值都是零。

求解一般矩阵的特征值主要有两种:(1)特征多项式(2)Power Method(包括 Arnoldi方法和QR分解)。 对于实对称正定矩阵,可以使用Choleskey分解。这些方法可以求出所有的特征值。而对于我们特别关心的实对称(Hermitain/Symmetric)正定(Positive Definite)矩阵,比如协方差矩阵,这里我们关心的是其最大(或者最小)的特征值。这里常用的方法的方法叫做Linczos方法,可以看做是Arnoldi方法的特例。Linczos方法计算A的特征值是(1)先随机生成一个向量v,(2)再利用一系列向量: v, A*v, A*A*v, …, A^(n-1) * v 来得到一个特殊矩阵T,这个矩阵的大小是m by m,这里m往往远大于n。(3)得到T之后,可以以O(m^2)的复杂度来的到特征值。

实现Lanczos方法时会使用Multiple restart,最有名的实现是ARPACK。ARPACK是用Fortran实现的,对Fortran不熟悉就不容易直接调用ARPACK的函数。不过很多其他软件都对ARPACK的功能进行了封装。在R里面使用ARPACK,可以用igraph包,R-blogger有一个这方面的例子(链接)。其他软件中一般也都有用ARPACK来实现求一部分特征值的功能:Matlab的eigs,Scipy的eigs和eighs(分别对应实对称矩阵和Hermitian矩阵, 链接),Julia的eigs(链接)。

后话:现在开源的科学计算软件很多,当我们想寻找合适高效的数值计算工具时,可以从源代码出发。比如我们知道Matlab有eigs函数,就可以用“julia eigs”来看julia语言里是怎么实现eigs的。

 

 

Sep 052013
 

Bioinformatics Tools (build for Windows)

1. Tabix (Download)

2. Bgzip (Download)

These tools are build using MingW, with necessary changes to source codes:

  • Use _WIN32 definition to make sure fopen() use “rb” to open files, not “r”. (fopen is equivalent to _bgzf_open)
  • Link with zlib-2.1.8
  • Tested with example data.