对于克里金法,我需要计算长度为20000
的大型mesh
数组。下面的代码运行良好,特别是对于较小的mesh
长度(<;100),但是,对于如此大的mesh
的计算时间非常长(约45分钟)。data
的长度在300到800之间。以下是工作代码:
import numpy as np
from scipy.spatial.distance import pdist, squareform
icovfct = [36, 6524.62, 1383.13, 2]
imesh = np.array([[632230, 632090, 632110, 632130, 632150, 632170, 632190, 632210, 632230, 632250, 632270, 632290, 632310, 632070],
[3045160, 3045180, 3045180, 3045180, 3045180, 3045180, 3045180, 3045180, 3045180, 3045180, 3045180, 3045180, 3045180, 3045200]], np.float64)
imesh = imesh.T
idata = np.array([[634026.049, 633901.182, 634001.365, 634007.132, 633893.706, 633802.327, 634144.246, 634015.993, 633897.326, 633779.479],
[3048117.579, 3048201.031, 3048191.922, 3047891.355, 3047994.462, 3048084.562, 3047633.421, 3047719.845, 3047818.914, 3047902.179],
[256.550, 236.317, 249.458, 281.889, 262.321, 239.495, 303.144, 295.319, 281.270, 261.083]], np.float64)
idata = idata.T
def ordinary(covfct, data, mesh):
prediction = []
for i, dummy_val in enumerate(mesh):
# distance between data and each data point in mesh
d = np.sqrt((data[:, 0]-mesh[i, 0])**2.0 + (data[:, 1]-mesh[i, 1])**2.0)
# add these distances to P
P = np.vstack((data.T, d)).T
# apply the covariance model to the distances
k = (covfct[0] + covfct[1]*(1-np.exp(-3*(P[:,3]/covfct[2])**covfct[3])))
# cast as a matrix
k = np.matrix(k).T
# form a matrix of distances between existing data points
K1 = squareform(pdist(P[:,:2]))
# apply the covariance model to these distances
K = (covfct[0] + covfct[1]*(1-np.exp(-3*(K1.ravel()/covfct[2])**covfct[3])))
# re-cast as a NumPy array
K = np.array(K)
# reshape into an array
K = K.reshape(len(P), len(P))
# if kriging is exact K diag = 0
K[K1 == 0] = 0
# cast as a matrix
K = np.matrix(K)
# add a column and row of ones to Ks,
# with a zero in the bottom, right hand corner (Lagrangian)
K2 = np.matrix(np.ones((len(P)+1, len(P)+1)))
K2[:len(P), :len(P)] = K
K2[-1, -1] = 0.0
# add a one to the end of ks
k3 = np.matrix(np.ones((len(P)+1, 1)))
k3[:len(P)] = k
# calculate the kriging weights
weights = np.linalg.lstsq(K2, k3, rcond = 1)
weights = weights[:-3][0][:-1]
weights = np.array(weights)
# calculate the residuals
residuals = P[:, 2]
# calculate the estimation
prediction.append(np.dot(weights.T, residuals))
return prediction
interpolate = ordinary(icovfct, idata, imesh)
有没有办法优化代码从而减少计算时间
您可以使用GSTools作为代码的替代品,其中克里格求和是在Cython中实现的:
通常,您可以通过将一个进程拆分为多个进程来减少时间
在python中,我们有python process lib可以做到这一点, 因此,如果长度为20000,则可以同时运行代码两次 每道工序处理长度为10000
我希望你这个答案能解决你的问题
相关问题 更多 >
编程相关推荐