Author 作者: Patrick Mineault，原链接 Rewriting Matlab code in Python，https://xcorr.net/2020/03/04/rewriting-matlab-code-in-python/
Translator 译者: Wenyin Wei

用 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');
end

你可以通过返回结果而不是保存结果，将其转换为无副作用的函数：

function [data] = mycomplexop()
    data = zeros(100, 100);
    % Complex things happen in the middle...
    return
end

我更倾向于注释掉绘图代码，而不是尝试完全重现它。原因还是那个理由，我们要聚焦在数字运算上。

确定测试用例 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);
end

我们可以重写此函数，以使其随机性源于该函数之外：

function [dist] = rangedist(A)
    dist = max(A, 1) - min(A, 1);
end

然后，我们可以在脚本中生成变量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]
A(:)

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 代码中找到了细微的错误！将来，你会感谢自己对细节的重视。