用 Python 重写 Matlab 代码 Rewriting Matlab code in Python

你已经天天都在用 Python 啦，很棒！过去你可能有一堆很重要的 Matlab 脚本，那时的 Python 可能没法干净利落地实现它们。你现在怎样才能把这份代码重写成 Python ，同时又不引入很多麻烦？我尝试了几次这样的工作，并且总结了一种基于测试驱动的（test-driven）方法，可以避免代码转换中的最大问题：

决定你要转写的内容

根据现有代码制定测试计划

用 Python 编写测试

用 Python 编写代码

不断迭代直到测试通过

为了测试该方法，我迁移了在 Susillo和Abbott（2009）中找到的脚本。它模拟了混沌神经网络，并使用递归最小二乘法对其进行训练。在实践该方法进行转写代码时，我想确保会遇到的问题（陷阱！），而我确实遇到了问题！我概述的方法能够尽早发现了转写的 bug 并解决它们。

这是我从 Matlab 过渡到 Python 的指南的第 3 部分。阅读第 1 部分：学习 Python 的课程；第 2 部分：选择环境，IDE 和软件包。

选择要转写的代码 Pick your battles

进行矩阵处理的算法程序通常是值得转写为 Python 的不错选择。重要的是，你将要转写的代码目前需要能够在 Matlab 中运行，因为我们将得运行代码并记录它当前的行为。有些代码不必转写，比如将 Mex 代码稍加改动可以很直接地用到 Python 扩展中，保持核心 C numerics 不被改变。

尽量避免依赖 Simulink，GUIDE 或艰深难懂的工具箱直接转写代码。特别是对于 GUIDE，我发现 GUIDE 的局限性强制了某种编码风格-许多全局变量和复杂的逻辑。与直接转写相比，从头重写 GUI 代码可能更好。

并非每条 Matlab 代码都需要转写成 Python。如果你的目标是存档（archive）现有工作流，以便尽管弃用仍可继续运行，那么使用 Docker 容器化 Matlab 是一个不错的选择。你还可以使用 MathWorks 官方的 Python Matlab 引擎直接从 Python 调用旧版 Matlab 代码，或通过 Shell 调用。但是，在某些情况下你确实需要转写 Matlab 代码，例如在云上运行时。

最后一点，请 首先考虑针对较小的项目 来给你小试身手。转写并不是那么容易！通过一次较小规模的练习，你能学到很多东西。我建议你为项目创建一个 git repo 并经常性地提交，因为你可能需要进行很多次迭代才能使一切运行正常。

用基于测试的策略来重写代码 Creating a test-based strategy for rewriting the code

我们现在确有一个合适的脚本或函数需要重写。我们应该马上开始编写 Python 吗？不！用另一种语言重写代码很有可能会出现这样一种情况，代码似乎能够正常工作，但却存在细微的 bugs。你可能需要花数天或数周的时间来解决这些错误，这会让你抓狂！因此，我们将采用基于测试的策略：

记录（capture）现有代码在一系列的输入下表现出的行为。

用 Python 重写代码，并确保它确实完成了我们希望它执行的操作。也就是说，对于相同的输入，它给我们相同的输出。

我们将创建单元测试以验证我们代码的行为正确，这意味着它与旧代码完全匹配。此时，你可能会说：“等等，我不只是想转移代码，我想改善它”！我们的目标是记录代码的逻辑部分，也就是数字运算。这不妨碍我们以多种方式改进参考代码，包括：

使常数可变

使代码更模块化

使用 Python 式的数据结构，例如


    pandas

用于存储时间序列。

无论如何，我们都必须在坚实的基础上进行改进。一旦复现了原先的行为，并将其 commit 到 git 库后，我们可以在以后的 commit 中对其进行改进。

确定输入和输出 Identify input and outputs

确定待重写的脚本（或函数）的 输入和输出 。Matlab 的一件非常好的事情是，函数往往是无状态的（stateless），并且不会变动函数参数（除非输入是引用对象类型，这种情况很少见）。因此，我们可以很容易地推断出函数的输入和输出：

输入是函数参数。

输出是函数的


    return

。

函数可能会有副作用（side-effect），例如绘图或保存文件。这里我们说函数的副作用是指除了将值（主要效果）返回给操作的调用者之外的可观察到的作用。如果需要转写的代码有副作用，那应该移除它们。

例如，考虑一个函数输入数据并且什么也不返回，将结果保存到磁盘，如下所示：

function [] = mycomplexop(outfile)
    data = zeros(100, 100);
    % Complex things happen in the middle...
    save(outfile, 'data');
你可以通过返回结果而不是保存结果，将其转换为无副作用的函数：
function [data] = mycomplexop()
    data = zeros(100, 100);
    % Complex things happen in the middle...
    return
我更倾向于注释掉绘图代码，而不是尝试完全重现它。原因还是那个理由，我们要聚焦在数字运算上。
确定测试用例 Identify test cases
许多 Matlab 脚本都是这样一种样式，先初始化函数的输入，调用一个或一组函数，然后用输出来作图。我们需要记录此功能集合的输入以及输出。一个简单的方法是使用 save：
% Set up the inputs. 初始化输入
t = 0:1000;
dt = .1;
g = 10;
save inputs.mat -hdf5
% Call the function 调用函数
complexoutput = mycomplexfun(t, dt, g);
save outputs.mat -hdf5
% Plotting the output 给输出绘图
plot(t, complexoutput);
我在这里使用了 Octave，但是如果你使用 Matlab，可以使用-v7.3而不是-hdfd5保存。mat以这些格式保存的文件实际上是 Python 可以读取的 hdf5 文件。
要特别注意随机化。很难完全用 Python 复制 Matlab 的随机数生成。理想情况下，函数应与随机化无关，因为具有随机输出的函数很难测试正确性。如果函数需要随机变量来实现算法，则可以通过为函数提供所需的随机数据来规避此问题。例如，假设我们有一个函数可以计算随机正态变量的分布范围：
function [dist] = rangedist()
    A = randn(100, 1000);
    dist = max(A, 1) - min(A, 1);
我们可以重写此函数，以使其随机性源于该函数之外：
function [dist] = rangedist(A)
    dist = max(A, 1) - min(A, 1);
然后，我们可以在脚本中生成变量A，并保留该特定随机抽奖。
选择不超过几秒钟即可运行完成的测试用例。你会运行很多次测试，而这些时间加起来就很多了！玩具级的数据就可以了。在简单的输入上测试我们的代码也是一个好习惯，简单的输入能够清楚地看到正确的输出，例如全零、全一等情况。请保存尽可能多的、必要的输入-输出对，以记录所需功能的行为范围。
除了端到端地（end-to-end）测试脚本外，我们还可以选择测试脚本中涉及到的中间计算。测试更多细粒的计算可以帮助你调试代码。如果你的函数已经被编写为函数，那它可能还有私有函数，这就使事情变得有些棘手。你可能需要重构一些 Matlab 代码来呈现私有函数，以便记录它们的输入和输出。
制作测试的脚手架 Writing a test scaffold
单元测试，即对给定一组输入的情况下，检查函数结果是否符合预期。使用 Python 进行测试的范围可以从简单的内联 assert 语句到复杂的自动化解决方案，该解决方案在每次将代码推送到存储库时都会检查测试是否通过。在我们的案例中，我们将使用一种轻量级的手动测试方法：将测试用例隐藏在 __name__ == '__main__' 后面。
我们将创建一个要导入的 python 模块，例如通过 import mymodule。但是，当我们在命令行上运行文件时，通过python module.py，判断语句将运行后面隐藏的代码：
import numpy as np
def my_complex_fun(A):
    # Complicated things occur. 复杂的事情放在这里
    return A * 2
def _run_tests():
    # Load inputs. 载入输入
    with h5py.File('inputs.mat', 'r') as f:
        A = np.array(f['A/value'])
    with h5py.File('outputs.mat', 'r') as f:
        B_ref = np.array(f['B/value'])
    # Call the function. 唤起函数
    B = my_complex_fun(A)
    # Check shapes are good. 检查形状是好习惯
    assert A.shape == B.shape
    # Check values are similar. 检查值是否相似
    np.testing.assert_all_close(B, B_ref)
if __name__ == '__main__':
    # Run the tests. 运行测试
    _run_tests()
测试在_run_tests()函数内部。这样，模块命名空间将不会被函数内部其他变量污染。
运行代码时，如果输出与预期不符，则测试会引发错误并停止执行。在命令行上会看到一个很大的错误，然后你知道必须修正代码。你可以通过以下方式引发错误：
使用assert语句。每当语句为 False 时，就引发错误。
使用np.testing中的方法，比如，检查 numpy 数组的元素是否全都接近参考（在容许误差内）。
使用raise语句手动引发错误。
这就是测试代码所需的全部内容。Python 中存在许多单元测试库，但是它们的核心想法是相同的：创建代码，当输出结果超出预期时，引发错误。
如何组织 Picking  Organization
在单元测试框架中，主作用（main effect)的测试比副作用要容易得多。因此，你经常会发现自己编写的是无副作用的、模块化的、函数式的代码。那很棒！我们倾向于构造较小的函数，以便测试个体组件。你的代码将更有条理。
你编写的 Python 代码可能是面向对象的。尽管 Matlab 现在对类提供了出色的支持，但是许多 Matlab 代码（尤其是较旧的代码）都避免使用类。如果你对编写面向对象的代码感到不舒服，就不用勉力为之，Pythonic 的代码不一定要使用类。当某些对象需要记住和管理自己的状态时，类才有意义。如果你感兴趣，这是有关 Python 类的教程。
编写代码—常见陷阱  Writing the Code - common gotchas
现在我们有了测试，并知道了如何组织代码，我们可以开始写代码了。Matlab 的矩阵运算自然可以转换为 Numpy。请注意以下常见陷阱：
reshape 操作很麻烦。Matlab 使用 Fortran 风格的顺序，Python 默认使用 C 风格的。这意味着在 Matlab 中：
A = [1, 2, 3; 4, 5, 6]
ans = 
    [1, 4, 2, 5, 3, 6]
而在Python中：
>>> A = np.array([[1, 2, 3], [4, 5, 6]])
>>> print(A.ravel())
[1, 2, 3, 4, 5, 6]
由于 C 和 Fortran 风格的顺序之间的区别，在 Matlab 中保存到 hdf5 而后在 Python 中加载会使轴反转。即，你得从前到后交换轴。例如，如果你在 Matlab 中保存了一个三维张量 A，然后通过 hdf5 又将其在 Python 载入，你需要调用 A.swapaxes((2, 1, 0)) 来恢复原始张量。
小心一维向量！在 Python 中，如果你具有a的一维向量，形状为(N, )，a.dot(a.T)会给你形状为(N, )的向量。它既不是标量，也不是矩阵。如果你想计算a本身的外积，请使用 a.reshape((-1, 1)).dot(a.reshape((1, -1))) ，或者最好 np.outer(a, a)。
隐藏的点积可能会带来麻烦！*操作符适用于标量和在 Matlab 载体和语义取决于操作数的尺寸是不同的。A.dot(B)取代 Matlab 的 A * B。你可能还会看到 A @ B，但它太艰深难懂了。另外，Python 中的对象按引用传递，numpy 中的一些方法直接进行原地（in plce）修改。要注意，有时，你需要显式地复制 copy() 数据。
Python 使用从 0 开始的索引而 Matlab 从 1 开始。小心！
盯着你的数据！别瞎搞！
如果你编写了一段可能会导致错误的代码，则可以使用内联 assert 语句来确保其正常工作。我就喜欢用它来检查中间数组的尺寸。
通过在适当的地方使用 pandas 的 dataframe，可以解决直接索引操作的一些问题。xarray 对张量使用命名维度，这可以帮助避免因交换 axes 和 reshape 而引起的许多错误。np.einsum 可以用易于阅读、不易出错的方式表达张量复杂的和（sum）、积 （product）操作。
你可能想要将部分数字运算直接转换为 PyTorch 或 Tensorflow，尤其是当你想利用自动微分功能的时候。可以使用 jax 或 numba 对性能敏感的代码进行 JIT 处理，或用 Cython 重写。从小处着手，并根据需要进行优化。
应用这种方法到真实的脚本上 Using the method on  a real script
我重写的示例脚本大约有165行代码。它以教学风格编写，其中散落着各种解释、逻辑和作图。从这种意义上讲，它类似于 Jupyter Notebook。我着重转写主循环的逻辑，该逻辑对网络进行正向仿真，并使用 FORCE 算法对其进行训练。
为了使它在 Octave 中运行，我不得不修改一些与绘图有关的代码，但是我能够在大约 10 分钟内使其运行。对于十年的代码来说还不错！代码以这种方式组织：
初始化变量。
运行计算主体（一个大循环）。这个主体有副作用（绘图）。
以包含计算结果的最终状态结束。
尽管代码没有被打包为函数，但是将输入挑出来是很容易的，它们是for循环开始便需要的：
x0，z0，wf，dt，M，N，ft
输出是在循环中计算的变量，即：
zt 和 w0_len
我记录了这些输入和输出作为脚本的测试用例。我注释了在主循环之后对已训练网络进行正向仿真的部分，因为它是代码中另外的逻辑上独立的组件，我可以之后再处理它。
除了该脚本端到端的结果，我还记录了网络状态第一次迭代前向传递的结果，以实现更精细的测试，另外也检查了数组尺寸。
我决定将代码写成一个类，并加上方法函数  train 和  simulate。最后代码看起来不像 Matlab 代码。但是，我已经有了测试的脚手架了，因此可以放心，它会与旧的 Matlab 代码运行得一样好。
转写很难！我遇到了很多非常精细的问题，转写了我正在工作的相当短的代码（真实逻辑也许只有20行）。这些问题都过于细节—对一维矢量的错误处理、循环中错开了一位等等。有了基于参考代码的功能测试后，我可以定位这些问题，修复它们并准确地还原一开始的功能。这花了大约 3 个小时，确实很费时间！转写代码时，我小心地设立了定义明确的目标，所以我知道什么时候能够完成，这使我得以持续。
此处提供了转写后的代码。我将代码导入到 Notebook 中并从中生成了一些图片——它真的可以跑起来！
使用测试驱动的开发，可以系统性地、一步步地将旧的 Matlab 代码逐步转移到 Python 环境中。放心，这份代码真的能够运行，因为你已经定义了什么是正常工作（because you will have defined what working means）。你的新 Python 代码还有了一个覆盖全面的测试套件。在此过程中，你甚至可能在原来的 Matlab 代码中找到了细微的错误！将来，你会感谢自己对细节的重视。