当你用数字电网运行,并希望真正好的表现,你应该考虑使用Numpy。令人惊讶的是,使用起来非常简单,并且可以让您根据网格运行而不是网格上的循环进行思考。性能来自这样的事实,即操作随后通过优化的SSE代码在整个网格上运行。
例如下面是一些使用我编写的代码,可以对由弹簧连接的带电粒子进行蛮力数值模拟的代码。该代码计算具有100个节点的3d系统的时间步长和31ms内的99个边沿。这比我能想出的最好的纯Python代码快10倍以上。
from numpy import array, sqrt, float32, newaxis
def evolve(points, velocities, edges, timestep=0.01, charge=0.1, mass=1., edgelen=0.5, dampen=0.95):
"""Evolve a n body system of electrostatically repulsive nodes connected by
springs by one timestep."""
velocities *= dampen
# calculate matrix of distance vectors between all points and their lengths squared
dists = array([[p2 - p1 for p2 in points] for p1 in points])
l_2 = (dists*dists).sum(axis=2)
# make the diagonal 1's to avoid division by zero
for i in xrange(points.shape[0]):
l_2[i,i] = 1
l_2_inv = 1/l_2
l_3_inv = l_2_inv*sqrt(l_2_inv)
# repulsive force: distance vectors divided by length cubed, summed and multiplied by scale
scale = timestep*charge*charge/mass
velocities -= scale*(l_3_inv[:,:,newaxis].repeat(points.shape[1], axis=2)*dists).sum(axis=1)
# calculate spring contributions for each point
for idx, (point, outedges) in enumerate(izip(points, edges)):
edgevecs = point - points.take(outedges, axis=0)
edgevec_lens = sqrt((edgevecs*edgevecs).sum(axis=1))
scale = timestep/mass
velocities[idx] += (edgevecs*((((edgelen*scale)/edgevec_lens - scale))[:,newaxis].repeat(points.shape[1],axis=1))).sum(axis=0)
# move points to new positions
points += velocities*timestep
这不是一回事,在第二个范围内取了alist [i]的len,为什么你删除了那个索引? – 2008-10-09 21:08:09
我绝对喜欢这个。这可能不是最好的答案,但我很可能会使用它。 – 2008-10-09 21:57:33
这只是Eugene写的东西。如果电网足够大,恐怕这台发电机可能会很慢。毕竟,每个收益率仍然需要四个索引查找。 – DzinX 2008-10-09 22:14:10