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??
```

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.
```

Create dynamics DNS

Create dynamics DNS

```#!/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)
result = urllib2.urlopen(request)

print data
```

安装Theano

Install Theano

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

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)]')
```

```[blas]
ldflags = -lf77blas -latlas -lgfortran
```

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

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

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

2. 保留临时文件

```[global]
nocleanup = False
```

给作者编号

Author Number

http://zhanxw.com/author

1. 学习Django

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

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

2. 学习jQuery

3. 一点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;
}

4. 一点Emacs web-mode经验

5. 一般性经验

PyCon2013 有意思的幻灯

Interesting slides from PyCon 2013

1. BeautifulCode
Raymond Hettinger是一个善用Python的高手，他的code非常简洁，有Python的风味。

David Beazley是一个Python作家，对Python GIL有过详细的介绍，以前写过如何把Generator组合成一个workflow的幻灯。

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

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

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

用pip安装Matplotlib (Mac)

Install matplotlib using pip (Mac)

``` port install py27-pygtk
```

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

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

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

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

发布Python程序

Distribute Python script

1. 使用Freeze之类的工具

2. 重写代码

```# 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
```

https://github.com/ghewgill/pyqver

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

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

How to mix C/C++ code in Python

1. 包装接口 C/C++ wrap functions up
2. 打包成共享库 Compiling C/C++ code and pack it to shared library
3. Python中导入共享库 Python imports shared library

```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 ................ */
}
```

```#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;
if (!gs->setReferenceName(str)){
gs->open();
} else {
}
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];
};
};
}
```

```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
```

```from ctypes import cdll

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"
```

如何扩展Python的容器类

How to extend Python container class (using some idiom)

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

```
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])
```

【1】Python Data Model:
http://docs.python.org/reference/datamodel.html
【2】Python in a nut shell: