F2Py avec des tableaux de formes attribuables et supposés

18

Je voudrais utiliser f2pyavec le Fortran moderne. En particulier, j'essaie de faire fonctionner l'exemple de base suivant. Ceci est le plus petit exemple utile que j'ai pu générer.

! alloc_test.f90
subroutine f(x, z)
  implicit none

! Argument Declarations !
  real*8, intent(in) ::  x(:)
  real*8, intent(out) :: z(:)

! Variable Declarations !
  real*8, allocatable :: y(:)
  integer :: n

! Variable Initializations !
  n = size(x)
  allocate(y(n))

! Statements !
  y(:) = 1.0
  z = x + y

  deallocate(y)
  return
end subroutine f

Notez que cela nest déduit de la forme du paramètre d'entrée x. Notez que yest alloué et désalloué dans le corps du sous-programme.

Quand je compile cela avec f2py

f2py -c alloc_test.f90 -m alloc

Et puis exécutez en Python

from alloc import f
from numpy import ones
x = ones(5)
print f(x)

J'obtiens l'erreur suivante

ValueError: failed to create intent(cache|hide)|optional array-- must have defined dimensions but got (-1,)

Je vais donc créer et éditer le pyffichier manuellement

f2py -h alloc_test.pyf -m alloc alloc_test.f90

Original

python module alloc ! in 
    interface  ! in :alloc
        subroutine f(x,z) ! in :alloc:alloc_test.f90
            real*8 dimension(:),intent(in) :: x
            real*8 dimension(:),intent(out) :: z
        end subroutine f
    end interface 
end python module alloc

Modifié

python module alloc ! in 
    interface  ! in :alloc
        subroutine f(x,z,n) ! in :alloc:alloc_test.f90
            integer, intent(in) :: n
            real*8 dimension(n),intent(in) :: x
            real*8 dimension(n),intent(out) :: z
        end subroutine f
    end interface 
end python module alloc

Maintenant, il s'exécute mais les valeurs de la sortie zsont toujours 0. Une impression de débogage révèle que cela na de la valeur 0dans le sous-programme f. Je suppose qu'il me manque une f2pymagie d'en-tête pour gérer correctement cette situation.

Plus généralement, quelle est la meilleure façon de lier le sous-programme ci-dessus à Python? Je préférerais fortement ne pas avoir à modifier le sous-programme lui-même.

MRocklin
la source
Matt, connaissez-vous le guide des meilleures pratiques d'Ondrej Certik, en particulier la section Interfaçage avec Python ? Nous avons discuté d'un problème d'interface similaire pour PyClaw et ne l'avons pas encore résolu à ce stade non plus :)
Aron Ahmadia

Réponses:

23

Je ne suis pas très familier avec les composants internes de f2py, mais je suis très familier avec l'emballage de Fortran. F2py automatise simplement tout ou partie des éléments ci-dessous.

  1. Vous devez d'abord exporter vers C à l'aide du module iso_c_binding, comme décrit par exemple ici:

    http://fortran90.org/src/best-practices.html#interfacing-with-c

    Avertissement: je suis l'auteur principal des pages fortran90.org. C'est le seul moyen indépendant de la plate-forme et du compilateur d'appeler Fortran à partir de C. Il s'agit de F2003, donc ces jours-ci, il n'y a aucune raison d'utiliser un autre moyen.

  2. Vous pouvez uniquement exporter / appeler des tableaux avec toute la longueur spécifiée (forme explicite), c'est-à-dire:

    integer(c_int), intent(in) :: N
    real(c_double), intent(out) :: mesh(N)

    mais pas prendre la forme:

    real(c_double), intent(out) :: mesh(:)

    En effet, le langage C ne prend pas lui-même en charge ces tableaux. Il est question d'inclure une telle prise en charge dans F2008 ou version ultérieure (je ne suis pas sûr), et la façon dont cela fonctionnerait consiste à utiliser certaines structures de données C, car vous devez transporter des informations de forme sur le tableau.

    Dans Fortran, vous devez principalement utiliser la forme supposée, uniquement dans des cas particuliers, vous devez utiliser une forme explicite, comme décrit ici:

    http://fortran90.org/src/best-practices.html#arrays

    Cela signifie que vous devez écrire un simple wrapper autour de votre sous-programme de forme supposée, qui enveloppera les choses dans des tableaux de formes explicites, selon mon premier lien ci-dessus.

  3. Une fois que vous avez une signature C, appelez-la simplement à partir de Python comme vous le souhaitez, j'utilise Cython, mais vous pouvez utiliser ctype ou C / API à la main.

  4. La deallocate(y)déclaration n'est pas nécessaire, Fortran se désalloue automatiquement.

    http://fortran90.org/src/best-practices.html#allocatable-arrays

  5. real*8ne doit pas être utilisé, mais plutôt real(dp):

    http://fortran90.org/src/best-practices.html#floating-point-numbers

  6. L'instruction y(:) = 1.0attribue 1.0 en simple précision, donc le reste des chiffres sera aléatoire! Ceci est un piège courant:

    http://fortran90.org/src/gotchas.html#floating-point-numbers

    Vous devez utiliser y(:) = 1.0_dp.

  7. Au lieu d'écrire y(:) = 1.0_dp, vous pouvez simplement écrire y = 1, c'est tout. Vous pouvez attribuer un entier à un nombre à virgule flottante sans perdre en précision, et vous n'avez pas besoin de mettre le redondant (:)à l'intérieur. Beaucoup plus simple.

  8. Au lieu de

    y = 1
    z = x + y

    il suffit d'utiliser

    z = x + 1

    et ne vous embêtez pas du tout avec le ytableau.

  9. Vous n'avez pas besoin de l'instruction "return" à la fin du sous-programme.

  10. Enfin, vous devriez probablement utiliser des modules, et simplement mettre implicit noneau niveau du module et vous n'avez pas besoin de le répéter dans chaque sous-programme.

    Sinon, ça me semble bien. Voici le code suivant les suggestions 1-10 ci-dessus:

    module test
    use iso_c_binding, only: c_double, c_int
    implicit none
    integer, parameter :: dp=kind(0.d0)
    
    contains
    
    subroutine f(x, z)
    real(dp), intent(in) ::  x(:)
    real(dp), intent(out) :: z(:)
    z = x + 1
    end subroutine
    
    subroutine c_f(n, x, z) bind(c)
    integer(c_int), intent(in) :: n
    real(c_double), intent(in) ::  x(n)
    real(c_double), intent(out) :: z(n)
    call f(x, z)
    end subroutine
    
    end module

    Il montre le sous-programme simplifié ainsi qu'un wrapper C.

    En ce qui concerne f2py, il essaie probablement d'écrire ce wrapper pour vous et échoue. Je ne sais pas non plus s'il utilise le iso_c_bindingmodule. Donc pour toutes ces raisons, je préfère envelopper les choses à la main. Il est alors parfaitement clair ce qui se passe.

Ondřej Čertík
la source
Pour autant que je sache, f2py ne s'appuie pas sur les liaisons ISO C (sa cible principale est le code Fortran 77 et Fortran 90).
Aron Ahmadia
Je savais que j'étais un peu stupide, ymais je voulais que quelque chose soit alloué (mon code réel a des allocations non triviales). Je ne connaissais cependant pas beaucoup d'autres points. On dirait que je devrais aller chercher une sorte de guide des meilleures pratiques Fortran90 .... Merci pour la réponse complète!
MRocklin
Notez qu'en utilisant les compilateurs Fortran d'aujourd'hui, vous encapsulez F77 exactement de la même manière --- en écrivant un simple wrapper iso_c_binding et en appelant le sous-programme hérité f77.
Ondřej Čertík
6

Tout ce que vous avez à faire est le suivant:

!alloc_test.f90
subroutine f(x, z, n)
  implicit none

! Argument Declarations !
  integer :: n
  real*8, intent(in) ::  x(n)
  real*8, intent(out) :: z(n)

! Variable Declarations !
  real*8, allocatable :: y(:)

! Variable Initializations !
  allocate(y(n))

! Statements !
  y(:) = 1.0
  z = x + y

  deallocate(y)
  return
end subroutine f

Bien que la taille des tableaux x et z soit désormais passée comme argument explicite, f2py rend l'argument n facultatif. Voici la docstring de la fonction telle qu'elle apparaît en python:

Type:       fortran
String Form:<fortran object>
Docstring:
f - Function signature:
  z = f(x,[n])
Required arguments:
  x : input rank-1 array('d') with bounds (n)
Optional arguments:
  n := len(x) input int
Return objects:
  z : rank-1 array('d') with bounds (n)

Importation et appel depuis python:

from alloc import f
from numpy import ones
x = ones(5)
print f(x)

donne la sortie suivante:

[ 2.  2.  2.  2.  2.]
Prométhée
la source
Existe-t-il un moyen d'utiliser une expression non triviale comme taille? Par exemple, je passe net je veux obtenir un tableau de taille 2 ** n. Jusqu'à présent, je dois également passer 2 ** n comme argument séparé.
Alleo