编写一个更快的Python物理模拟器

16
我一直在使用Python编写自己的物理引擎,作为物理和编程练习。我最初是按照这里的教程进行的。那个很顺利,但后来我发现了thomas jakobsen撰写的“高级角色物理”文章,其中涵盖了使用Verlet积分进行模拟,我觉得很有趣。
我一直在尝试使用Verlet积分编写自己的基本物理模拟器,但事实证明它比我最初预期的要稍微困难一些。我正在浏览示例程序以阅读,并偶然发现了这个用Python编写的程序,还发现了使用Processing的这个教程
令我印象深刻的是Processing版本的运行速度非常快。仅布料就有2400个不同的点在进行模拟,这还不包括物体。
Python示例仅使用256个粒子用于布料,并且以大约30帧每秒的速度运行。我尝试将粒子数增加到2401(该程序必须为正方形),它以大约3帧每秒的速度运行。
这两个示例均通过在列表中存储一个粒子对象的实例,然后迭代该列表,调用每个粒子的“更新位置”方法来工作。例如,这是Processing示例代码的一部分,用于计算每个粒子的新位置:
for (int i = 0; i < pointmasses.size(); i++) {
    PointMass pointmass = (PointMass) pointmasses.get(i);
    pointmass.updateInteractions();
    pointmass.updatePhysics(fixedDeltaTimeSeconds);
}

编辑:以下是我之前提供的 Python 版本代码:

"""
verletCloth01.py
Eric Pavey - 2010-07-03 - www.akeric.com

Riding on the shoulders of giants.
I wanted to learn now to do 'verlet cloth' in Python\Pygame.  I first ran across
this post \ source:
http://forums.overclockers.com.au/showthread.php?t=870396
http://dl.dropbox.com/u/3240460/cloth5.py

Which pointed to some good reference, that was a dead link.  After some searching,
I found it here:
http://www.gpgstudy.com/gpgiki/GDC%202001%3A%20Advanced%20Character%20Physics
Which is a 2001 SIGGRAPH paper by Thomas Jakobsen called:
"GDC 2001: Advanced Characer Physics".

This code is a Python\Pygame interpretation of that 2001 Siggraph paper.  I did
borrow some code from 'domlebo's source code, it was a great starting point.  But
I'd like to think I put my own flavor on it.
"""

#--------------
# Imports & Initis
import sys
from math import sqrt

# Vec2D comes from here: http://pygame.org/wiki/2DVectorClass
from vec2d import Vec2d
import pygame
from pygame.locals import *
pygame.init()

#--------------
# Constants
TITLE = "verletCloth01"
WIDTH = 600
HEIGHT = 600
FRAMERATE = 60
# How many iterations to run on our constraints per frame?
# This will 'tighten' the cloth, but slow the sim.
ITERATE = 2
GRAVITY = Vec2d(0.0,0.05)
TSTEP = 2.8

# How many pixels to position between each particle?
PSTEP = int(WIDTH*.03)
# Offset in pixels from the top left of screen to position grid:
OFFSET = int(.25*WIDTH)

#-------------
# Define helper functions, classes

class Particle(object):
    """
    Stores position, previous position, and where it is in the grid.
    """
    def __init__(self, screen, currentPos, gridIndex):
        # Current Position : m_x
        self.currentPos = Vec2d(currentPos)
        # Index [x][y] of Where it lives in the grid
        self.gridIndex = gridIndex
        # Previous Position : m_oldx
        self.oldPos = Vec2d(currentPos)
        # Force accumulators : m_a
        self.forces = GRAVITY
        # Should the particle be locked at its current position?
        self.locked = False
        self.followMouse = False

        self.colorUnlocked = Color('white')
        self.colorLocked = Color('green')
        self.screen = screen

    def __str__(self):
        return "Particle <%s, %s>"%(self.gridIndex[0], self.gridIndex[1])

    def draw(self):
        # Draw a circle at the given Particle.
        screenPos = (self.currentPos[0], self.currentPos[1])
        if self.locked:
            pygame.draw.circle(self.screen, self.colorLocked, (int(screenPos[0]),
                                                         int(screenPos[1])), 4, 0)
        else:
            pygame.draw.circle(self.screen, self.colorUnlocked, (int(screenPos[0]),
                                                         int(screenPos[1])), 1, 0)

class Constraint(object):
    """
    Stores 'constraint' data between two Particle objects.  Stores this data
    before the sim runs, to speed sim and draw operations.
    """
    def __init__(self, screen, particles):
        self.particles = sorted(particles)
        # Calculate restlength as the initial distance between the two particles:
        self.restLength = sqrt(abs(pow(self.particles[1].currentPos.x -
                                       self.particles[0].currentPos.x, 2) +
                                   pow(self.particles[1].currentPos.y -
                                       self.particles[0].currentPos.y, 2)))
        self.screen = screen
        self.color = Color('red')

    def __str__(self):
        return "Constraint <%s, %s>"%(self.particles[0], self.particles[1])

    def draw(self):
        # Draw line between the two particles.
        p1 = self.particles[0]
        p2 = self.particles[1]
        p1pos = (p1.currentPos[0],
                 p1.currentPos[1])
        p2pos = (p2.currentPos[0],
                 p2.currentPos[1])
        pygame.draw.aaline(self.screen, self.color,
                           (p1pos[0], p1pos[1]), (p2pos[0], p2pos[1]), 1)

class Grid(object):
    """
    Stores a grid of Particle objects.  Emulates a 2d container object.  Particle
    objects can be indexed by position:
        grid = Grid()
        particle = g[2][4]
    """
    def __init__(self, screen, rows, columns, step, offset):

        self.screen = screen
        self.rows = rows
        self.columns = columns
        self.step = step
        self.offset = offset

        # Make our internal grid:
        # _grid is a list of sublists.
        #    Each sublist is a 'column'.
        #        Each column holds a particle object per row:
        # _grid =
        # [[p00, [p10, [etc,
        #   p01,  p11,
        #   etc], etc],     ]]
        self._grid = []
        for x in range(columns):
            self._grid.append([])
            for y in range(rows):
                currentPos = (x*self.step+self.offset, y*self.step+self.offset)
                self._grid[x].append(Particle(self.screen, currentPos, (x,y)))

    def getNeighbors(self, gridIndex):
        """
        return a list of all neighbor particles to the particle at the given gridIndex:

        gridIndex = [x,x] : The particle index we're polling
        """
        possNeighbors = []
        possNeighbors.append([gridIndex[0]-1, gridIndex[1]])
        possNeighbors.append([gridIndex[0], gridIndex[1]-1])
        possNeighbors.append([gridIndex[0]+1, gridIndex[1]])
        possNeighbors.append([gridIndex[0], gridIndex[1]+1])

        neigh = []
        for coord in possNeighbors:
            if (coord[0] < 0) | (coord[0] > self.rows-1):
                pass
            elif (coord[1] < 0) | (coord[1] > self.columns-1):
                pass
            else:
                neigh.append(coord)

        finalNeighbors = []
        for point in neigh:
            finalNeighbors.append((point[0], point[1]))

        return finalNeighbors

    #--------------------------
    # Implement Container Type:

    def __len__(self):
        return len(self.rows * self.columns)

    def __getitem__(self, key):
        return self._grid[key]

    def __setitem__(self, key, value):
        self._grid[key] = value

    #def __delitem__(self, key):
        #del(self._grid[key])

    def __iter__(self):
        for x in self._grid:
            for y in x:
                yield y

    def __contains__(self, item):
        for x in self._grid:
            for y in x:
                if y is item:
                    return True
        return False


class ParticleSystem(Grid):
    """
    Implements the verlet particles physics on the encapsulated Grid object.
    """

    def __init__(self, screen, rows=49, columns=49, step=PSTEP, offset=OFFSET):
        super(ParticleSystem, self).__init__(screen, rows, columns, step, offset)

        # Generate our list of Constraint objects.  One is generated between
        # every particle connection.
        self.constraints = []
        for p in self:
            neighborIndices = self.getNeighbors(p.gridIndex)
            for ni in neighborIndices:
                # Get the neighbor Particle from the index:
                n = self[ni[0]][ni[1]]
                # Let's not add duplicate Constraints, which would be easy to do!
                new = True
                for con in self.constraints:
                    if n in con.particles and p in con.particles:
                        new = False
                if new:
                    self.constraints.append( Constraint(self.screen, (p,n)) )

        # Lock our top left and right particles by default:
        self[0][0].locked = True
        self[1][0].locked = True
        self[-2][0].locked = True
        self[-1][0].locked = True

    def verlet(self):
        # Verlet integration step:
        for p in self:
            if not p.locked:
                # make a copy of our current position
                temp = Vec2d(p.currentPos)
                p.currentPos += p.currentPos - p.oldPos + p.forces * TSTEP**2
                p.oldPos = temp
            elif p.followMouse:
                temp = Vec2d(p.currentPos)
                p.currentPos = Vec2d(pygame.mouse.get_pos())
                p.oldPos = temp

    def satisfyConstraints(self):
        # Keep particles together:
        for c in self.constraints:
            delta =  c.particles[0].currentPos - c.particles[1].currentPos
            deltaLength = sqrt(delta.dot(delta))
            try:
                # You can get a ZeroDivisionError here once, so let's catch it.
                # I think it's when particles sit on top of one another due to
                # being locked.
                diff = (deltaLength-c.restLength)/deltaLength
                if not c.particles[0].locked:
                    c.particles[0].currentPos -= delta*0.5*diff
                if not c.particles[1].locked:
                    c.particles[1].currentPos += delta*0.5*diff
            except ZeroDivisionError:
                pass

    def accumulateForces(self):
        # This doesn't do much right now, other than constantly reset the
        # particles 'forces' to be 'gravity'.  But this is where you'd implement
        # other things, like drag, wind, etc.
        for p in self:
            p.forces = GRAVITY

    def timeStep(self):
        # This executes the whole shebang:
        self.accumulateForces()
        self.verlet()
        for i in range(ITERATE):
            self.satisfyConstraints()

    def draw(self):
        """
        Draw constraint connections, and particle positions:
        """
        for c in self.constraints:
            c.draw()
        #for p in self:
        #    p.draw()

    def lockParticle(self):
        """
        If the mouse LMB is pressed for the first time on a particle, the particle
        will assume the mouse motion.  When it is pressed again, it will lock
        the particle in space.
        """
        mousePos = Vec2d(pygame.mouse.get_pos())
        for p in self:
            dist2mouse = sqrt(abs(pow(p.currentPos.x -
                                      mousePos.x, 2) +
                                  pow(p.currentPos.y -
                                      mousePos.y, 2)))
            if dist2mouse < 10:
                if not p.followMouse:
                    p.locked = True
                    p.followMouse = True
                    p.oldPos = Vec2d(p.currentPos)
                else:
                    p.followMouse = False

    def unlockParticle(self):
        """
        If the RMB is pressed on a particle, if the particle is currently
        locked or being moved by the mouse, it will be 'unlocked'/stop following
        the mouse.
        """
        mousePos = Vec2d(pygame.mouse.get_pos())
        for p in self:
            dist2mouse = sqrt(abs(pow(p.currentPos.x -
                                      mousePos.x, 2) +
                                  pow(p.currentPos.y -
                                      mousePos.y, 2)))
            if dist2mouse < 5:
                p.locked = False

#------------
# Main Program
def main():
    # Screen Setup
    screen = pygame.display.set_mode((WIDTH, HEIGHT))
    clock = pygame.time.Clock()

    # Create our grid of particles:
    particleSystem = ParticleSystem(screen)
    backgroundCol = Color('black')

    # main loop
    looping = True
    while looping:
        clock.tick(FRAMERATE)
        pygame.display.set_caption("%s -- www.AKEric.com -- LMB: move\lock - RMB: unlock - fps: %.2f"%(TITLE, clock.get_fps()) )
        screen.fill(backgroundCol)

        # Detect for events
        for event in pygame.event.get():
            if event.type == pygame.QUIT:
                looping = False
            elif event.type == MOUSEBUTTONDOWN:
                if event.button == 1:
                    # See if we can make a particle follow the mouse and lock
                    # its position when done.
                    particleSystem.lockParticle()
                if event.button == 3:
                    # Try to unlock the current particles position:
                    particleSystem.unlockParticle()

        # Do stuff!
        particleSystem.timeStep()
        particleSystem.draw()

        # update our display:
        pygame.display.update()

#------------
# Execution from shell\icon:
if __name__ == "__main__":
    print "Running Python version:", sys.version
    print "Running PyGame version:", pygame.ver
    print "Running %s.py"%TITLE
    sys.exit(main())

由于这两个程序的工作方式大致相同,但Python版本要慢得多,这让我想知道:
- 这种性能差异是Python本质的一部分吗? - 如果我想从自己的Python程序中获得更好的性能,我应该如何与上述不同?例如,将所有粒子的属性存储在数组中而不是使用单个对象等。
编辑:已回答!!
@Mr E在评论中提供的PyCon讲座和@A. Rosa的答案以及链接的资源都极大地帮助了我更好地理解如何编写高效的Python代码。我现在将此页面加为书签以备将来参考:D

2
一个普遍的观点。有一个关于过度使用类的不错的Pycon视频。演讲者不断指出一些只有“两个方法,其中一个是__init__”的类的例子,并说它们最好表示为方法(这里忽略了__str__函数)。你可以很容易地用namedtuple替换你的粒子,或者使用一个draw_particle函数。 - YXD
哦,我也看到这不是你的代码,所以可能不相关... - YXD
@MrE我喜欢那个演示文稿!但我发现人们不太喜欢被指向它,这种反应让我想起了这个链接(http://abstrusegoose.com/381)。 - Jaime
哈哈!是的,你必须仔细挑选时机。 - YXD
4个回答

8
在Python Wiki的性能提示部分中,有一篇Guido van Rossum的文章。在结论中,您可以阅读以下句子:“如果您需要速度,请选择内置函数-无法击败用C编写的循环。”
该文章继续列出了有关循环优化的指南。我推荐这两个资源,因为它们提供了有关优化Python代码的具体和实用建议。
还有一个众所周知的基准测试组benchmarksgame.alioth.debian.org,在其中您可以找到不同程序和语言在不同机器上的比较。如可见,有许多变量在起作用,使得不可能声明像“Java比Python更快”这样广泛的内容。这通常总结为“语言没有速度;实现才有”。
在您的代码中可以使用内置函数应用更Pythonic和更快的替代方案。例如,有几个嵌套循环(其中一些不需要处理整个列表),可以使用imaplist comprehensions进行重写。PyPy也是另一个提高性能的有趣选择。我不是Python优化方面的专家,但有很多非常有用的技巧(请注意不要用Java写Python就是其中之一!)。 资源和其他相关问题,请参见 Stack Overflow:

5
如果你写Python的方式和你写Java的方式一样,那么显然它会变得更慢,惯用的Java代码翻译成惯用的Python代码可能效果会不尽人意。

Python的性能差异是否与生俱来? 如果我想从自己的Python程序中获得更好的性能,我应该如何做呢?例如,将所有粒子的属性存储在数组中而不是使用单独的对象等。

没有看到你的代码很难说。

以下是Python和Java之间的一些不完整的差异列表,这些差异有时会影响性能:

  1. 处理使用立即模式画布,如果你想在Python中获得可比的性能,你也需要使用立即模式画布。 大多数GUI框架(包括Tkinter画布)都是保留模式,这样更容易使用,但固有地比立即模式慢。 你需要使用提供的pygame,SDL或Pyglet之类的立即模式画布。

  2. Python是动态语言,这意味着实例成员访问,模块成员访问和全局变量访问在运行时解析。 在Python中,实例成员访问,模块成员访问和全局变量访问实际上是字典访问。 在Java中,它们在编译时解析,并且本质上更快。 将经常访问的全局变量,模块变量和属性缓存到本地变量中。

  3. 在Python 2.x中,range()生成一个具体的列表,在Python中,使用迭代器进行迭代,for item in list通常比使用迭代变量for n in range(len(list))进行迭代要快。 你几乎总是应该直接使用迭代器进行迭代,而不是使用range(len(...))进行迭代。

  4. Python的数字是不可变的,这意味着任何算术计算都会分配一个新对象。 这是纯Python不太适合低级别计算的原因之一; 大多数想要能够编写低级别计算而不必使用C扩展的人通常使用cython、psyco或numpy。 但这通常只有在有数百万个计算时才会成为问题。

这只是部分的,非常不完整的列表,将Java翻译成Python可能会产生次优代码的原因有很多其他。 没有看到你的代码,无法确定你需要做什么不同。 优化的Python代码通常看起来与优化的Java代码非常不同。


我已经添加了来自链接的Python程序代码,因为它几乎与我的代码完全相同,除了它实际上是可以工作的。 - mooglinux

5
我建议你阅读其他物理引擎的相关内容。有一些开源引擎使用各种方法计算“物理”。
- Newton Game Dynamics - Chipmunk - Bullet - Box2D - ODE(Open Dynamics Engine)
这些引擎大多都有移植版本:
- Pymunk - PyBullet - PyBox2D - PyODE
如果你阅读这些引擎的文档,你会发现它们通常会声明自己针对速度进行了优化(30fps - 60fps)。但是,如果你认为它们能够在计算“真实”物理时做到这一点,那么你就错了。大多数引擎计算物理的程度已经足够让普通用户无法区分“真实”物理行为和“模拟”的物理行为。然而,如果你仔细研究误差,你会发现在编写游戏时可以忽略不计。但是,如果你想做物理模拟,所有这些引擎对你来说都没有用处。
这就是为什么我要说,如果你正在进行真正的物理模拟,你的速度肯定比这些引擎慢,并且你永远也无法超过另一个物理引擎。

0

基于粒子的物理模拟可以轻松地转化为线性代数运算,即矩阵运算。Numpy提供了这样的操作,其在Fortran/C/C++下实现。良好编写的Python/Numpy代码(充分利用语言和库)可以编写出相当快速的代码。


网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接