Ottimizzazione del calcolo della gravitazione per particelle a gravità zero spazio 2d

Question

Ho creato una piccola visualizzazione di particelle in python. Sono in grado di caclulare il movimento di particelle in uno spazio 2D senza gravità. Come ogni particella attira tutte le altre particelle in base alla massa e alla distanza delle particelle.

Ho realizzato una visualizzazione in pygame e tutto funziona come da piano (con la calucazione), tuttavia ho bisogno di ottimizzare il calcolo in modo estrattivo. Oggi il sistema può calcolare circa 100-150 particelle in un framerate deacent. Ho messo tutto il calcolo in un thread separato che mi ha dato un po 'di più ma non quasi quello che voglio.

Ho guardato scipy e insensibile, ma dal momento che non sono uno scienziato o un matemogruppo mi sento solo confuso. Sembra che io sia sulla strada giusta ma non ho idea di come.

Ho bisogno di calcolare tutta l'attrazione su tutte le particelle che ho a un ciclo in un ciclo. E dal momento che ho bisogno di trovare se qualcuno si è scontrato, devo fare lo stesso tutto da capo.

Mi si spezza il cuore a scrivere quel tipo di codice ....

Numpy ha la capacità di calcolare array con array, ma non ho trovato alcuna cosa per calcolare tutti gli elementi in array con tutte le voci dallo stesso/altro array. Ce n'è uno? Se è così ho potuto creare e un paio di matrici e calcolare molto più veloce e ci deve essere una funzione per ottenere indice da a 2 array dove il loro valori partita (Collitiondetect IOW)

Qui è oggi calcolo di attrazione/collsion:

class Particle: 
    def __init__(self): 
     self.x = random.randint(10,790) 
     self.y = random.randint(10,590) 
     self.speedx = 0.0 
     self.speedy = 0.0 
     self.mass = 4 

#Attraction  
for p in Particles: 
    for p2 in Particles: 
     if p != p2: 
      xdiff = P.x - P2.x 
      ydiff = P.y - P2.y 
      dist = math.sqrt((xdiff**2)+(ydiff**2)) 
      force = 0.125*(p.mass*p2.mass)/(dist**2) 
      acceleration = force/p.mass 
      xc = xdiff/dist 
      yc = ydiff/dist 
      P.speedx -= acceleration * xc 
      P.speedy -= acceleration * yc 
for p in Particles: 
    p.x += p.speedx 
    p.y += p.speedy 

#Collision 
for P in Particles: 
    for P2 in Particles: 
     if p != P2: 
      Distance = math.sqrt( ((p.x-P2.x)**2) + ((p.y-P2.y)**2) ) 
      if Distance < (p.radius+P2.radius): 
       p.speedx = ((p.mass*p.speedx)+(P2.mass*P2.speedx))/(p.mass+P2.mass) 
       p.speedy = ((p.mass*p.speedy)+(P2.mass*P2.speedy))/(p.mass+P2.mass) 
       p.x = ((p.mass*p.x)+(P2.mass*P2.x))/(p.mass+P2.mass) 
       p.y = ((p.mass*p.y)+(P2.mass*P2.y))/(p.mass+P2.mass) 
       p.mass += P2.mass 
       p.radius = math.sqrt(p.mass) 
       Particles.remove(P2) 

Answer 1

(Questo può dovrebbe andare in un commento, ma non ho la reputazione necessaria per farlo)

io non vedo come si fa il tempo di fare un passo. Hai

P.speedx -= acceleration * xc 
P.speedy -= acceleration * yc 

ma per ottenere la nuova velocità al tempo t + delta_t fareste

P.speedx -= acceleration * xc * delta_t 
P.speedy -= acceleration * yc * delta_t 

e quindi aggiornare la posizione in questo modo:

P.x = P.x + P.speedx * delta_t 
P.y = P.y + P.speedy * delta_t 

Quindi per la vostra velocità di preoccupazione . Forse sarebbe meglio memorizzare le informazioni sulle particelle non in una classe, ma in array numpy? Ma non penso che tu possa evitare i loop.

Inoltre, avete guardato allo wikipedia, qui vengono descritti alcuni metodi per accelerare il calcolo.

(a cura dovuta al commento di Mike)

Answer 2

Ho lavorato su questo in precedenza, e una delle cose che ho visto in passato per accelerare i calcoli collisione è quello di archiviare una lista di particelle vicine.

In sostanza, l'idea è all'interno del calcolo gravità si fa qualcosa di simile:

for (int i = 0; i < n; i++) 
{ 
    for (int j = i + 1; j < n; j++) 
    { 
     DoGravity(Particle[i], Particle[j]); 
     if (IsClose(Particle[i], Particle[j])) 
     { 
      Particle[i].AddNeighbor(Particle[j]); 
      Particle[j].AddNeighbor(Particle[i]); 
     } 
    } 
} 

Poi, è sufficiente passare sopra tutte le particelle e si fa rilevare la collisione di ognuno su a turno.Questo è solitamente qualcosa come O(n) nel migliore dei casi, ma può facilmente degradarsi a O(n^2) nel peggiore dei casi.

Un'altra alternativa è provare a collocare le particelle all'interno di un Octree. Costruire uno è qualcosa come O(n), quindi è possibile interrogarlo per vedere se qualcosa è vicino l'un l'altro. A quel punto dovresti solo rilevare le collisioni sulle coppie. Faccio questo è O(n log n) credo.

Non solo, ma è possibile utilizzare l'Octree per accelerare anche il calcolo della gravità. Invece del comportamento O(n^2), scende anche a O(n log n). La maggior parte delle implementazioni di Octree include un "parametro di apertura" che controlla la velocità rispetto all'accuratezza che farai. Quindi gli Octies tendono ad essere meno precisi di un calcolo diretto a coppie e complicati da codificare, ma rendono anche possibili simulazioni su larga scala.

Se si utilizza l'Octree in questo modo, si farà ciò che è noto come Barnes-Hut Simulation.

Nota: poiché si lavora in 2D, l'analogico 2D a Octree è noto come Quadtree. Vedere il seguente articolo di Wikipedia per ulteriori informazioni: http://en.wikipedia.org/wiki/Quadtree

Answer 3

per eseguire il calcolo veloce, è necessario memorizzare x, y, veloce, veloce, m in array numpy. Ad esempio:

import numpy as np 

p = np.array([ 
    (0,0), 
    (1,0), 
    (0,1), 
    (1,1), 
    (2,2), 
], dtype = np.float) 

p è un array 5x2 che memorizza la posizione x, y delle particelle. Per calcolare la distanza tra ogni coppia, è possibile utilizzare:

print np.sqrt(np.sum((p[:, np.newaxis] - p[np.newaxis, :])**2, axis=-1)) 

l'output è:

[[ 0.   1.   1.   1.41421356 2.82842712] 
[ 1.   0.   1.41421356 1.   2.23606798] 
[ 1.   1.41421356 0.   1.   2.23606798] 
[ 1.41421356 1.   1.   0.   1.41421356] 
[ 2.82842712 2.23606798 2.23606798 1.41421356 0.  ]] 

oppure è possibile utilizzare cdist da SciPy:

from scipy.spatial.distance import cdist 
print cdist(p, p) 

Answer 4

È possibile provare prima a lavorare con i numeri complessi : le formule relative alla gravitazione e alla dinamica sono molto semplici in questo formalismo e possono anche essere abbastanza veloci (perché NumPy può fare il calcolo internamente, invece di gestire separatamente le coordinate xey). Ad esempio, la forza fra due particules a Z e Z 'è semplicemente:

(z-z')/abs(z-z')**3 

NumPy può calcolare una quantità molto rapidamente, per ogni z/z' coppie. Ad esempio, la matrice di tutti i valori z-z 'viene semplicemente ottenuta dall'array 1D Z di coordinate come Z-Z[:, numpy.newaxis] (i termini diagonali [z = z'] richiedono una particolare attenzione nel calcolo di 1/abs(z-z')**3: devono essere impostati su zero).

Per quanto riguarda l'evoluzione nel tempo, si può certamente utilizzare veloce equazione differenziale di SciPy routines: sono molto più preciso rispetto al passo dopo passo l'integrazione di Eulero.

In ogni caso, scavare in NumPy sarebbe molto utile, specialmente se si prevede di fare calcoli scientifici, dato che NumPy è molto veloce.

Answer 5

Non sono sicuro se questo ti aiuterà molto, ma fa parte di una soluzione su cui ho lavorato per lo stesso problema. Non ho notato un enorme guadagno in termini di prestazioni in questo modo, ancora iniziando a fermarmi a circa 200 particelle, ma forse ti darà qualche idea.

modulo

C++ per il calcolo del componenti xey di attrazione gravitazionale su un piano 2D:

#include <Python.h> 
#include <math.h> 

double _acceleration(double &Vxa, double &Vya, double &Vxb, double &Vyb, double xa, double ya, double xb, double yb, double massa, double massb) 
{ 
    double xdiff = xa - xb; 
    double ydiff = ya - yb; 
    double distance = sqrt(xdiff*xdiff + ydiff*ydiff) * pow(10, 5); 

    if (distance <= 0) 
     distance = 1; 

    double force = (6.674 * pow(10, -11))*(massa*massb)/(distance*distance); 

    double acca = force/massa; 
    double accb = force/massb; 
    double xcomponent = xdiff/distance; 
    double ycomponent = ydiff/distance; 

    Vxa -= acca * xcomponent; 
    Vya -= acca * ycomponent; 
    Vxb += accb * xcomponent; 
    Vyb += accb * ycomponent; 

    return distance; 
} 

static PyObject* gforces(PyObject* self, PyObject* args) 
{ 
    double Vxa, Vya, Vxb, Vyb, xa, ya, xb, yb, massa, massb, distance; 

    if (!PyArg_ParseTuple(args, "dddddddddd", &Vxa, &Vya, &Vxb, &Vyb, &xa, &ya, &xb, &yb, &massa, &massb)) 
     return NULL; 

    distance = _acceleration(Vxa, Vya, Vxb, Vyb, xa, ya, xb, yb, massa, massb); 

    return Py_BuildValue("ddddd", Vxa, Vya, Vxb, Vyb, distance); 
} 

static PyMethodDef GForcesMethods[] = { 
{"gforces", gforces, METH_VARARGS, "Calculate the x and y acceleration of two masses and the distance between the two."}, 
{NULL, NULL, 0, NULL} 
}; 

PyMODINIT_FUNC 
initgforces(void) 
{ 
(void) Py_InitModule("gforces", GForcesMethods); 
} 

Se compilate come file PYD si dovrebbe ottenere un oggetto Python che è possibile importare. Devi però avere tutte le opzioni del compilatore e del linker corrette. Sto usando dev-C++ e le mie opzioni del compilatore sono impostate su -shared -o gforces.pyd e linker sono impostate su -lpython27 (assicurati di usare la stessa versione che hai installato) e aggiungi il tuo percorso di directory python agli include e alle librerie schede.

L'oggetto accetta gli argomenti (p1.speedx, p1.speedy, p2.speedx, p2.speedy, p1.x, p1.y, p2.x, p2.y, p1.mass, p2.mass) e restituisce il nuovo p1.speedx, p1.speedy, p2.speedx, p2.speedy e la distanza tra p1 p2.

Utilizzando il modulo di cui sopra, ho anche cercato di tagliare fuori un paio di passi per il rilevamento delle collisioni confrontando la distanza restituita alla somma dei raggi delle particelle in quanto tali:

def updateForces(self):   #part of a handler class for the particles 
    prevDone = [] 
    for i in self.planetGroup: #planetGroup is a sprite group from pygame 
     prevDone.append(i.ID) 
     for j in self.planetGroup: 
      if not (j.ID in prevDone):    #my particles have unique IDs 
       distance = i.calcGForce(j)   #calcGForce just calls the above 
       if distance <= i.radius + j.radius: #object and assigns the returned 
        #collision handling goes here #values for the x and y speed 
                #components to the particles 

Spero che questo aiuti un po. Qualsiasi ulteriore consiglio o segnalazione di errori grossolani da parte mia è il benvenuto, anch'io sono nuovo a questo.