Leistungsvergleich Fortran, Numpy, Cython und numexpr

Question

Ich habe folgende Funktion:

def get_denom(n_comp,qs,x,cp,cs): 
''' 
len(n_comp) = 1 # number of proteins 
len(cp) = n_comp # protein concentration 
len(qp) = n_comp # protein capacity 
len(x) = 3*n_comp + 1 # fit parameters 
len(cs) = 1 

''' 
    k = x[0:n_comp] 
    sigma = x[n_comp:2*n_comp] 
    z = x[2*n_comp:3*n_comp] 

    a = (sigma + z)*(k*(qs/cs)**(z-1))*cp 
    denom = np.sum(a) + cs 
    return denom 

Ich vergleiche es gegen eine Fortran-Implementierung (Meine erste Fortran-Funktion überhaupt):

subroutine get_denom (qs,x,cp,cs,n_comp,denom) 

! Calculates the denominator in the SMA model (Brooks and Cramer 1992) 
! The function is called at a specific salt concentration and isotherm point 
! I loops over the number of components 

implicit none 

! declaration of input variables 
integer, intent(in) :: n_comp ! number of components 
double precision, intent(in) :: cs,qs ! salt concentration, free ligand concentration 
double precision, dimension(n_comp), INTENT(IN) ::cp ! protein concentration 
double precision, dimension(3*n_comp + 1), INTENT(IN) :: x ! parameters 

! declaration of local variables 
double precision, dimension(n_comp) :: k,sigma,z 
double precision :: a 
integer :: i 

! declaration of outpur variables 
double precision, intent(out) :: denom 

k = x(1:n_comp) ! equlibrium constant 
sigma = x(n_comp+1:2*n_comp) ! steric hindrance factor 
z = x(2*n_comp+1:3*n_comp) ! charge of protein 

a = 0. 
do i = 1,n_comp 
    a = a + (sigma(i) + z(i))*(k(i)*(qs/cs)**(z(i)-1.))*cp(i) 
end do 

denom = a + cs 

end subroutine get_denom 

ich die kompilierte. F95-Datei durch die Verwendung:

1) f2py -c -m get_denom get_denom.f95 --fcompiler=gfortran

2) f2py -c -m get_denom_vec get_denom.f95 --fcompiler=gfortran --f90flags='-msse2' (Die letzte Option auf Auto-Vektorisierung drehen sollte)

ich die Funktionen testen:

import numpy as np 
import get_denom as fort_denom 
import get_denom_vec as fort_denom_vec 
from matplotlib import pyplot as plt 
%matplotlib inline 

def get_denom(n_comp,qs,x,cp,cs): 
    k = x[0:n_comp] 
    sigma = x[n_comp:2*n_comp] 
    z = x[2*n_comp:3*n_comp] 
    # calculates the denominator in Equ 14a - 14c (Brooks & Cramer 1992) 
    a = (sigma + z)*(k*(qs/cs)**(z-1))*cp 
    denom = np.sum(a) + cs 
    return denom 

n_comp = 100 
cp = np.tile(1.243,n_comp) 
cs = 100. 
qs = np.tile(1100.,n_comp) 
x= np.random.rand(3*n_comp+1) 
denom = np.empty(1) 
%timeit get_denom(n_comp,qs,x,cp,cs) 
%timeit fort_denom.get_denom(qs,x,cp,cs,n_comp) 
%timeit fort_denom_vec.get_denom(qs,x,cp,cs,n_comp) 

Ich habe folgende Cython Code:

import cython 
# import both numpy and the Cython declarations for numpy 
import numpy as np 
cimport numpy as np 

@cython.boundscheck(False) 
@cython.wraparound(False) 
def get_denom(int n_comp,np.ndarray[double, ndim=1, mode='c'] qs, np.ndarray[double, ndim=1, mode='c'] x,np.ndarray[double, ndim=1, mode='c'] cp, double cs): 

    cdef int i 
    cdef double a 
    cdef double denom 
    cdef double[:] k = x[0:n_comp] 
    cdef double[:] sigma = x[n_comp:2*n_comp] 
    cdef double[:] z = x[2*n_comp:3*n_comp] 
    # calculates the denominator in Equ 14a - 14c (Brooks & Cramer 1992) 
    a = 0. 
    for i in range(n_comp): 
    #a += (sigma[i] + z[i])*(pow(k[i]*(qs[i]/cs), (z[i]-1)))*cp[i] 
     a += (sigma[i] + z[i])*(k[i]*(qs[i]/cs)**(z[i]-1))*cp[i] 

    denom = a + cs 

    return denom 

EDIT:

Added numexpr mit einem Gewinde:

def get_denom_numexp(n_comp,qs,x,cp,cs): 
    k = x[0:n_comp] 
    sigma = x[n_comp:2*n_comp] 
    z = x[2*n_comp:3*n_comp] 
    # calculates the denominator in Equ 14a - 14c (Brooks & Cramer 1992) 
    a = ne.evaluate('(sigma + z)*(k*(qs/cs)**(z-1))*cp') 
    return cs + np.sum(a) 

ne.set_num_threads(1) # using just 1 thread 
%timeit get_denom_numexp(n_comp,qs,x,cp,cs) 

Das Ergebnis ist (kleiner ist besser):

Warum ist die Geschwindigkeit des Fortran näher an Numpy bekommen Größe der Arrays mit zunehmender? Und wie könnte ich Cython beschleunigen? Zeiger verwenden?

Answer 1

Sussed es.

OK, endlich, wir durften Numpy etc auf einer unserer Boxen installieren, und das hat erlaubt, was eine umfassende Erklärung Ihres ursprünglichen Posts sein kann.

Die kurze Antwort ist, dass Ihre ursprünglichen Fragen in gewisser Weise "die falsche Frage" sind. Darüber hinaus gab es viele ärgerliche Beschimpfungen und Fehlinformationen durch einen der Mitwirkenden, und diese Fehler und Fälschungen verdienen Aufmerksamkeit, damit niemand den Fehler macht, ihnen zu glauben, und zu ihren Kosten.

Auch ich habe beschlossen, diese Antwort als separate Antwort, anstatt meine Antwort vom 14. April aus Gründen, siehe unten, und Anstand zu bearbeiten.

Teil A: Die Antwort auf die OP

Das Wichtigste zuerst, sich mit der Frage in der ursprünglichen Nachricht zu tun: Sie erinnern sich vielleicht, ich nur WRT auf die Fortran Seite kommentieren könnte, da unsere Richtlinien streng sind über das, was Software kann installiert werden und wo auf unseren Maschinen, und wir hatten Python usw. nicht zur Hand (bis jetzt). Ich habe auch wiederholt gesagt, dass der Charakter Ihres Ergebnisses interessant ist in Bezug auf das, was wir seinen gekrümmten Charakter oder vielleicht "Konkavität" nennen können.

Darüber hinaus arbeitet rein mit „relative“ Ergebnisse (wie man nicht die absoluten Timings haben Post, und ich habe nicht Numpy zu der Zeit zur Hand haben), hatte ich ein paar Mal darauf hingewiesen, dass einige wichtige Informationen sein kann lauern darin.

Das ist genau der Fall.

Zuerst wollte ich sicher sein, dass ich Ihre Ergebnisse reproduzieren konnte, da wir Python/F2py normalerweise nicht verwenden, es war nicht offensichtlich, welche Compilereinstellungen usw. in Ihren Ergebnissen enthalten sind, also führte ich eine Vielzahl von Tests durch seien Sie sicher, dass wir Äpfel-zu-Äpfel sprechen (meine Antwort vom 14. April zeigte, dass Debug vs Release/O2 einen großen Unterschied macht).

Abbildung 1 zeigt meine Python-Ergebnisse für nur die drei Fälle von: Ihrem ursprünglichen internen Python/Numpy-Unterprogramm (nennen Sie dieses P, ich habe gerade Ihr Original ausgeschnitten/eingefügt), Ihr ursprüngliches Do based Fortran s/r gepostet (rufen Sie diesen FDo auf, ich kopiere einfach das Original wie zuvor), und eine der Varianten, die ich zuvor vorgeschlagen hatte, auf Array-Abschnitte angewiesen ist und Sum() benötigt (nennen Sie das FSAS, erstellt durch Bearbeiten Ihres ursprünglichen FDO) . Abbildung 1 zeigt die absoluten Timings über Zeit.

Abbildung 2 zeigt die relativen Ergebnisse basierend auf Ihrer relativen Strategie der Division durch Python/Numpy (P) -Timing. Nur die zwei (relativen) Fortran-Varianten werden gezeigt.

Natürlich reproduzieren diese das Zeichen in Ihrer ursprünglichen Handlung, und wir können darauf vertrauen, dass meine Tests mit Ihren Tests übereinstimmen.

Nun war Ihre ursprüngliche Frage "Warum ist die Geschwindigkeit von Fortran, die mit zunehmender Größe der Arrays Numpy näher kommt?".

In der Tat ist es nicht. Es ist rein ein Artefakt/eine Verzerrung, sich rein auf "relative" Timings zu verlassen, die diesen Eindruck vermitteln können.

Betrachtet man Abbildung 1, mit den absoluten Timings, ist klar, dass die Numpy und Fortran divergieren. In der Tat bewegen sich die Fortran-Ergebnisse von Numpy oder umgekehrt, wenn Sie möchten.

Eine bessere Frage, und eine, die wiederholt in meinen früheren Kommentaren auftauchte, ist, warum sind diese aufwärts kurven in erster Linie (zB linear, zum Beispiel)? Meine bisherigen Fortran-Ergebnisse zeigten ein "überwiegend" flaches relatives Performance-Verhältnis (gelbe Linien in meinem Apr 14 Chart/Antwort), und so hatte ich spekuliert, dass auf der Python-Seite etwas Interessantes passierte, und schlug vor, dies zu überprüfen.

Eine Möglichkeit, dies zu zeigen, ist eine andere Art von relativen Maßen. Ich habe jede (absolute) Reihe mit ihrem eigenen höchsten Wert (dh bei n_comp = 10k) geteilt, um zu sehen, wie sich diese "interne relative" Leistung entfaltet (diese werden als die 10k-Werte bezeichnet, die die Zeiten für n_comp = 10.000 darstellen) . Abbildung 3 zeigt diese Ergebnisse für P, FDo und FSAS als P/P10k, FDo/FDo10k und FSAS/FSAS10k. Zur Klarheit hat die y-Achse eine logarithmische Skala. Es ist klar, dass die Fortran-Varianten mit abnehmendem n_comp.cf relativ viel besser vorformuliert sind. die P-Ergebnisse (z. B. der mit dem roten Kreis versehene Abschnitt).

Mit anderen Worten, Fortran ist sehr sehr (nicht linear) effizienter für die Verringerung der Array-Größe. Nicht sicher, warum Python mit sinkendem n_comp so viel schlimmer machen würde ... aber da ist es und kann ein Problem mit internem Overhead/Setup usw. sein, und die Unterschiede zwischen Interpretern und Compilern usw.

Es ist also nicht so, dass Fortran mit Python "aufholt", ganz im Gegenteil, es distanziert sich weiterhin von Python (siehe Abbildung 1). Die Unterschiede in den Nichtlinearitäten zwischen Python und Fortran mit abnehmendem n_comp erzeugen jedoch "relative" Timings mit scheinbar kontraintuitivem und nicht-linearem Charakter.

Somit kann, wie n_comp zunimmt und jede Methode „stabilisiert“ zu einem mehr oder weniger linearen Modus, die Kurven ihre Differenzen linear zu zeigen, abflachen, die wachsen, und die relativen Verhältnisse absetzen zu einem ungefähren konstanten (memory Conten ignorierend, smp Probleme, usw.) ... das ist leichter zu sehen, wenn n_comp> 10k erlaubt ist, aber die gelbe Linie in meiner Antwort vom 14. April zeigt dies bereits für die Fortran-only s/r's.

Beiseite: Meine Präferenz ist es, meine eigenen Timing-Routinen/Funktionen zu erstellen. Zeit scheint bequem, aber in dieser "Black Box" ist viel los. Das Setzen eigener Schleifen und Strukturen und das Sicherstellen der Eigenschaften/Auflösung Ihrer Timing-Funktionen ist wichtig für eine korrekte Bewertung.

Answer 2

Es gibt nicht genügend Informationen im Anhang, aber einige der folgenden Punkte können helfen:

1) Fortran hat intrinsische Funktionen optimiert wie „Sum()“ und „Dot_Product“, die Sie gerne ausleihen können anstelle der Do-Schleife für Summierungen usw.

In einigen Fällen (nicht unbedingt hier), kann es besser sein, ForAll oder was auch immer zu verwenden, um "Meta" -Arrays zu summieren und dann die Summierung anzuwenden auf den "Meta" -Arrays.

Fortran erlaubt jedoch Array-Abschnitte, sodass Sie die automatischen/intermediären Arrays sigma, k und z und den zugehörigen Overhead nicht erstellen müssen. Statt etwas wie

n_compP1 = n_comp+1 
n_compT2 = n_comp*2 
a = Sum(x(1:n_comp)+2*x(n_compP1,n_compT2)) ! ... just for example 

2) Manchmal haben kann (abhängig von den Compiler, Maschine, etc.), kann es sein „Gedächtnis Behauptungen“, wenn die Array-Größen sind nicht auf bestimmte „binary Intervallen“ (zB 1024 vs. 1000) etc.

Sie können Ihre Experimente an einigen weiteren Punkten in Ihrem Diagramm wiederholen (dh an verschiedenen anderen "n_comps") und insbesondere in der Nähe solcher "Grenzen".

3) Kann nicht sagen, ob Sie die vollständigen Compiler-Optimierungen (Flags) zum Kompilieren Ihres Fortran-Codes verwenden. Vielleicht möchten Sie die verschiedenen "-o" -Flags usw. nachschlagen.

4) Möglicherweise möchten Sie die OpemMP-Anweisung einbeziehen (oder zumindest openmp in Ihre Flags usw. aufnehmen). Dies kann manchmal bestimmte Overhead-Probleme verbessern, auch wenn Sie nicht explizit auf OpenMP-Anweisungen in Ihren Schleifen usw. angewiesen sind.

5) Allgemeines: Dies würde zu jedem Ihrer Methoden wahrscheinlich anwendbar, wenn Schleifen

a) „constant operations“ in dem „Aufsummierung“ -Formel verwendet werden, kann außerhalb der Schleife durchgeführt werden, zum Beispiel. Erstellen Sie etwas wie qsDcs = qs/cs und verwenden Sie qsDcs in der Schleife.

b) In ähnlicher Weise ist es manchmal nützlich, etwas wie zM1 (:) = z (:) - 1 zu erstellen und stattdessen zM1 (:) in der Schleife zu verwenden.

Answer 3

Weiter zu meiner vorherigen Antwort und Vladimirs schwacher Spekulation, habe ich zwei s/rs eingerichtet: eins als das Original gegeben, und eins, das Array-Abschnitte und Sum() verwendet. Ich wollte auch zeigen, dass Wladimirs Ausführungen zur Do-Loop-Optimierung so schwach sind.

Auch ein Punkt, den ich normalerweise für das Benchmarking mache, ist die Größe von n_comp in dem oben gezeigten Beispiel zu klein. Die Ergebnisse unten setzen jede der "ursprünglichen" und der "besseren" SumArraySection (SAS) Variation in Schleifen, die 1000 Mal innerhalb der Timing-Aufrufe wiederholt werden, so dass die Ergebnisse für 1000 calcs jedes s/r sind. Wenn Ihre Zeiten Bruchteile von Sekunden sind, sind sie wahrscheinlich unzuverlässig.

Es gibt eine Reihe von Varianten, die in Betracht gezogen werden, keine mit expliziten Zeigern. Die eine Variante für diese Abbildungen verwendet wird, ist

subroutine get_denomSAS (qs,x,cp,cs,n_comp,denom) 

! Calculates the denominator in the SMA model (Brooks and Cramer 1992) 
! The function is called at a specific salt concentration and isotherm point 
! I loops over the number of components 

implicit none 

! declaration of input variables 
integer, intent(in) :: n_comp ! number of components 
double precision, intent(in) :: cs,qs ! salt concentration, free ligand concentration 
double precision, Intent(In)   :: cp(:) 
double precision, Intent(In)   :: x(:) 

! declaration of local variables 
integer :: i 

! declaration of outpur variables 
double precision, intent(out) :: denom 
! 
! 
double precision      :: qsDcs 
! 
! 
qsDcs = qs/cs 
! 
denom = Sum((x(n_comp+1:2*n_comp) + x(2*n_comp+1:3*n_comp))*(x(1:n_comp)*(qsDcs) & 
              **(x(2*n_comp+1:3*n_comp)-1))*cp(1:n_comp)) + cs 
! 
! 
end subroutine get_denomSAS 

Die Hauptunterschiede sind:

a) übergeben Arrays sind (:) b) keine Array-Zuweisungen in s/r, anstelle Feldabschnitte verwenden (äquivalent zu "effektive" Zeiger). c) Verwenden Sie Sum() anstelle von Do

Dann versuchen Sie auch zwei verschiedene Compiler-Optimierungen, um Auswirkungen zu demonstrieren.

Wie die beiden Diagramme zeigen, ist der Origcode (blaue Rauten) wesentlich langsamer. SAS (rote Quadrate) mit niedriger Optimierung. SAS ist immer noch besser mit hoher Optimierung, aber sie nähern sich. Dies wird teilweise dadurch erklärt, dass Sum() "besser optimiert" ist, wenn eine niedrige Compiler-Optimierung verwendet wird.

Die gelben Linien zeigen das Verhältnis zwischen den beiden S/Timings des r. Ignoriere den gelben Linienwert bei "0" im oberen Bild (n_comp zu klein verursacht eines der Timings wackelig)

Da ich die ursprünglichen Daten des Benutzers nicht im Verhältnis zu Numpy habe, kann ich nur das machen Aussage, dass die SAS-Kurve auf seinem Chart unter seinen aktuellen Fortran-Ergebnissen liegen sollte und möglicherweise flacher oder sogar fallend sein sollte.

Mit anderen Worten, es gibt möglicherweise nicht die Divergenz in der ursprünglichen Buchung, oder zumindest nicht in diesem Ausmaß.

... obwohl weitere Experimente hilfreich sein könnten, um auch die anderen bereits gegebenen Kommentare/Antworten zu demonstrieren.

Lieber Moritz: Hoppla, ich habe vergessen zu erwähnen, und zu Ihrer Frage über Zeiger. Wie bereits erwähnt, besteht ein Hauptgrund für die Verbesserung mit der SAS-Variante darin, dass sie "effektive Zeiger" besser ausnutzt, da sie die Notwendigkeit, das Array x() in drei neue lokale Arrays neu zuzuweisen, überflüssig macht (dh weil x durch ref übergeben wird) , Verwenden von Array-Abschnitten ist eine Art von Pointer-Ansatz, der in Fortran integriert ist und daher keine expliziten Zeiger benötigt), erfordert aber dann Sum() oder Dot_Product() oder was auch immer.

Stattdessen können Sie Do beibehalten und etwas Ähnliches erreichen, indem Sie x entweder in ein n_compx3-2D-Array ändern oder die drei expliziten 1D-Arrays der Ordnung n_comp direkt übergeben. Diese Entscheidung hängt wahrscheinlich von der Größe und Komplexität Ihres Codes ab, da Sie die aufrufenden/sr-Anweisungen usw. neu schreiben müssen und wo auch immer sonst x() verwendet wird. Einige unserer Projekte sind> 300.000 Zeilen Code, also ist es viel viel billiger, den Code lokal zu ändern, zB zu SAS etc.

Ich warte immer noch auf die Erlaubnis, Numpy auf einem der zu installieren unsere Boxen. Wie bereits erwähnt, ist es von Interesse, warum Ihre relativen Timings implizieren, dass sich Numpy mit steigendem n_comp ...

Natürlich die Kommentare über "richtige" Benchmarking usw., sowie die Frage, welche Compiler-Schalter beteiligt sind durch Ihre Verwendung von fpy, immer noch anwenden, da diese den Charakter Ihrer Ergebnisse stark verändern können.

Ich wäre daran interessiert, Ihre Ergebnisse zu sehen, wenn sie für diese Permutationen aktualisiert wurden.

Answer 4

in der anderen Antwort genannt zu werden, muss ich antworten.

Ich weiß, dass dies nicht wirklich die ursprüngliche Frage beantwortet, aber das ursprüngliche Plakat ermutigt, diese Richtung in seinen Kommentaren zu verfolgen.

Meine Punkte sind diese:

1. Ich glaube nicht, das Array intrinsische sind besser in irgendeiner Weise optimiert. Wenn man Glück hat, werden sie in den gleichen Schleifencode wie die manuellen Schleifen übersetzt. Wenn nicht, können Leistungsprobleme auftreten. Deshalb muss man vorsichtig sein. Es besteht die Möglichkeit, temporäre Arrays auszulösen.

Ich übersetzte die angebotenen SAS-Arrays zu üblichen Do-Schleife. Ich nenne es DOS. Ich zeige, dass die DO-Schleifen in keiner Weise langsamer sind, beide Unterprogramme führen in diesem Fall mehr oder weniger zum selben Code.

qsDcs = qs/cs 

denom = 0 
do j = 1, n_comp 
    denom = denom + (x(n_comp+j) + x(2*n_comp+j)) * (x(j)*(qsDcs)**(x(2*n_comp+j)-1))*cp(j) 
end do 

denom = denom + cs 

Es ist wichtig zu sagen, dass ich glaube nicht, diese Version weniger lesbar ist, nur weil es eine oder zwei weitere Linien. Es ist eigentlich ziemlich einfach zu sehen, was hier passiert.

Jetzt sind die Zeiten für diese

f2py -c -m sas sas.f90 --opt='-Ofast' 
f2py -c -m dos dos.f90 --opt='-Ofast' 


In [24]: %timeit test_sas(10000) 
1000 loops, best of 3: 796 µs per loop 

In [25]: %timeit test_sas(10000) 
1000 loops, best of 3: 793 µs per loop 

In [26]: %timeit test_dos(10000) 
1000 loops, best of 3: 795 µs per loop 

In [27]: %timeit test_dos(10000) 
1000 loops, best of 3: 797 µs per loop 

Sie sind genau das gleiche. Es gibt keine versteckte Leistungszauberei in der array intrinsics und Array-Ausdruckarithmetik. In diesem Fall werden sie nur in Schleifen unter der Haube übersetzt.

Wenn Sie den erzeugten GIMPLE Code überprüfen, sowohl die SAS und DOS auf den gleichen Hauptblock von optimierten Code übersetzt werden, keine magische Version von SUM() hier genannt:

<bb 8>: 
    # val.8_59 = PHI <val.8_49(9), 0.0(7)> 
    # ivtmp.18_123 = PHI <ivtmp.18_122(9), 0(7)> 
    # ivtmp.25_121 = PHI <ivtmp.25_120(9), ivtmp.25_117(7)> 
    # ivtmp.28_116 = PHI <ivtmp.28_115(9), ivtmp.28_112(7)> 
    _111 = (void *) ivtmp.25_121; 
    _32 = MEM[base: _111, index: _106, step: 8, offset: 0B]; 
    _36 = MEM[base: _111, index: _99, step: 8, offset: 0B]; 
    _37 = _36 + _32; 
    _40 = MEM[base: _111, offset: 0B]; 
    _41 = _36 - 1.0e+0; 
    _42 = __builtin_pow (qsdcs_18, _41); 
    _97 = (void *) ivtmp.28_116; 
    _47 = MEM[base: _97, offset: 0B]; 
    _43 = _40 * _47; 
    _44 = _43 * _42; 
    _48 = _44 * _37; 
    val.8_49 = val.8_59 + _48; 
    ivtmp.18_122 = ivtmp.18_123 + 1; 
    ivtmp.25_120 = ivtmp.25_121 + _118; 
    ivtmp.28_115 = ivtmp.28_116 + _113; 
    if (ivtmp.18_122 == _96) 
    goto <bb 10>; 
    else 
    goto <bb 9>; 

    <bb 9>: 
    goto <bb 8>; 

    <bb 10>: 
    # val.8_13 = PHI <val.8_49(8), 0.0(6)> 
    _51 = val.8_13 + _17; 
    *denom_52(D) = _51; 

der Code ist funktional identisch mit die Do-Schleife Version, nur der Name der Variablen sind unterschiedlich.

---

2. Sie nahmen Argumente Form Array (:) ein Potenzial haben, die Leistung zu beeinträchtigen.Während das im angenommenen Größenargument (*) oder explizite Größe (n) erhaltene Argument immer einfach zusammenhängend ist, muss die angenommene Form theoretisch nicht sein und der Compiler muss darauf vorbereitet sein. Daher empfehle ich immer, das Attribut contiguous für Ihre angenommenen Formargumente zu verwenden, wo immer Sie wissen, dass Sie es mit zusammenhängenden Arrays aufrufen.

In der anderen Antwort war die Änderung ziemlich sinnlos, weil sie keinen der Vorteile angenommener Formargumente verwendet hat. Nämlich, dass Sie die Argumente mit den Array-Größen nicht übergeben müssen, und Sie können die intrinsics wie size() und shape() verwenden.

---

Hier sind die Ergebnisse eines umfassenden Vergleichs. Ich habe es so fair wie möglich gemacht. Fortran-Dateien werden mit -Ofast zusammengestellt, wie oben gezeigt:

import numpy as np 
import dos as dos 
import sas as sas 
from matplotlib import pyplot as plt 
import timeit 
import numexpr as ne 

#%matplotlib inline 



ne.set_num_threads(1) 

def test_n(n_comp): 

    cp = np.tile(1.243,n_comp) 
    cs = 100. 
    qs = np.tile(1100.,n_comp) 
    x= np.random.rand(3*n_comp+1) 

    def test_dos(): 
     denom = np.empty(1) 
     dos.get_denomsas(qs,x,cp,cs,n_comp) 


    def test_sas(): 
     denom = np.empty(1) 
     sas.get_denomsas(qs,x,cp,cs,n_comp) 

    def get_denom(): 
     k = x[0:n_comp] 
     sigma = x[n_comp:2*n_comp] 
     z = x[2*n_comp:3*n_comp] 
     # calculates the denominator in Equ 14a - 14c (Brooks & Cramer 1992) 
     a = (sigma + z)*(k*(qs/cs)**(z-1))*cp 
     denom = np.sum(a) + cs 
     return denom 

    def get_denom_numexp(): 
     k = x[0:n_comp] 
     sigma = x[n_comp:2*n_comp] 
     z = x[2*n_comp:3*n_comp] 
     loc_cp = cp 
     loc_cs = cs 
     loc_qs = qs 
     # calculates the denominator in Equ 14a - 14c (Brooks & Cramer 1992) 
     a = ne.evaluate('(sigma + z)*(k*(loc_qs/loc_cs)**(z-1))*loc_cp') 
     return cs + np.sum(a) 

    print 'py', timeit.Timer(get_denom).timeit(1000000/n_comp) 
    print 'dos', timeit.Timer(test_dos).timeit(1000000/n_comp) 
    print 'sas', timeit.Timer(test_sas).timeit(1000000/n_comp) 
    print 'ne', timeit.Timer(get_denom_numexp).timeit(1000000/n_comp) 


def test(): 
    for n in [10,100,1000,10000,100000,1000000]: 
     print "-----" 
     print n 
     test_n(n) 

Ergebnisse:

  py    dos    sas    numexpr 
10   11.2188110352 1.8704519272 1.8659651279 28.6881871223 
100   1.6688809395 0.6675260067 0.667083025  3.4943861961 
1000  0.7014708519 0.5406000614 0.5441288948 0.9069931507 
10000  0.5825948715 0.5269498825 0.5309231281 0.6178650856 
100000  0.5736029148 0.526198864  0.5304090977 0.5886831284 
1000000  0.6355218887 0.5294830799 0.5366530418 0.5983200073 
10000000 0.7903120518 0.5301260948 0.5367569923 0.6030929089 

Sie können sehen, dass es sehr kleiner Unterschied zwischen den beiden Versionen Fortran ist. Die Array-Syntax ist marginal langsamer, aber wirklich nichts zu sagen.

Schlussfolgerung 1: In diesem Vergleich Aufwand für alle sollte fair und Sie sehen, dass für ideale Vektorlänge Numpy und numexpr CAN fast Fortran Leistung erreichen, aber wenn der Vektor zu klein ist oder vielleicht sogar zu groß die Overhead der Python-Lösungen herrscht vor.

Fazit 2: Die leistungsfähigere SAS-Version im anderen Vergleich wird durch den Vergleich mit der ursprünglichen OP-Version verursacht, die nicht äquivalent ist. Die äquivalente optimierte DO-Loop-Version ist oben in meiner Antwort enthalten.