Feb 192015
 

2015 春节
2015 Chinese New Year

爆竹声中一岁除,春风送暖入屠苏。
千门万户曈曈日,总把新桃换旧符。

——王安石的《元日》

除夕夜里,稍稍反思一下2014年经历的投稿。

背景:我们写的软件速度快,功能全。比如在文章中,我们强调该软件可以处理很大规模的数据,这是已知其他软件达不到的。

审稿:除去离题万里的审稿意见以及错误的审稿意见,我们得到的反馈是对自己的软件过分称赞,对有竞争关系的其他软件过分贬低。同时这些意见往往对我们有强烈的负面情绪。

反思:比较软件从某种意义上就像是比较孩子。说自己的孩子好没问题,但要是说自己的孩子比别人的好,就易引起强烈的情绪反扑。即使是数字上的比较(比如身高说自己的孩子高于其他小孩XX厘米),有意见的人也可能会说:身高不是重要的标准,你比较的方式不对等等。这个经历告诉我指出其他人/事在科研上的不足应该非常谨慎。或许寻找“和气”、容易接受的比较方式是有必要的。

2015年我们会继续努力,写有价值的文章,做好科研,也希望能通过温和的方式,得到其他科研同行的肯定。此外,希望看到这里的各位 新春快乐!

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