LEVMpy: a python wrapper for Fortran LEVM code
Many excellent scientific programs and routines have been written in Fortran over the years since the language was first developed in the late 1950s. Even today it is often the language of choice in the scientific community for intensive numerical and high-performance computing. NASA's recent challenge to improve their FUN3D fluid dynamics is a good example of how important Fortran still is.
Wrapping Fortran code with more user friendly languages such as Python is therefore a very valuable task with libraries such as Python's scipy and numpy being prime examples. I decided to do this for the complex non-linear least squares fitting program, LEVM, that is designed for fitting impedance spectroscopy data to a vast selection of different circuit models. The code was developed by J Ross Macdonald and is popular in the field, however, it only exists as a DOS command line or basic Windows GUI program. Building a Python library would allow me to integrate it much more easily into my standard data analysis and processing scripts.
The tools to do this are relatively straightforward, in particular I used f2py to compile the Fortran code. However, there were definitely some learning points that I picked up along the way.
First, I had very little experience with Fortran coding (Fortran 77 in this case), so I had to get my head around column formatting of the code (a relic of punch cards), COMMON blocks, endless GOTO statements, and strange subroutine behaviour. But in general, it is an easy language to learn.
The next step was to separate all of the subroutines and functions from the source code and put each into an individual file. This is vital for your sanity for any reasonably sized project! Checking that all of these compiled correctly took some time since incompatible PAUSE statements had to be removed, many variables hadn't been declared explicitly, and column formatting had to be fixed. I also removed several WRITE statements since it is impossible to catch these in Python with f2py, so they only remain as info and warning statements printed to the terminal.
The code had a lot global variables stored in 25 COMMON blocks, which is not ideal by modern coding practices! Although it would be better to have these as function arguments, it was the easiest way to input and output variables without significantly changing the source code. Therefore I made sure a separate Fortran file existed with all the COMMON block variables so they could be accessed from the Python module once f2py is run.
There were also several subroutines and functions that had input variables with assumed size arrays, e.g. DIMENSION N(*). In order to wrap to numpy arrays correctly, I changed those that would be used such that they were of fixed dimension. This usually meant using the global maximum array sizes set in the SIZE.INC file.
Most importantly, there is the issue of SAVE statements. In the original code several subroutines relied on static storage and so SAVE statements were not included (for counter variables for example). This caused me a significant headache at first but once I realised the problem it was a relatively easy fix to add them in!
The module only really needed to wrap 5 subroutines out of a total of 73 and access the COMMON block variables. Therefore I generated a raw .pyf file with:
f2py *.f -m LEVMpyFortran -h levmpyfortran_raw.pyf
I then edited this file to remove the subroutines and functions I did not need, added intent(in), intent(out) or intent(in,out) statements, and added depend() statements for array dimensions. For example here is an array x of size nfrei, where nfrei is an input, and x an input and output of the subroutine:
integer intent(in) :: nfreireal*8 dimension(nfrei),depend(nfrei),intent(in,out) :: x
Once done I could compile and link the Fortran with:
f2py -c --f77flags='-Wall -Wno-tabs -O3 -ffast-math' levmpyfortran_edited.pyf *.f
Finally, I wrote a Python class that could read an input data file, set all the flags and COMMON variables, and call the fitting function. The code can be found here, and here are the resulting fits on some of my impedance data for zirconia thin films...
Wrapping Fortran code with more user friendly languages such as Python is therefore a very valuable task with libraries such as Python's scipy and numpy being prime examples. I decided to do this for the complex non-linear least squares fitting program, LEVM, that is designed for fitting impedance spectroscopy data to a vast selection of different circuit models. The code was developed by J Ross Macdonald and is popular in the field, however, it only exists as a DOS command line or basic Windows GUI program. Building a Python library would allow me to integrate it much more easily into my standard data analysis and processing scripts.
The tools to do this are relatively straightforward, in particular I used f2py to compile the Fortran code. However, there were definitely some learning points that I picked up along the way.
First, I had very little experience with Fortran coding (Fortran 77 in this case), so I had to get my head around column formatting of the code (a relic of punch cards), COMMON blocks, endless GOTO statements, and strange subroutine behaviour. But in general, it is an easy language to learn.
The next step was to separate all of the subroutines and functions from the source code and put each into an individual file. This is vital for your sanity for any reasonably sized project! Checking that all of these compiled correctly took some time since incompatible PAUSE statements had to be removed, many variables hadn't been declared explicitly, and column formatting had to be fixed. I also removed several WRITE statements since it is impossible to catch these in Python with f2py, so they only remain as info and warning statements printed to the terminal.
The code had a lot global variables stored in 25 COMMON blocks, which is not ideal by modern coding practices! Although it would be better to have these as function arguments, it was the easiest way to input and output variables without significantly changing the source code. Therefore I made sure a separate Fortran file existed with all the COMMON block variables so they could be accessed from the Python module once f2py is run.
There were also several subroutines and functions that had input variables with assumed size arrays, e.g. DIMENSION N(*). In order to wrap to numpy arrays correctly, I changed those that would be used such that they were of fixed dimension. This usually meant using the global maximum array sizes set in the SIZE.INC file.
Most importantly, there is the issue of SAVE statements. In the original code several subroutines relied on static storage and so SAVE statements were not included (for counter variables for example). This caused me a significant headache at first but once I realised the problem it was a relatively easy fix to add them in!
The module only really needed to wrap 5 subroutines out of a total of 73 and access the COMMON block variables. Therefore I generated a raw .pyf file with:
f2py *.f -m LEVMpyFortran -h levmpyfortran_raw.pyf
I then edited this file to remove the subroutines and functions I did not need, added intent(in), intent(out) or intent(in,out) statements, and added depend() statements for array dimensions. For example here is an array x of size nfrei, where nfrei is an input, and x an input and output of the subroutine:
integer intent(in) :: nfreireal*8 dimension(nfrei),depend(nfrei),intent(in,out) :: x
Once done I could compile and link the Fortran with:
f2py -c --f77flags='-Wall -Wno-tabs -O3 -ffast-math' levmpyfortran_edited.pyf *.f
Finally, I wrote a Python class that could read an input data file, set all the flags and COMMON variables, and call the fitting function. The code can be found here, and here are the resulting fits on some of my impedance data for zirconia thin films...
Comments
Post a Comment