Code optimieren

[Absolvent1993]

Guten Tag,

ich habe nummerische Simulation geschrieben, die sehr rechenintensiv ist. Die Rechendauer an meinem Coputer betraegt fuer den Parameter T = 10^4 (Anzahl der Schleifendurchlaeufe) fast 42 Stunden.

Um ein exaktes Resultat zu erhalten, muss der Parameter jedoch auf T = 10^6 gesetzt werden ---> Rechnendauer ca. 4000 std.

Kann sich bitte jemand von euch mal meinen Code anschauen und vielleicht eine Moeglichkeit finden um die Rechendauer zu verringern.

Code: Alles auswählen

from numpy import cos, sin, arctan, sqrt, pi, array, random, mean, linspace,log, zeros, shape, array, savetxt, dot, std, fill_diagonal
import matplotlib.pyplot as plt
from numpy.linalg import norm
import time
import numpy as np

# Initial Condition

#---------------------------------------------------------------
T  = 10**4                         # time steps 
#---------------------------------------------------------------

Teq = 10**3                        # transition time

N  = 50                            # Number of Oscillators
dt = 0.01                          # time increment

Tepsilon  = 33                     # number of steps
depsilon = 0.05                    # increment of coupling constant

d = zeros([N,N])
fill_diagonal(d,1)
L_values = []
epsilon_points = [] 

def field( phi ) :                              
    real = 1 / N * sum( cos( phi ))
    imag = 1 / N * sum( sin( phi ))
    R   = sqrt( real**2 + imag**2 )
    Phi = 2 * arctan( imag / ( real+R ))
    
    return R,Phi

def phi_dot( phi, epsilon, R, Phi ) :               # time evolutions of phases
    return  omega + epsilon * R * sin( Phi - phi ) 


def jacobian(phi):                 # compute jacobi-matrix of phi_dot
    A = zeros([N,N])
    for i in range(N):
        for j in range(N):
            if i != j :
                A[i,j] = cos(phi[j] - phi[i])  # non diagonal elements
            else:
                mask = np.ones(N, bool)
                mask[i] = False
                phi_new = phi[mask]                 # phi[i] is excluded from phi
    
                A[i,i] = sum(-cos(phi_new - phi[i])) # diagonal elements of the jacobian
        
    return (epsilon / N) * A

# The tangent space (linearized) flow 
def g(d, phi):
    M = jacobian(phi)
    return dot(d,M.T)

# Runge-Kutta integrator
def RK4(phi, d):
    k1  = dt * phi_dot(phi, epsilon, R, Phi)                 
    k11 = dt * g(d, phi)
    k2  = dt * phi_dot(phi + 0.5 * k1, epsilon, R, Phi)
    k22 = dt * g(d + 0.5 * k11, phi + 0.5 * k1)
    k3  = dt * phi_dot(phi + 0.5 * k2, epsilon, R, Phi)
    k33 = dt * g(d + 0.5 * k22, phi + 0.5 * k2)
    k4  = dt * phi_dot(phi + k3, epsilon, R, Phi)
    k44 = dt * g(d + k33, phi + k3)
    
    phi  += (k1  + 2 * k2  + 2 * k3  + k4)  / 6
    d    += (k11 + 2 * k22 + 2 * k33 + k44) / 6
    
    return phi, d

# Gram-Schmidt - Scheme
def gram_schmidt(A, method = 'transient'):
    n = A.shape[1]
    if method == 'transient':
        for i in range(n):
            for j in range(i):
                A[i] -= dot(d[i], d[j]) * d[j]
                
            A[i]  = A[i] / norm(A[i])   
        
        return A
    
    if method == 'computation':
        for i in range(n):
            for j in range(i):
                A[i] -= dot(d[i], d[j]) * d[j]
        
            L[i] += log(norm(A[i]))
            A[i]  = A[i] / norm(A[i])   
        
        return A

#---------------------------------------------------------------------------------------
#-------------------------------Main Loop-----------------------------------------------   
#---------------------------------------------------------------------------------------
    
start_time = time.time()

for i in range( Tepsilon ) :
    
    epsilon = i * depsilon
    epsilon_points.append(epsilon)
    omega = linspace(-1,1,N)                          # regular sampling from [-1,1]
    phi = random.uniform(low=0.0, high=2*pi, size=N)  # reset initial condition for every new epsilon
    
    L = zeros( N )                                    # ...
    
    # Transient
    for j in range( Teq ) :
        R,Phi = field(phi)
        phi, d = RK4(phi, d)
        d = gram_schmidt(d, 'transient')
        
    # Computation    
    for k in range(T):
        R,Phi = field(phi)
        r, d = RK4(phi, d)
        d = gram_schmidt(d, 'computation')
    
    L_values.append( list( L ))            # contains avg freq of oscl for a given epsilon

    
print("--- %s seconds ---" % (time.time() - start_time))

[__deets__]

Du kannst mal PyPy probieren. Das kompiliert sowas wie deine RK Implementierung in effizienten Maschinencode.

[Sirius3]

Wenn Du Schleifen bei numerischen Berechnungen hast, dann machst Du etwas falsch.
Die Berechnung in `jacobian` könnte z.B. so aussehen:

Code: Alles auswählen

def jacobian(phi):
    """ compute jacobi-matrix of phi_dot """
    A = np.cos(phi[:, None] - phi[None, :])
    A[np.diag_indices(len(phi))] = 1 - A.sum(axis=0)
    return (epsilon / len(phi)) * A

Dass da viele Namen direkt aus numpy importiert werden, macht es schwierig zu sehen, ob man auch wirklich immer die richtige Funktion benutzt. `sum` ist z.B. das eingebaute sum aus Python und nicht das effektive aus numpy.

Globale Variablen sollte man nicht benutzten, und sprechendere Namen als griechische Groß- und Kleinbuchstaben wären gut.

[einfachTobi]

Das Runge-Kutta-Verfahren ist auch bereits in scipy implementiert und deutlich schneller. Die Berechnung einer Jacobi-Matrix ebenfalls. Gram-Schmidt 'transient' kannst du ebenfalls du die im anderen Thread genannte Methode ersetzen, weil du dort deine L-Matrix ja nicht berechnest.
Was berechnest du da überhaupt? Ggf. gibt es dafür ja auch andere Wege.
Wenn du die einzelnen Algorithmen geändert hast, könntest du noch über Parallelisierung (falls nicht voneinander abhängig, weil dann ists einfach) nachdenken und dir numba mit dem @jit(nopython=True) Decorator ansehen.

[Sirius3]

Nachdem ich mir jetzt den Code genauer angeschaut habe, da ist in `RK4` einiges im Argen. Du änderst die übergebenen Arrays und gibts sie aber auch wieder zurück. In der ersten Schleife werden ja die Arrays wieder an die selben Namen gebunden, im zweiten aber nicht. Da weiß ich jetzt nicht, ob das ändern von phi ein Fehler ist oder beabsichtigt. Dann sind einige Arrays nicht Argumente sondern globale Variablen, was die Benutzung unten völlig unleserlich macht. Warum sollte man R ausrechen, wenn es scheinbar nicht benutzt wird?
Sich ändernde Arrays sollte man nicht zurückgeben. Klarer wird es, wenn man nur die Deltas zurückgibt und die Addition in der Hauptschleife macht.

`gram_schmidt` habe ich noch mit numba beschleunigt, der Rest ist ja jetzt schon vektoriell, da ist also nicht mehr viel rauszuholen.

Code: Alles auswählen

import time
from numpy.linalg import norm
import numpy as np
from numba import njit
import matplotlib.pyplot as plt

T  = 10**4                         # time steps 
Teq = 10**3                        # transition time

N  = 50                            # Number of Oscillators
dt = 0.01                          # time increment

Tepsilon  = 33                     # number of steps
depsilon = 0.05                    # increment of coupling constant
omega = np.linspace(-1,1,N)                          # regular sampling from [-1,1]

def field(phi):                              
    real = np.cos(phi).mean()
    imag = np.sin(phi).mean()
    r   = (real**2 + imag**2) ** 0.5
    phi = 2 * np.arctan( imag / (real + r))
    return r, phi

def phi_dot( phi, epsilon, R, Phi ) :               # time evolutions of phases
    return omega + epsilon * R * np.sin(Phi - phi) 

def jacobian(phi, epsilon):
    """ compute jacobi-matrix of phi_dot """
    A = np.cos(phi[:, None] - phi[None, :])
    A[np.diag_indices(len(phi))] = 1 - A.sum(axis=0)
    return (epsilon / len(phi)) * A

# The tangent space (linearized) flow 
def g(d, phi, epsilon):
    M = jacobian(phi, epsilon)
    return d.dot(M.T)

# Runge-Kutta integrator
def RK4(R, Phi, phi, d, epsilon):
    k1  = dt * phi_dot(phi, epsilon, R, Phi)                 
    k11 = dt * g(d, phi, epsilon)
    k2  = dt * phi_dot(phi + 0.5 * k1, epsilon, R, Phi)
    k22 = dt * g(d + 0.5 * k11, phi + 0.5 * k1, epsilon)
    k3  = dt * phi_dot(phi + 0.5 * k2, epsilon, R, Phi)
    k33 = dt * g(d + 0.5 * k22, phi + 0.5 * k2, epsilon)
    k4  = dt * phi_dot(phi + k3, epsilon, R, Phi)
    k44 = dt * g(d + k33, phi + k3, epsilon)
    
    dphi = (k1  + 2 * k2  + 2 * k3  + k4)  / 6
    dd = (k11 + 2 * k22 + 2 * k33 + k44) / 6
    return dphi, dd

# Gram-Schmidt - Scheme
@njit
def gram_schmidt(L, A):
    n = A.shape[1]
    for i in range(n):
        A[i] -= (np.dot(A[i], A[:i].T)*A[:i].T).sum(axis=1)
        m = norm(A[i])
        L[i] += np.log(m)
        A[i] /= m

def main():
    start_time = time.time()
    d = np.eye(N)
    L_values = []
    epsilon_points = np.linspace(0, (Tepsilon-1)*depsilon, Tepsilon)
    for epsilon in epsilon_points:
        # reset initial condition for every new epsilon
        phi = np.random.uniform(low=0.0, high=2*np.pi, size=d.shape[0])

        # Transient
        L = np.zeros_like(phi)
        for j in range(Teq):
            R, Phi = field(phi)
            dphi, dd = RK4(R, Phi, phi, d, epsilon)
            phi += dphi
            d += dd
            gram_schmidt(L, d)
            
        # Computation    
        L = np.zeros_like(phi)
        for k in range(T):
            R,Phi = field(phi)
            dphi, dd = RK4(R, Phi, phi, d, epsilon)
            phi += dphi
            d += dd
            gram_schmidt(L,d)
        
        L_values.append(L.tolist())            # contains avg freq of oscl for a given epsilon
        print("--- %s seconds ---" % (time.time() - start_time))

if __name__ == '__main__':
    main()

[__blackjack__]

Was man sowieso machen sollte ist das Hauptprogramm nicht auf Modulebene zu schreiben sondern auch in eine Funktion zu stecken. Nicht aus Performance-Gründen sondern weil es sicherer ist das man sich nicht aus versehen globalen Zustand einfängt. Aber es hat den Nebeneffekt das CPython lokale Namen geringfügig schneller auflöst als globale Namen.

Da fällt dann auf dass das nicht sauber gelöst ist, denn Deine Funktionen verändern globale Variablen – was nicht sein sollte. Funktionen/Methoden bekommen alles was sie benötigen als Argument(e) übergeben. Schön wäre dann noch wenn sie nichts verändern was sie übergeben bekommen, also ”reine” Funktionen sind, denn da kann man dann überlegen ob und was man eventuell parallelisieren kann.

`d` sollte übergeben werden und bis auf in `RK4()` wird das auch nicht verändert — dort ist es trivial die Veränderung zu vermeiden. Die API von `RK4()` ist an der Stelle auch falsch `d` wird als Argument übergeben, in der Funktion verändert *und* als Rückgabewert an den Aufrufer zurückgegeben. Wenn man ein übergebenes Objekt verändert macht es keinen Sinn das als Ergebnis zurück zu geben, denn die Veränderung hat der Aufrufer ja auch schon so, da er das Objekt bereits haben muss um es übergeben zu können. Das gleiche gilt für das andere Argument der Funktion (`phi`).

An `gram_schmidt()` müsste man das `d` als Argument übergeben. Wenn man das macht, dann sehen die beiden Aufrufe der Funktion so aus:

Code: Alles auswählen

            ...
            d = gram_schmidt(d, d, "transient")
            ...
            d = gram_schmidt(d, d, "computation")
            ...

Das heisst die Argumente `A` und `d` sind immer das *selbe* Objekt, was wenn man jetzt noch bedenkt das `d` vorher in der Funktion sowohl über den Namen `A` als auch global über `d` angesprochen wurde ziemlich schräg/verwirrend ist. In `gram_schmidt()` kann man also alle `d` in `A` ändern, denn das ist faktisch das selbe Objekt über zwei Namen erreichbar. Und `gram_schmidt()` hat auch wieder das Problem, dass `A`/`d` übergeben, verändert, und als Ergebnis zurückgegeben wird.

Zudem würde ein Tippfehler im `method`-Argument dazu führen das die Funktion nichts macht und `None` zurück gibt. Es ist auch nicht sinnvoll zweimal ``if`` mit `method` zu verwenden wenn sich beide Fälle gegenseitig ausschliessen. Dafür gibt es ``elif``. Und dann noch ein ``else`` dahinter das den Fall behandelt wenn `method` einen falschen Wert hat. Einen Tippfehler bei den Zeichenketten kann man auch früher ausschliessen wenn man dafür Konstanten definiert. Vielleicht auch als `enum.Enum` weil man dann mit ``is`` statt mit ``==`` testen kann. Andererseits gibt es ja nur zwei Methoden wobei eine etwas zusätzlich macht was die andere nicht macht, man könnte da also auch einfach einen Wahrheitswert übergeben der darüber entscheidet ob der zusätzliche Schritt gemacht wird oder nicht. Dann kommt man mit einem ``if``/``else`` aus, ohne einen Vergleichsoperator.

Das da in beiden ``if``-Zweigen fast der gleiche Code steht ist unschön. Diese Duplikation kann man vermeiden wenn man den Test weiter nach innen verschiebt.

Das `L` in der Funktion ist auch eine globale Variable. Da sollte eigentlich ein neuer Wert berechnet und zurückgegeben werden.

Wobei das eventuell auch einfach sehr ungünstig aufgeteilt ist. Man bräuchte `method` oder ein ein Flag gar nicht wenn diese Funktion nur den Teil machen würde der beiden Methoden gleich ist und den Rest ausserhalb der Funktion. Ein `L` braucht die Funktion dann gar nicht kennen.

Dann ist da noch `epsilon` das global ist und ausser bei `phi_dot()` nirgends als Argument übergeben wird wo es verwendet wird.

Zuguterletzt dann noch `Phi` mit grossem P und `R` bei `RK4()`.

`epsilon_points` kann man doch eigentlich *vor* der Hauptschleife schneller und speicherplatzsparender mit `numpy.arange()` und Multiplikation erstellen und dann darüber iterieren.

`omega` wird in der Haupt-``for``-Schleife immer wieder neu definiert, aber immer mit dem selben Wert der auch nirgends verändert wird.

Namen die man nur aus syntaktischen Gründen braucht, die man aber nicht benutzt, beispielsweise Laufvariablen in Schleifen die man nur wegen der Wiederholung des Schleifenkörpers hat, werden per Konvention `_` genannt. Dann weiss der Leser sofort das der Wert nicht verwendet wird, und dass das auch Absicht ist, und man braucht nicht darauf achten das keine andere Laufvariable im gleichen Namensraum damit kollidieren könnte.

Bei der zweiten inneren ``for``-Schleife ändert sich `phi` nie, weil das Ergebnis von `RK4()` an den Namen `r` gebunden wird statt an `phi`, das `r` dann aber auch nirgends verwendet wird. Soll das so sein? Falls ja, sollte man auch `r` in `_` umbenennen, damit der Leser leicht erkennt, dass das so gewollt und kein Fehler ist.

Warum wird `L` am Ende in eine Liste umgewandelt? Speichereffizienter wäre ein Array. Schneller vielleicht auch.

Für Laufzeitmessungen würde ich `time.monotonic()` verwenden.

``%`` für Zeichenkettenformatierung würde ich nicht mal mehr in Python 2.7 verwenden.

Ungetestet:

Code: Alles auswählen

#!/usr/bin/env python3
import time

import numpy as np
from numpy import (
    arange,
    arctan,
    cos,
    dot,
    fill_diagonal,
    linspace,
    log,
    pi,
    random,
    sin,
    sqrt,
    zeros,
)
from numpy.linalg import norm


T = 10 ** 4  # time steps

Teq = 10 ** 3  # transition time

N = 50  # Number of Oscillators
dt = 0.01  # time increment

Tepsilon = 33  # number of steps
depsilon = 0.05  # increment of coupling constant

omega = linspace(-1, 1, N)  # regular sampling from [-1,1]


def field(phi):
    real = 1 / N * sum(cos(phi))
    imag = 1 / N * sum(sin(phi))
    R = sqrt(real ** 2 + imag ** 2)
    Phi = 2 * arctan(imag / (real + R))

    return R, Phi


def phi_dot(phi, epsilon, R, Phi):  # time evolutions of phases
    return omega + epsilon * R * sin(Phi - phi)


def jacobian(phi, epsilon):  # compute jacobi-matrix of phi_dot
    A = zeros([N, N])
    for i in range(N):
        for j in range(N):
            if i != j:
                A[i, j] = cos(phi[j] - phi[i])  # non diagonal elements
            else:
                mask = np.ones(N, bool)
                mask[i] = False
                phi_new = phi[mask]  # phi[i] is excluded from phi
                #
                # diagonal elements of the jacobian
                #
                A[i, i] = sum(-cos(phi_new - phi[i]))

    return (epsilon / N) * A


# The tangent space (linearized) flow
def g(d, phi, epsilon):
    return dot(d, jacobian(phi, epsilon).T)


# Runge-Kutta integrator
def RK4(phi, d, epsilon, R, Phi):
    k1 = dt * phi_dot(phi, epsilon, R, Phi)
    k11 = dt * g(d, phi, epsilon)
    k2 = dt * phi_dot(phi + 0.5 * k1, epsilon, R, Phi)
    k22 = dt * g(d + 0.5 * k11, phi + 0.5 * k1, epsilon)
    k3 = dt * phi_dot(phi + 0.5 * k2, epsilon, R, Phi)
    k33 = dt * g(d + 0.5 * k22, phi + 0.5 * k2, epsilon)
    k4 = dt * phi_dot(phi + k3, epsilon, R, Phi)
    k44 = dt * g(d + k33, phi + k3, epsilon)

    phi = phi + (k1 + 2 * k2 + 2 * k3 + k4) / 6
    d = d + (k11 + 2 * k22 + 2 * k33 + k44) / 6

    return phi, d


# Gram-Schmidt - Scheme
def gram_schmidt(A):
    norm_of_A = zeros(A.shape)

    for i in range(A.shape[1]):
        for j in range(i):
            A[i] -= dot(A[i], A[j]) * A[j]
        norm_of_A[i] = norm(A[i])

    return norm_of_A


def main():
    d = zeros([N, N])
    fill_diagonal(d, 1)

    L_values = []
    epsilon_points = arange(Tepsilon) * depsilon

    start_time = time.monotonic()

    for epsilon in epsilon_points:
        #
        # reset initial condition for every new epsilon
        #
        phi = random.uniform(low=0.0, high=2 * pi, size=N)

        # Transient
        for _ in range(Teq):
            R, Phi = field(phi)
            phi, d = RK4(phi, d, epsilon, R, Phi)
            d /= gram_schmidt(d)

        # Computation
        L = zeros(N)
        for _ in range(T):
            R, Phi = field(phi)
            #
            # TODO Soll `_phi` wirklich nirgends verwendet werden?
            #
            _phi, d = RK4(phi, d, epsilon, R, Phi)
            d_new = gram_schmidt(d)
            L += log(d_new)
            d /= d_new
        #
        # contains avg freq of oscl for a given epsilon
        #
        L_values.append(L.copy())

    print(f"--- {(time.monotonic() - start_time)} seconds ---")


if __name__ == "__main__":
    main()

Wenn man etwas schneller machen will, wäre der erste Schritt messen was eigentlich wie viel Zeit verbraucht. Der nächste Schritt die langsamen Stellen versuchen mit besseren Algorithmen schneller zu machen. Als letzter Schritt kämen ”Micro-Optimierungen”, die würde ich mir aber gut überlegen wenn sie den Code schlechter machen. Wie das kopieren fast einer ganzen Funktion wie bei Deinem Quelltext. Da würde ich eher zu Cython oder so greifen statt unschöne Code-Duplikation ins Spiel zu bringen.