Call a Fortran function in Python

longior, non legi (tl;dr)

In this (third) post about interfacing compiled code with different programming languages, I focus on how to access Fortran functions with Python.

To this end, I divide my post in four sections. In the first section, I explain why knowing how to call shared library functions from (some might say: yet another programming language) Python may prove useful.¹ In the subsequent section, I briefly discuss the library, data, and model we will use in our Python project. The third section describes the Python code required to run the shared library’s function. The last section concludes.

About Python

Python differs in several ways from the languages discussed before.² Unlike C, Python code does not require compilation to be run. Instead, the Python interpreter takes the users’ code (or interactive commands) and executes it accordingly. Variable types are inferred by the interpreter at runtime, hence freeing the programmer from declaring them explicitly. Although interpretation and type inference substantially reduce programmers’ burden writing code, these very characteristics also lead to a slower startup and execution speed. For this reason, packages for scientific computing such as NumPy actually wrap various C and Fortran routines.

In contrast to C and Java, Python ships with quite a comprehensive standard library, including many functionalities ranging from IO, working with lists, decoding/encoding JSON, to retrieving data from URLs.³ Moreover, Python distributions include a package-manager (called pip), which makes the installation and use of additional packages very straightforward.⁴ Inter alia due to these characteristics, Python has become a rather popular programming language⁵ being used in various areas such as scientific computing/machine learning, web frameworks, and application development.

Summarizing the above, Python has been adopted in a multitude of areas and fields, yet some algorithms may be substantially slower than the same algorithms written in C, Fortran, or Java. To alleviate this issue and to get the best of both worlds – ease of code development and computational speed – I believe it to be most helpful to know how to interface compiled (Fortran) code with Python.

Project requirements

As in the previous posts, we need our Fortran linear regression library liblinreg.so⁶, the R attitude data saved as a text file, and a Python script file.

The first row of the space separated text file contains the variable names, the second and third the number of observations and variables respectively. The respondents’ response values are stored in the subsequent rows (see below).

rating complaints privileges learning raises critical advance
30 30 30 30 30 30 30
7 7 7 7 7 7 7
43 51 30 39 61 92 45
63 64 51 54 63 73 47
71 70 68 69 76 86 48
61 63 45 47 54 84 35
81 78 56 66 71 83 47

The linear model we want to estimate remains unchanged:⁷

with denoting the dependent variable rating and the design matrix (containing a column of 1 for the intercept as well as the values of the remaining variables complaints, privileges, etc.).

Python code

Similarly to previous posts, the main tasks of the Python script consist of:

1. loading and parsing the data
2. transforming the data into suitable inputs for linreg_
3. calling linreg_
4. printing results

I implemented these steps in the Python script call_fortran.py (see below). Lines 4-17 display how to import the space delimited text data into lists, where each list element corresponds to one observation in the attitudes data. Given that Fortran requires two-dimensional arrays to be layed out in column-major order (i.e. blocks of contiguous memory correspond to variables, not observations), we further need to convert the data into the appropriate format (see lines 20-28).

Still, Python’s lists are incompatible with linreg_ (which requires arrays of doubles and integer scalars). We have seen in Interfacing a Fortran library function with C that we can easily interface Fortran with C, as the latter possesses a set of widely interoperable data types. Most conveniently, Python ships with a builtin foreign function library – ctypes – which provides C – and Fortran – compatible data types. Lines 30-34 display one way to convert the data into appropriate inputs for linreg_.

The next step consists of calling linreg_ with the arguments defined above. Therefore, we first need to dynamically load our Fortran library liblinreg with ctype’s CDLL class (see line 37). This imports the library functions as methods of a CDLL class instance (called flib in the code below). To run the linear regression function, simply call the linreg_ method with the required arguments passed by reference (see lines 40-41).

Printing the regression estimates is carried out in lines 44-52.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52

from ctypes import byref, c_int, c_double, CDLL

# ---- import data ----
f = open("attitude.txt", "r")

# variable names, number of observations, number of variables
varnames = f.readline().rstrip("\n").split(" ")
varnames[0] = "intercept"  # for printing regression results
n = int(f.readline().rstrip("\n").split(" ")[0])
k = int(f.readline().rstrip("\n").split(" ")[0])

# survey data
data = [None] * n
for i in range(0, n):
  data[i] = [
    int(v) for v in f.readline().rstrip("\n").split(" ")
  ]

# ---- prepare objects to be passed to flib.linreg_() ----
fdata = [None] * n * k
for i in range(0, n):
  for j in range(0, k):
    fdata[j * n + i] = data[i][j]

y = [r[0] for r in data]
X = fdata
for i in range(0, n):
  X[i] = 1

c_y = (c_double * n)(*y)
c_X = (c_double * (n * k))(*X)
c_beta = (c_double * k)(*[0])
c_vcov = (c_double * (k * k))(*[0])
c_se = (c_double * k)(*[0])

# ---- load library ----
flib = CDLL("liblinreg.so")

# ---- run function ----
s = flib.linreg_(c_y, c_X, byref(c_int(n)), byref(c_int(k)),
		 c_beta, c_vcov, c_se)

# ---- display results ----
print("Results:")
print("=======================================")
print("           Estimate  Std. Err  z-Value ")
print("=======================================")
for i in range(len(c_beta)):
  print("%-10s   %6.3f    %6.3f   %6.3f" %
	(varnames[i], c_beta[i], c_se[i],
	 (c_beta[i] / c_se[i])))
print("=======================================")

Code execution and results

The shared library, data, and script being set, invoke the Python interpreter to run the code.

python3 call_fortran.py

Executing the script yields – quite unsurprisingly ;) – the same results that we already witnessed in the previous blog posts. For more information about the interpretation and meaning of the results, please see this post.

Results:
=======================================
	   Estimate  Std. Err  z-Value
=======================================
intercept    10.787    11.589    0.931
complaints    0.613     0.161    3.809
privileges   -0.073     0.136   -0.538
learning      0.320     0.169    1.901
raises        0.082     0.221    0.369
critical      0.038     0.147    0.261
advance      -0.217     0.178   -1.218
=======================================

Conclusion

In this third post of my series on interfacing shared libraries – which I shall finally call “How to Interface a Fortran Library with Yet Another Programming Language” – we have seen how to use the linreg Fortran function in Python. The code resembles very much what we’ve already seen with C and Java: (1) load the data, (2) convert it to input types appropriate for linreg, (3) load and run linreg, and (4) display the results. Similarly to Java, Python allows to import the space delimited test data with fewer lines of code than C. In contrast to Java however, there is no need to fetch the external JNA library,⁸ as Python’s standard library already ships with an appropriate foreign function library – ctypes.

Note: For the Python code used in this blog post, see https://git.staudtlex.de/blog/call-fortran-from-python.