NumPy C 程式碼說明#

狂熱主義在於當你忘記目標時，加倍努力。— 喬治·桑塔亞納

權威人士是指比你真正想知道的更多地告訴你某件事的人。— 佚名

此頁面嘗試解釋一些新程式碼背後的邏輯。這些說明的目的是使人們能夠比僅僅盯著程式碼更容易理解實作背後的想法。或許透過這種方式，演算法可以被更多人改進、借鑒和/或優化。

記憶體模型#

ndarray 的一個基本方面是，陣列被視為從某個位置開始的「記憶體區塊」。此記憶體的解釋取決於 stride 資訊。對於 \(N\) 維陣列中的每個維度，一個整數 (stride) 指示必須跳過多少位元組才能到達該維度中的下一個元素。除非您有一個單段陣列，否則在遍歷陣列時必須查閱此 stride 資訊。編寫接受 stride 的程式碼並不困難，您只需使用 char* 指標，因為 stride 的單位是位元組。另請記住，stride 不必是元素大小的單位倍數。此外，請記住，如果陣列的維度數為 0（有時稱為 rank-0 陣列），則 strides 和 dimensions 變數為 NULL。

除了 PyArrayObject 的 strides 和 dimensions 成員中包含的結構資訊外，flags 還包含有關如何存取資料的重要資訊。特別是，當記憶體位於根據資料類型陣列的合適邊界上時，會設定 NPY_ARRAY_ALIGNED 旗標。即使您有一個 contiguous 的記憶體區塊，您也不能僅僅假設對元素取消引用資料類型特定的指標是安全的。只有當 NPY_ARRAY_ALIGNED 旗標設定時，這才是安全的操作。在某些平台上它可以運作，但在其他平台（如 Solaris）上，它會導致匯流排錯誤。NPY_ARRAY_WRITEABLE 也應確保您計劃寫入陣列的記憶體區域。也可以取得指向不可寫記憶體區域的指標。有時，當未設定 NPY_ARRAY_WRITEABLE 旗標時，寫入記憶體區域只會是不禮貌的行為。其他時候，它可能會導致程式崩潰（例如，作為唯讀記憶體映射檔案的資料區域）。

資料類型封裝#

另請參閱

資料類型物件 (dtype)

資料類型 是 ndarray 的重要抽象概念。運算將查看資料類型，以提供運算陣列所需的關鍵功能。此功能在 PyArray_Descr 結構的 f 成員指向的函數指標列表中提供。透過這種方式，只需在 f 成員中提供具有合適函數指標的 PyArray_Descr 結構，即可簡單地擴展資料類型的數量。對於內建類型，有一些最佳化會繞過此機制，但資料類型抽象的重點是允許新增新的資料類型。

內建資料類型之一，void 資料類型允許任意 結構化類型，其中包含 1 個或多個欄位作為陣列的元素。欄位 只是另一個資料類型物件，以及當前結構化類型中的偏移量。為了支援任意巢狀欄位，為 void 類型實作了多個資料類型存取的遞迴實作。常見的慣用語是循環遍歷字典的元素，並根據給定偏移量處儲存的資料類型物件執行特定操作。這些偏移量可以是任意數字。因此，如果需要，必須識別並考慮遇到未對齊資料的可能性。

N 維迭代器#

另請參閱

迭代陣列

在許多 NumPy 程式碼中，一個非常常見的操作是需要迭代一般、跨步、N 維陣列的所有元素。通用 N 維迴圈的此操作在迭代器物件的概念中被抽象化。若要編寫 N 維迴圈，您只需從 ndarray 建立迭代器物件，使用迭代器物件結構的 dataptr 成員，並在迭代器物件上呼叫巨集 PyArray_ITER_NEXT 以移動到下一個元素。next 元素始終以 C 連續順序排列。巨集的工作方式是首先特殊處理 C 連續、1 維和 2 維的情況，這些情況非常簡單。

對於一般情況，迭代的工作方式是在迭代器物件中追蹤座標計數器列表。在每次迭代中，最後一個座標計數器都會增加（從 0 開始）。如果此計數器小於該維度中陣列大小少 1（預先計算並儲存的值），則計數器會增加，並且 dataptr 成員會按該維度中的 strides 增加，並且巨集結束。如果達到維度的末尾，則最後一個維度的計數器會重設為零，並且 dataptr 會透過減去 strides 值乘以該維度中元素數少 1 的值（這也會預先計算並儲存在迭代器物件的 backstrides 成員中）移回該維度的開頭。在這種情況下，巨集不會結束，但會遞減本機維度計數器，以便倒數第二個維度取代最後一個維度所扮演的角色，並且先前描述的測試會在倒數第二個維度上再次執行。透過這種方式，dataptr 會針對任意 strides 進行適當調整。

coordinates 成員會維護目前的 N 維計數器，除非底層陣列是 C 連續的，在這種情況下，座標計數會被繞過。index 成員會追蹤迭代器目前的平面索引。它由 PyArray_ITER_NEXT 巨集更新。

廣播#

另請參閱

廣播

在 NumPy 的前身 Numeric 中，廣播是在埋藏在 ufuncobject.c 深處的幾行程式碼中實作的。在 NumPy 中，廣播的概念已被抽象化，因此可以在多個位置執行。廣播由函數 PyArray_Broadcast 處理。此函數需要傳入 PyArrayMultiIterObject（或二進位等效物）。PyArrayMultiIterObject 會追蹤廣播維度的數量和每個維度的大小，以及廣播結果的總大小。它還會追蹤正在廣播的陣列數量，以及指向每個正在廣播的陣列的迭代器的指標。

PyArray_Broadcast 函數採用已定義的迭代器，並使用它們來判斷每個維度中的廣播形狀（若要在廣播發生時同時建立迭代器，請使用 PyArray_MultiIterNew 函數）。然後，調整迭代器，使每個迭代器都認為它正在迭代具有廣播大小的陣列。這是透過調整迭代器的維度數量和每個維度中的 形狀 來完成的。這之所以有效，是因為迭代器 strides 也會被調整。廣播僅調整（或新增）長度為 1 的維度。對於這些維度，strides 變數會簡單地設定為 0，以便該陣列的迭代器在廣播操作在擴展維度上執行時不會移動資料指標。

廣播始終在 Numeric 中使用針對擴展維度的 0 值 strides 實作。它在 NumPy 中的執行方式完全相同。最大的不同在於，現在 strides 陣列在 PyArrayIterObject 中追蹤，廣播結果中涉及的迭代器在 PyArrayMultiIterObject 中追蹤，而 PyArray_Broadcast 呼叫實作 一般廣播規則。

陣列純量#

另請參閱

純量

陣列純量提供了一個 Python 類型層次結構，允許在陣列中儲存的資料類型與從陣列中提取元素時傳回的 Python 類型之間建立一對一的對應關係。物件陣列是此規則的例外。物件陣列是任意 Python 物件的異質集合。當您從物件陣列中選取項目時，您會取回原始 Python 物件（而不是物件陣列純量，後者確實存在，但很少用於實際目的）。

陣列純量也提供與陣列相同的方法和屬性，目的是可以使用相同的程式碼來支援任意維度（包括 0 維度）。陣列純量是唯讀的（不可變的），但 void 純量除外，後者也可以寫入，以便結構化陣列欄位設定可以更自然地運作（a[0]['f1'] = value）。

索引#

另請參閱

ndarray 的索引、索引常式

所有 Python 索引運算 arr[index] 都透過首先準備索引並找到索引類型來組織。支援的索引類型有

  • 整數

  • newaxis

  • 切片

  • 省略號

  • 整數陣列/類陣列（進階）

  • 布林值（單個布林值陣列）；如果有多個布林值陣列作為索引，或形狀不完全匹配，則布林值陣列將轉換為整數陣列。

  • 0 維布林值（以及整數）；0 維布林值陣列是一種特殊情況，必須在進階索引程式碼中處理。它們表示 0 維布林值陣列必須解釋為整數陣列。

以及純量陣列特殊情況，表示整數陣列被解釋為整數索引，這很重要，因為整數陣列索引強制複製，但如果傳回純量，則會被忽略（完整整數索引）。準備好的索引保證有效，但超出邊界的值和進階索引的廣播錯誤除外。這包括為不完整的索引新增 Ellipsis，例如，當使用單個整數為二維陣列建立索引時。

下一步取決於找到的索引類型。如果所有維度都使用整數建立索引，則會傳回或設定純量。單個布林值索引陣列將呼叫專門的布林值函數。包含 Ellipsis 或 切片 但沒有進階索引的索引將始終透過計算新的 strides 和記憶體偏移量來建立舊陣列的檢視。然後，可以傳回此檢視，或者對於賦值，可以使用 PyArray_CopyObject 填寫。請注意，PyArray_CopyObject 也可能在其他分支中的臨時陣列上呼叫，以支援當陣列的物件 dtype 時的複雜賦值。

進階索引#

到目前為止，最複雜的情況是進階索引，它可能與也可能不與典型的基於檢視的索引結合使用。此處的整數索引被解釋為基於檢視。在嘗試理解這一點之前，您可能需要先熟悉它的微妙之處。進階索引程式碼有三個不同的分支和一個特殊情況

  • 有一個索引陣列，並且索引陣列和賦值陣列都可以輕鬆迭代。例如，它們可能是連續的。此外，索引陣列必須是 intp 類型，並且賦值中的值陣列應為正確的類型。這純粹是一條快速路徑。

  • 只有整數陣列索引，因此不存在子陣列。

  • 基於檢視的索引和進階索引混合在一起。在這種情況下，基於檢視的索引定義了一組子陣列，這些子陣列由進階索引組合在一起。例如，arr[[1, 2, 3], :] 是透過垂直堆疊子陣列 arr[1, :]、arr[2, :] 和 arr[3, :] 來建立的。

  • 存在子陣列，但它只有一個元素。這種情況可以像沒有子陣列一樣處理，但在設定期間需要小心。

判斷適用哪種情況、檢查廣播以及判斷所需的轉置種類都在 PyArray_MapIterNew 中完成。設定完成後，有兩種情況。如果沒有子陣列或只有一個元素，則不需要子陣列迭代，並且準備一個迭代器，該迭代器迭代所有索引陣列以及結果或值陣列。如果存在子陣列，則會準備三個迭代器。一個用於索引陣列，一個用於結果或值陣列（減去其子陣列），一個用於原始陣列和結果/賦值陣列的子陣列。前兩個迭代器給出（或允許計算）子陣列起點的指標，然後允許重新啟動子陣列迭代。

當進階索引彼此相鄰時，可能需要轉置。所有必要的轉置都由 PyArray_MapIterSwapAxes 處理，並且必須由呼叫者處理，除非要求 PyArray_MapIterNew 分配結果。

準備完成後，取得和設定相對簡單，儘管需要考慮不同的迭代模式。除非在項目取得期間只有一個索引陣列，否則會預先檢查索引的有效性。否則，為了最佳化，它會在內部迴圈本身中處理。

通用函數#

另請參閱

通用函數 (ufunc)、通用函數 (ufunc) 基礎知識

通用函數是可呼叫物件，它透過將逐元素運作的基本一維迴圈封裝到易於使用的完整函數中，來取得 \(N\) 個輸入並產生 \(M\) 個輸出，這些函數無縫實作廣播、類型檢查、緩衝強制型轉 和 輸出引數處理。新的通用函數通常在 C 中建立，儘管有一種從 Python 函數建立 ufunc 的機制 (frompyfunc)。使用者必須提供一個一維迴圈，該迴圈實作基本函數，該函數接受輸入純量值並將產生的純量放置到適當的輸出槽中，如實作中所述。

設定#

每個 ufunc 計算都涉及一些與設定計算相關的額外負荷。此額外負荷的實際意義在於，即使 ufunc 的實際計算非常快，您也將能夠編寫陣列和類型特定的程式碼，這些程式碼對於小型陣列的運作速度會比 ufunc 更快。特別是，使用 ufunc 對 0 維陣列執行許多計算會比其他基於 Python 的解決方案慢（靜默匯入的 scalarmath 模組的存在正是為了讓陣列純量具有基於 ufunc 計算的外觀和感覺，並顯著減少了額外負荷）。

當呼叫 ufunc 時，必須完成許多事情。從這些設定操作收集的資訊儲存在迴圈物件中。此迴圈物件是一個 C 結構（可能會變成 Python 物件，但未初始化為此類物件，因為它僅在內部使用）。此迴圈物件具有與 PyArray_Broadcast 一起使用的佈局，以便可以像在程式碼的其他部分中處理廣播一樣處理廣播。

首先完成的事情是在執行緒特定的全域字典中查閱緩衝區大小、錯誤遮罩和關聯的錯誤物件的目前值。錯誤遮罩的狀態控制在找到錯誤條件時會發生什麼情況。應注意，硬體錯誤旗標的檢查僅在每個一維迴圈執行後執行。這表示，如果輸入和輸出陣列是連續的且類型正確，以便執行單個一維迴圈，則可能要到陣列的所有元素都已計算完畢後才會檢查旗標。在執行緒特定的字典中查閱這些值需要時間，對於除非常小的陣列之外的所有陣列，這段時間都很容易忽略。

檢查執行緒特定的全域變數後，會評估輸入以判斷 ufunc 應如何繼續，並在必要時建構輸入和輸出陣列。任何不是陣列的輸入都會轉換為陣列（必要時使用上下文）。會記錄哪些輸入是純量（因此轉換為 0 維陣列）。

接下來，根據輸入陣列類型，從 ufunc 可用的一維迴圈中選取適當的一維迴圈。此一維迴圈是透過嘗試將輸入的資料類型簽章與可用的簽章進行比對來選取的。與內建類型對應的簽章儲存在 ufunc 結構的 ufunc.types 成員中。與使用者定義類型對應的簽章儲存在函數資訊的連結列表中，頭部元素儲存為 CObject 在 userloops 字典中，索引鍵為資料類型編號（引數列表中第一個使用者定義的類型用作索引鍵）。搜尋簽章，直到找到一個簽章，輸入陣列都可以安全地強制型轉為該簽章（忽略任何不允許判斷結果類型的純量引數）。此搜尋程序的含義是，當儲存簽章時，「較小的類型」應放置在「較大的類型」下方。如果找不到一維迴圈，則會報告錯誤。否則，argument_list 會使用儲存的簽章更新 — 以防需要強制型轉，並修正一維迴圈假設的輸出類型。

如果 ufunc 有 2 個輸入和 1 個輸出，且第二個輸入是 Object 陣列，則會執行特殊情況檢查，以便在第二個輸入不是 ndarray、具有 __array_priority__ 屬性且具有 __r{op}__ 特殊方法時，傳回 NotImplemented。透過這種方式，Python 會收到訊號，讓另一個物件有機會完成運算，而不是使用泛型物件陣列計算。這允許（例如）稀疏矩陣覆寫乘法運算子一維迴圈。

對於小於指定緩衝區大小的輸入陣列，會複製所有非連續、未對齊或位元組順序錯誤的陣列，以確保對於小型陣列，使用單個迴圈。然後，為所有輸入陣列建立陣列迭代器，並將產生的迭代器集合廣播為單個形狀。

然後處理輸出引數（如果有），並建構任何遺失的傳回陣列。如果任何提供的輸出陣列沒有正確的類型（或未對齊）且小於緩衝區大小，則會建構一個新的輸出陣列，並設定特殊的 NPY_ARRAY_WRITEBACKIFCOPY 旗標。在函數結束時，會呼叫 PyArray_ResolveWritebackIfCopy，以便其內容將複製回輸出陣列。然後處理輸出引數的迭代器。

最後，會決定如何執行迴圈機制，以確保輸入陣列的所有元素都組合起來以產生正確類型的輸出陣列。迴圈執行的選項包括單迴圈（適用於 :term:`contiguous`、對齊且資料類型正確）、跨步迴圈（適用於非連續但仍對齊且資料類型正確）和緩衝迴圈（適用於未對齊或資料類型錯誤的情況）。根據呼叫哪種執行方法，然後設定並計算迴圈。

函數呼叫#

本節說明如何針對三種不同的執行種類中的每一種設定和執行基本通用函數計算迴圈。如果在編譯期間定義了 NPY_ALLOW_THREADS，則只要不涉及物件陣列，Python 全域解譯器鎖定 (GIL) 就會在呼叫迴圈之前釋放。如果需要處理錯誤條件，則會重新取得它。硬體錯誤旗標僅在一維迴圈完成後檢查。

單迴圈#

這是所有情況中最簡單的案例。ufunc 的執行方式是精確呼叫底層的 1 維迴圈一次。只有在以下情況才有可能：輸入和輸出的對齊資料具有正確的型別（包括位元組順序），且所有陣列都具有一致的步幅（無論是連續的、0 維或 1 維）。在這種情況下，1 維計算迴圈會被呼叫一次，以計算整個陣列的結果。請注意，硬體錯誤標 flags 僅在整個計算完成後才會檢查。

跨步迴圈#

當輸入和輸出陣列對齊且型別正確，但步幅不一致（非連續且為 2 維或更高維度）時，則會採用第二種迴圈結構進行計算。此方法會轉換輸入和輸出引數的所有迭代器，使其迭代除了最大維度之外的所有維度。然後，內部迴圈由底層的 1 維計算迴圈處理。外部迴圈是轉換後迭代器上的標準迭代器迴圈。硬體錯誤 flags 會在每次 1 維迴圈完成後檢查。

緩衝迴圈#

當輸入和/或輸出陣列與底層 1 維迴圈預期的資料未對齊或資料型別錯誤（包括位元組交換）時，會使用此程式碼處理這種情況。陣列也被假設為非連續的。此程式碼的工作方式非常類似於跨步迴圈，不同之處在於內部 1 維迴圈經過修改，以便在輸入上執行預處理，並在輸出上以 bufsize 區塊（其中 bufsize 是使用者可設定的參數）執行後處理。底層 1 維計算迴圈會在複製過來的資料上呼叫（如果需要複製）。在這種情況下，設定程式碼和迴圈程式碼會複雜得多，因為它必須處理

  • 暫時緩衝區的記憶體配置

  • 決定是否在輸入和輸出資料上使用緩衝區（未對齊和/或資料型別錯誤）

  • 複製以及可能轉換任何需要緩衝區的輸入或輸出的資料。

  • 特殊處理 Object 陣列，以便在需要複製和/或轉換時正確處理參考計數。

  • 將內部 1 維迴圈分解為 bufsize 區塊（可能會有餘數）。

同樣，硬體錯誤 flags 會在每個 1 維迴圈結束時檢查。

最終輸出操作#

Ufuncs 允許其他類似陣列的類別無縫地通過介面傳遞，特定類別的輸入將導致輸出為相同的類別。其運作機制如下。如果任何輸入不是 ndarrays 並且定義了 __array_wrap__ 方法，則具有最大 __array_priority__ 屬性的類別決定所有輸出的型別（傳入的任何輸出陣列除外）。輸入陣列的 __array_wrap__ 方法將會被呼叫，並將從 ufunc 返回的 ndarray 作為其輸入。支援 __array_wrap__ 函數的兩種呼叫樣式。第一種樣式將 ndarray 作為第一個引數，將 “context” tuple 作為第二個引數。context 為 (ufunc、arguments、輸出引數編號)。這是第一次嘗試呼叫。如果發生 TypeError，則只會使用 ndarray 作為第一個引數來呼叫該函數。

方法#

ufunc 的三種方法需要類似於通用 ufuncs 的計算。這些方法是 ufunc.reduce、ufunc.accumulate 和 ufunc.reduceat。這些方法中的每一種都需要一個設定命令，後跟一個迴圈。這些方法有四種可能的迴圈樣式，分別對應於無元素、單元素、跨步迴圈和緩衝迴圈。這些迴圈樣式與為通用函數呼叫實作的基本迴圈樣式相同，只是無元素和單元素的情況是特殊情況，分別在輸入陣列物件有 0 個和 1 個元素時發生。

設定#

所有三種方法的設定函數都是 construct_reduce。此函數會建立一個歸約迴圈物件，並填入完成迴圈所需的參數。所有方法僅適用於接受 2 個輸入並傳回 1 個輸出的 ufuncs。因此，底層 1 維迴圈的選擇假設簽名為 [otype, otype, otype]，其中 otype 是請求的歸約資料型別。然後從（每個執行緒）全域儲存空間檢索緩衝區大小和錯誤處理。對於未對齊或資料型別不正確的小型陣列，會建立副本，以便使用未緩衝的程式碼區段。然後，選擇迴圈策略。如果陣列中有 1 個或 0 個元素，則選擇簡單的迴圈方法。如果陣列未未對齊且具有正確的資料型別，則選擇跨步迴圈。否則，必須執行緩衝迴圈。然後建立迴圈參數，並建構傳回陣列。輸出陣列的形狀不同，具體取決於方法是 reduce、accumulate 還是 reduceat。如果已提供輸出陣列，則會檢查其形狀。如果輸出陣列不是 C 連續的、已對齊且具有正確的資料型別，則會使用設定了 NPY_ARRAY_WRITEBACKIFCOPY 標 flags 建立暫時副本。這樣，這些方法將能夠使用良好行為的輸出陣列，但當在函數完成時呼叫 PyArray_ResolveWritebackIfCopy 時，結果將會複製回真正的輸出陣列。最後，迭代器設定為在正確的軸上迴圈（取決於提供給方法的軸的值），並且設定常式返回到實際的計算常式。

Reduce#

所有 ufunc 方法都使用相同的底層 1 維計算迴圈，並調整輸入和輸出引數，以便進行適當的歸約。例如，reduce 運作的關鍵在於，1 維迴圈被呼叫時，輸出和第二個輸入指向記憶體中的相同位置，並且兩者都具有 0 的步幅。第一個輸入指向輸入陣列，其步幅由所選軸的適當步幅給定。透過這種方式，執行的操作是

\begin{align*} o & = & i[0] \\ o & = & i[k]\textrm{<op>}o\quad k=1\ldots N \end{align*}

其中 \(N+1\) 是輸入 \(i\) 中的元素數量，\(o\) 是輸出，而 \(i[k]\) 是沿所選軸的 \(i\) 的第 \(k^{\textrm{th}}\) 個元素。對於維度大於 1 的陣列，此基本操作會重複執行，以便沿所選軸對每個 1 維子陣列進行歸約。已移除所選維度的迭代器會處理此迴圈。

對於緩衝迴圈，必須注意在呼叫迴圈函數之前複製和轉換資料，因為底層迴圈期望對齊且資料型別正確（包括位元組順序）的資料。緩衝迴圈必須在以不大於使用者指定的 bufsize 的區塊呼叫迴圈函數之前，處理此複製和轉換。

Accumulate#

accumulate 方法與 reduce 方法非常相似，因為輸出和第二個輸入都指向輸出。不同之處在於，第二個輸入指向記憶體，其步幅比目前的輸出指標落後一個步幅。因此，執行的操作是

\begin{align*} o[0] & = & i[0] \\ o[k] & = & i[k]\textrm{<op>}o[k-1]\quad k=1\ldots N. \end{align*}

輸出具有與輸入相同的形狀，並且當所選軸的形狀為 \(N+1\) 時，每個 1 維迴圈會對 \(N\) 個元素進行操作。同樣，緩衝迴圈會注意在呼叫底層 1 維計算迴圈之前複製和轉換資料。

Reduceat#

reduceat 函數是 reduce 和 accumulate 函數的泛化。它實作了在索引指定的輸入陣列範圍上的 reduce。在迴圈計算發生之前，會檢查額外的索引引數，以確保每個輸入對於沿所選維度的輸入陣列而言不會太大。迴圈實作是使用與 reduce 程式碼非常相似的程式碼處理的，該程式碼重複的次數與索引輸入中的元素數量相同。特別是：傳遞到底層 1 維計算迴圈的第一個輸入指標指向索引陣列指示的正確位置的輸入陣列。此外，傳遞到底層 1 維迴圈的輸出指標和第二個輸入指標指向記憶體中的相同位置。1 維計算迴圈的大小固定為目前索引與下一個索引之間的差值（當目前索引是最後一個索引時，則下一個索引假設為沿所選維度的陣列長度）。透過這種方式，1 維迴圈將在指定的索引上實作 reduce。

未對齊或迴圈資料型別與輸入和/或輸出資料型別不符的情況，會使用緩衝程式碼處理，其中資料會複製到暫時緩衝區，並在呼叫底層 1 維函數之前轉換為正確的資料型別（如果必要）。暫時緩衝區是以不超過使用者可設定的緩衝區大小值的（元素）大小建立的。因此，迴圈必須足夠靈活，才能呼叫底層 1 維計算迴圈足夠的次數，以在不超過緩衝區大小的區塊中完成總計算。