न्यूप्टीयन बलों की अधिकांश प्रदर्शन गणना numpy/scipy

Question

विश्वविद्यालय में एक अभ्यास के लिए, हमें पाइथन में सटीक न्यूटनियन बलों के साथ Leapfrog integrator लागू करना पड़ा। पाठ्यक्रम खत्म हो गया है और हमारे समाधान काफी अच्छे थे, लेकिन मुझे आश्चर्य हुआ कि अगर कोई बल गणना के प्रदर्शन में और भी सुधार कर सकता है।

टोंटी (त्वरण उर्फ) सभी बलों की गणना करने के लिए है:

एक बड़े के लिए

(1000 और बड़ा) कणों एन (i, j < एन) की संख्या।

यहाँ r और आर कणों की स्थिति आकार की एक ndarray में संग्रहीत की 3 आयामी वैक्टर हैं (एन, 3) और ग्राम कण द्रव्यमान बार गुरुत्वाकर्षण स्थिरांक जो मैं बचाया है आकार (एन) के एक अंडाकार में।

def a(self): 
    sep = self.r[np.newaxis, :] - self.r[:, np.newaxis] 
    # Calculate the distances between all particles with cdist 
    # this is much faster than by hand 
    dists = cdist(self.r, self.r) 
    scale =dists*dists*dists 
    # set diagonal elements of dist to something != 0, to avoid division by 0 
    np.fill_diagonal(scale,1) 
    Fsum = (sep/scale.reshape(self.particlenr,self.particlenr,1))*self.Gm[:,None] 
    return np.add.reduce(Fsum,axis=1) 

यह कीड़े मुझे कि यह शायद सबसे तेजी से संस्करण नहीं है लेकिन:

सबसे तेजी से संस्करण मैं अब तक पाया निम्नलिखित है। सीडीआईटी की तुलना करते समय पहली पंक्ति बहुत धीमी प्रतीत होती है जो अनिवार्य रूप से वही गणना कर रही है। साथ ही, यह समाधान समस्या में आर और आर स्विचिंग की समरूपता का उपयोग नहीं करता है और सभी तत्वों की गणना दो बार करता है।

क्या आप कोई प्रदर्शन सुधार जानते हैं (कुछ अनुमान के लिए बल-गणना को बदलने या प्रोग्रामिंग भाषा को बदलने के बिना)?

Answer 1

मुझे संदेह है कि numpy वास्तव में दूरी को कंप्यूटिंग कर रहा है (क्योंकि यह हमेशा सममित होगा)। यह शायद एक गणना कर रहा है और दो स्थानों में एक ही मूल्य असाइन कर रहा है।

विचारों में से कुछ, मेरे लिए होते हैं, हालांकि:

1. आप numpy स्रोत का पालन करें और cdist का अनुकूलित संस्करण लिख सकते हैं। ऐसा हो सकता है कि कार्यक्रम प्रत्येक बार कई बार को पार्स कर रहा है। ज्यादा नहीं, लेकिन शायद यह आपको एक छोटा प्रतिशत टक्कर दे सकता है।
2. प्री-आवंटित करें। प्रत्येक बार जब आप एक() चलाते हैं, तो आप संभावित रूप से सभी मध्यवर्ती मैट्रिक्स मानों के लिए स्मृति को पुन: आवंटित कर रहे हैं। क्या आप इन लगातार मात्रा बना सकते हैं?

मैंने अभी तक गणनाओं के माध्यम से भाग नहीं लिया है, लेकिन अगर आप अनावश्यक रूप से अनावश्यक सममित गणना को कम कर सकते हैं तो मुझे आश्चर्य नहीं होगा।

Answer 2

मैं इसे आजमाइए हूँ: मैं एक दिनचर्या, जो निर्धारित करता है कार्यान्वित एक भी a_i:

import numpy as np 

GM = .01 # article mass times the gravitation 

def calc_a_i(rr, i): 
    """ Calculate one a_i """ 
    drr = rr - rr[i, :] # r_j - r_i 
    dr3 = np.linalg.norm(drr, axis=1)**3 # |r_j - r_i|**3 
    dr3[i] = 1 # case i==j: drr = [0, 0, 0] 
    # this would be more robust (elimnate small denominators): 
    #dr3 = np.where(np.abs(dr3) > 1e-12, dr3, 1) 
    return np.sum(drr.T/dr3, axis=1) 

n = 4000 # number of particles 
rr = np.random.randn(n, 3) # generate some particles 

# Calculate each a_i separately: 
aa = np.array([calc_a_i(rr, i) for i in range(n)]) * GM # all a_i 

यह परीक्षण करने के लिए, मैं भाग गया:

In [1]: %timeit aa = np.array([calc_a_i(rr, i) for i in range(n)]) 
1 loops, best of 3: 2.93 s per loop 

तेजी लाने के लिए सबसे आसान तरीका

import numexpr as ne 
ne.set_num_threads(1) # multithreading causes to much overhead 

def ne_calc_a_i(i): 
    """ Use numexpr - here rr is global for easier parallelization""" 
    dr1, dr2, dr3 = (rr - rr[i, :]).T # r_j - r_i 
    drrp3 = ne.evaluate("sqrt(dr1**2 + dr2**2 + dr3**2)**3") 
    drrp3[i] = 1 
    return np.sum(np.vstack([dr1, dr2, dr3])/drrp3, axis=1) 

# Calculate each a_i separately: 
aa_ne = np.array([ne_calc_a_i(i) for i in range(n)]) * GM # all a_i  

थी: इस तरह के कोड, सरणी भाव के तेजी से मूल्यांकन के लिए numexpr उपयोग करने के लिए है रों 2 का एक पहलू से गति को बेहतर बनाता है:

In [2]: %timeit aa_ne = np.array([ne_calc_a_i(i) for i in range(n)]) 
    1 loops, best of 3: 1.29 s per loop 

आगे कोड में तेजी लाने के लिए, के लिए यह एक IPython Cluster पर चलाते हैं:

# Start local cluster with 4 clients in a shell with: 
# ipcluster start -n 4 

rc = Client() # clients of cluster 
dview = rc[:] # view of clusters 

dview.execute("import numpy as np") # import libraries on clients 
dview.execute("import numexpr as ne") 
dview.execute("ne.set_num_threads(1)") 

def para_calc_a(dview, rr): 
    """ Only in function for %timeit """ 
    # send rr and ne_calc_a_i() to clients: 
    dview.push(dict(rr=rr, ne_calc_a_i=ne_calc_a_i), block=True) 
    return np.array(dview.map_sync(ne_calc_a_i, range(n)))*GM 

speedup चार का एक पहलू से भी अधिक है:

यह अपने CPU वास्तुकला पर निर्भर करता है कि स्मृति बस या चल बिन्दु इकाई सीमित कारक है, जो च कहता है:

In[3] %timeit aa_p = para_calc_a(dview, rr) 
1 loops, best of 3: 612 ms per loop 

@mathdan पहले से ही बताया गया है, यह कैसे ऐसी समस्या का अनुकूलन करने के स्पष्ट नहीं है या विभिन्न तकनीकें।

अधिक लाभ के लिए आप Theano पर देख सकते हैं: यह गतिशील रूप से पाइथन से कोड जीपीयू कोड उत्पन्न कर सकता है।

Answer 3

निम्नलिखित में थोड़ा और अधिक इष्टतम है:

import numpy as np 
from scipy.spatial.distance import pdist, squareform  

def a6(r, Gm): 
    dists = pdist(r) 
    dists *= dists*dists 
    dists = squareform(dists) 
    np.fill_diagonal(dists, 1.) 
    sep = r[np.newaxis, :] - r[:, np.newaxis] 
    return np.einsum('ijk,ij->ik', sep, Gm/dists) 

गति हासिल einsum लाइन की वजह से मुख्य रूप से है, pdist और squareform का उपयोग करके cdist के साथ मूल तरीके से केवल मामूली तेज़ है।

आप इसे थोड़ा और आगे ले सकते हैं, उदा। थ्रेडिंग और नुंबा (संस्करण 0.17.0 आवश्यक) का उपयोग करके। हालांकि नीचे दिया गया कोड बहुत बदसूरत है और निश्चित रूप से बहुत सुधार किया जा सकता है, यह बहुत तेज़ है।

import numpy as np 
import math 
from numba import jit 
from threading import Thread 
NUM_THREADS = 2 # choose wisely 

def a_numba_par(r, Gm): 
    a = np.zeros_like(r) 
    N = r.shape[0] 

    offset = range(0, N+1, N//NUM_THREADS) 
    chunks = zip(offset, offset[1:]) 
    threads = [Thread(target=_numba_loop, args=(r,Gm,a)+c) for c in chunks] 

    for thread in threads: 
     thread.start() 
    for thread in threads: 
     thread.join() 

    return a 

@jit(nopython=True, nogil=True) 
def _numba_loop(r, Gm, a, i1, i2): 
    N = r.shape[0] 
    for i in range(i1, i2): 
     _helper(r, Gm, i, 0 , i, a[i,:]) 
     _helper(r, Gm, i, i+1, N, a[i,:]) 
    return a 

@jit(nopython=True, nogil=True) 
def _helper(r, Gm, i, j1, j2, a): 
    for j in range(j1, j2): 
     dx = r[j,0] - r[i,0] 
     dy = r[j,1] - r[i,1] 
     dz = r[j,2] - r[i,2] 

     sqeuc = dx*dx + dy*dy + dz*dz 
     scale = Gm[j]/(sqeuc * math.sqrt(sqeuc)) 

     a[0] += scale * dx 
     a[1] += scale * dy 
     a[2] += scale * dz