2016-01-03 15 views
5

Ho una subroutine Fortran che vorrei usare in Python.errore f2py con array allocatable

subroutine create_hist(a, n, dr, bins, hist) 
    integer, intent(in) :: n 
    real(8), intent(in) :: a(n) 
    real(8), intent(in) :: dr 
    integer, intent(out), allocatable :: hist(:) 
    real(8), intent(out), allocatable :: bins(:) 

    n_b = n_bins(a, n, dr) ! a function calculating the number of bins 
    allocate(bins(n_b+1)) 
    allocate(hist(n_b)) 
    hist = 0 

    !... Create the histogram hist by putting elems of a into bins 
end subroutine 

Si tratta di un semplice programma per ottenere una matrice di numeri a e creare un istogramma in base alle dimensioni proposta bin dr. Innanzitutto, ottiene il numero di contenitori utilizzando la funzione n_bins e quindi alloca di conseguenza lo spazio per gli array bins e hist.

Mentre gfortran compila questo codice bene, f2py genera un errore:

/tmp/tmpY5badz/src.linux-x86_64-2.6/mat_ops-f2pywrappers2.f90:32.37: 
     call create_hist_gnu(a, n, dr, bins, hist) 
Error: Actual argument for 'bins' must be ALLOCATABLE at (1) 

Da quanto ho capito, f2py non tollera l'allocazione dello spazio per gli array in fase di esecuzione (idea del perché, in quanto questo sembra essere un naturale bisogno scientifico).

Qualcuno potrebbe suggerire un modo per allocare gli array Fortran in fase di esecuzione in modo che f2py sia soddisfatto?

risposta

1

Per quanto ne so f2py non supporta gli argomenti fittizi (talvolta noti come parametri di funzione) che hanno l'attributo ALLOCATABLE. Un'alternativa sarebbe quella di renderlo un array di forme esplicite abbastanza grande. Dovresti quindi utilizzare solo una determinata parte dell'array.

subroutine create_hist(a, n, dr, n_b_max, bins, hist) 
    integer, intent(in) :: n 
    real(8), intent(in) :: a(n) 
    real(8), intent(in) :: dr 
    integer, intent(in) :: n_b_max 
    integer, intent(out) :: hist(n_b_max) 
    real(8), intent(out) :: bins(n_b_max+1) 

Tra l'altro, utilizzando real(8) non è un modo portabile come specificare a 64 bit di precisione e fallirà con alcuni compilatori. Anche il buon vecchio double precision è migliore in quanto verrà sempre compilato a qualcosa e spesso a 64 bit.

+0

Un argomento abbastanza grande a forma assunta funziona anche? E dovrebbe esserci anche un argomento di "dimensione restituita" in modo che il chiamante sappia quanta parte dell'array è valida (dipende dal caso d'uso, ne sono sicuro)? [Chiedere all'ignoranza delle cose più di ogni altra cosa.] – francescalus

+0

Sì, anche la forma assunta funzionerebbe. Hai anche ragione sulle dimensioni restituite. Altrimenti dovrebbe essere dedotto dai valori di 'bin '(potrebbero essere inizializzati su un numero negativo, o qualcosa di simile. –

5

Esistono alcune opzioni valide/round di lavoro che è possibile utilizzare.

1. Probabilmente il più semplice sarebbe quella di provvedere alla pitone chiamare la funzione n_bins, e quindi utilizzare il risultato da quello a chiamare (una versione leggermente modificata) la funzione create_hist con array di uscita non allocatable della taglie corrette

cioè in Python:

n_b = n_bins(a,dr) # don't need to pass n - inferred from a 
bins,hist = create_hist(a,dr,n_b) # bins and hist are auto-allocated 

con l'interfaccia Fortran per create_hist ora definito come

subroutine create_hist(a, n, dr, n_b, bins, hist) 
    integer, intent(in) :: n 
    real(8), intent(in) :: a(n) 
    real(8), intent(in) :: dr 
    integer, intent(in) :: n_b 
    integer, intent(out) :: hist(n_b) 
    real(8), intent(out) :: bins(n_b+1) 

    ! code follows 

! you also need to specify the function n_bins... 

che funziona solo nei casi in cui è possibile chiamare n_bins a buon mercato e dall'esterno create_hist. Conosco 2 metodi per simulare array allocabili per i casi in cui ciò non si applica (ad esempio il codice per calcolare le dimensioni dell'array è costoso e non può essere facilmente separato).

2. Il primo è utilizzare gli array allocabili a livello di modulo (descritti nella documentazione here). Questa è essenzialmente una variabile globale allocabile: si chiama la funzione, si salvano i dati nella variabile globale e quindi si accede a tale da Python. Lo svantaggio è che non è thread-safe (ad es.è male se si chiama create_hist simultaneamente in parallelo)

module something 
    real(8), allocatable :: bins(:) 
    integer, allocatable :: hist(:) 
contains 
    subroutine create_hist(a,n,dr) 
     integer, intent(in) :: n 
     real(8), intent(in) :: a(n) 
     real(8), intent(in) :: dr 
     integer :: n_b 

     n_b = n_bins(a,n,dr) 

     allocate(bins(n_b+1)) 
     allocate(hist(n_b)) 
     ! code follows 
    end subroutine 
end module 

La chiamata Python sembra quindi come

something.create_hist(a,n,dr) 
bins = something.bins # or possible bins = copy.copy(something.bins) 
hist = something.hist # or possible hist = copy.copy(something.hist) 

3. L'altro modo che mi piace molto è quello di allocare solo gli array all'interno della funzione (cioè non passarli dentro/fuori come parametri). Tuttavia, ciò che si passa è una funzione di callback Python che viene chiamata alla fine e salva gli array. Questo è thread-safe (credo).

Il codice FORTRAN quindi sembra qualcosa di simile

subroutine create_hist(a,n,dr,callback) 
    integer, intent(in) :: n 
    real(8), intent(in) :: a(n) 
    real(8), intent(in) :: dr 
    external callable ! note: better to specify the type with an interface block (see http://www.fortran90.org/src/best-practices.html#callbacks) 
    integer :: n_b 
    real(8), allocatable :: bins(:) 
    integer, allocatable :: hist(:) 

    n_b = n_bins(a,n,dr) 

    allocate(bins(n_b+1)) 
    allocate(hist(n_b)) 
    ! code follows 
    call callable(bins,hist,n_b) 
end subroutine 

purtroppo è quindi un po 'più coinvolti. È necessario creare un file di firma con il comando f2py -m fortran_module -h fortran_module.pyf my_fortran_file.f90 (questo è quello di costruire un modulo Python chiamato fortran_module - cambiare i nomi a seconda dei casi), e quindi modificare le linee rilevanti di esso per chiarire le dimensioni della funzione di callback:

python module create_hist__user__routines 
    interface create_hist_user_interface 
     subroutine callable(bins,hist,n_b) ! in :f:my_fortran_file.f90:stuff:create_hist:unknown_interface 
      real(kind=8), dimension(n_b+1) :: bins 
      integer, dimension(n_b) :: hist 
      integer :: n_b 
     end subroutine callable 
    end interface create_hist_user_interface 
end python module create_hist__user__routines 

compilare che con f2py -c fortran_module.pyf my_fortran_file.f90

e poi una breve involucro Python (per darvi una semplice interfaccia) che assomiglia

def create_hist(a,dr): 
    class SaveArraysCallable(object): 
     def __call__(self,bins,hist): 
      self.bins = bins.copy() 
      self.hist = hist.copy() 

    f = SaveArrayCallable() 
    fortran_module.create_hist(a,dr,f) 
    return f.bins, f.hist 

Sommario

Per molti casi l'opzione 1 è probabilmente la migliore. Altrimenti opzione 2 o opzione 3. Preferisco l'opzione 3 perché non ci sono problemi di multithreading da evitare (ma in realtà questo è un caso d'angolo che probabilmente non vedrai mai). L'opzione 2 è più facile da implementare.

+0

C'è qualche documentazione per l'opzione 3? L'hai testata? La dimensione allocatable, (n_b + 1 '' è chiaramente contro le regole di Fortran ma forse f2py lo consente per le sue interfacce di callback.Tuttavia, 'callable' ha ancora un'interfaccia implicita in 'create_hist'. –

+0

Ho testato l'opzione 3 e funziona. È un orrendo trucco della mia invenzione quindi non c'è documentazione. Penso che tu abbia ragione - dovrei eliminare il "allocatable" nelle specifiche (modifico la risposta per farlo) – DavidW

+0

"Tuttavia, callable ha ancora un'interfaccia implicita all'interno di create_hist" - sì, questo è quasi certamente vero. Non conosco abbastanza bene Fortran per sistemarlo, ma accetterei con gratitudine un suggerimento (o modifica) che lo faccia. – DavidW

Problemi correlati