R的诡异错误

zhanxw Leave a comment

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要好一点。

Create dynamics DNS

zhanxw Leave a comment

自动更新IP
Create dynamics DNS

办公室电脑的IP经常变动。为了能远程登录，我需要知道它的IP是什么。
我的解决方法是写一个小程序，这个程序把办公室电脑的IP绑定到域名 zhanxw.x64.me
然后，我让办公室电脑每小时运行这个脚本一次，这样zhanxw.x64.me就能指向办公室电脑了。

这里x64.me域名是由DNSDynamic (www.dnsdynamic.org)提供的.

下面是脚本，注意使用的时候要替换username和password。

#!/usr/bin/env python
import socket
ip=str(socket.gethostbyname(socket.gethostname()))
print ip

import urllib2, base64
request = urllib2.Request("http://www.dnsdynamic.org/api/?hostname=zhanxw.x64.me&myip=%s" % ip)
username = "######"
password = "######"
base64string = base64.encodestring('%s:%s' % (username, password)).replace('\n', '')
request.add_header("Authorization", "Basic %s" % base64string)   
result = urllib2.urlopen(request)

data = result.read()
print data

对于Windows主机，可以用计划任务每小时运行脚本一次。对于Linux主机，可以用crontab。

gcc参数顺序

zhanxw 2 Comments

gcc参数顺序
Ordering of gcc parameters matters

2013年最后一个月最后一个星期谈一个奇怪的bug，我在尝试pinfo的时候，configure脚本会提示ncurses有问题：

checking location of curses.h file... /usr/include/ncurses.h
checking if curses is usable... no
configure: error: Curses not found. You need curses to compile pinfo

打开config.log文件，发现这个提示：

configure:12563: checking if curses is usable
configure:12587: gcc -o conftest -g -O2 -I/usr/include    -L/usr/lib -lncurses conftest.c  >&5
/tmp/ccBLkSqd.o: In function `main':
/net/fantasia/home/zhanxw/software/pinfo-0.6.10/conftest.c:38: undefined reference to `initscr'
/net/fantasia/home/zhanxw/software/pinfo-0.6.10/conftest.c:39: undefined reference to `printw'
/net/fantasia/home/zhanxw/software/pinfo-0.6.10/conftest.c:40: undefined reference to `stdscr'
/net/fantasia/home/zhanxw/software/pinfo-0.6.10/conftest.c:40: undefined reference to `wrefresh'
/net/fantasia/home/zhanxw/software/pinfo-0.6.10/conftest.c:41: undefined reference to `stdscr'
/net/fantasia/home/zhanxw/software/pinfo-0.6.10/conftest.c:41: undefined reference to `wgetch'
/net/fantasia/home/zhanxw/software/pinfo-0.6.10/conftest.c:42: undefined reference to `endwin'
collect2: ld returned 1 exit status

从这个命令行来看，curses库的头文件目录和库文件目录都对，但编译的时候（链接/tmp/ccBLkSqd.o）却找不到几个函数的定义（比如：initscr，printw等等）。为了重现这个看起来不该出现的问题，我把conftest.c创建好，然后用相同的gcc命令行编译代码，的确是可以重现错误。

更令我不解的是如果把conftest.c放到命令行的末尾，这个编译就能通过了：

# This does not work
gcc -o conftest -g -O2 -I/usr/include    -L/usr/lib -lncurses conftest.c 
# This works
gcc -o conftest conftest.c -g -O2 -I/usr/include    -L/usr/lib -lncurses

为什么把conftest.c移到“-L/usr/lib -lncurses”前面就能解决问题？
看到StackOverFlow上的一个帖子：

这个帖子讲到gcc编译和链接一个源程序时，会保留参数的顺序，也就是说在链接时，把conftest.c放到前面和后面，效果不同。具体来讲：

conftest.c 放到后面

gcc -o conftest -g -O2 -I/usr/include    -L/usr/lib -lncurses conftest.c 
<=> 等价于
gcc -o conftest -g -O2 -L/usr/lib -lncurses tmporary_object_file.o

conftest.c 放到前面

gcc -o conftest -g -O2 -I/usr/include conftest.c -L/usr/lib -lncurses 
<=> 等价于
gcc -o conftest -g -O2 tmporary_object_file.o -L/usr/lib -lncurses

也就是说，gcc在链接目标文件时，会依照命令行的顺序，对于每个在前面目标文件中没有定义的函数，都试着在其后面的目标文件找找，如果最后一个目标文件也找不到，就报错。因此在我们的例子里，tmporary_object_file.o一定要放到-L/usr/lib -lncurses前面才能正确编译。

回到本文最开始，我们如何让configure把conftest.c放到“-L”前面？答案不是改动configure，而是使用LIBS环境变量：

LIBS="-L/usr/lib -lncurses" ./configure

与LIBS不同，LDFLAGS会放到conftest.c前面，这样configure还是没法正确的使用configure。

在本文最后，回想我们本来的问题，就是说gcc的参数有一定的顺序要求。也许gcc的手册里写了顺序的要求，但是很少有人能够看完gcc的手册。如果见到类似问题，应该怎么办呢？我觉得可以有两个思路，第一个就是想想写gcc软件的人是怎么想的。比如现在gcc的选择就是一个很方便的方法，如果当前没某个函数没有定义，就在他后面制定的目标文件里找找。如果这个要求反过来，或者任何顺序都支持，那么链接的函数就会复杂（因为链接程序需要存储所有定义过的函数，参见链接）；第二个就是想想传统的（或者标准的）写法是什么，把已知的经验和这个错误联系起来。比如configure的错误看起来很奇怪，但是先按照传统的写法先编译后链接，找到一个能凑合工作的方案（把conftest.c放到前面），再逐步改变到出错时的命令行，这个逐步的过程可以帮助查错也能巩固一下已有的知识。

安装Theano

zhanxw Leave a comment

安装Theano
Install Theano

Theano是一个Machine Learning的工具箱，因为和Deep Learning相关，现在非常流行。

这篇文章简单记录一下安装Theano是遇到的问题。

1. import theano出错

>>> import theano
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/net/fantasia/home/zhanxw/python/lib/python2.7/site-packages/Theano-0.6.0rc3-py2.7.egg/theano/__init__.py", line 76, in <module>
    import theano.tests
  File "/net/fantasia/home/zhanxw/python/lib/python2.7/site-packages/Theano-0.6.0rc3-py2.7.egg/theano/tests/__init__.py", line 6, in <module>
    import unittest_tools
  File "/net/fantasia/home/zhanxw/python/lib/python2.7/site-packages/Theano-0.6.0rc3-py2.7.egg/theano/tests/unittest_tools.py", line 8, in <module>
    import theano.tensor as T
  File "/net/fantasia/home/zhanxw/python/lib/python2.7/site-packages/Theano-0.6.0rc3-py2.7.egg/theano/tensor/__init__.py", line 10, in <module>
    import blas
  File "/net/fantasia/home/zhanxw/python/lib/python2.7/site-packages/Theano-0.6.0rc3-py2.7.egg/theano/tensor/blas.py", line 393, in <module>
    StrParam(default_blas_ldflags()))
  File "/net/fantasia/home/zhanxw/python/lib/python2.7/site-packages/Theano-0.6.0rc3-py2.7.egg/theano/tensor/blas.py", line 332, in default_blas_ldflags
    blas_info = numpy.distutils.__config__.blas_opt_info
AttributeError: 'module' object has no attribute '__config__'

加入下面几行可以解决问题：

import numpy
import numpy.distutils
import numpy.distutils.__config__

2.Theano无法编译

运行示例程序是遇到编译错误：

Traceback (most recent call last):
  File "logistic_sgd.py", line 373, in <module>
    sgd_optimization_mnist()
  File "logistic_sgd.py", line 271, in sgd_optimization_mnist
    y: test_set_y[index * batch_size: (index + 1) * batch_size]})
...
  File "/net/fantasia/home/zhanxw/python/lib/python2.7/site-packages/Theano-0.6.0rc3-py2.7.egg/theano/gof/cmodule.py", line 1641, in compile_str
    return dlimport(lib_filename)
  File "/net/fantasia/home/zhanxw/python/lib/python2.7/site-packages/Theano-0.6.0rc3-py2.7.egg/theano/gof/cmodule.py", line 263, in dlimport
    rval = __import__(module_name, {}, {}, [module_name])
ImportError: ('/net/fantasia/home/zhanxw/.theano/compiledir_Linux-3.4.61-x86_64-with-debian-wheezy-sid-x86_64-2.7.5/tmphSYif1/baf8a6fdee34b135ea9110785e356489.so: undefined symbol: ATL_dptsyrk', '[_dot22(<TensorType(float64, matrix)>, W)]')

这里的提示信息是找不到ATL_dptsyrk函数。在没有管理员权限的服务器上，Theano找不到合适的ATLAS函数库。
解决方法是加入以下几行到 $HOME/.theanorc

[blas]
ldflags = -lf77blas -latlas -lgfortran

对于其他错误来讲，还有一些更一般的解决方法。

1. 可以打开Theano的日志信息，比如在程序中加入这两行：

import logging                                                                  |
#logging.getLogger("theano.gof.cmodule").setLevel(logging.DEBUG)

Theano在运行是可以输出更多的信息。

2. 保留临时文件
可以在 $HOME/.theanorc 中加入以下两行

[global]
nocleanup = False

这样可以保留Theano编译过程中的临时文件，便于调试。

题外话

在有管理员权限的Ubuntu系统， 可以用OpenBLAS来代替ATLAS，因为OpenBLAS是多线程的，而且速度更快。具体可以参考Theano 的手册。

现在Theano可以正常使用，希望接下来几篇blog有空能写写学习Deep Learning的心得。

R语言如何用索引

zhanxw Leave a comment

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

Ridge regression

zhanxw Leave a comment

Ridge regression 相关的资料

以前学Ridge regression，大概记住了如何求回归系数 \(\beta^{ridge}\)，但有些细节不了解。我找到下面两个slide感觉又学到了新东西。

 

第一个讲义 【Link or RidgeRegression – Ryan Tibshirani】

http://www.stat.cmu.edu/~ryantibs/datamining/lectures/16-modr1.pdf

Slide 6：Shrink之后的回归系数看起来什么样？

Slide 7：为什么要center和normalize predictors？

Slide 13：当真正的回归系数很大的时候，Ridge regression还有好的表现么？或者说当true coefficient很大的时候，应该向什么方向调整\(\lambda\)可以减少prediction error？

 

第二个讲义 【Link or RidgeRegression-Rudyregularization】

http://www-stat.stanford.edu/~owen/courses/305/Rudyregularization.pdf

Slide 17：Ridge regression的回归系数的bias具体是多少？

Slide 23-24：Ridge regression和Principal component analysis （PCA）的关系是什么？

Slide 37：选择Ridge regression的tuning parameter时一般用cross validation。如果用GCV（Generalized cross validation）的话，可以有close form.

矩阵最大特征值的求法

zhanxw Leave a comment

矩阵最大特征值的求法

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的。

 

 

对RAII的思考

zhanxw Leave a comment

对RAII的思考

Some random thoughts about RAII

写程序时常常需要申请系统资源，比如打开文件，申请一块内存。申请到这些资源后，在程序退出或者资源使用完毕后，应当正确的释放。如果不能正确释放，会造成一系列问题。比如申请的内存没有释放造成内存泄露Memory Leak，申请的进程锁没有被解锁Unlock，造成进程间的死锁DeadLock。在C++语言里，解决这类资源管理问题的管用手法是RAII (Resource Acuiquistion Is Initialization)

这篇笔记是对RAII的一点思考。

 

1. 什么是RAII

简单来讲，把获取资源的代码放到类的构造函数里，把释放资源的代码放到析构函数里。比如用ofstream file(“output.txt”) 可以打开文件，当file变量不起作用是，文件会被自动关闭。

比如下面这张图（from：The RAII Programming Idiom），看看这里面有多少地方需要写释放资源的代码。如果使用RAII，这些其实地方都不用留代码。

RAII Example
RAII Example

 

 

2. RAII的优缺点

RAII的好处是利用C++语言优势安全、正确的管理资源。同时RAII是C++建议的资源获取方式，这种代码可以被广大C++用户理解。

不方便之处是，使用RAII有一些陷阱。比如不要用RAII一次获取多个资源。

 

3.为什么C++有RAII

C++语言保证了一个类构造之后，析构函数会被自动调用。这个使用方式与资源管理的方式相似。因此可以用类的生命周期来管理资源。

 

4. 为什么C/Java/Python没有RAII

C语言没有原生的构造和析构函数，获取的资源不能有任何自动机制来释放。

Java/Python有语言中的支持，即Dispose Pattern。举例来说就是 try…catch..finally语句。使用者只要把释放资源的语句写到finally，资源就会被释放。

 

4. RAII 和Exception的关系

RAII和Exception紧密相关，更确切的说，构造函数和异常这两个特性在某种程度上互相依赖。

对于构造函数来说（获取资源的语句在构造函数里），构造函数没有返回值，因此想知道资源是否成功获取是不能从函数返回值来判断的，唯一可以用的手法是在资源获取失败时抛出异常。也就是说构造函数需要使用异常。

另一方面，使用异常之后，需要用构造函数来管理资源。因为异常抛出以后，很可能处理异常的代码和异常发生的代码不在一个层次（异常在Call Stack上逐层向上）。为了实现异常安全（Exception Safe），应该使用构造函数（另一个选择是智能指针，但智能指针有智能指针的问题，详见C++FAQ的讨论）。

对已有的C++代码来说，实现或检查代码是不是异常安全不是一个的简单人物。这种情况下，异常这个特性往往会被禁用（比如Google C++ style guide）。如果异常被禁用了，我们就没法从构造函数本身获知资源是否成功获取，那是不是说我们没法使用RAII特性呢？

答案是否定的。我们可以在获取资源后，用其他的类函数来检查资源获取是否成功。比如ostream::is_open()就可以检查文件是否被正常打开。

 

5. 怎么绕开RAII

在C等不提供RAII支持的语言里，可以直接绕开RAII，即保证获取资源后，程序的每一个出口都有释放资源的语句。

这种方法有可能造成多处重复的资源释放代码，或者使用goto语句把所有程序跳转到一处资源释放代码。

 

6. 实践中怎么用RAII

实践中除了把资源获取的语句写到构造函数，把资源释放的语句写到析构函数，还应当注意：

1）获取多个资源时，可以写在多个类的构造函数里，使得每一个类的构造函数对应一个资源。这样在任何资源获取失败时，已经获得的资源会得到释放

2）有时候获取资源失败等于程序失败（Fatal Condition），这种情况下可以直接退出（exit），把清理资源的任务留给操作系统。

 

给作者编号

zhanxw Leave a comment

给作者编号

Author Number

http://zhanxw.com/author

 

最近写了一个网页，主要解决写文章给作者编号的问题。这个问题怎么来的呢？现在的文章需要提供作者，作者工作单位和作者的贡献，但是提供的方式是给每个作者单位按照作者出现的顺序编号。比如我的文章有79个作者，现在突然要加一个新的合作者，把他拍到第50位，那么从第51位到最后一位作者的编号可能都需要改变。这个过程很容易出错。我想写个软件给所有人用，顺便凑个热闹，学习一下jQuery + Django，最后的作品放到我的网页(link)。在这个过程中，还学习了Nginx, Emacs web-mode, 这篇笔记主要是把我的经验记录下来。

 

1. 学习Django

Django是Python语言实现的Web架构，它最开始是用于展示新闻（Newroom），国内的豆瓣也用Django。作为Python的“粉丝”，我学一下Django，希望以后可以把有用的工具放到网上。

想入门，最好最省时间的方法是Django tutorial，就在Django的官方网页。这个Tutorial一共有六个部分，介绍的常用Django的功能，包括： 模型（Model），模板（Template），视图（View），静态文件（Static），管理界面（Admin），网址分发（URL Dispatching）等等。学习之后，我感到Django的强项是利用Python语言做到来简化数据库使用。对于一个简单的App，基本开发的流程是： （1）设定project 的settings.py 和urls.py，注册新的app ；（2）设定urls.py来确定网址和视图（View）的关系；（3）用HTML开发模板（/template/app/XXX.html）同时实现各种需要的视图（views.py）

架设Django，一般简易使用uWSGI。我使用的是Nginx处理静态页面，同时用uWSGI处理Django相关界面。简易先在本机用uWSGI调好程序，再放到服务器上并配置Nginx和uWSGI。

Djano默认使用Unicode，而我们一般都用str()，这是ASCII编码。两种编码对于字符串处理来讲（包括join, index）没有什么区别，但在print的时候，Unicode的字符应该先转成UTF-8字符，比如： print u”aaa”.encode(‘UTF-8’)

 

2. 学习jQuery

用jQuery的主要目的是用Ajax来更新网页数据，这样用户用起来有行云流水的感觉。jQuery有.get 和 .post两个方法，使用起来和访问网页很相似。不过这里面有几个陷阱。第一个是GET方式适合小数据，对于Django来讲，默认的大小是4096个字符，大一点的数据应该用POST；第二个陷阱就是POST，为了安全（CSRF: cross-site request forgery），Django要求POST的数据必须有csrftoken，一般的HTML表格Form必须有csrf的标记才会被Django接收。

为了调试Ajax的输入输出，简易用Firefox下的FireBug或者最新版本的Firefox，它们都可以显示Ajax请求的参数和返回值。最新版本的Firefox还可以给Javascript设置断点，这就更加降低了调试的难度。

jQuery的功能很多，这个网页把jQuery 1.9以及之前版本的功能用不同颜色区分出来，很方便查询（link）。

在这个网页中，我还用了handsontable，这是一个在jQuery的基础上开发的实用的javascript spreadsheet，语法简单，使用方便。

3. 一点Nginx经验

对Nginx来讲，我们的配置是用location语句来把特定的网址传给uWSGI进程。注意Nginx设置里如果有多个location语句，Nginx并不一定选择最先匹配的网址，而是选择最长的网址。比如：

location ~* /author/static/.+.(jpg|jpeg|gif|css|png|js|ico|xml)$ {
root /var/django/zhanxw/; # STATIC_ROOT
expires 30d;
}

location ~* ^.+.(jpg|jpeg|gif|css|png|js|ico|xml)$ {
access_log        off;
expires           30d;
root /var/www;
}

如果第一个location写成： “location /author/static”，那么.jpg之类的文件会被第二个location处理，这就不是Django static文件的正确处理方式。

 

4. 一点Emacs web-mode经验

用Emacs写Django的模板HTML，最好用的不是django-html-mode，而是web-mode。只有这个mode可以识别 “{% static ‘polls/index.html’ %}”这样的记号，并正确缩进。

此外web-mode可以自动补全HTML tab，比如你在<p>后面打</，web-mode会帮你补全</p>。还可以用C-c C-e b 和C-c C-e e跳到一对tag的最前面（beginning）和最后面（end）。

 

5. 一般性经验

最开始设计应该以最少功能，最小实现为好，不要一下子把界面设计复杂。可以想几个用例（Use Case），保证最基本最重要的功能，其他功能应该越少越好。

网页的布局应该少用<br/>这种硬回车。在有Bootstrap的情况下，完全可以用<div>和<p>来用更少的HTML tag来实现更自然、自动的布局。