在numpy中執行許多手段

早上好，我在實現Numpy中的距離加權平均值時使用Cressman過濾器..我使用Ball Tree implimentation（感謝Jake VanderPlas）返回一個位於每個點的位置列表請求數組..查詢數組（q）是形狀[n，3]，並且在每個點處都有x，y，z點，我想要對存儲在樹中的點進行加權平均。樹返回一定距離內的點，所以我得到一個可變長度數組的數組..我使用一個在哪裏找到非空的條目（即在影響半徑內至少有一些點的位置）創建isgood數組。 ..在numpy中執行許多手段

然後我遍歷所有查詢點以返回加權平均值self.z（請注意，這可以是dims = 1或dims = 2以允許多個協同網格）

因此，使用map或其他更快速的方法完成的事情是長度的非均勻性self.distances和self.locations內的陣列......我仍然相當綠色numpy的/ Python的，但我不能想辦法來明智做陣列（即不恢復到環路）

self.locations, self.distances = self.tree.query_radius(q, r, return_distance=True) 
t2=time() 
if debug: print "Removing voids" 
isgood=np.where(np.array([len(x) for x in self.locations])!=0)[0] 
interpol = np.zeros((len(self.locations),) + np.shape(self.z[0])) 
interpol.fill(np.nan) 
for dist, ix, posn, roi in zip(self.distances[isgood], self.locations[isgood], isgood, r[isgood]): 
    interpol[isgood[jinterpol]] = np.average(self.z[ix], weights=(roi**2-dist**2)/(roi**2 + dist**2), axis=0) 
    jinterpol += 1

所以...任何提示如何加快循環？

對於一個典型的映射，適用於從一個範圍，方位角，高程網格映射到天氣雷達數據的cartesi一個網格，我有240x240x34點和4個變量需要99s查詢樹（由傑克在C和cython書面..這是艱難的一步，因爲你需要搜索數據！）和100秒來做計算...在我看來，這是慢？我的開銷在哪裏？是np.mean高效率還是因爲它被稱爲數百萬次在這裏獲得加速嗎？我會通過使用float32而不是默認的64 ...或甚至縮放到整數（這將很難避免在加權中纏繞...感謝任何提示！

來源

2011-06-03 Scott Collis

您是否試圖在您的代碼上運行分析器？嘗試在循環外預先計算roi ** 2和dist ** 2，並使用乘法（r [isgood] * r [isgood]）而不是**。使用np.empty而不是np.zeros。 – 2011-06-03 15:44:43

只是簡單的時間測試time.time（）和timeit ... – 2011-06-03 16:45:26

你可以找到關於所述克雷斯曼方案的相對優點VS在使用高斯加權函數：

http://www.flame.org/~cdoswell/publications/radar_oa_00.pdf

的關鍵是匹配平滑參數的數據（我建議使用接近的數據點之間的平均間距的值）。一旦你知道了平滑參數，你可以設置一個「影響半徑」等於權重函數下降到0.01（或其他）的半徑。

速度有多重要？如果你願意，而不是調用指數函數來確定重量，你可以爲一些固定數量的半徑增量組成一個離散的權重表，這大大加快了計算速度。理想情況下，您應該有網格邊界外的數據，可用於映射網格點周圍的值（即使在網格的邊界點上）。注意這不是一個真正的插值方案 - 它不會精確地返回數據點處的觀測值。像Cressman計劃一樣，它是一個低通濾波器。

來源

2011-06-03 15:37:15

是的..我已經引用了你的OBAN（DWA vrs線性，即Vaghan和Mohr）論文:) ...改變巴恩斯是一個改變權重的函數形式......編寫這段代碼的要點是REORDER的靈活性......最後，我希望將它與一個校準雷達網絡一起使用......這將與任意分散的點和任意值一組查詢點..例如在一個網站上獲取列是非常快速和容易... – 2011-06-03 16:12:29

在numpy中執行許多手段

回答

相關問題