R的诡异错误

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

自动更新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。

安装Theano

安装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的心得。

给作者编号

给作者编号

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来实现更自然、自动的布局。

 

PyCon2013 有意思的幻灯

Interesting slides from PyCon 2013
今年PyCon在加州的Santa Clara召开。我虽然没去,但一如既往的关心。
挑出和我相关的一些有意思的幻灯,在此分享。

1. BeautifulCode
Raymond Hettinger是一个善用Python的高手,他的code非常简洁,有Python的风味。
这个幻灯里,他介绍很多Python中常用的简洁的写法,包括怎么用iterator, list comprehension.
我感觉有意思的是defaultDict (不需要再用dict.get()), decorator(修饰方法,类比CSS)和context(干净的获取和释放资源)

原始链接

2. Python: A “Toy” Language

David Beazley是一个Python作家,对Python GIL有过详细的介绍,以前写过如何把Generator组合成一个workflow的幻灯。
现在在芝加哥教Python。这个人擅长Python的教学,并能给出有趣的例子。这次也不例外,他介绍了如何用Python和Shapeoko (包括Arduino)来组装并使用CNC (数控机床?)。这个例子告诉我们Python可以做计算机程序之外很有用的 应用。

原始链接

3. Awesome Big Data Algorithms

作者是MSU的老师。这个Blog的背景是土壤Genetics。因为土壤中的genetics比单纯人的DNA更复杂,数据量更大,因此需要Big Data Algorithm。这个幻灯介绍三种算法:SkipList, HyperLogLog, Bloom filter(CountMin Sketch)。

SkipList是一种基于链表的数据结构,相对羽平衡二叉树,这个算法的优点是更好的支持并发操作。本质上,SkipList是一个分层次的链表。在最底层,链表元素按顺序排列。在更高的层次,(按照概率)简历一部分低层的数据的索引。这种数据结构在查找时非常有效:从高层开始查找,直到最底层顺序查找,整个查找是log(N)

HyperLogLog是应用于大数据的算法,用来计算一个很大集合的基数(即合理总共有多少不相同的元素)。大致思路是用一组相互独立的哈希函数依次处理输入,然后对哈希值分块计数:对高位统计有多少连续的0;用低位的值当做数据块。比如:011000|01,就是高位有3个连续的0,低位是1,就表示第一个数据块。因为连续观测的三个0的概率大约是1/8,所以对数据块1来讲,可以把计数乘8,作为集合基数的估计。因为低位可能有0,1,2,3这四种数值,总基数可以取上述4中技术的几何平均数。在HyperLogLog中,具体的还有一些系数可以调整,使得估计更准确。
这片Blog详细介绍了HyperLogLog算法,图文并茂。

BloomFilter或CountMin Sketch是两个不同的算法,但又有紧密联系。相似之处是两个算法都需要一族独立的哈希函数。不同之处是处理的问题不同。对BloomFilter,在预处理阶段对每一个特定的输入算出所有哈希函数的值,并在这些值上做出标记。最后,当查找一个特定的输入是否出现过,只需查找这一系列的哈希函数对应值上有没有标记。对于BloomFilter,可能有False Positive,但不可能有False Negative。此外,BloomFilter可看做查找一个数据有或者没有的数据结构(数据的频率是否大于1)。CountMin Sketch在BloomFilter的基础上更进一步,它可用来估算某一个输入的频率(不局限于大于1)。具体思路是对哈希函数值对应的sketch上计数(对BloomFilter则只是标记是或否)。最后估计频率的时候,对每个估算出的频率取最小值。

原始链接

4. Why you should use Python 3 for text processing

这个讲座关注的Python3,而不是现在我使用的Python2.7.3。但在Python社区,有时好的功能会从版本3移植回版本2。
这个讲座介绍了Python3里面的新功能,例如ChainMap,startswith(tuple作为参数),unicode支持,textwrap模块(可以方便的排版)和email模块。

原始链接

用pip安装Matplotlib (Mac)

Install matplotlib using pip (Mac)

在Mac OSX下安装matplotlib并不简单 (比如Linux可以用apt-get,Windows可以直接下载Binary build),因为matplotlib需要一些底层库(例如freetype,pygtk,而pygtk又需要gobject, gobject需要cairo)。

这里总结一下我在Mac OSX下安装Matplotlib的经验。

首先安装MacPorts,具体参见Install MacPorts.

然后安装py27-pygtk

 port install py27-pygtk

然后把pygtk的路径添加到PYTHONPATH里:

export PYTHONPATH=$PYTHONPATH:/opt/local//Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/

最后用pip安装

alias pipInstall='pip install  --install-option="--prefix=/Users/zhanxw/python" '
export CC=/opt/local/bin/g++; export CFLAGS=; export LDFLAGS= ; pipInstall -U matplotlib

这里要指定GCC编译器的路径,否则系统默认用gcc 4.2版编译。
还需要清除CFLAG和LDFLAGS,不然安装的时候不用pkg-config,也就得不到gtk库文件的位置。

额外两个tip:

【1】 安装pip需要先安装setuptools
【2】 一个命令升级所有pip安装的软件:

  pip freeze --local | grep -v '^\-e' | cut -d = -f 1  | xargs pip install -U

此外如何在Mac OSX下用源码安装numpy和scipy?
这里有官网的说明
特别注意的是,编译numpy和scipy要用gcc 4.0版。
另外scipy还需要cython (pip install cython)

发布Python程序

发布Python程序

Distribute Python script

最近需要把一个Python程序发布给其他学校的Collaborators使用,发现最主要的问题是版本不兼容:我使用的是Python 2.7版本,但其他学校可能使用2.5或者更老的2.4版本。因为我使用了大量新版的特有功能,直接把代码发布给对方是没法让对方使用的。因为老版本的Python不支持2.7版里的函数比如:os.SEEK_SET, all(), str.format() 。解决方法有两个:

1. 使用Freeze之类的工具

在Python source code里有一个Tools文件夹,里面有freeze这个工具。他可以把Python代码编译成C语言中的Python。这种方式可以解决简单的脚本,但复杂的脚本这种方法可能会带来其他的错误。

2. 重写代码

这是最根本的解决方法,把新版本特有的函数重新定义,或者用另一种方法实现。比如all()

# all is a keyword since Python 2.7                                                                                                                                                  
try:
    all
except:
    def all(iterable):
        for element in iterable:
            if not element:
                return False
            return True

有一个工具可以检查当前的Python的向前兼容性:

https://github.com/ghewgill/pyqver

检查checkVCF.py的结果如下所示:

/net/fantasia/home/zhanxw/rvtests/rvMeta/sftp.clean/checkVCF/checkVCF.py
2.5 all
2.4 set, generator expression
2.3 logging, sum, enumerate

就是说我使用了2.5特有的all()函数,如果想在2.4版的Python上运行,需要重新定义all()函数。
类似的,set(), generator expression是2.4版才加入的特性,如果要支持2.3版本,就必须改写这两个地方。

 

 

如何在Python中调用C/C++代码

如何在Python中调用C/C++代码
How to mix C/C++ code in Python

本文介绍一种手动的、简单的在Python中使用C/C++代码的方式。这个方法主要使用了ctypes模块。其他的混合Python,C/C++编程的方法还有Swig Boost.Python。前一种方法需要写一个接口文件(interface),而后一种需要使用庞大、深奥的boost类库,后两者适合可能适合更复杂的情况,这里只介绍第一种方法。

混合C/C++代码需要这几步:
1. 包装接口 C/C++ wrap functions up
2. 打包成共享库 Compiling C/C++ code and pack it to shared library
3. Python中导入共享库 Python imports shared library

先介绍一下北京,这里我的C++类GenomeSequence使用了模板(Template)和Memorymap,这是一个访问基因序列的类,比如如果一个生物序列是GAGTTTTATCGCTTCCATGACGCAGAAGTTAACACT… 我们的类是gs,那么gs[0] = ‘G’, gs[1]=’A’ …. 摘录相关的函数如下:

class GenomeSequence : public genomeSequenceArray
{
public:
    /// Simple constructor - no implicit file open
    GenomeSequence();
    /// set the reference name that will be used in open()
    /// \param referenceFilename the name of the reference fasta file to open
    /// \return false for success, true otherwise
    ///
    /// \sa open()
    bool setReferenceName(std::string referenceFilename);
    /// return the number of bases represented in this reference
    /// \return count of bases
    genomeIndex_t   getNumberBases() const
    {
        return getElementCount();
    }
    inline char operator[](genomeIndex_t index) const
    {
        uint8_t val;
        if (index < getNumberBases())
        {
            if ((index&1)==0)
            {
                val = ((uint8_t *) data)[index>>1] & 0xf;
            }
            else
            {
                val = (((uint8_t *) data)[index>>1] & 0xf0) >> 4;
            }
        }
        else
        {
            val = baseNIndex;
        }
        val = isColorSpace() ? int2colorSpace[val] : int2base[val];
        return val;
    }
    /* ........... more codes omitted ................ */
}

但实际上这些细节并不重要,重要是如何包装,我们编写GenomeSequence_wrap.cpp文件,包括对上述4个函数的封装,源码如下:

#include "GenomeSequence.h"
#include <string>

extern "C"{
    GenomeSequence* GenomeSequence_new(){ return new GenomeSequence();}
    bool GenomeSequence_setReferenceName(GenomeSequence* gs, char* s) { 
        if (!gs) return false;
        std::string str = s;
        //printf("Loading %s ...\n", s);
        if (!gs->setReferenceName(str)){
            gs->open();
        } else {
            printf("Loading FAIL\n");
        }
        return (gs->setReferenceName(str));
    }
    void GenomeSequence_close(GenomeSequence* gs) {if (gs) gs->close();};
    int GenomeSequence_getNumBase(GenomeSequence* gs) {
        if (!gs) {
            printf("invalid gs\n");
            return -1;
        }
        return (gs->getNumberBases());
    }
    char GenomeSequence_getBase(GenomeSequence* gs, unsigned int i) { 
        if (gs) {
            return (*gs)[i];
        };
    };
}

第二步是编译,记住单个C/C++文件编译时使用-fPIC参数,最后打包的时候编译成共享库,摘录Makefile文件中片段如下:

lib:
	g++ -c -fPIC -I./lib GenomeSequence_wrap.c
	g++ -shared -Wl,-soname,libstatgen.so -o libstatgen.so  lib/*.o lib/samtools/k*.o lib/samtools/bgzf.o *.o

最后一步是在Python中写一个封装类,注意前两行引入ctypes库,之后就用这个库调用包装函数就行。
注意:我在GenomeSequence类的__getitem__中使用了如何扩展Python的容器类一文中介绍的一些技巧,这样可以更灵活的使用下标来访问数组中的元素。

from ctypes import cdll
lib = cdll.LoadLibrary("./libstatgen.so")

class GenomeSequence:
    def __init__ (self):
        self.obj = lib.GenomeSequence_new()
    def open(self, filename):
        lib.GenomeSequence_setReferenceName(self.obj, filename)
    def __len__ (self):
        return lib.GenomeSequence_getNumBase(self.obj)
    def __getitem__(self, key):
        if isinstance(key, int):
            return chr(lib.GenomeSequence_getBase(self.obj, key))
        elif isinstance(key, slice):
            return ''.join([self[x] for x in xrange(*key.indices(len(self)))])
        elif isinstance(key, tuple):
            return ''.join([self[i] for i in key])

    def at(self, i):
        return chr(lib.GenomeSequence_getBase(self.obj, i))
    def close(self):
        lib.GenomeSequence_close(self.obj)
    
if __name__ == '__main__':
    gs = GenomeSequence ()
    gs.open("/home/zhanxw/statgen/src/karma/test/phiX.fa");
    print len(gs)
    seq = [(gs.at(i)) for i in xrange(60)]
    print ''.join(seq)
    print gs[0:10],gs[20:30]
    print gs[0:10, 20:30]
    print gs[-10:]
    gs.close()
    print "DONE"

本文主要参考【1】。这里的方法基本重复了【1】中的步骤。写出本文中的代码在于进一步验证ctypes库可以灵活的处理C/C++和Python中的简单数据类型int, char*。

【1】Calling C/C++ from python?

如何扩展Python的容器类

如何扩展Python的容器类
How to extend Python container class (using some idiom)

本文假设已经有一个C++语言写的array类型的数据结构,可以用v.getBase(unsigned int i) 来得到v 数组在下标i的数值。我们想利用Python灵活的slice功能,比如1:10, 1:10:2, -10:-5等方式来指定不同的下标。这种灵活的下标在Python中可以有三种形式:

1. 整数: v[1]
2. slice 对象: v[1:10]
3. tuple 对象: v[1:10, 20:30]

这三种对象都会被传到__getitem__(self, key)的key参数中。通过参考【1】,【2】,我发现下面的代码可以简洁的处理上述所有情况:

 
class ContainerClass:
    def __getitem__(self, key):
        if isinstance(key, int):
            return chr(v.getBase(self.obj, key))
        elif isinstance(key, slice):
            return ''.join([self[x] for x in xrange(*key.indices(len(self)))])
        elif isinstance(key, tuple):
            return ''.join([self[i] for i in key])

注意:
这里只是代码片段。全部代码见另一片Blog:如何在Python中调用C/C++代码

【1】Python Data Model:
http://docs.python.org/reference/datamodel.html
【2】Python in a nut shell:
http://books.google.com/books?id=JnR9hQA3SncC&pg=PA110&lpg=PA110&dq=python+slice+object+idiom&source=bl&ots=Jb1XIv_71t&sig=-_NHkwycfC8yipkc4Tl_e4sruKc&hl=en&ei=uRXfTcr5Jsro0QHa4sG5Cg&sa=X&oi=book_result&ct=result&resnum=10&ved=0CF0Q6AEwCQ#v=onepage&q=python%20slice%20object%20idiom&f=false