比较gzip压缩和解压缩性能

Benchmark gzip compression and decompression performance

gzip是很多测序数据的默认压缩格式。在学校HPC环境中,几个GB的文件的压缩和解压缩往往需要几分钟到几十分钟的时间。如果处理的数据较大,gzip的压缩和解压缩会占用不少时间。因此我想看看有没有比Linux默认的gzip更有效率的软件。这里主要比较:

  1. libdeflate 这是samtools的选择
  2. igzip 这包含在Intel ISA-L软件中
  3. pigz 这是gzip的多线程实现

具体的测试分为压缩和解压缩两部分。用的是一个FASTQ文件,压缩后大约2G,原始文件6G。

测试结果如下:

Decompression
User Time (s)Real Time (s)Memory
gzip39.73u1.47s41.24r876kB
igzip7.60u1.81s11.31r2844kB
libdeflate20.45u6.53s28.05r7909292kB
pigz30.18u4.86s23.49r1000kB
Decompression Benchmark
Compression
User Time (s)Real Time (s)MemoryCompressed File Size
gzip620.17u2.54s623.23r1116kB1.81G
igzip15.61u1.92s17.69r2628kB2.03G
libdeflate152.88u10.29s163.43r7879072kB1.78G
pigz995.38u4.15s18.18r38352kB1.82G
Compression Benchmark

结论:Intel igzip在耗时指标中完胜。在默认设置下,igzip解压缩时间是系统提供gzip的1/5,压缩时间是1/40。唯一不足的是压缩后的文件偏大。

讨论:

这个结果和samtools 的benchmark差别很大。我发现samtools用的GitHub仓库(https://github.com/jtkukunas/zlib)不同于Intel Tuning Guide里提供的仓库(https://github.com/intel/isa-l/blob/master/README.md)。 因此性能的差别可能来源于samtools没有用官方的实现,当然也有可能是测试方法不同:samtools测量的是bam文件的读写速度,而我只是用FASTQ文件作为测试。

其他相关的gzip工具还有zlib-ng 和cloudflare优化的zlib。这些性能应该在igzip和系统自带的gzip之间,这里就不再测试了。

节省编译时间的工具

节省编译时间的工具
Tools save compiling time

最近发现编译C++程序很花时间。
比如说仅仅编译一个用了模版库的程序需要花10秒。
为了加快开发速度,我想找一个能提高单机编译速度的程序。结果如下:

Benchmark
ProgramUserTime (s)SysTime (s)RealTime (s)Memory (kB)
g++-4.9.4137.35u10.11s56.44r1392268kB
clang++-4.093.24u6.56s36.21r949524kB
zapcc++-1.01.81u0.57s33.95r39888kB
ccache-g++ (1st)109.62u9.01s49.43r1392392kB
ccache-g++ (2nd)1.47u0.87s1.67r55832kB

结论:
从编译速度来看,首次编译速度最快的是zapcc(不免费)和clang(免费),耗时~35秒左右,比g++快40%。
重复编译的时候ccache最快,提速10倍以上。

注:表中zapcc的用户时间很短,内容占用少,这是因为它是起动其他进程去真正编译,其他程序的用户时间和内存使用没有被计入。

卡方检验的公式(2×2)

卡方检验的公式(2×2)
A short-cut formula for Chi-square test a 2 by 2 table

最近发现一个快速计算卡方检验的公式。
对于2×2的表格

\(
\begin{array}{|c|c|c|}
\hline
a & b & (a+b)\\ \hline
c & d & (c+d)\\ \hline
(a+c) & (b+d) & n\\ \hline
\end{array}
\)
上式中\(n=a+b+c+d\)

卡方检验的统计量是:
\(
\chi^2_1 = \frac{n(ad-bc)^2}{(a+b)(c+d)(a+c)(b+d)}
\)

这个方法和R计算的结果(不要用continuity correction, e.g. chisq(…, correct = FALSE))是一致的。

-- February 17, 2017更新 --

刚刚发现Matthew’s Correlation Coefficient (MCC)的定义和2×2 table的卡方检验有如下联系
\(
|MCC| = \sqrt{\frac{\chi^2}{n}}
\)
因为MCC的范围是-1到1,因此\(\chi^2 \)的范围是0到\(n\)

Makefile 笔记

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)

解决一个奇怪的网络连接错误

解决一个奇怪的网络连接错误
Solve a strange network problem

我的服务器bunny最近出现了一个奇怪的网络连接错误:只有台式机和服务器能连接到bunny,笔记本就连不上。
假设网络IP如下:

A. 服务器bunny:198.215.54.48 10G 网络
B. 服务器: 198.215.54.5 10G 网络
C. 服务器bronco: 129.112.7.169 1G 网络
D. 台式机:129.112.185.246 局域网
E. 笔记本:172.17.157.121 无线网

现在的问题是可以从B, C, D 连接到A,但不能从E连接到A。

解决思路如下:

1. 怀疑E->A的路由有问题

通过traceroute,发现E到A的最后一跳是*,状态是Host unavailable.
但进一步发现E->B,A和B在同一个子网,因此E->A的路由不大可能有问题

2. 检查A的防火墙设置

先关闭所有防火墙,用:

ufw disable

但E仍然不能连A

3. 检查A的路由

> route -n
Kernel IP routing table
Destination Gateway Genmask Flags Metric Ref Use Iface
0.0.0.0 198.215.54.254 0.0.0.0 UG 0 0 0 em1
172.17.0.0 0.0.0.0 255.255.0.0 U 0 0 0 docker0
192.168.1.0 0.0.0.0 255.255.255.0 U 0 0 0 em4
198.215.54.0 0.0.0.0 255.255.255.0 U 0 0 0 em1

> ip route --list
ip rule list
0: from all lookup local
32766: from all lookup main
32767: from all lookup default

仔细检查发现172.17.0.0这个路由有问题:因为E在172.17.0.0网络,从E发出的ICMP包到达A,而从A返回的ICMP包没有通过em1接口而是docker0发出,这样E->A就显示了Host unavailable。

值得注意的是在路由IP包时,它的终点并不是第一个可以匹配的路由,而是所有路由里的第一个匹配最好的路由。
原文:Any entry whose first field matches the destination IP address completely(a host) or partially (a network) would signal the IP address of the next router. (link)

找到原因后,我们可以让docker0的网络IP地址避开172.17.0.0子网,这就可以解决E->A的连接问题。
这个解决方案里的重要命令是 (参考链接


ip link set dev docker0 down
ip addr del 172.17.42.1/16 dev docker0 # 172.17.42.1是docker0的ip
ip addr add 10.0.0.10/24 dev docker0
ip link set dev docker0 up
ip addr show docker0 #检查docker0的ip

此外要在/etc/default/docker里加上DOCKER_OPTS (–bip参数见Docker文档Customize the docker0 bridge)。这样docker0的IP不会在服务器重启后改变。


DOCKER_OPTS="--bip=10.0.0.10/24"

题外话

1. 为什么docker和无线网都用172.17.0.0子网?

根据RFC 1918 3. Private Address Space, 保留IP除了常见的192.168.0.0外, 还有:

10.0.0.0 - 10.255.255.255 (10/8 prefix)
172.16.0.0 - 172.31.255.255 (172.16/12 prefix)
192.168.0.0 - 192.168.255.255 (192.168/16 prefix)

因为172.17.0.0是保留地址,正巧docker和无线网同时都使用了这个子网,因此造成了我遇到的网络问题。

2. 怎么监测服务器是否接收到IP包?

可以使用iptables

例如参考 http://askubuntu.com/questions/348439/iptables-log-file-and-how-change-it, 记录192.168.11.0/24发来的包:


iptables -A INPUT -s 192.168.11.0/24 -j LOG --log-prefix='[netfilter] '

可以在/var/log/kern.log 找含有“[netfilter]”的日志。

如果想监测ICMP包,可以用:


iptables -A INPUT -p icmp -j LOG --log-prefix='[netfilter] '

其他常用命令(参考这个链接)包括:


sudo iptables -L # 打印iptables内容
sudo iptables -L INPUT -v # 检查INPUT表收到了多少数据
sudo iptables -Z # 清零
sudo iptables -Z INPUT # 清零INPUT表数据清零
sudo iptables -A INPUT -p icmp -j LOG --log-prefix='[netfilter] ' # 记录ICMP包
sudo iptables -D INPUT -p icmp -j LOG --log-prefix='[netfilter] ' # 删除防火墙规则(记录ICMP包)
sudo iptables -F #防火墙规则生效

以下命令重置防火墙

sudo iptables -P INPUT ACCEPT
sudo iptables -P FORWARD ACCEPT
sudo iptables -P OUTPUT ACCEPT
sudo iptables -t nat -F
sudo iptables -t mangle -F
sudo iptables -F
sudo iptables -X

iptable是个复杂的防火墙,它使用了table,chain等概念,比如下图:

iptables Overview
iptables Overview

(来自链接
更详细的介绍可以参考[2]。

除iptables外,Ubuntu还提供了UFW,这个工具建立在iptables上,提供比iptables更简易的配置。比如:


ufw enable # 启用防火墙
ufw disable # 关闭防火墙
sudo ufw status # 显示防火墙设置
ufw deny 80/tcp # 封锁tcp的80端口
sudo ufw delete deny 80/tcp # 删除上一条(封锁tcp的80端口)

更多的信息可以参考Ubuntu UFW文档 (链接)。

 

参考

[1]7 Linux Route Command Examples (How to Add Route in Linux)

[2]Linux Firewall Tutorial: IPTables Tables, Chains, Rules Fundamentals

建立SFTP共享账户

建立SFTP共享账户
Create SFTP share account

做学术的时候往往需要交流数据。目前最安全简便的解决方案是用SFTP传输文件协议。
为此,需要配置Linux服务器,并为学术伙伴建立一个不能登录、只能用SFTP的账号。

下面是步骤:

1. 建立用户组

addgroup --system filetransfer

2. 设定用户目录和相应权限
以用户名username为例:

mkdir /home/username
mkdir /home/username/data
chown root:root /home/username
chown username:filetransfer /home/username/data
chmod 755 /home/username
chmod 700 /home/username/data

权限755可以允许sftp用户登录到其HOME目录/home/username
权限700可以防止上传的数据被其他用户访问

3. 建立用户

useradd -d /home/username -s /usr/lib/sftp-server -M -N -g filetransfer username

-d 指定HOME目录
-s 指定登录shell
-M 不创建HOME目录
-N 不创建以username为名的 用户组
-g 指定user group

4. 配置ssh

打开 /etc/ssh/sshd_config

见到这行之后:Subsystem sftp /usr/lib/openssh/sftp-server
改称:Subsystem sftp internal-sftp
然后在末尾加上:

Match group filetransfer
ChrootDirectory %h
X11Forwarding no
AllowTcpForwarding no
ForceCommand internal-sftp

5. 配置登录shell

在/etc/shells末尾加上一行:

/usr/lib/sftp-server

6. 最后重启ssh 服务器就好了

service ssh restart

R的诡异错误(3)

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)

深度学习MNIST数据

深度学习MNIST数据

Deep learning on MNIST data

最近试了一下Deep Learning算法,在Kaggle的Digit Recognizer排行榜上目前排19名。

简单说说经验。

从工具上讲,建议先用Keras来做原型(代码简单,很容易构建复杂的网络,缺点是太占内存/显存),再用Caffe做更细致的调节(计算快,占内存/显存少,但用起来复杂,比如不支持RMSProp,需要手动把数据分成Training/Testing,要手动用proto buffer来构建神经网络)。此外,用GPU加速是必须的,我观察到GPU速度至少是CPU速度的60倍。

从模型上讲,基本上建立不同的Convolution Network。套路似乎都是重复这个结构:Convolution Layer + ReLU Layer + MaxPooling Layer。模型可以更Deep (层数更多),比如VGG网络,但是代价是参数数量多,计算量太大。对于MNIST这个整理好(Pre-processed)的数据,这些网络的效果很相似,一般默认设置都能达到98%正确率以上(用42k数据训练,28k数据测试)。如果要更好的成绩,基本上需要更多的输入数据(比如把原始的图像上下左右移动),更深的模型(比如增加层数),更多的模型(比如重复5遍,再Ensemble结果),这样一般能稳定达到99%以上的正确率。

目前能达到99.486%正确性的模型是:把已有的42k数据按照80%, 20%的比例分成training,testing数据,然后用Lenet(120 convolution layers [kernel = 5] + ReLU + MaxPooling + 200 convolution layer [kernel = 3] + ReLU + MaxPooling + InnerProduct [param = 200] + InnerProduct [param=10]), 迭代60,000次。然后这个过程重复15次,然后选多数Ensemble。

最后放几个神经网络识别起来费劲的图案供欣赏,还有一个介绍性质的Slides。

 

DigitRecognizerRank
DigitRecognizerRank

probably 4 looks like 9 3 or 2 7 or 2

不要用C++11特性开发R包

不要用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++代码吧。

安装GCC 4.9

安装GCC
Install GCC

这是我找到的最简单的安装gcc,并且不需要root权限的方法 (https://gcc.gnu.org/wiki/InstallingGCC

这个方法适用与GCC的多个版本(4.6,4.7,4.8,4.9)

tar xzf gcc-4.6.2.tar.gz
cd gcc-4.6.2
./contrib/download_prerequisites
cd ..
mkdir objdir
cd objdir
$PWD/../gcc-4.6.2/configure --prefix=$HOME/gcc-4.6.2
make
make install

第七行:如果只需要64位编译器,不需要32位的,可以加上–disable-multilib
第八行:可以用-j N参数指定N个进程同时编译

    class='wp-pagenavi' role='navigation'>
  • Page 1 of 9
  • 1
  • 2
  • 3
  • 4
  • 5
  • ...
  • Last »