Large data sets in R

近來正在學習如何處理 large data sets，非常頭痕。

何謂 large data set，不同人有不同定義。在某些應用範疇，要起碼以 TB 來計算的才稱得上是真正的大筆資料；在另外一些範疇，幾個 GB 已算，總之 RAM 不夠裝載的，就是大筆。以我的情況來說，資料只有數百 MB，完全可以讀入記憶體，就其他人的標準而言，只是很細筆資料。然而我用的只是單機 PC，幾百 MB 到兩三 GB 的資料量，不但讀入資料須時，還佔去大量記憶體。完全讀入記憶體，並不划算。

要提高效率，首先要將資料轉為二進位檔，令讀取速度較讀取文字檔為快。然而此舉措不能省回多少讀入整筆資料的時間。由於我們每次使用資料，往往只須整筆資料的其中一小部份，因此更重要的，是要找方法 on-demand 地直接從硬碟讀取所需部份。

長遠來說，也許將資料匯入某些 SQL database 最徹底，不過暫時我只希望找到一些 quick and dirty 的方法，可以使 R/Python/Octave/C++ 都可以快速使用同一筆資料，其中，R 最令我頭痕。

舉例說，我有一個 temp.csv 檔，大約 500MB，當中每橫行均有四欄資料

n, t1, t2, s

其中 n 是整數、t1 與 t2 是一個 2008 年或以後的 UTC time strings（即 2012-01-16 00:03:59 之類），而 s 是一列長短不一，但介乎 0 與 280 的整數（例如 "1 5 236 247 248" 或 "16 135 277" 等等）。若單純將此 csv file 存為 R 的 rda format，一來每次仍須讀入整筆資料，二來，出乎意料地，讀入整筆資料的時間依然很長（在我的電腦上要大約兩分鐘），原因除了資料量較大以外，好像還因為 R 要花時間解讀 s 這個字串。因此，我首先要做的，是將 s 轉化為一個 fixed width 的欄位。280 bits 大約有 9 bytes，所以我們可用 9 個 32-bit integers 作儲存。每當 s 包含某數字 m ── 例如 m=236，我們可以將第 $\lfloor m/32\rfloor$ 個 byte 的第 $m \textrm{mod} 32$ 個 bit 設為 1。這可用 Python 輕鬆搞掂：

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:
        headers = line.split(",")
    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 讀回數據。

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())

這種寫法看起來有點蠢，不過它富彈性的地方，在於每個欄位的資料種類也可不同，此處我每一欄都是一個 32-bit integer，但若第四戙是一個雙精密度的浮點數字，你可以將第四個 int32() 改成 double()，餘此類推。定義格式之後，就可以將 temp.bin 映射到一個 mmap object:

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

這樣，就可以讀寫 temp.bin 這筆數據。mmap 的原理，是有需要時，才將硬碟中所需的數據映射到記憶體之中，因此毋須讀入整筆資料。你可以用 m[i,j] 的方式提取數據：

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

此處即可見到 mmap 與 data frame 的分別。若 m 是一個 data frame，當鍵入 m[1,2] 之後，並不會有以上 "[[1]]" 那一行。事實上，若要顯示 m 的第一行的頭三個數字，mmap 與 data frame 會有資料結構上的分別。

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

[[2]]
[1] 83017033

[[3]]
[1] 83017098

若換了是 data frame，應會有如下結果：

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

出現此分別的原因，在於 mmap 將一個 table 視為 a list of columns，而每個 column 又是一個 list，因此，要讀取一個 column，所用的語法就與 data.frame 有點不同。例如要讀取第二個 column，你就應該用以下其中一種方法：

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

這些過多的方括號很煩人，但起碼能夠 get things done。

最後，當用完一筆資料時，記得釋放數據，否則不知會否鎖死資料，令其他程式不能讀取。

> munmap(m)

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