yablog

python tips: set up IPython Notebook on Ubuntu12.04LTS

noreply@blogger.com (Anonymous) — Sun, 04 Nov 2012 11:58:00 +0000

set up pythonbrew, which makes it easier to maintain multiple python binaries.

# install curl for downloading pythonbrew package
$ sudo apt-get install curl

# install packages for compiling python2.7.3 with pythonbrew
$ sudo apt-get install build-essential, libbz2-dev, libsqlite3-dev, zlib1g-dev, \
                       libxml2-dev, libxslt1-dev, libreadline5, libreadline-dev, \
                       libgdbm-dev, libgdb-dev, libxml2, libssl-dev, tk-dev, \
                       libgdbm-dev, libexpat1-dev, libncursesw5-dev

# install pythonbrew
$ curl -kLO https://github.com/utahta/pythonbrew/raw/master/pythonbrew-install
$ chmod +x pythonbrew-install
$ ./pythonbrew-install

# edit .bashrc then, 
$ source ${HOME}/.pythonbrew/etc/bashrc

# setup python-2.7.3 with pythonbrew
$ pythonbrew install 2.7.3

# change default python binary
$ pybrew switch 2.7.3

set up virtualenv, which provides virtual pythons' executable environment

# install virtualenv
$ pip virtualenv virtualenvwrapper
$ echo "VIRTUALENVWRAPPER_PYTHON=${HOME}/.pythonbrew/pythons/Python-2.7.3/bin/python" >> ~/.bashrc
$ source ${HOME}/.pythonbrew/pythons/Python-2.7.3/bin/virtualenvwrapper.sh

set up ipython notebook.

# get package for pyqt compile
$ sudo apt-get install libqtcore4
$ sudo apt-get install libqt4-dev

# download SIP package from http://www.riverbankcomputing.co.uk/software/sip/download
# install SIP following README
 
# download PyQt package from http://www.riverbankcomputing.co.uk/software/pyqt/download 
# install SIP following README

# install numpy and scipy
$ sudo apt-get install liblapck3gf liblapack-dev libblas3gf libblas-dev gfortran
$ pip install numpy
$ pip install scipy

# install matplotlib
$ sudo apt-get install libpng12-dev libjpeg8-dev libfreetype6-dev
$ pip install matplotlib

# install pyzmq
$ pip install pyzmq

# install tornado
$ pip install tornado

# install pygments
$ pip install pygmenta

# install pandas
$ pip install pandas

# install sympy
$ pip install sympy

# install nose
$ pip install nose

# install ipython
$ pip install ipython
$ easy_install ipython[test]
$ easy_install ipython[zmq,qtconsole,notebook,test]
$ iptest

# install packages for ipython extension
$ sudo apt-get install r-base-core ruby
$ pip install rpy2
$ pip install rtools

This is just a test of nbconvert.py

noreply@blogger.com (Anonymous) — Tue, 30 Oct 2012 06:14:00 +0000

In [1]:

%pylab inline
x = linspace(0, 2*pi)

Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.zmq.pylab.backend_inline].
For more information, type 'help(pylab)'.

In [2]:

plot(x, sin(x), label=r'$\sin(x)$')
plot(x, cos(x), 'ro', label=r'$\cos(x)$')
title(r'Two familiar function')

Out [2]:

<matplotlib.text.Text at 0x105507cd0>

In [3]:

%load_ext sympyprinting

/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/nose/plugins/manager.py:418: UserWarning: Module dateutil was already imported from /Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/dateutil/__init__.pyc, but /Users/yaboo/src/matplotlib/lib is being added to sys.path
  import pkg_resources
/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/nose/plugins/manager.py:418: UserWarning: Module six was already imported from /Users/yaboo/src/matplotlib/lib/six.pyc, but /Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages is being added to sys.path
  import pkg_resources
/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/IPython/extensions/sympyprinting.py:119: UserWarning: The sympyprinting extension in IPython is deprecated, use sympy.interactive.ipythonprinting
  warnings.warn("The sympyprinting extension in IPython is deprecated, "

In [4]:

import sympy as sym
from sympy import *
x, y, z = sym.symbols("x y z")

In [5]:

Rational(3,2)*pi + exp(I*x)/(x**2 + y)

Out [5]:

$$\frac{3}{2} \pi + \frac{e^{\mathbf{\imath} x}}{x^{2} + y}$$

In [6]:

eq = ((x+y)**2 * (x+1))
eq

Out [6]:

$$\left(x + 1\right) \left(x + y\right)^{2}$$

In [7]:

expand(eq)

Out [7]:

$$x^{3} + 2 x^{2} y + x^{2} + x y^{2} + 2 x y + y^{2}$$

In [8]:

diff(cos(x**2)**2 / (1 + x), x)

Out [8]:

$$- 4 \frac{x \sin{\left (x^{2} \right )} \cos{\left (x^{2} \right )}}{x + 1} - \frac{\cos^{2}{\left (x^{2} \right )}}{\left(x + 1\right)^{2}}$$

You can italicize, boldface

build

list

and embed code meant for illustration instead of execution in Python:

def f(x):
    """a docstring"""
    return x**2

or other languages:

In [9]:

from IPython.display import Image
Image(filename=u'/Users/yaboo/Downloads/hoge.jpeg')

Out [9]:

In [10]:

from IPython.display import YouTubeVideo
YouTubeVideo('iwVvqwLDsJo')

Out [10]:

In [11]:

%%ruby
puts "Hello from Ruby #{RUBY_VERSION}"

Hello from Ruby 1.8.7

In [12]:

%%bash
echo "hello from $BASH"

hello from /bin/bash

In [13]:

%load_ext rmagic
X = np.array([0,1,2,3,4])
Y = np.array([3,5,4,6,7])

In [14]:

%%R -i X,Y -o XYcoef
XYlm = lm(Y~X)
XYcoef = coef(XYlm)
print(summary(XYlm))
par(mfrow=c(2,2))
plot(XYlm)

Call:
lm(formula = Y ~ X)

Residuals:
   1    2    3    4    5 
-0.2  0.9 -1.0  0.1  0.2 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)   3.2000     0.6164   5.191   0.0139 *
X             0.9000     0.2517   3.576   0.0374 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.7958 on 3 degrees of freedom
Multiple R-squared:  0.81, Adjusted R-squared: 0.7467 
F-statistic: 12.79 on 1 and 3 DF,  p-value: 0.03739

In [42]:

XYcoef

Out [42]:

[ 3.2  0.9]

In [16]:

load_ext cythonmagic

In [17]:

%%cython -lm
from libc.math cimport sin
print 'sin(1)=', sin(1)

sin(1)= 0.841470984808

R tips: read and visualize shapefile with R

noreply@blogger.com (Anonymous) — Thu, 25 Oct 2012 15:36:00 +0000

Below sample code shows how to read shapefiles with R and plot them.

library(maptools)

# downloaded via http://www.globalmap.org/download/data/shp/1_1/shape_jp_1_1.zip
boundp <- readShapePoly("shape_jp_1_1/bnda_1_1.shp")
transl <- readShapeLines("shape_jp_1_1/transl_1_1.shp")

# Plot shapefile
plot(boundp)
plot(transl)

# plot shapefile within specified boundary
plot(boundp, xlim=c(135, 136), ylim=c(34,35))
lines(transl, xlim=c(135, 136), ylim=c(34,35), col='red')

# plot polygons which name equals 'OSAKA'
boundp_hyogo <- shp[(shp$nam == 'HYOGO'),]
x_range <- bbox(shp_tokyo)[1,]
y_range <- bbox(shp_tokyo)[2,]
plot(boundp_hyogo, xlim=x_range, ylim=y_range)
lines(transl, xlim=x_range, ylim=y_range, col='green')

R tips: operating probability density function in R

noreply@blogger.com (Anonymous) — Wed, 22 Aug 2012 17:02:00 +0000

There are several probability density functions (PDFs) defined in R. You can find defined PDFs with "library(help='stats')" in R-console. Now, let's plot Gaussian distribution N(0,1) with *norm function.

# Probability density function
> curve(dnorm(x, mean=0, sd=1, log=FALSE), -3, 3, n=1001, col = "blue")

# Cumulative probability density function
> curve(pnorm(x, mean=0, sd=1, log=FALSE), -3, 3, n=1001, col = "blue")

# Probabilistic generative model
> rnorm(n=100, mean=0, sd=1) 
  [1]  0.76291389  0.32249468  1.43329768  0.39455640 -0.96082434 -0.03498441
  [7] -0.98776053 -0.69866914  1.31514002 -1.42647303 -0.46658770 -1.32506735
 [13] -1.52074716  0.45169290 -1.11993007 -0.27547321  0.86232800 -2.36687435
 [19] -1.34228574 -0.72800613 -0.41979553 -2.24104445  1.86935264 -0.88038071
 [25] -0.71816560 -0.53788946  0.50974550 -1.23530999  0.77035374  0.06483997
 [31] -0.62251296  1.60636131  0.36695288  0.76654543  0.86672801 -0.75831284
 [37] -1.48413034 -0.37661399  1.34613939 -1.20912506 -0.22231654  2.66835415
 [43]  0.49528498  0.26435606 -0.74862141 -0.41349885 -0.94850251  0.83038086
 [49]  0.29994578 -0.15778865 -1.16049616  2.24109521  0.86357121  1.28825845
 [55]  0.50549591  0.08668917 -0.34392882  2.08620620 -0.02514741 -0.38974372
 [61] -0.07943300  1.01109304  1.26749487  1.65848787  1.40626520  0.13476383
 [67] -0.81541581 -0.23397894 -1.54089193  2.69271063  1.02584790  0.34583390
 [73]  0.35808556 -2.23486925 -0.97331144 -0.75742258 -1.49065873 -0.24702416
 [79]  1.46008300  0.26804260  1.25457656  0.55169794 -0.34509144 -1.24602872
 [85]  0.65070530 -0.69068100  0.98802706 -1.76179226 -0.30577751  0.30197949
 [91] -0.18236187 -0.86313716  1.36986087 -1.23400546 -0.78780711 -2.14942716
 [97]  1.17743473 -0.95787888 -0.27079632  1.79351280

geometry and python tips: automatic generation of Japanese standardized regional meshes

noreply@blogger.com (Anonymous) — Tue, 21 Aug 2012 17:11:00 +0000

Huge number of geometric object takes time to load (Especially, loading standardized regional meshes are time consuming). Therefore, I implemented standardized regional mesh generator which generates meshes according to input boundary box.

#!/usr/bin/env python
#coding: utf-8

import sys
import datetime

from rtree import index
from rtree.index import Rtree
from shapely.geometry import Polygon

class SrmGenerator():
    def get_mesh1st_to_polygon(self, bbox):
        return dict([(mesh1st, Polygon(self._get_rect(*bbox)))
                     for mesh1st, bbox in self.generate_mesh1st_to_bbox(bbox).items()])

    def generate_mesh1st_to_bbox(self, bbox):
        d_mesh1st_to_bbox = self._get_d_mesh1st_to_bbox()
        rindex = self._create_rtree_index(d_mesh1st_to_bbox)
        return dict([(mesh1st, d_mesh1st_to_bbox[mesh1st])
                     for mesh1st in list(rindex.intersection(bbox))])

    def get_mesh2nd_to_polygon(self, bbox):
        return dict([(mesh2nd, Polygon(self._get_rect(*bbox)))
                     for mesh2nd, bbox in self.generate_mesh2nd_to_bbox(bbox).items()])

    def generate_mesh2nd_to_bbox(self, bbox):
        d_mesh1st_to_bbox = self.generate_mesh1st_to_bbox(bbox)
        d_mesh2nd_to_bbox = self._get_d_mesh2nd_to_bbox(d_mesh1st_to_bbox)
        rindex = self._create_rtree_index(d_mesh2nd_to_bbox)
        return dict([(mesh2nd, d_mesh2nd_to_bbox[mesh2nd])
                     for mesh2nd in list(rindex.intersection(bbox))])

    def get_mesh3rd_to_polygon(self, bbox):
        return dict([(mesh3rd, Polygon(self._get_rect(*bbox)))
                     for mesh3rd, bbox in self.generate_mesh3rd_to_bbox(bbox).items()])

    def generate_mesh3rd_to_bbox(self, bbox):
        d_mesh2nd_to_bbox = self.generate_mesh2nd_to_bbox(bbox)
        d_mesh3rd_to_bbox = self._get_d_mesh3rd_to_bbox(d_mesh2nd_to_bbox)
        rindex = self._create_rtree_index(d_mesh3rd_to_bbox)
        return dict([(mesh3rd, d_mesh3rd_to_bbox[mesh3rd])
                     for mesh3rd in list(rindex.intersection(bbox))])

    def get_mesh4th_to_polygon(self, bbox):
        return dict([(mesh4th, Polygon(self._get_rect(*bbox)))
                     for mesh4th, bbox in self.generate_mesh4th_to_bbox(bbox).items()])

    def generate_mesh4th_to_bbox(self, bbox):
        d_mesh3rd_to_bbox = self.generate_mesh3rd_to_bbox(bbox)
        d_mesh4th_to_bbox = self._get_d_mesh4th_to_bbox(d_mesh3rd_to_bbox)
        rindex = self._create_rtree_index(d_mesh4th_to_bbox)
        return dict([(mesh4th, d_mesh4th_to_bbox[mesh4th])
                     for mesh4th in list(rindex.intersection(bbox))])

    def _create_rtree_index(self, d_meshid_to_bbox):
        p = index.Property()
        rindex = index.Index(properties=p)
        [rindex.add(meshid, bbox) for meshid, bbox in d_meshid_to_bbox.items()]
        return rindex

    def _get_d_mesh1st_to_bbox(self):
        l_x = [str(x).zfill(2) for x in range(20, 50)]
        l_y = [str(y).zfill(2) for y in range(30, 70)]
        return  dict([self._generate_mesh1st_bbox_pair(x, y)
                      for x, y in self._get_combinations(l_x, l_y)])
        
    def _generate_mesh1st_bbox_pair(self, x, y):
        id_mesh1st = int("{0}{1}".format(y, x))
        left = float("1{0}".format(x))
        bottom = float(y) / 1.5
        right = float("1{0}".format(str(int(x) + 1)))
        top = float(int(y) + 1) / 1.5
        return (id_mesh1st, (left, bottom, right, top))

    def _get_d_mesh2nd_to_bbox(self, d_mesh1st_to_bbox, num_split=8):
        fn = self._get_meshid_for_2nd_3rd
        return self._get_d_meshid_to_bbox(d_mesh1st_to_bbox, num_split, fn)

    def _get_d_mesh3rd_to_bbox(self, d_mesh2nd_to_bbox, num_split=10):
        fn = self._get_meshid_for_2nd_3rd
        return self._get_d_meshid_to_bbox(d_mesh2nd_to_bbox, num_split, fn)

    def _get_d_mesh4th_to_bbox(self, d_mesh3rd_to_bbox, num_split=2):
        fn = self._get_meshid_for_4th
        return self._get_d_meshid_to_bbox(d_mesh3rd_to_bbox, num_split, fn)

    def _get_d_meshid_to_bbox(self, d_meshcode_to_bbox, num_split, fn):
        D = {}
        for meshcode, bbox in d_meshcode_to_bbox.items():
            left, bottom, right, top = bbox
            x_delta = (right - left) / num_split
            y_delta = (top - bottom) / num_split
            x_coord = [left] + [left+x_delta*i for i in range(1, num_split)] + [right]
            y_coord = [bottom] + [bottom+y_delta*i for i in range(1, num_split)] + [top]
            for i, j in self._get_combinations(range(len(y_coord)-1), range(len(x_coord)-1)):
                new_meshcode = fn(meshcode, i, j)
                bbox = (x_coord[j], y_coord[i], x_coord[j+1], y_coord[i+1])
                D.update({new_meshcode: bbox})
        return D

    def _get_meshid_for_2nd_3rd(self, meshcode, i, j):
        return int(str(meshcode) + str(i) + str(j))

    def _get_meshid_for_4th(self, meshcode, i, j):
        return int(str(meshcode) + self._get_mesh4th_last_code(i, j))

    def _get_mesh4th_last_code(self, y, x):
        if (y == 0 and x == 0):
            return '1'
        elif (y == 0 and x == 1):
            return '2'
        elif (y == 1 and x == 0):
            return '3'
        elif (y == 1 and x == 1):
            return '4'
        else:
            print "unexpected error in _get_mesh4th_last_code"
            sys.exit(0)

    def _get_combinations(self, L1, L2):
        ret = []
        for e1 in L1 :
            for e2 in L2:
                ret.append([e1, e2])
        return ret

    def _get_rect(self, left, bottom, right, top):
        return ((left, bottom),
                (left, top),
                (right, top),
                (right, bottom),
                (left, bottom))
                
if __name__ == "__main__":
    L = 130.0
    R = 140.0
    B = 34.0
    T = 34.2
    bbox = (L, B, R, T)

    from datetime import datetime
    s = datetime.now()
    print len(SrmGenerator().get_mesh4th_to_polygon(bbox))
    e = datetime.now()
    print(e-s)

$ python generate_area_srg.py
80100
0:00:07.602852

geometry and python tips: How to load shapefile correctly using python?

noreply@blogger.com (Anonymous) — Tue, 21 Aug 2012 15:19:00 +0000

Shapely provides an easy way to operate geometrical objects. However, shapely does not have good functionality for loading ESRI shapefile directly and formatting topologically complex geometry for making shapely's method arguments is tired work. Below code will solve this issue: combination of shapely.wkb.loads and osgeo library. To write this code, I referred here.

#!/usr/bin/env python
#coding: utf-8

import sys
import cPickle as pickle
import collections

from osgeo import ogr
from shapely.wkb import loads

def main():
    path_to_shp = argvs[1]
    path_to_out = argvs[2]

    d_fid2pol = get_d_fid2pol(path_to_shp)
    dump_object_to_file(d_fid2pol, path_to_out)

def get_d_fid2pol(path_to_shp):
    """
    :param path_to_shp: path to *.shp
    :type path_to_shp: string

    :returns: {fid: MultiPolygon(), ...}
    :rtype: dictionary
    """
    data = ogr.Open(path_to_shp)
    elements = data.GetLayer()
    return dict([get_fid_and_poly(element) for element in elements])

def get_fid_and_poly(element):
    """
    :param element: element of shapefile
    :type element: object

    :returns: {fid: MultiPolygon()}
    :rtype: element of dictionary
    """
    fid = element.GetField(0)
    poly = loads(element.GetGeometryRef().ExportToWkb())
    return (fid, poly)
    
def dump_object_to_file(obj, path_to_out):
    with open(path_to_out, 'wb') as F:
        pickle.dump(obj, F, -1)

if __name__ == "__main__":
    main()

R tips: first contact with linear pattern recognition in R.

noreply@blogger.com (Anonymous) — Mon, 20 Aug 2012 15:54:00 +0000

This is just test code for linear pattern recognition. Below code is introduced in http://www.slideshare.net/sleipnir002/05-12739580.

> index <- sample(nrow(data), nrow(data)*0.2) #sample data
> index
 [1]  83  47 146  59 111  94 121  14  15 122  33  22  87  78  66 133  18  79  19
[20]   9  88  32 116  60   7 104  67  20  86   1
> data[index,]
    Sepal.Length Sepal.Width Petal.Length Petal.Width    Species Species2
83           5.8         2.7          3.9         1.2 versicolor   vervir
47           5.1         3.8          1.6         0.2     setosa   setosa
146          6.7         3.0          5.2         2.3  virginica   vervir
59           6.6         2.9          4.6         1.3 versicolor   vervir
111          6.5         3.2          5.1         2.0  virginica   vervir
94           5.0         2.3          3.3         1.0 versicolor   vervir
121          6.9         3.2          5.7         2.3  virginica   vervir
14           4.3         3.0          1.1         0.1     setosa   setosa
15           5.8         4.0          1.2         0.2     setosa   setosa
122          5.6         2.8          4.9         2.0  virginica   vervir
33           5.2         4.1          1.5         0.1     setosa   setosa
22           5.1         3.7          1.5         0.4     setosa   setosa
87           6.7         3.1          4.7         1.5 versicolor   vervir
78           6.7         3.0          5.0         1.7 versicolor   vervir
66           6.7         3.1          4.4         1.4 versicolor   vervir
133          6.4         2.8          5.6         2.2  virginica   vervir
18           5.1         3.5          1.4         0.3     setosa   setosa
79           6.0         2.9          4.5         1.5 versicolor   vervir
19           5.7         3.8          1.7         0.3     setosa   setosa
9            4.4         2.9          1.4         0.2     setosa   setosa
88           6.3         2.3          4.4         1.3 versicolor   vervir
32           5.4         3.4          1.5         0.4     setosa   setosa
116          6.4         3.2          5.3         2.3  virginica   vervir
60           5.2         2.7          3.9         1.4 versicolor   vervir
7            4.6         3.4          1.4         0.3     setosa   setosa
104          6.3         2.9          5.6         1.8  virginica   vervir
67           5.6         3.0          4.5         1.5 versicolor   vervir
20           5.1         3.8          1.5         0.3     setosa   setosa
86           6.0         3.4          4.5         1.6 versicolor   vervir
1            5.1         3.5          1.4         0.2     setosa   setosa
> res <- glm(Species~.,iris[index,], family=binomial) #linear regression
 警告メッセージ： 
 glm.fit: 数値的に 0 か 1 である確率が生じました  
> res$coefficient
 (Intercept) Sepal.Length  Sepal.Width Petal.Length  Petal.Width 
    6.587948    -7.617771   -13.925058    36.005598   -31.809827 
> logLik(res) #log-scale likelyhood
'log Lik.' -1.785923e-10 (df=5)
> predict(res,iris,type="response") # predictions from fitting functions
           1            2            3            4            5            6 
2.220446e-16 4.328295e-13 2.220446e-16 3.870520e-11 2.220446e-16 2.220446e-16 
           7            8            9           10           11           12 
2.220446e-16 2.220446e-16 7.856922e-11 9.478634e-11 2.220446e-16 4.737275e-12 
          13           14           15           16           17           18 
2.231699e-11 2.220446e-16 2.220446e-16 2.220446e-16 2.220446e-16 2.220446e-16 
          19           20           21           22           23           24 
2.220446e-16 2.220446e-16 1.795636e-12 2.220446e-16 2.220446e-16 2.220446e-16 
          25           26           27           28           29           30 
2.326153e-07 2.709491e-10 2.220446e-16 2.220446e-16 2.220446e-16 1.643915e-10 
          31           32           33           34           35           36 
3.088882e-10 2.220446e-16 2.220446e-16 2.220446e-16 3.937881e-12 2.220446e-16 
          37           38           39           40           41           42 
2.220446e-16 2.220446e-16 5.330791e-13 2.220446e-16 2.220446e-16 1.769196e-10 
          43           44           45           46           47           48 
2.220446e-16 2.220446e-16 1.556449e-13 2.220446e-16 2.220446e-16 2.626083e-13 
          49           50           51           52           53           54 
2.220446e-16 2.220446e-16 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 
          55           56           57           58           59           60 
1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 
          61           62           63           64           65           66 
1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 9.999926e-01 1.000000e+00 
          67           68           69           70           71           72 
1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 
          73           74           75           76           77           78 
1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 
          79           80           81           82           83           84 
1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 
          85           86           87           88           89           90 
1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 
          91           92           93           94           95           96 
1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 
          97           98           99          100          101          102 
1.000000e+00 1.000000e+00 9.974026e-01 1.000000e+00 1.000000e+00 1.000000e+00 
         103          104          105          106          107          108 
1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 
         109          110          111          112          113          114 
1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 
         115          116          117          118          119          120 
1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 
         121          122          123          124          125          126 
1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 
         127          128          129          130          131          132 
1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 
         133          134          135          136          137          138 
1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 
         139          140          141          142          143          144 
1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 
         145          146          147          148          149          150 
1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 
> levels(data$Species2)
[1] "setosa" "vervir"
> yhat <- levels(data$Species2)[(predict(res,iris,type="response")>0.5)+1]
> yhat
  [1] "setosa" "setosa" "setosa" "setosa" "setosa" "setosa" "setosa" "setosa"
  [9] "setosa" "setosa" "setosa" "setosa" "setosa" "setosa" "setosa" "setosa"
 [17] "setosa" "setosa" "setosa" "setosa" "setosa" "setosa" "setosa" "setosa"
 [25] "setosa" "setosa" "setosa" "setosa" "setosa" "setosa" "setosa" "setosa"
 [33] "setosa" "setosa" "setosa" "setosa" "setosa" "setosa" "setosa" "setosa"
 [41] "setosa" "setosa" "setosa" "setosa" "setosa" "setosa" "setosa" "setosa"
 [49] "setosa" "setosa" "vervir" "vervir" "vervir" "vervir" "vervir" "vervir"
 [57] "vervir" "vervir" "vervir" "vervir" "vervir" "vervir" "vervir" "vervir"
 [65] "vervir" "vervir" "vervir" "vervir" "vervir" "vervir" "vervir" "vervir"
 [73] "vervir" "vervir" "vervir" "vervir" "vervir" "vervir" "vervir" "vervir"
 [81] "vervir" "vervir" "vervir" "vervir" "vervir" "vervir" "vervir" "vervir"
 [89] "vervir" "vervir" "vervir" "vervir" "vervir" "vervir" "vervir" "vervir"
 [97] "vervir" "vervir" "vervir" "vervir" "vervir" "vervir" "vervir" "vervir"
[105] "vervir" "vervir" "vervir" "vervir" "vervir" "vervir" "vervir" "vervir"
[113] "vervir" "vervir" "vervir" "vervir" "vervir" "vervir" "vervir" "vervir"
[121] "vervir" "vervir" "vervir" "vervir" "vervir" "vervir" "vervir" "vervir"
[129] "vervir" "vervir" "vervir" "vervir" "vervir" "vervir" "vervir" "vervir"
[137] "vervir" "vervir" "vervir" "vervir" "vervir" "vervir" "vervir" "vervir"
[145] "vervir" "vervir" "vervir" "vervir" "vervir" "vervir"
> mean(yhat[index] != data$Species2[index]) #training error
[1] 0
> mean(yhat[-index] != data$Species2[-index]) #estimation error
[1] 0
> table(true=data$Species2,prediction=yhat) #show result
        prediction
true     setosa vervir
  setosa     50      0
  vervir      0    100

R tips: inserting factor value into data frame

noreply@blogger.com (Anonymous) — Mon, 20 Aug 2012 15:20:00 +0000

Here is sample R script for inserting factor value into data frame.

> data <- data.frame(iris)
> cind <- length(data) + 1
> data[data$Species == 'setosa', cind] <- c('setosa')
> data[data$Species != 'setosa', cind] <- c('vervir')
> data[,cind] <- as.factor(data[,cind])
> names(data)
[1] "Sepal.Length" "Sepal.Width"  "Petal.Length" "Petal.Width"  "Species"     
[6] "V6"          
> names(data)[6]
[1] "V6"
> names(data)[6] <- c('Species2')

geometry tips: generates identifier of standard grid cell which covers Japan region from latitude and longitude.

noreply@blogger.com (Anonymous) — Wed, 15 Aug 2012 01:44:00 +0000

In geometric operation with python, sometimes we map point (which is latitude and longitude pair) to corresponding standard grid cell. Here is sample code for mapping point to identifier of the cell.

#!/usr/bin/env python
#coding: utf-8

def get_meshcode_from_latlon(lat, lon):
    lat, lon = lat * 3600, lon * 3600

    l_code = []

    l_code += str(int(lat/2400+100))[1:3]
    l_code += str(int(lon/3600+100))[1:3]
    t_lat, t_lon = int(lat%2400), int(lon%3600)

    l_code += str(t_lat/300)
    l_code += str(t_lon/450)
    t_lat, t_lon = t_lat % 300, t_lon % 450

    l_code += str(t_lat/30)
    l_code += str(t_lon/45)
    t_lat, t_lon = t_lat % 30, t_lon % 45

    if (t_lat < 15):
        if (t_lon < 23):
            l_code += '1'
        else:
            l_code += '2'
    else:
        if (t_lat < 23):
            l_code += '3'
        else:
            l_code += '4'

    return ('').join(l_code)

if __name__ == "__main__":
    lat = 34.12345
    lon = 139.12345
    print get_meshcode_from_latlon(lat, lon)

    lat = 30.12345
    lon = 130.12345
    print get_meshcode_from_latlon(lat, lon)

python tips: usage of scipy.spatial.KDTree

noreply@blogger.com (Anonymous) — Thu, 26 Jul 2012 16:41:00 +0000

kd-tree is a well-known algorithm for searching spatially distributed points. Below code shows how to use kd-tree implemented in scipy.spatial.

#!/usr/bin/env python                                                                                                                                                                                          
#coding:utf-8                                                                                                                                                                                                  
                                                                                                                                                                                                               
import numpy as np                                                                                                                                                                                             
import scipy.spatial as ss                                                                                                                                                                                     
from itertools import combinations
from itertools import permutations

def main():
    x, y = np.mgrid[0:100, 0:100]
    points = zip(x.ravel(), y.ravel())
    tree = ss.KDTree(points)
    # get index of each point which distance from (x=0, y=0) is under 1
    a = tree.query_ball_point([0, 0], 1)
    print [points[i] for i in a]

if __name__ == "__main__":
    main()

output of above code is like:

~$ python above_code.py 
[(0, 0), (0, 1), (1, 0)]

python tips: usage of nkf.python2

noreply@blogger.com (Anonymous) — Tue, 24 Jul 2012 14:57:00 +0000

NKF.python2 enables us to convert Japanese text encoding without specifying input text coding. This article explains how to install and use NKF.python2. First, install NKF.python2 on your linux system:

# download nkf using git
$ git clone git://git.sourceforge.jp/gitroot/nkf/nkf.git
Cloning into 'nkf'...
remote: Counting objects: 1378, done.
remote: Compressing objects: 100% (432/432), done.
remote: Total 1378 (delta 945), reused 1378 (delta 945)
Receiving objects: 100% (1378/1378), 503.08 KiB, done.
Resolving deltas: 100% (945/945), done.

$ cd nkf
$ sudo make install #if you want
$ cd NKF.python2
$ python setup.py build
$ python setup.py install

Second, check NKF.python2 is installed correctly:

$ python
Python 2.7.3 (default, Apr 20 2012, 22:39:59) 
[GCC 4.6.3] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> import nkf
>>>

Finally, test NKF.python2 (path_to_file points japanese enceded text file):

with open(path_to_file, 'r') as F:
    data = [nkf("-w -Lu -d", row) for row in F]

Above code converts encoding of each row to utf-8 *without* specifying types of encoding of input file. This feature is very useful if you do not know which encode is used for input file or you need to read several input files which have different encoding each simultaneously.

python tips: calculate distance between N points and M points using numpy array

noreply@blogger.com (Anonymous) — Mon, 23 Jul 2012 07:38:00 +0000

This is memory consuming way of calculate distance between N points and M points. Let l_lat and l_lon are vectors for N points and v_lat and v_lon are for M points. Using numpy array, distance calculation can be written as follows.

#!/usr/bin/env python                                                                                                                                                                                        
#coding:utf-8                                                                                                                                                                                                

import numpy

M = 3
N = 5

# points1                                                                                                                                                                                                    
l_lat = numpy.arange(N)
l_lon = numpy.arange(N)
m_l_lat = numpy.tile(l_lat, (M, 1))
m_l_lon = numpy.tile(l_lon, (M, 1))

# points2                                                                                                                                                                                                    
v_lat = numpy.arange(M)
v_lon = numpy.arange(M)
m_v_lat_t = numpy.tile(v_lat, (N, 1)).T
m_v_lon_t = numpy.tile(v_lon, (N, 1)).T

distance = numpy.sqrt((m_v_lat_t - m_l_lat) ** 2 + (m_v_lon_t - m_l_lat) ** 2)
print distance

Although code generates lat and lon with numpy.arange(), in real case you should fill each vector with appropriate lat and lon values. Below is an output of previous code.

~$ python test_calc_distance.py 
[[ 0.          1.41421356  2.82842712  4.24264069  5.65685425]
 [ 1.41421356  0.          1.41421356  2.82842712  4.24264069]
 [ 2.82842712  1.41421356  0.          1.41421356  2.82842712]]

python tips: effective sparse matrix operation ?

noreply@blogger.com (Anonymous) — Thu, 31 May 2012 12:56:00 +0000

#!/usr/bin/env python
#coding:utf-8
 
import numpy as np
import scipy as sp
from scipy.sparse import *
from numpy.random import *
from collections import OrderedDict
from datetime import datetime
 
# performance logger
import cProfile as profile
import pstats
 
from itertools import groupby
from operator import itemgetter

def gen_mat(M, N, R, format="csc"):
    """generating lil_matrix with random valuables"""
    mat = sp.sparse.rand(M, N, R, format=format, dtype=np.dtype('float64'))
    mat = mat * 250
    print("{0}x{1} matrix (sparsity:{2}, # of data:{3}) is generated...".format(M, N, R, mat.nnz))
    return mat
 
def join_columns1(spamat, target_col):
    """join columns using getcol() which requires object copying"""
    return reduce(lambda x,y: x+y, [spamat[:,index] for index in target_col])

def join_columns2(lilmat, target_col):
    """join columns using lil_matrix's getrowview() which does not require object copying """
    lilt = lilmat.T
    return reduce(lambda x,y: x+y, (lilt.getrowview(index) for index in target_col))

def join_columns3(coomat, N):
    """join columns using coo indices, numpy's fancy index and pydict"""
    dic = {}

    def dic_update(key, val):
        if (dic.has_key(key)):
            dic[key] += val
        else:
            dic[key] = val

    mask = (coomat.col % 26 == 0)
    row = coomat.row[mask]
    data = coomat.data[mask]
    n = len(row)

    [dic_update(row[index], data[index]) for index in range(n)]
    return dic

def get_elapsed_time(func, listargs):
    """wrapper for calculating consuming time of function"""
    s = datetime.now()
    ret = func(*listargs)
    e = datetime.now()
    print("{0}\t{1}".format((e-s), func.__name__))
    return ret
     
def bench(M, N, R, L):
    column_indexes = frozenset([i for i in range(0, N, 26)])

    print("# lil_matrix")
    lilmat = gen_mat(M, N, R, format="lil")
    # get_elapsed_time(join_columns1, (lilmat, column_indexes)) # very slow!
    get_elapsed_time(join_columns2, (lilmat, column_indexes)) # bit slow

    print("# csc_matrix")
    cscmat = lilmat.tocsc()
    get_elapsed_time(join_columns1, (cscmat, column_indexes)) # bit fast

    print("# csr_matrix")
    csrmat = cscmat.tocsr()
    get_elapsed_time(join_columns1, (cscmat, column_indexes)) # bit fast

    print("# coo_matrix")
    coomat = cscmat.tocoo()
    get_elapsed_time(join_columns1, (cscmat, column_indexes)) # bit fast
    get_elapsed_time(join_columns3, (coomat, N))              # very fast

def main():
    M = 1400000      # n_row
    N = 50000        # n_col
    R = 0.0005       # sparsity
    L = int(N * 0.2) # column_indexes
 
    print("-------------------------------------")
    print("M:{0}, N:{1}, R:{2}, L:{3}".format(M, N, R, L))
    print("-------------------------------------")
    bench(M, N, R, L)
     
if __name__ == "__main__":
   main()

-------------------------------------
M:1400, N:5000, R:0.001, L:1000
-------------------------------------
1400x5000 matrix (sparsity:0.001, # of data:7000) is generated...
# lil_matrix
0:00:01.576205 join_columns1
0:00:00.044499 join_columns2
# csc_matrix
0:00:00.047233 join_columns1
# csr_matrix
0:00:00.047329 join_columns1
# coo_matrix
0:00:00.047277 join_columns1
0:00:00.000282 join_columns3

-------------------------------------
M:140000, N:50000, R:0.001, L:10000
-------------------------------------
140000x50000 matrix (sparsity:0.001, # of data:7000000) is generated...
# lil_matrix
0:00:03.003079 join_columns2
# csc_matrix
0:00:02.640237 join_columns1
# csr_matrix
0:00:02.640100 join_columns1
# coo_matrix
0:00:02.647828 join_columns1
0:00:00.300576 join_columns3

statistics tips: my covariance matrix implementation on python

noreply@blogger.com (Anonymous) — Wed, 30 May 2012 16:44:00 +0000

#!/usr/bin/env python                                                                                                                
#coding: utf-8                                                                                                                       
                                                                                                                                     
import numpy as np                                                                                                                   
from numpy.random import *                                                                                                           
                                                                                                                                     
def get_covariance_matrix(A):                                                                                                        
    N, M = A.shape                                                                                                                   
    ret = np.reshape(np.zeros(N*M), (N,M))                                                                                           
                                                                                                                                     
    for n in range(0, N):                                                                                                            
        for m in range(0, M):                                                                                                        
            ret[n][m] = calc_covariance(A, n, m)                                                                                     
                                                                                                                                     
    return ret                                                                                                                       
                                                                                                                                     
def calc_covariance(A, n, m):                                                                                                        
    nominator = np.sum((A[n] - np.average(A[n])) * (A[m] - np.average(A[m])))                                                        
    denominator = A.shape[0] - 1                                                                                                     
    return nominator/denominator                                                                                                     
                                                                                                                                     
def test():                                                                                                                          
    D = 4                                                                                                                            
    A = np.reshape([randint(100) for i in range(D**2)], (D, D))                                                                      
                                                                                                                                     
    print("# numpy.cov() output:")                                                                                                   
    print(np.cov(A))                                                                                                                 
                                                                                                                                     
    print("# my covariance matrix:")                                                                                                 
    print(get_covariance_matrix(A))                                                                                                  
                                                                                                                                     
    print("# my precision matrix:")                                                                                                  
    print(get_covariance_matrix(A) ** -1)  
                                                                                                                                     
if __name__ == "__main__":                                                                                                           
    test()

yaboo@overtones:~$ python kyoubunsan.py
# numpy.cov() output:
[[  226.91666667   -26.83333333   -39.33333333   499.75      ]
 [  -26.83333333   328.33333333   629.          -257.83333333]
 [  -39.33333333   629.          1432.66666667  -882.        ]
 [  499.75        -257.83333333  -882.          1976.25      ]]                                                                      
# my covariance matrix:
[[  226.91666667   -26.83333333   -39.33333333   499.75      ]
 [  -26.83333333   328.33333333   629.          -257.83333333]
 [  -39.33333333   629.          1432.66666667  -882.        ]
 [  499.75        -257.83333333  -882.          1976.25      ]]
# my precision matrix:                                                                                                               
[[ 0.0044069  -0.03726708 -0.02542373  0.002001  ]
 [-0.03726708  0.00304569  0.00158983 -0.00387847]
 [-0.02542373  0.00158983  0.000698   -0.00113379]
 [ 0.002001   -0.00387847 -0.00113379  0.00050601]]

algorithm tips: Sparse Matrix Multiplication Package (SMMP)

noreply@blogger.com (Anonymous) — Tue, 22 May 2012 11:33:00 +0000

UNDERWRITING... Matrix operation is very useful for improving computing speed of our codes. For this operation, python has very popular packages such as numpy and scipy. In this article, I show sample code to calculate dot product using sparse matrix and vector with numpy and scipy.

import scipy
import numpy
from scipy.sparse import *
from numpy.random import *

vec = create_vec(N)
mat = create_mat(N,M)
ans = mat.dot(vec)

In above code, mat.dot(vec) calculates dot product. This method implements sparse matrix multiplication pacakage called SMMP which computational complexity is : O(n_row*K^2 + max(n_row,n_col)) where n_row is # of matrix rows, n_col is # of matrix columns, K is the maximum index for nonzero value in a row of A and column of B.

python tips: unpickling object is very slow because of garbage collection

noreply@blogger.com (Anonymous) — Tue, 22 May 2012 10:28:00 +0000

In python, unpickling gigantic object from file is very much time consuming because of python's garbage collection. Python's garbage collection starts when cPickle loads gigantic object contains a lot of container object and number of container objects in the object exceed threshold of collect() method, python's garbage collection starts.
Now, assume you have 2GB pickle object which contains huge dictionary object in your file system and your python code load it. I compare loading time with/without garbage collection in python 2.7.2 on Ubuntu12.04.

#!/usr/local/bin/python                                                         
#coding: utf-8                                                                  
                                                                                
import cPickle as pickle                                                        
from datetime import datetime                                                   
import gc
                                                                                
F = file('gigantic_data.pickle', 'rb')                                          
s = datetime.now()                                                              
myobj = pickle.load(F)                                                          
e = datetime.now()                                                              
F.close()
                                                                                
print("load_pickle takes is {0}".(e-s))

# disabling cyclic garbage collection
gc.disable()

F = file('gigantic_data.pickle', 'rb')                                          
s = datetime.now()                                                              
myobj = pickle.load(F)                                                          
e = datetime.now()                                                              
F.close()

print("load_pickle takes is {0}".(e-s))

> /usr/local/bin/python above_script.py
load_pickle takes 0:15:16.530199
load_pickle takes 0:00:14.195409

python tips: multiprocessing and shared memory

noreply@blogger.com (Anonymous) — Fri, 18 May 2012 11:23:00 +0000

With multiprocessing package, it is extremely easy to make programs which share tasks with multiple workers(processes). This post is memo for how to share object which initialized by parent process with child process.

General way for sharing object between parent process and child processes is using Queue or Pipe. However, it's the easiest way for me to use global valuable. Below code is just a example of how to share object between processes using global valuable.

#!/usr/bin/env python
# coding:utf-8

import os
import multiprocessing

x = 10
y = 20
z = []
a = None
b = [1,2,3]
c = None

c = 3

def main():
    a = 1
    print("###############################")
    print("x:: {0}".format(print_value(x)))
    print("y:: {0}".format(print_value(y)))
    print("z:: {0}".format(print_value(z)))
    print("a:: {0}".format(print_value(a)))
    print("b:: {0}".format(print_value(b)))
    print("c:: {0}".format(print_value(c)))
    print("###############################")
    pool = multiprocessing.Pool(3)
    pool.map(func, [1,2,3], chunksize=1)
    

def print_value(val):
    return "{0}, id:{1}".format(val, id(val))
    
def func(val):
    print("pid: {0}, process_number: {1}".format(os.getpid(), val))
    print("x:: {0}".format(print_value(x)))
    print("y:: {0}".format(print_value(y)))
    print("z:: {0}".format(print_value(z)))    
    print("a:: {0}".format(print_value(a)))
    print("b:: {0}".format(print_value(b)))
    print("c:: {0}".format(print_value(c)))    

if __name__ == "__main__":
    main()

Here is output from above code. 'id' means memory-address pointed by valuable.

yaboo@ubuntu:~$ python shared_memory.py 
###############################
x:: 10, id:24606432
y:: 20, id:24606192
z:: [], id:139931380327848
a:: 1, id:24606648
b:: [1, 2, 3], id:139931380327920
c:: 4, id:24606576
###############################
pid: 3747, process_number: 1
x:: 10, id:24606432
pid: 3748, process_number: 2
y:: 20, id:24606192
x:: 10, id:24606432
z:: [], id:139931380327848
pid: 3749, process_number: 3
y:: 20, id:24606192
a:: None, id:9071744
x:: 10, id:24606432
z:: [], id:139931380327848
y:: 20, id:24606192
b:: [1, 2, 3], id:139931380327920
a:: None, id:9071744
c:: 3, id:24606600
z:: [], id:139931380327848
a:: None, id:9071744
b:: [1, 2, 3], id:139931380327920
b:: [1, 2, 3], id:139931380327920
c:: 3, id:24606600
c:: 3, id:24606600
yaboo@ubuntu:~$

tips: port-forwarding with ssh

noreply@blogger.com (Anonymous) — Tue, 01 May 2012 01:54:00 +0000

These are memos for usage of port forwarding with ssh command. First of all, you need to enable ssh port forwarding with below command.

ssh -f -N -L [#port]:[dest_hostname]:22 [user_name_for_proxy_host]@[proxy_hostname]

If you want to transfer some files use below commands.

scp -P [#port] [user_name_for_dest_host]@localhost:[filename] ./[filename]
scp -P [#port]  ./[filename] [user_name_for_dest_host]@localhost:[filename]

If you want to login with ssh, use below commands.

ssh -p [#port] [user_name_for_dest_host]@localhost

python tips: array computation with numpy II

noreply@blogger.com (Anonymous) — Wed, 25 Jan 2012 14:31:00 +0000

As previous post shows, for-loop calculation in array computations needs more processing time than vectorized calculation. There are some techniques for fasten calculation over arrays. Sample code shown below is one of them.

>>> from numpy import *
>>> 
>>> A = array([[1,2,3], [4,5,6], [7,8,9]])
>>> 
>>> def somefunc(x):
...     x1 = x
...     x2 = x - 3
...     return where(x < 7, x1, x2)
... 
>>> somefunc(A)
array([[1, 2, 3],
       [4, 5, 6],
       [4, 5, 6]])
>>>

python tips: array computation with numpy

noreply@blogger.com (Anonymous) — Wed, 25 Jan 2012 13:55:00 +0000

I compared performance of array computations. Listed below are results and codes of some experiments which introduced in the book named 'Python Scripting for Computational Science (3rd edition)'.

>>> from numpy import *
>>> import time
>>> 
>>> a = linspace(0, 1, 1E+07)
>>> t0 = time.clock()
>>> b = 3*a - 1
>>> 
>>> t1 = time.clock()
>>> 
>>> for i in xrange(a.size):
...     b[i] = 3*a[i] -1
... 
>>> t2 = time.clock()
>>> 
>>> print '3*a-1: %g sec, loop: %g sec' %(t1-t0, t2-t1)
3*a-1: 0.05 sec, loop: 15.46 sec

python tips: create shapefile with python

noreply@blogger.com (Anonymous) — Tue, 24 Jan 2012 17:20:00 +0000

Long time has been passed since last post.

Here, I want to show the way to create Shapefile, which is standard data format for geometry data, using python.

Among several python tools such as geoscript, pyshp and pyproj which support operations of Shapefiles, I introduce how to use pyshp.

Pyshp is open source, provides shapefile library written all in pure python and distributed under MIT lisence.

You can download it as shapfile.py from the official web site or just install via easy_install command:

# easy_install pyshp

However most of available operations of pyshp are introduced in the web site, I write some code using pyshp here just for my memo.

If you want to know more, I recommend looking over the web site and make some sample code according to tutorials.

#!/usr/bin/env python

import shapefile

def create_square(x=0., y=0., i=1.):
    return [[[x, y], [x+i, y], [x+i, y+i], [x, y+i], [x, y]]]

def define_geometry(w):
    # 1st record's geometry
    w.poly(parts=create_square(0, 0, 1))
    # 2nd record's geometry
    w.poly(parts=create_square(1, 1, 1))
    # 3rd record's geometry
    w.poly(parts=create_square(0.5, 0.5, 1))
    return w

def define_fields(w):
    w.field('FLD_1st')
    w.field('FLD_2nd', 'N', '10')
    return w

def fill_attributes(w):
    w.record('hoge', 'hogehogeho')
    w.record('mage', 'mage1')
    w.record('moge', 'mogemoge')
    return w

def create_proj_file(filename):
    # this function creates .proj file which defines projection.
    prj = open("%s.prj" % filename, "w")
    epsg = 'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]]'
    prj.write(epsg)
    prj.close()

def main():
    filename = 'test'
    w = shapefile.Writer(shapefile.POLYGON)
    w = define_geometry(w)
    w = define_fields(w)
    w = fill_attributes(w)
    w.save('./' + filename)
    create_proj_file(filename)

if __name__ == "__main__":
    main()

sphinx tips: add subsubsection into pdf compiled by make all-pdf-ja

noreply@blogger.com (Anonymous) — Sun, 11 Sep 2011 04:51:00 +0000

When sphinx compiles latex files, it chooses cls file which written in conf.py.

# Grouping the document tree into LaTeX files. List of tuples
# (source start file, target name, title, author, documentclass [howto/manual]).
latex_documents =
   ('index', 'sample_document.tex', u'This is sample',
     u'Author Name', 'manual'),
]

If 'manual' is specified in conf.py, sphinx read sphinxmanual.cls when it compiles latex file. So, if you want to enable subsubsection, you should edit sphinxmanual.cls or other specified cls file. In my environment, original sphinxmanual.cls is located in /usr/local/lib/python2.6/dist-packages/Sphinx-1.0.7-py2.6.egg/sphinx/texinputs/sphinxmanual.cls. And you can enable subsubsection with editing this cls file like below.

before:
\setcounter{secnumdepth}{2}

after:
\setcounter{secnumdepth}{3}

enjoy.

sphinx tips: generate center-aligned table in pdf document

noreply@blogger.com (Anonymous) — Fri, 02 Sep 2011 08:21:00 +0000

When we type 'make latex' command in sphinx root directory, sphinx generates set of latex files from rst files based on '/usr/share/pyshared/sphinx/writers/latex.py' which defines rules for generating latex files. If you do not like default behavior of latex generator, you can change it by modifying the code.

Sometimes I was annoyed by left-aligned table which is automatically created by using csv-table directive when creating a pdf document. So, I changed some part of latex.py. Below is difference between original latex.py and bit-of modified latex.py which can generate center-aligned table by using csv-table.

( ... just added '\begin{center}' and '\end{center}' command )

617c617
<             self.body.append('\n\\begin{threeparttable}\n'
---
>             self.body.append('\n\\begin{center}\n\\begin{threeparttable}\n'
664c664
<             self.body.append('\\end{threeparttable}\n\n')
---
>             self.body.append('\\end{threeparttable}\n\\end{center}\n\n')

postgres tips: enable plpgsql for implementing stored procedure

noreply@blogger.com (Anonymous) — Sun, 21 Aug 2011 03:01:00 +0000

This entry is just a memo for me 8^)
If we want to implement stored procedure to any database, it is necessary to type below command.

createlang -d dbname plpgsql

python tips: draw beta distribution with matplotlib

noreply@blogger.com (Anonymous) — Tue, 16 Aug 2011 10:25:00 +0000

import numpy as np
import numpy as np
import pylab as pl
import scipy.special as ss

def beta(a, b, mew):
    e1 = ss.gamma(a + b)
    e2 = ss.gamma(a)
    e3 = ss.gamma(b)
    e4 = mew ** (a - 1)
    e5 = (1 - mew) ** (b - 1)
    return (e1/(e2*e3)) * e4 * e5

def plot_beta(a, b):
    Ly = []
    Lx = []
    mews = np.mgrid[0:1:100j]
    for mew in mews:
        Lx.append(mew)
        Ly.append(beta(a, b, mew))
    pl.plot(Lx, Ly, label="a=%f, b=%f" %(a,b))

def main():
    plot_beta(0.1, 0.1)
    plot_beta(1, 1)
    plot_beta(2, 3)
    plot_beta(8, 4)
    pl.xlim(0.0, 1.0)
    pl.ylim(0.0, 3.0)
    pl.legend()
    pl.show()

if __name__ == "__main__":
    main()