完全用Python工作——F2PY简明使用指南

前言

本文是计划中的“完全使用Python工作”系列文章中的第一篇，这个系列的目的是将Python这一强大的工具引入海洋科学领域的绘图和计算中。先前打算本学期的海洋要素计算作业都基于F2PY工具来实现，但是相关的中文资料不多，文档也有些过时。于是在完成了第一份比较复杂点的作业之后，发篇博客总结一下用法以备日后查阅。后续会不定期更新，如果感兴趣的话请移步本人博客查看。

F2PY这个东东是一个绝佳的将能够处理高性能计算的Fortran程序与功能强大且灵活的Python语言融合为一体的胶水模块。我们知道MATLAB本质上就是一个包含已经编译好的高性能底层数学计算模块的工具箱，搭配上用自行开发的脚本语言写成的各种实现上层功能的脚本，从而提供了一套用于工程和科学的数值运算解决方案。

在科学计算领域，Python走的也是这个路线。Python的Numpy模块提供了一种通用的多维数组对象，并为之内建了诸多操作数据的方法。在此基础上，Scipy提供了高性能底层数值计算模块，它们都是采用了一些久经考验的开源库，比如线性代数方面是用LAPACK库等等，效率也不低。

相比于完全闭源的MATLAB，Python的最大优势是它的所有模块都可以被任意地二次开发。所以像Array对象这样设计完善并且内建了常用方法的模块完全可以拿来用作通用的存储容器，成为沟通各个程序的桥梁，不管它是用编译型语言开发的计算密集模块还是普通的Python代码。这也是我的Fortran作业要用到Array对象的原因——它能够显著地提升程序的通用性，更好地让自己开发的模块与已有的模块进行功能整合。

效果示例

一般而言，应用于科学计算的Fortran代码主要都是写成子程序的形式，接受一些数组参数，然后输出也是数组参数。F2PY的作用就是把Fortran的数组和Python的Array对象进行结合，使得在Python代码中可以传递Python的Array对象作为输入参数来调用Fortran的subroutine，同时返回的也是Array对象。至于数据长度、内存布局之类的转换，都可以自动完成。这也就是说，我们不需要再为Fortran的subroutine编写主程序了，因为调用过程是在Python运行时环境中完成的。而当数据由Fortran模块处理完成并返回为Array对象之后，我们还可以继续用Numpy/Scipy中的函数对其进行进一步的处理。

1. 基本的数组输入

   subroutine dprod(x, y, n)  
   integer, intent(in) :: n  
   real(kind=8), intent(in) :: x(n)  
   real(kind=8), intent(out) :: y  
   y = 1.0  
   do i=1, n  
   y = y * x(i)  
   end do  
   end
   

   上述代码实现了输入一个Array对象数组并计算其中元素乘积的功能，并且不需要在输入时指定数组维数或者大小，所以使用起来非常方便。
   在Python中的调用方法大致如下：

   import test #导入模块  
   a = linspace(1,10,10) #生成一个行向量  
   # 输出：array([  1.,   2.,   3.,   4.,   5.,   6.,   7.,   8.,   9.,  10.])  
   test.dprod(a) #计算其乘积  
   # 输出：3628800.0
   

2. 操作多维数组：

   SUBROUTINE FOO(A,N,M)  
   INTEGER N,M,I,J  
   REAL(kind=8) A(N,M)  
   !f2py intent(in,out,copy) a  
   !f2py integer,intent(hide),depend(a) :: n=shape(a,0), m=shape(a,1)  
   DO J=1,M  
   A(1,J) = A(1,J) + 1D0  
   ENDDO  
   DO I=1,N  
   A(I,1) = A(I,1) - 1D0  
   ENDDO  
   END
   

   添加上述的编译指导语句并按照下文所述的方法编译后，得到了一个可以接受二维Array对象的Python模块，调用过程就不再展示了。

3. 基于F2PY编写的中期潮汐观测资料调和分析工具以及整个编写过程。

编译命令

F2PY模块在Linux下的使用应该不会有什么问题，但是在Windows下用就很令人纠结。首先要使用合适的编译器，我现在仅仅能够使用Mingw中的gfortran编译器对代码进行编译，这就意味着它在环境变量里的位置一定要比cygwin之类编译器的位置靠前！

亲测可用的正确编译参数：

f2py -c --fcompiler=gnu95 --compiler=mingw32 -m test test.f90

其中-m后面的test为编译出的模块名称，文件名可以指定多个，用以将多个文件中的不同子程序（subroutine）编译到一个模块中。编译完成之后的使用方法和正常的Python模块是一样的，只要import就可以随便调用了，只不过不能乱传递参数不然也会引发类型异常。当然，如果一个文件中引用了其他文件中定义的subroutine，就得在编译的时候把那个文件也包含进来，否则就会出现无法解析的外部名称了。

如果报错提示找不到某个bat批处理文件，需要设置环境变量：

VS90COMNTOOLS => %VS110COMNTOOLS%（对应VS2012）
VS90COMNTOOLS => %VS120COMNTOOLS%（对应VS2013）

指导语句

F2PY需要解析输入变量之间的依赖关系，比如哪些变量用于输入，哪些变量用于输出，哪些变量取决于输入的数组（指Array对象，下同）大小——这也是它的方便之处。但是对于变量关系复杂一点的程序，F2PY自己显然是做不到自动处理的，所以这时就需要人为地对其进行指定。

方法和区别

F2PY提供了两种方式来实现这个目的，一种是在代码中用编译指导语句进行指定；另一种则是生成一个Signature File，可以理解为是利用F2PY来parse一下Fortran的语义，生成对应的信息文件，其中包含了类似Fortran语法的接口说明。我们可以人为地修改这个文件，以使得各个变量的依赖关系与自己设想的相符。使用方法如下：

首先利用下述命令生成对应的签名文件
f2py -m test -h test.pyf test.f90
然后修改其中的变量IO属性，再执行编译
f2py -c --fcompiler=gnu95 --compiler=mingw32 test.pyf test.f90

个人比较偏好前一种做法，其优点在于代码可以直接拿来编译而不用多出一步。后一种方法可以作为对前一种方法的检验，如果代码编译出来的模块行为与预期的不符，就可以对代码生成一个签名文件，检查里面F2PY具体生成了怎样的变量依赖关系。此外，对于一个已经编译好的模块，可以直接查看它的输入输出参数，和查看普通Python函数帮助的方法是一样的，函数名后面加个问号即可。

具体语法

F2PY编译指导语句允许在Fortran77/90源码中使用F2PY签名文件的扩展属性来描述变量属性，这个功能使得我们可以跳过生成签名文件，直接对Fortran源码应用F2PY编译成Python模块。F2PY指令格式如下：
<comment char>f2py ...
其中，固定格式的fortran代码的注释字符是“cC!#”，自由格式则是“!”。对固定格式的代码而言，<comment char>必须出现在第一列，而对于自由格式，F2PY指令可以出现在文件的任意地方。
编译器会忽略掉<comment char>f2py后的所有东东，但F2PY会像普通代码行一样读取。当F2PY发现某行含有F2PY指令时，首先用5个空格替换掉指令，然后重读该行，这也就意味着注释字符后面绝对不能有空格*！

指令中的C表达式（即下文提及的<init_expr>）可以包含：

• 标准的C结构；
• math.h和Python.h中定义的函数；
• 利用给出的依赖关系，从参数列表中计算出并初始化的变量；
• 下列C++宏：
  • rank(<name>) 返回数组的维数
  • shape(<name>,<n>) 返回数组的第n维大小，n从0开始
  • len(<name>) 返回数组长度
  • size(<name>) 返回数组大小
  • slen(<name>) 返回字符串长度

扩展属性及其语法

在f90中，扩展的变量修饰属性有以下这些，它们用于调整F2PY的具体行为。

• 可选（optional）
  相应的参数被移到可选参数列表的末尾。可选参数的缺省值由<init_expr>指定。缺省值必须是有效的C语言表达式，在使用<init_expr>时，F2PY会自动将变量设为可选属性。
  可选数组参数的所有维数都必须是有界的。
• 必须（required）
  相应的参数被视为必须的参数，这是默认属性。只有当使用了<init_expr>，而又需要禁用自动optional设置时，才需要指定required属性。
• dimension(<arrayspec>)
  相应的变量视为数组，其维数由<arrayspec>指定。
  <arrayspec>是用逗号分隔的维数边界列表。
• intent(<intentspec>)
  这个参数指定变量的IO属性，是一个逗号分割并包含下列属性的表达式：
  • in
    将参数作为输入变量，并且函数不能改变参数的值。
  • inout
    这个属性指明将参数作为一个既用于输入又用于输出的变量，换句话说，一个就地输出的变量。这里的变量必须是连续的数值数组。
    通常不推荐使用intent(inout)，而代之以intent(in,out)，或者inplace属性。
  • inplace
    与上一个属性类似，但是如果数组的类型不完全匹配，或者数组非连续，则进行就地自动转换以使得类型匹配。
    通常也不推荐使用intent(inplace)，因为这个属性直接修改输入变量。但是如果输入参数仅仅是临时变量，比如它是一个数组的Slice，则当这部分内存被释放之后会导致访问无效内存。
  • out
    将变量作为返回值，并自动附加到<returned variables>这个List的末尾。这一属性会自动设置intent(hide)属性，除非显式设定了其他IO属性如in或者inout。
    缺省情况下，返回多维数组的默认内存布局是列优先的，即Fortran类型。
  • hide
    这一属性将变量从必须或者可选参数的列表中移除。一般情况下，intent(hide)仅仅在使用了intent(out)属性时，或者<init_expr>可以完全决定变量的取值时才会用到。例如：
    integer intent(hide),depend(a) :: n = len(a)
real intent(in),dimension(n) :: a
  • copy
    保证带有intent(in)属性的变量的原始内容一定会被保留。常用于与intent(in,out)属性的连接。比如这种用法：
    !f2py intent(in,out,copy) a
  • overwrite
    指明带有intent(in)属性的变量的原始内容可能会被函数所修改。
• check([<C-booleanexpr>])
  对<C-booleanexpr>求值，从而进行参数变量的一致性检验。如果返回False则抛出异常。如果没有显式使用这个属性，则F2PY会自动生成一些标准检验语句，如检验数组的大小是否符合等等。
• depend([<names>])
  声明具有该属性的参数变量依赖于<names>这个list中变量的取值，即指定变量的依赖关系，例如用于描述数组的大小依赖于某个参数，或者某个标量依赖于输入数组的大小。例如，<init_expr>可能会用到其他参数的数值，因此使用depend属性所给出的信息，F2PY可以保证所有的参数都按照正确的顺序被初始化。如果没有显式指明这一属性，F2PY会自动生成。
  如果需要手动修改F2PY自动生成的depend属性描述，注意不要破坏任何相关的关系，并且不要生成循环依赖，否则会报错。

此外，Python语言中只有函数但是没有子过程（即Fortran中的subroutine），而且所有非数组的参数均为传值调用。所以当Python调用Fortran的子过程时，修饰为“intent(out)”的变量将作为函数返回值被返回；而当多于一个变量具有这个属性时，返回的是一个tuple，包含所有的返回变量。

其他细节

双精度输入

由于某些未知原因，我使用的F2PY版本只能将以real(kind=8)声明的变量作为双精度处理，而不接受其他方式。详细的分析可以参考StackOverflow上面的这个问答。

数组大小的自动推导

同样由于某种未知原因，F2PY会强制将能够由数组大小推导出来的输入参数自动转换为隐藏参数（即不在输入参数列表中）。如果不需要这样的自动优化，理论上可以通过修改签名文件后重新进行编译来解决。