用於計算分子生物學的 Python 工具
用於計算分子生物學的 Python 工具
BioPython 是生物資訊學和計算生物學中非常流行的函式庫。Mocapy++ 是一個機器學習工具包，用於動態貝氏網路 (DBN) 中的參數學習和推論PaluszewskiHamelryck2010，它編碼了領域中隨機變數之間的機率關係。Mocapy++ 在 GNU 通用公共許可證 (GPL) 下免費提供，可從 SourceForge 取得。該函式庫支援廣泛的 DBN 架構和機率分佈，包括方向統計學的分佈。值得注意的是，球體上的肯特分佈和環面上的雙變數馮·米塞斯分佈，已被證明在制定蛋白質和 RNA 結構的機率模型中很有用。
這樣一個非常有用且功能強大的函式庫，已成功應用於 TorusDBNBMTFKH2008、BasiliskHBPFJH2010、FB5HMMHKK2006PaluszewskiWinter2008 等專案中，是長期努力的成果。Mocapy 的原始實作可以追溯到 2004 年，從那時起，該函式庫已用 C++ 重寫。然而，C++ 是一種靜態類型和編譯的程式語言，不利於快速原型設計。因此，目前的 Mocapy++ 沒有動態載入自訂節點類型的功能，並且需要一個機制來插入新的節點類型，而無需修改和重新編譯函式庫，這是一個值得關注的問題。這樣一個外掛介面將有助於快速原型設計，允許快速實作和測試新的機率分佈，進而可以大幅減少開發時間和精力；使用者將能夠在不修改和後續重新編譯的情況下擴展 Mocapy++。意識到這種需求，T. Hamelryck 提出了該專案 (以下稱為 MocapyEXT)，旨在改進目前的 Mocapy++ 節點類型擴充機制。
MocapyEXT 專案主要是一個工程努力，旨在為 Mocapy++ 帶來一個透明的 Python 外掛介面，其中內建和動態載入的節點類型可以以統一的方式使用。此外，外部實作和動態載入的節點可以由使用者修改，並且這些變更不需要重新編譯客戶端程式，也不需要重新編譯隨附的 Mocapy++ 函式庫。這將有助於快速原型設計，簡化現有程式碼的調整，並提高軟體互通性，同時對現有的 Mocapy++ 介面進行最小的變更，從而有助於順利接受 MocapyEXT 引入的變更。
Justinas V. Daugmaudis (vygis.d@gmail.com)
指導老師
Thomas Hamelryck
Eric Talevich
'''了解 S-EM 和方向統計'''
'''研究 Mocapy++ 的使用案例'''
'''研究 Mocapy++ 內部結構和程式碼'''
'''設計 Mocapy++ 外掛介面'''
'''實作 Mocapy++ 外掛模組'''
'''實驗模組化的 Mocapy++ 架構'''
'''第 1-2 週 [5 月 23 日 – 6 月 5 日]'''
介面策略設計：對於 Mocapy++ 使用者來說，外掛 API 感覺「自然」的重要性再怎麼強調也不為過。大部分時間將用於介面設計。
'''第 3-5 週 [6 月 6 日 – 6 月 19 日]'''
實作外掛模組。
'''第 6-7 週 [6 月 20 日 – 6 月 30 日]'''
實作一些範例，以展示 Mocapy++ 外掛系統的實際應用；例如，如何使用外部實作的邏輯分佈。
'''第 7-8 週（期中）[7 月 1 日 – 7 月 10 日]'''
'''第 9-10 週 [7 月 11 日 – 7 月 24 日]'''
範例應用程式。應著重於外掛模組的類型反映能力。
'''第 11-12 週 [7 月 25 日 – 8 月 7 日]'''
更新文件以反映新功能。同時記錄範例，進行任何程式碼清理工作等。計畫的「停止撰寫」日期。
'''第 13-14 週 [8 月 8 日 – 8 月 21 日]'''
向更廣泛的受眾介紹綁定，收集社群中的意見和評論。
'''第 15 週 [8 月 22 日]'''
專案結束。
託管於 gSoC11 Mocapy 分支
對於 Mocapy++ 術語不熟悉的人來說，必須了解什麼是 ESS 電腦和密度，因為這些術語將在本文中通篇使用。
在下表中，X 表示 ESS 電腦類別，而 u 是 X 的可變值。
| 表示式 | 傳回類型 | 前提條件/後置條件 |
|---|---|---|
| u.construct(parent_sizes) | void | 定義了 ESS 的適當形狀 |
| u.estimate(ess) | void | 將樣本點新增至 ESS |
類別 mocapy::BippoESS 是 ESS 電腦的範例模型。
在下表中，X 表示密度電腦類別，v 是 X 的常數值，而 u 是 X 的可變值。請注意，ptv 代表「父節點和此節點的值」：父節點的值以及該方法所屬節點的值。
| 表示式 | 傳回類型 | 前提條件/後置條件 | |
|---|---|---|---|
| u.construct(parent_sizes) | void | 初始化參數（平均值、共變異數、CPD 等） | |
| u.estimate(ess) | void | 根據給定的 ESS 估計節點的參數 | |
| u.sample(ptv) | std::vector |
根據指示的父節點值傳回樣本。請注意，後續呼叫 sample 成員函式可能會傳回不同的值，因此此操作是可變的 | |
| v.get_lik(ptv, log) | double | 傳回似然值，P(子節點 | 父節點) |
| v.get_parameters() | std::vector<MDArray |
傳回分佈參數 | |
| v.get_node_size() | unsigned int | 傳回節點大小 | |
| v.get_output_size() | unsigned int | 傳回節點的輸出大小 |
類別 mocapy::BippoDensities 是密度電腦的範例模型。
為了讓客戶端程式能夠執行 Python 程式碼，必須初始化指令碼環境。這是通過調用 Py_Initialize() 函式來完成的。解譯器通過調用 Py_Finalize() 函式來釋放。
然而，必須注意，Py_Initialize() 和 Py_Finalize() 呼叫之間的任何陳述都可能拋出例外。如果拋出例外，則必須在 try/catch 區塊中處理，否則必須終止程式。考慮到這一點，可以通過將 Py_Initialize() 和 Py_Finalize() 之間的陳述包裝在 try/catch 區塊中，使先前的範例程式更加強大。Py_Finalize() 的安全性再怎麼強調也不為過。目前，Boost.Python 有多個全域（或函式靜態）物件，它們的存在使得引用計數無法降至零，直到卸載 Boost.Python 共用物件為止。這可能會導致崩潰，因為當引用計數歸零時，沒有解譯器。這提出了一個問題，即這種初始化 Python 解譯器的方法是否可以被認為是「易於使用」、「安全」甚至是「非侵入式」的，而這些都是 MocapyEXT 的主要設計宗旨。
這個問題的解決方案是易於使用、安全且非侵入式的，令人驚訝的是它非常優雅，並展示了 RAII（資源獲取即初始化）的慣用法。
Python 解譯器在進入 main 函式之前初始化，並在退出 main 函式之後釋放。不需要額外的終端使用者努力，只需包含定義必要 Python 外掛封裝函式的標頭即可。
Python 解譯器應連結到客戶端程式（而不是 Mocapy++ 函式庫）。MocapyEXT 的終端使用者應負責額外的包含和/或程式庫路徑，以指向成功編譯和連結客戶端程式所必需的 Python 標頭和程式庫。
MocapyEXT 也以這樣的方式實例化其靜態資料成員，即它們僅實例化一次，並且在程式庫的編譯單元中實例化。C++ 標準中明確指出
「…特別是，除非靜態資料成員本身以需要靜態資料成員定義存在的方式使用，否則不會發生靜態資料成員的初始化（以及任何相關的副作用）。」[cpp_std2003, temp.inst]
使用者不需要管理實例化的靜態資料成員。
MocapyEXT 外掛具有以下執行緒安全性保證
對傳入的參數容器進行同時可變存取需要同步處理，例如透過讀寫鎖。可變存取包括修改容器內的值、呼叫會使迭代器失效的成員函式，以及透過移動建構函式移動容器。
當匯入包含 ESS 或密度 (Densities) 定義的模組 foo 時，內嵌的 Python 解譯器會在環境變數 PYTHONPATH 和目前路徑所指定的目錄列表中搜尋名為 foo.py 的檔案。
在繼續匯入模組之前，會將目前路徑明確附加到 sys.path 列表。如果使用者依賴腳本存在於客戶端程式以外的其他目錄中，則應在 PYTHONPATH 環境變數中列出這些路徑。
請注意，由於解譯器也會在與安裝相關的預設路徑中搜尋，因此包含 ESS 或密度定義的檔案名稱不應與標準模組的名稱相同，這一點非常重要。
最小介面匯入範例顯示了三個外掛程式類別的用法：densities_plugin、ess_plugin 和 aggregate_plugin。densities_plugin 和 ess_plugin 的預期用法是在 ESS 和密度類型於不同檔案中實作的情況下，即每個檔案一個類別。每個類別都會在其自己的解譯器環境中載入。aggregate_plugin 的目的是簡化載入類型的管理，並最佳化為內嵌 Python 解譯器配置的資源，因為 aggregate_plugin 只會為 ESS 和密度類型建立一個內嵌解譯器的執行個體。
#include <mocapy/plugin/ess_plugin.hpp> #include <mocapy/plugin/densities_plugin.hpp> #include <mocapy/plugin/aggregate_plugin.hpp> int main() { using namespace mocapy::ext; densities_plugin dens_pl("plugin_tests", "DensitiesPython"); ess_plugin ess_pl("plugin_tests", "ESSPython"); aggregate_plugin aggr_pl("plugin_tests", "ESSPython", "DensitiesPython"); return 0; }基本上，使用者會提供 Python 程式庫的名稱、類別的名稱，以及建構類別模型所需的實際引數清單。然後，MocapyExt 程式庫會傳回新類別的執行個體。
densities_plugin p("test_module", "DensitiesClassName"); densities_adapter n = p.densities(arg1, arg2, ..., argN);給定的引數列表（最多支援 N = 6 個引數）將參數化密度的物件初始化。參數會被轉送到 Python API，並使用 Boost.Python 程式庫提供的類型轉換功能進行轉換。這個機制是通用的，使用者可以將任何任意類型 T 轉換為 Python 中對應的類型。
參數會透過 const-reference 轉送。這個方法接受並轉送任意類型的引數，但代價是始終將引數視為 const。這個解決方案通常用於建構函式的引數。這種方法的一個深奧問題是，無法形成函式類型的 const reference，但這個缺陷（正在核心問題 295 中處理），實際上對我們有利，因為這可以防止惡意嘗試將函式傳遞給 Python 解譯器。
任何進一步呼叫 ess 和 dens 成員函式的行為，都會自動呼叫 Python 中對應的類別方法。
ESS 計算機和密度物件 (分別為 ess 和 dens) 的生命週期可能會超過其各自的工廠外掛程式的生命週期。物件 ess 和 dens 會保留對 Python 解譯器的參考計數指標，因此，有效的 Python 解譯器執行個體會一直存在，直到物件 ess 和 dens 被銷毀為止。
這是一個簡單的範例，展示了兩種不同的方法來初始化外掛程式節點。在這個特定案例中，程式中使用的模組 plugin_tests 實作了一個固定長度為 1 的虛擬離散節點，並且始終傳回 [0,] 作為取樣值。
#include <mocapy/plugin/ess_plugin.hpp> #include <mocapy/plugin/densities_plugin.hpp> #include <mocapy/plugin/aggregate_plugin.hpp> #include <mocapy/plugin/plugin_node.hpp> int main() { using namespace mocapy::ext; // 從兩個單獨載入的模組初始化節點 { ess_plugin ess("plugin_tests", "ESSPython"); densities_plugin dens("plugin_tests", "DensitiesPython"); plugin_node_type node(ess.ess(arg1, arg2, ..., argN), dens.densities(arg1, arg2, ..., argN)); } // 從單一載入的模組初始化節點 { aggregate_plugin pl("plugin_tests", "ESSPython", "DensitiesPython"); plugin_node_type node(pl.ess(arg1, arg2, ..., argN), pl.densities(arg1, arg2, ..., argN)); } return 0; }plugin_node_type 是 ChildNode 類別範本的 typedef。在此我們也注意到，節點的生命週期與外掛程式工廠的生命週期無關。
外掛程式節點是可串流的，也就是說，它可以透過運算子輸出到 std::ostream：
輸出的內容與 repr(Z) 相同，其中 Z 是 ESS 計算機或密度物件。
MocapyEXT 保留了舊的 Mocapy++ ChildNode 類別範本序列化行為，這表示對於不需要序列化的 ESS 計算機（包括 kentess、gaussianess、poissoness 等），不會發生序列化。但是，MocapyEXT 必須序列化實作 ESS 計算機的 Python 外掛程式，才能保留模組和類別的名稱，以便稍後可以將 ESS 計算機反序列化。
還必須注意的是，反序列化後，即使 ESS 計算機和密度在同一個模組中實作，它們也不會共用同一個 Python 解譯器執行個體。
測量外掛程式介面與 ESS 和密度計算機本身的實作相比，成本有多高，這很有趣。
我們測量名為 N、S 和 A 的測試的相對效能。所有測試都會呼叫 ESS 和密度計算機的成員函式，儘管沒有執行任何特定的計算。
N 測試代表 N(ative)，這表示在測試中使用了 ESS 和密度計算機的純 C++ 實作。
S 測試代表 S(eparate)；Python 類別會在個別的 Python 解譯器執行個體中載入。
A 測試代表 A(ggregate)；ESS 和密度計算機會在同一個 Python 解譯器執行個體中載入。
信賴區間。 4.53e+04 5.03e+04 5.53e+04
基於 Wilcoxon 配對信賴區間的加速百分比。
函式 A vs 函式 B 最小值 中位數 最大值 (快 %) (快 %) (快 %) N vs S 717 719 724 N vs A 718 725 730 A vs S -0.673 -0.251 0.0774
測試結果顯示，使用 MocapyEXT 外掛程式介面會產生一些效能損失。但是必須注意的是，測試也執行了重複建構 ESS 和密度執行個體的行為。
信賴區間。 1.94e+04 2.15e+04 2.37e+04
基於 Wilcoxon 配對信賴區間的加速百分比。
函式 A vs 函式 B 最小值 中位數 最大值 (快 %) (快 %) (快 %) N vs S 240 242 243 N vs A 239 243 245 A vs S -1.24 -0.51 -0.25
在沒有建構 ESS 和密度物件的情況下重複測試顯示，透過 MocapyEXT 介面呼叫方法只比呼叫本機實作的 ESS 和密度物件的成員函式慢約 2.5 倍。顯然，在實際情況中，與方法本身的呼叫相比，節點建構、取樣、概似計算等更有可能構成執行時間的主要部分。基本上，這些測試顯示，成員函式呼叫所花費的時間，甚至比節點建構這種「輕量級」的操作還要少得多。原則上，現在可以編寫相當通用的演算法。請考慮這個範例
template inline void do_something_with_densities(Densities dens) // 通用函式 { std::vector ps(1, 2); dens.construct(ps); /* 在此執行某些操作 */ } // 透過 Python 介面使用 Bippo 密度 densities_plugin dens("mocapy.bippo", "BippoDensities"); do_something_with_densities( dens.densities(lambdas, thetas, ns) ); // 直接使用 Bippo 密度 mocapy::BippoDensities bd(lambdas, thetas, ns); do_something_with_densities(bd); 我們可以合理地確定，使用 MocapyEXT 的效能損失雖然不能忽略，但並不會超過通用性所帶來的附加價值。而在純 C++ 的情況下，它不會產生任何效能成本。 ### 未來規劃 MocapyEXT 是一個正在進行中的專案。在目前狀態下，它是可用的，但是，必須強調的是，所描述的功能絕非完整，並且未來可能會以不同的方式運作。為了支援 MocapyEXT 基礎結構，Mocapy++ 中進行了許多變更。最值得注意的是，Mocapy++ 程式庫在開發分支中已部分修正為 const 正確。const 正確性的重要性再怎麼強調也不為過，因為它會影響從 Python 程式碼傳遞參數的協定。但是，這些變更是以透明的方式完成的。如果終端使用者不依賴現在不可變運算的副作用，則應該能夠以較小的修改來編譯和執行程式碼。毫無疑問，還有許多其他變更在等待。例如，在某些情況下，如果能夠透過將參數傳遞給建構函式來複製密度，會很有用。