F2Py mit zuweisbaren und angenommenen Formarrays

18

Ich möchte f2pymit modernen Fortran verwenden. Insbesondere versuche ich, das folgende grundlegende Beispiel zum Laufen zu bringen. Dies ist das kleinste nützliche Beispiel, das ich generieren konnte.

! 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

Beachten Sie, dass dies naus der Form des Eingabeparameters abgeleitet wird x. Beachten Sie, dass dies yim Hauptteil der Unterroutine zugewiesen und freigegeben wird.

Wenn ich das mit kompiliere f2py

f2py -c alloc_test.f90 -m alloc

Und dann in Python ausführen

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

Ich erhalte den folgenden Fehler

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

Also erstelle und bearbeite ich die pyfDatei manuell

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

Geändert

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

Jetzt läuft es aber die Werte der Ausgabe zsind immer 0. Bei einigen Debug-Druckvorgängen wird angezeigt, dass nder Wert 0innerhalb der Unterroutine liegt f. Ich gehe davon aus, dass mir etwas Kopfzauber fehlt f2py, um diese Situation richtig zu handhaben.

Allgemeiner gesagt, was ist der beste Weg, um die obige Subroutine mit Python zu verknüpfen? Ich würde es nachdrücklich vorziehen, das Unterprogramm selbst nicht ändern zu müssen.

MRocklin
quelle
Matt, kennen Sie die Best Practices von Ondrej Certik, insbesondere den Abschnitt Interfacing with Python ? Wir haben ein ähnliches Schnittstellenproblem für PyClaw besprochen und es auch noch nicht gelöst :)
Aron Ahmadia

Antworten:

23

Ich bin nicht besonders vertraut mit f2py-Interna, aber ich bin sehr vertraut damit, Fortran zu verpacken. F2py automatisiert nur einige oder alle der folgenden Dinge.

  1. Sie müssen zuerst mit dem Modul iso_c_binding nach C exportieren, wie hier beschrieben:

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

    Haftungsausschluss: Ich bin der Hauptautor der fortran90.org-Seiten. Dies ist die einzige plattform- und compilerunabhängige Methode, um Fortran von C aus aufzurufen. Dies ist F2003. Heutzutage gibt es keinen Grund, eine andere Methode zu verwenden.

  2. Sie können nur Arrays mit der angegebenen vollständigen Länge (explizite Form) exportieren / aufrufen, dh:

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

    aber keine Form annehmen:

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

    Dies liegt daran, dass die C-Sprache solche Arrays selbst nicht unterstützt. Es ist die Rede, eine solche Unterstützung in F2008 oder höher aufzunehmen (ich bin nicht sicher), und die Art und Weise, wie dies funktionieren würde, ist durch einige unterstützende C-Datenstrukturen, da Sie Forminformationen über das Array tragen müssen.

    In Fortran sollten Sie hauptsächlich die Form annehmen, nur in besonderen Fällen sollten Sie die explizite Form verwenden, wie hier beschrieben:

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

    Das bedeutet, dass Sie einen einfachen Wrapper um Ihre Subroutine "Form annehmen" schreiben müssen, der die Dinge in explizite Form-Arrays einschließt, wie in meinem ersten Link oben angegeben.

  3. Sobald Sie eine C-Signatur haben, rufen Sie sie einfach in Python auf, wie Sie möchten. Ich verwende Cython, aber Sie können ctype oder C / API von Hand verwenden.

  4. Die deallocate(y)Anweisung wird nicht benötigt, Fortran wird automatisch freigegeben.

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

  5. real*8sollte nicht verwendet werden, sondern real(dp):

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

  6. Die Anweisung y(:) = 1.0weist 1,0 mit einfacher Genauigkeit zu, der Rest der Ziffern ist also zufällig! Dies ist eine häufige Falle:

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

    Sie müssen verwenden y(:) = 1.0_dp.

  7. Anstatt zu schreiben y(:) = 1.0_dp, können Sie einfach schreiben y = 1, das war's. Sie können einer Gleitkommazahl eine Ganzzahl zuweisen, ohne an Genauigkeit zu verlieren, und Sie müssen die Redundanz nicht (:)dort einfügen . Viel einfacher

  8. Anstatt

    y = 1
    z = x + y

    benutz einfach

    z = x + 1

    und kümmere dich überhaupt nicht um das yArray.

  9. Sie brauchen die "return" -Anweisung am Ende des Unterprogramms nicht.

  10. Schließlich sollten Sie wahrscheinlich Module verwenden und nur implicit nonedie Modulebene aufrufen, und Sie müssen sie nicht in jeder Unterroutine wiederholen.

    Ansonsten sieht es für mich gut aus. Hier ist der Code, der den obigen Vorschlägen 1-10 folgt:

    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

    Es zeigt das vereinfachte Unterprogramm sowie einen C-Wrapper.

    In Bezug auf f2py wird wahrscheinlich versucht, diesen Wrapper für Sie zu schreiben, und dies schlägt fehl. Ich bin mir auch nicht sicher, ob es das iso_c_bindingModul verwendet. Aus all diesen Gründen bevorzuge ich es, Dinge von Hand zu wickeln. Dann ist genau klar, was passiert.

Ondřej Čertík
quelle
Soweit ich weiß, stützt sich f2py nicht auf ISO-C-Bindungen (sein Hauptziel ist Fortran-77- und Fortran-90-Code).
Aron Ahmadia
Ich wusste, dass ich ein bisschen dumm war, yaber ich wollte, dass etwas zugewiesen wurde (mein aktueller Code weist nicht-triviale Zuweisungen auf). Über viele der anderen Punkte wusste ich allerdings nichts. Es sieht so aus, als sollte ich mir eine Art Fortran90-Leitfaden für bewährte Methoden ansehen. Vielen Dank für die gründliche Antwort!
MRocklin
Beachten Sie, dass Sie mit den heutigen Fortran-Compilern F77 genauso umbrechen - indem Sie einen einfachen iso_c_binding-Wrapper schreiben und die alte f77-Subroutine aufrufen.
Ondřej Čertík
6

Sie müssen lediglich Folgendes tun:

!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

Obwohl die Größe von Array x und z jetzt als explizites Argument übergeben wird, macht f2py das Argument n optional. Im Folgenden sehen Sie die Dokumentzeichenfolge der Funktion, wie sie für Python angezeigt wird:

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)

Importieren und Aufrufen von Python:

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

gibt die folgende Ausgabe aus:

[ 2.  2.  2.  2.  2.]
Prometheus
quelle
Gibt es eine Möglichkeit, einen nicht trivialen Ausdruck als Größe zu verwenden? Zum Beispiel übergebe ich nund möchte ein Array von Größe erhalten 2 ** n. Bisher muss ich auch 2 ** n als separates Argument übergeben.
Alleo