matlab的numpy當量dummyvar

Question

爲了很好地處理類別變量，matlab的dummyvar函數最常用的pythonic等效函數是什麼？

下面是一個說明我的問題的例子，用N×M矩陣表示將N個數據點劃分成< = N個類別的M種不同方式。

​​

的任務是有效地計數的任何兩個數據點被分類爲同一類別的次數並且將結果存儲在一個NxN矩陣。在matlab中，這可以通過dummyvar實現，它爲每個分區的每個類別創建一個列變量。

>> dummyvar(partitions)*dummyvar(partitions)' 
ans = 
3  2  1  1  1  1  1  0  1  2 
2  3  2  0  2  0  2  1  2  1 
1  2  3  1  1  1  3  2  1  0 
1  0  1  3  1  3  1  1  0  2 
1  2  1  1  3  1  1  1  2  2 
1  0  1  3  1  3  1  1  0  2 
1  2  3  1  1  1  3  2  1  0 
0  1  2  1  1  1  2  3  2  0 
1  2  1  0  2  0  1  2  3  1 
2  1  0  2  2  2  0  0  1  3 

，我能想到的解決這個任務是寫一個O（N * M）環模仿dummyvar行爲的最有效方式。 （請注意，下面的代碼更喜歡partition.shape[0] < < partition.shape[1]，這通常可能是真實的，但假設是不安全的）。

dv=np.zeros((0,10)) 
for row in partitions: 
    for val in xrange(1,np.max(row)+1): 
    dv=np.vstack((dv,row==val)) 
np.dot(dv.T,dv) 

，當然還有因爲vstack在一個循環是非常低效的，這可以通過查找所需的大小和創建陣列與開始時得到改善，但我真的找了一個襯墊做一樣在matlab中。

編輯：關於我在做什麼只是添加上下文的更多信息。我正在編寫用於分析大腦網絡的庫的python（沒有python實現存在）的庫函數。現有的工作matlab源是可用的。由於特定於域的約束，輸入的最大大小約爲幾千個節點的網絡。但是，基本上我寫的所有功能都必須很好地適應大型輸入。

Answer 1

你可以做一個小廣播魔術快讓你的虛擬陣列：

>>> partitions = np.array([[1, 1, 2, 2, 1, 2, 2, 2, 1, 1], 
...      [1, 2, 2, 1, 2, 1, 2, 2, 2, 1], 
...      [1, 1, 1, 2, 2, 2, 1, 3, 3, 2]]) 
>>> n = np.max(partitions) 
>>> d = (partitions.T[:, None, :] == np.arange(1, n+1)[:, None]).astype(np.int) 
>>> d = d.reshape(partitions.shape[1], -1) 
>>> d.dot(d.T) 
array([[3, 2, 1, 1, 1, 1, 1, 0, 1, 2], 
     [2, 3, 2, 0, 2, 0, 2, 1, 2, 1], 
     [1, 2, 3, 1, 1, 1, 3, 2, 1, 0], 
     [1, 0, 1, 3, 1, 3, 1, 1, 0, 2], 
     [1, 2, 1, 1, 3, 1, 1, 1, 2, 2], 
     [1, 0, 1, 3, 1, 3, 1, 1, 0, 2], 
     [1, 2, 3, 1, 1, 1, 3, 2, 1, 0], 
     [0, 1, 2, 1, 1, 1, 2, 3, 2, 0], 
     [1, 2, 1, 0, 2, 0, 1, 2, 3, 1], 
     [2, 1, 0, 2, 2, 2, 0, 0, 1, 3]]) 

有明顯的缺點，即使行只有幾個不同的價值觀，我們正在創造的意志仿真陣列對於具有最多值的行，具有該行所需的那麼多列。但除非你有巨大的陣列，否則它可能會比其他方法更快。

---

好吧，如果你是一個可擴展的解決方案之後，你想用一個稀疏數組爲您的虛擬矩陣。

import scipy.sparse as sps 
def sparse_dummyvar(partitions): 
    num_rows = np.sum(np.max(partitions, axis=1)) 
    nnz = np.prod(partitions.shape) 
    as_part = np.argsort(partitions, axis=1) 
    # You could get s_part from the indices in as_part, left as 
    # an exercise for the reader... 
    s_part = np.sort(partitions, axis=1) 
    mask = np.hstack(([[True]]*len(items_per_row), 
         s_part[:, :-1] != s_part[:, 1:])) 
    indptr = np.where(mask.ravel())[0] 
    indptr = np.append(indptr, nnz) 

    return sps.csr_matrix((np.repeat([1], nnz), as_part.ravel(), indptr), 
          shape=(num_rows, partitions.shape[1],)) 

這返回的dummyvar(partitions)轉置：如果你不熟悉的CSR稀疏格式的細節下面的代碼可能很難效仿。你可以通過簡單調用csc_matrix而不是csr_matrix並交換形狀值來獲得陣列。但是，由於你只是在矩陣的乘積之後，並且scipy在乘以之前將所有內容都轉換爲CSR格式，所以它可能稍微快一點。你現在可以這樣做：

>>> dT = sparse_dummyvar(partitions) 
>>> dT.T.dot(dT) 
<10x10 sparse matrix of type '<type 'numpy.int32'>' 
    with 84 stored elements in Compressed Sparse Column format> 
>>> dT.T.dot(dT).A 
array([[3, 2, 1, 1, 1, 1, 1, 0, 1, 2], 
     [2, 3, 2, 0, 2, 0, 2, 1, 2, 1], 
     [1, 2, 3, 1, 1, 1, 3, 2, 1, 0], 
     [1, 0, 1, 3, 1, 3, 1, 1, 0, 2], 
     [1, 2, 1, 1, 3, 1, 1, 1, 2, 2], 
     [1, 0, 1, 3, 1, 3, 1, 1, 0, 2], 
     [1, 2, 3, 1, 1, 1, 3, 2, 1, 0], 
     [0, 1, 2, 1, 1, 1, 2, 3, 2, 0], 
     [1, 2, 1, 0, 2, 0, 1, 2, 3, 1], 
     [2, 1, 0, 2, 2, 2, 0, 0, 1, 3]])