2012年1月16日星期一

Large data sets in R

n, t1, t2, s

import struct
import time
import calendar

t2008 = calendar.timegm(time.strptime('1/1/2008', '%d/%m/%Y'))

def seconds_since_2008(time_str, fmt_str):
if time_str == 'NULL':
return -1
else:
return calendar.timegm(time.strptime(time_str, fmt_str)) - t2008

# The last column (s) of the csv file is a list of ordered distinct
# integers between 0 and 280.  We represent it a 32-bit unsigned
# integer array of length 9.  If the list contains an integer t, we
# will flag the (t%32)-th LEAST significant bit (i.e. the (t%32)-th
# bit from the RIGHT) of the (t//32)-th byte as 1.
s_array_size = 9

s_array = [0]*s_array_size  # an array that will contain the numbers in s

infile = open('c:/temp/temp.csv', 'r')
outfile = open('c:/temp/temp.bin', 'wb')

iteration=0
for line in infile:
if iteration==0:
else:
data = line.split(",")
s_array = [0]*s_array_size

n = int(data[0])

# The following three entries are represented as number
# of seconds since 2008-01-01 00:00:00. NULL values are
# represented by -1.
t1 = seconds_since_2008(data[1], '%Y-%m-%d %H:%M:%S')
t2 = seconds_since_2008(data[2], '%Y-%m-%d %H:%M:%S')

# s
s = data[3].split()
for item in s:
m = int(item)
s_array[(m//32)] |= (1<<(m%32))

outfile.write( struct.pack('<iii', n, t1, t2) )
for m in s_array:
outfile.write( struct.pack('<I', m) )
iteration +=1
infile.close()
outfile.close()


R 有三個處理大筆資料的 packages：ff, bigmemory, mmap（其實還有一個 lazy.frame，不過看來並不很成熟）。ff 功能最廣，但也最複雜；bigmemory 較側重如何用 shared memory 讀寫大量數據，亦支援 memory-mapped files ，資料格式亦似乎最有彈性。ff 和 bigmemory 均可與 biglm package 配合，用來做 linear 或 generalized linear regression。兩個套件皆曾經獲獎（亦見 bigmemory 作者的 presentation）。

mmap 的作者乃現職 quant，作者自稱動輒要應付以 TB 計的數據。以匯入資料的介面而言，三個 packages 之中，以 mmap 最直觀，但不知 mmap 能否配合 biglm 或其他套件使用。由於暫時我仍在非常前期的清理數據與基本分析階段，所以暫且先嘗試 mmap。前面我用 Python 將一個叫 temp.csv 的檔案轉為一個稱作 temp.bin 的二進位檔。以下將用 mmap 來讀它。首先定義 temp.bin 每行的格式。先前我用一個 32-bit integer 來表達 n, t1 與 t2，而 s 就用 9 個 32-bit integers 來表達，故此每行共有 12 個 32-bit integers：
> library(mmap)
> record.type <- struct(int32(), int32(), int32(), int32(), int32(), int32(), int32(), int32(), int32(), int32(), int32(), int32())


> m <- mmap(file="c:/temp/temp.bin", mode=record.type)


> m[1,2]
[[1]]
[1] 83017033


> m[1,1:3]
[[1]]
[1] 85818

[[2]]
[1] 83017033

[[3]]
[1] 83017098


> m[1,1:3]
V1       V2       V3
1 85818 83017033 83017098


> m[][[2]]
> colnames(m) <- ('a', 'b', 'c')
> m[]$b > m[][['b']]  若純粹要顯示第二個 column，你仍可以像顯示 data frame 般鍵入 m[,2]，但是若要做頭兩戙的 linear regression，你就不能說 > library(lm) > lm(m[,1] ~ 1 + m[,2]) Error in model.frame.default(formula = m[, 1] ~ 1 + m[, 2], drop.unused.levels = TRUE) : invalid type (list) for variable 'm[, 1]'  而應該說 > library(lm) > lm(m[][[1]] ~ 1 + m[][[2]])  用家必須小心留神。又例如要找出 m 第二戙有幾多個 0，我們不能用以下方式查詢 > length(m[,2][m[,2]==0]) > length(m$b[m$b==0]) > length(m['b'][m['b']==0])  而必須使用下列其中一種方法： > length(m[][[2]][m[][[2]][]==0]) > length(m[]$b[m[]\$b[]==0])
> length(m[][['b']][m[][['b']]==0])


> munmap(m)

mmap 套件還可配合另一個稱為 indexing 的套件（只在 R-forge，CRAN 好像還未有，而且暫時只能在 *nix 上使用），令功能更強，兼可配合資料庫使用。見套件作者夫子自道

…竟然是為了替普京助選