F2py example

From Interactive System for Ice sheet Simulation
Revision as of 09:29, 12 August 2009 by Mauro (Talk | contribs)

(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Jump to: navigation, search

High level languages (Python, R, Matlab, IDL, etc.) are fast to write, powerful and fun. However, to do serious number crunching they are too slow. Low level languages (Fortran, C, etc.) are very fast for computing but slow to write. One way around is to use the best of both worlds: write most your program in a high level language and only write the number crunching bits in a fast low level language. One way of making such a connection is f2py, connecting Python and Fortran. f2py turns a Fortran subroutine into a module which can be imported into Python. The function can then be called within a Python program but will still be executed at the speed of the Fortran code, as it is still Fortran code. I found the documentation for f2py [1] a bit hard to read, especially only having learned Fortran this week and there are not so many examples around ([2], [3], [4]) , thus I put this together.


f2py example

This example is derived from the Kees'_assignment exercise where we developed a flow line model of an alpine glacier in Fortran. I modified the code of the Team_3_Solution such that it can be compiled into a python module using f2py. The source code is here f2py_examples.tar.gz. When going through this example, unpack it (with tar -xvvzf f2py_example.tar.gz), cd into the directory and test with python run_test.py.

Source code

F2py_examples.tar.gz We used some 'advanced' Fortran in that code, if you struggle see [5].

Idea behind f2py and usage

f2py takes a Fortran subroutine and some additional instructions, compiles the Fortran code and builds a module which can then be imported into Python and used there like a normal function. To make this possible a few things need to be told to f2py which are stored either in commented bits in the Fortran code or in the signature file or both. So the steps are:

1. Changes to the code

    1. remove everything to do with allocation in the subroutine you want to turn into the Python module as this doesn't work too well
    2. make the program into a subroutine
    3. add command strings for f2py having the form !f2py ...: example

2. Automatically generate a signature file

f2py -m module_name -h sig_file.pyf list_of_fortran_files

3. Edit the signature file

in this case: emove/comment everything with the submodule. If this is not done, there will be an ImportError in python. This might be due to the dimension(:) statements, which don't always work in f2py.

4. Compile the Python module

f2py -c --fcompiler=gnu95 sig_file.pyf list_of_fortran_files
Note the --fcompiler=gnu95 might be necessary if you compile Fortran 90 files and f2py can't find gfortran

I will refer to above steps in the examples.

Basic example

Found in the directory f2py_example/basic.

This is just very basic and no parameters can be passed to the Fortran function (this will be handled later). Furthermore, the Fortran output is just printed to screen (standard output) and thus in the python script we need to do some tricks to actually get the Fortran program output.

1. Changes to the code

Original Fortran code is is here: original code

2. Remove allocation stuff
  !! unstaggered grid'''1. Changes to the code'''
  real(8), dimension(grid) :: xx         ! node positions
  real(8), dimension(grid) :: thick      ! ice thickness
  real(8), dimension(grid) :: dthick_dt  ! ice thickness time derivative
  real(8), dimension(grid) :: mass_b     ! mass balance
  real(8), dimension(grid) :: bed        ! bed elevation
  !! staggered grid
  real(8), dimension(grid-1) :: diffus     ! diffusion
2. Make the program into a subrotine
subroutine simple_ice_model
end subroutine simple_ice_model

2. Automatically generate a signature file

Run this command (in folder basic)
f2py -m simple_ice_model_basic -h simple_ice_model_basic.pyf ../simple_ice_model_modules.f90 simple_ice_model_basic.f90

3. Edit the signature file

Uncomment all the subroutines used by the actual model: (Note: the source code of the uncommented stuff is still highlighted, so don't get confused)

!    -*- f90 -*-
! Note: the context of this file is case sensitive.
python module simple_ice_model_basic ! in 
    interface  ! in :simple_ice_model
!         module simple_ice_model_modules ! in :simple_ice_model:simple_ice_model_modules.f90
!             function on_stagger(xx) ! in :simple_ice_model:simple_ice_model_modules.f90:simple_ice_model_modules
!                 real(kind=8) dimension(:),intent(in) :: xx
!                 real(kind=8) dimension(size(xx)-1) :: on_stagger
!             end function on_stagger
!             function calc_diffus(thick,bed,dx,const_bn,n_glen) ! in :simple_ice_model:simple_ice_model_modules.f90:simple_ice_model_modules
!                 real(kind=8) dimension(:),intent(in) :: thick
!                 real(kind=8) dimension(:),intent(in) :: bed
!                 real(kind=8) intent(in) :: dx
!                 real(kind=8) intent(in) :: const_bn
!                 real(kind=8) intent(in) :: n_glen
!                 real(kind=8) dimension(size(thick)-1) :: calc_diffus
!             end function calc_diffus
!             function diff_1(yy,dx) ! in :simple_ice_model:simple_ice_model_modules.f90:simple_ice_model_modules
!                 real(kind=8) dimension(:),intent(in) :: yy
!                 real(kind=8) intent(in) :: dx
!                 real(kind=8) dimension(size(yy)-1) :: diff_1
!             end function diff_1
!             function time_step(yy,dyy_dt,dt) ! in :simple_ice_model:simple_ice_model_modules.f90:simple_ice_model_modules
!                 real(kind=8) dimension(:),intent(in) :: yy
!                 real(kind=8) dimension(:),intent(in) :: dyy_dt
!                 real(kind=8) intent(in) :: dt
!                 real(kind=8) dimension(size(yy)) :: time_step
!             end function time_step
!             function pad(on_stag,left,right) ! in :simple_ice_model:simple_ice_model_modules.f90:simple_ice_model_modules
!                 real(kind=8) dimension(:),intent(in) :: on_stag
!                 real(kind=8) intent(in) :: left
!                 real(kind=8) intent(in) :: right
!                 real(kind=8) dimension(size(on_stag)+2) :: pad
!             end function pad
!         end module simple_ice_model_modules
        subroutine simple_ice_model ! in :simple_ice_model:simple_ice_model.f90
            use simple_ice_model_modules
        end subroutine simple_ice_model
    end interface 
end python module simple_ice_model_basic
! This file was auto-generated with f2py (version:2_4422).
! See http://cens.ioc.ee/projects/f2py2e/

4. Make the Python-module

f2py -c --fcompiler=gnu95 simple_ice_model_basic.pyf ../simple_ice_model_modules.f90 simple_ice_model_basic.f90

Use it in python

The complied module can now be used in python. Import it with import simple_ice_model_basic and the function can then be run with simple_ice_model_basic.simple_ice_model().

Here is an example script on how to run the model.

import numpy as N
import pylab as P
import sys
import os
import simple_ice_model_basic
import py_fortran_tools
def run():
    # as the Frotran modul just prints to stdout (the screen) we need to do some
    # magic to capture that output:
    # from http://stackoverflow.com/questions/977840/redirecting-fortran-called-via-f2py-output-in-python
    output_file = 'out.tmp.dat'
    # open outputfile
    outfile = os.open(output_file, os.O_RDWR|os.O_CREAT)
    # save the current file descriptor
    save = os.dup(1)
    # put outfile on 1
    os.dup2(outfile, 1)
    # end magic
    # model paramaters
    grid = 51
    dt = 0.1
    t_final = 1000
    # run the program
    # restore the standard output file descriptor
    os.dup2(save, 1)
    # close the output file
    # read it
    dict_ = py_fortran_tools.parse_output(output_file)
    # delete it
    # plot
    sol_ind = -1 # index of solution to plot
    P.plot(dict_['xx']/1e3, dict_['bed']+dict_['thick'], 'b')
    P.plot(dict_['xx']/1e3, dict_['bed'], 'r')
    P.xlabel('x (km)')
    P.ylabel('b+H (m)')
if __name__=='__main__':

Improving it a bit: intermediate

Found in the directory f2py_example/intermediate.

We want to be able to pass arguments to the subroutine and get the output returned:

  • Input: grid, dt, t_final
  • Output: thick

1. Changes to the code

subroutine simple_ice_model(grid, dt, t_final, thick)
end subroutine simple_ice_model
Note: both the input and output variables need to be passed

Add f2py commands to the Fortran file:

!f2py integer intent(in) :: grid
!f2py real(8) intent(in) :: dt, t_final
!f2py real(8) intent(out) :: thick, bed, xx

2. Generate a signature file (.pyf)

f2py -m simple_ice_model_intermediate -h simple_ice_model_intermediate.pyf ../simple_ice_model_modules.f90 simple_ice_model_intermediate.f90

3. Edit the signature file

Remove/comment everything with the submodule, like done in the above example. If this is not done, there will be an ImportError in python.

4. Make the Python-module

f2py -c --fcompiler=gnu95 simple_ice_model_intermediate.pyf ../simple_ice_model_modules.f90 simple_ice_model_intermediate.f90

Run it

ipython run_intermediate.py

More advanced example

If we want to run the Fortran module from python, we want to have control over all the model parameters, so we include some more function parameters. And as an example of a vector being passed into Fortran, the boundary conditions are changed from two separate floats to a vector of two floats.

1. Changes to the code

subroutine simple_ice_model(n_t_out, grid, dt, t_final, xl, xr, rho, g_grav, A_rate, &
   n_glen, dbed_dx, mass_b_0, mass_b_1, b_cond, &
   thick_out, bed, xx, times_out) ! output
  use simple_ice_model_modules
  implicit none
  ! local variables
  real(8), parameter    :: pi = 3.1415926536
  !! input xxxxxxxxxx
  integer :: grid, n_t_out                    ! Number of nodes, number of output times
  real(8) :: dt, t_final             ! time step, stopping time
  real(8) :: dbed_dx, mass_b_0, mass_b_1, xl, xr ! model parameters, bed slope, mass balanace params 2x, left and right extent of domain
  real(8) :: rho, g_grav, A_rate, n_glen ! physical constants: dens. of ice, ice rate factor, n
  !! BC
  real(8), dimension(2) :: b_cond          ! thickness at left and right edge
end subroutine simple_ice_model

Note: both the input and output variables need to be passed into the Fortran subroutine.

Add f2py commands to the Fortran file:

!! f2py declarations
!f2py integer optional, intent(in) :: n_t_out=30, grid=51
!f2py real(8) optional, intent(in) :: dt=0.1, t_final=1000., xl=0., xr=50.e3
!f2py real(8) optional, intent(in) :: dbed_dx=0.025, mass_b_0=4., mass_b_1=0.2e-3 
!f2py real(8) optional, intent(in) :: rho=920., g_grav=9.8, A_rate=1.e-16, n_glen=3.
!f2py real(8) optional, dimension(2), intent(in) :: b_cond=(0.,0.)
!f2py real(8) optional, intent(out) :: thick_out, bed, xx, times_out

2. Generate a signature file (.pyf)

f2py -m simple_ice_model_advanced -h simple_ice_model_advanced.pyf ../simple_ice_model_modules.f90 simple_ice_model_advanced.f90

3. Edit the signature file

Remove/comment everything with the submodule, like done in the above example. If this is not done, there will be an ImportError in python.

4. Make the Python-module

f2py -c --fcompiler=gnu95 simple_ice_model_advanced.pyf ../simple_ice_model_modules.f90 simple_ice_model_advanced.f90

Run it

ipython run_intermediate.py

Vector in and output

Here just a generic example of how to do vector input and output. Generally f2py can handle in and output of the size of the vectors/matrices automatically, if you tell it what to do. Find this in the folder vector_input_output.

Fortran code

subroutine vec_in_out(in_matrix, in_out_vec, n ,m, p, out_vector)
  implicit none
  integer :: n,m,p
  real(8) :: in_matrix(n,m)
  real(8) :: out_vector(n+m-1)
  real(8) :: in_out_vec(p)
!f2py real(8), intent(out), dimension(n+m-1) :: out_matix
!f2py real(8), intent(in), dimension(n,m) :: in_matrix
!f2py real(8), intent(inout), dimension(p) :: in_out_vec
!f2py integer, intent(in) :: n,m, p
  out_vector = (/ in_matrix(:,1), in_matrix(:,2) /)
  in_out_vec = in_out_vec*5
end subroutine vec_in_out

Signature file

The autogenerated signature file (python make_signature_file.py) has a bug and the intent(out) needs to be put there manually. Even though in the source code the line !f2py real(8), intent(out), dimension(n+m-1) :: out_matix is included.

After the edit it looks like this:

!    -*- f90 -*-
! Note: the context of this file is case sensitive.
python module vector_in_out ! in 
    interface  ! in :vector_in_out
        subroutine vec_in_out(in_matrix,in_out_vec,n,m,p,out_vector) ! in :vector_in_out:vector_in_out.f90
            real(kind=8) dimension(n,m),intent(in) :: in_matrix
            real(kind=8) dimension(p),intent(inout,out) :: in_out_vec
            integer optional,intent(in),check(shape(in_matrix,0)==n),depend(in_matrix) :: n=shape(in_matrix,0)
            integer optional,intent(in),check(shape(in_matrix,1)==m),depend(in_matrix) :: m=shape(in_matrix,1)
            integer optional,intent(in),check(len(in_out_vec)>=p),depend(in_out_vec) :: p=len(in_out_vec)
            real(kind=8) intent(out), dimension(n+m-1),depend(n,m) :: out_vector
        end subroutine vec_in_out
    end interface 
end python module vector_in_out
! This file was auto-generated with f2py (version:2_4422).
! See http://cens.ioc.ee/projects/f2py2e/


Make the module, similarly to above (or just do python make_module.py) and run it with python run.py.