Cannot reproduce simulation with same space and actions (impulse)

[ShaoxiongYao]

I am using Pymunk version 5.6.0 on ubunutu 18.04.
I write a simulation environment with a set of dynamics balls connected by damped springs. I apply a series of actions to these balls in the space and I recorded the actions and velocity changes of these balls. I also store the initial space variable into a pickle. However, when I load the space and apply the same actions to balls in the space, the velocity changes are different. The difference in the L2 norm is about 0.02 for 15 balls. But when I apply actions for 2000 steps, the error will accumulate to 10-20. I am wondering whether it is correct to expect the same change of velocities. Also, I am thinking whether space itself contains all the information I need to reproduce the result.

[viblo]

There are some hidden state in the underlying Chipmunk library which Pymunk uses for the physics calculations. Basically it has to do with caching impulses for resolving contact between steps.

So, if you create a new space, add some things and run simulation I think it should behave the same. If you dont rewind all the way to the beginning the cached impulses will be missing, and you will get different behavior.

Does this match the behavior you are seeing?

[ShaoxiongYao]

Thank you for your reply but an issue is that I need to resume not only the position but also velocity. I currently manually set these values using body.velocity = (v_x, v_y) after I add this object to space. I believe this does not create cached impulses that may exist in the underlying Chipmunk library. May I ask is there a possible way to resume the space with specific velocities?

[ShaoxiongYao]

I tried to use apply force instead of apply impulse, and the simulation can be reproduced. I think applying force does not change underlying cached impulse in this space, is this correct?

[viblo]

Both methods will activate (wake from sleep) the body.
Apply impulse will update velocity and angular velocity of body directly.
Apply force will set the force and torque (angular force). The force is then used in the update velocity function to update velocity once per space.step...

I should add that under to hood the damped spring works by applying impulses, so if the base apply impulse is unreliable so will damped springs.

Some things I can think of:

1. Just to double check, you have a fixed value for dt when you call space.step, right?
2. Could you make a (small) example with the behavior I could try myself?
3. You could try to set Space.collision_persistence to 0. Will probably make the simulation worse, but maybe it will make it repeatable.

[ShaoxiongYao]

1. Just double-checked dt is fixed value 0.02, is this value acceptable?
2. write an example

`

 import pymunk
 import pickle
 import numpy as np

def print_state(body_list):
    print("print state")
    for body in body_list:
        print(body.position, body.velocity)

def add_rope(space):
    balls = []

    # create the rope
    for i in range(15):
        moment = pymunk.moment_for_circle(1, 0, 5)
        body = pymunk.Body(1, moment)
        body.position = (10.5 * i, 0)
        shape = pymunk.Circle(body, 5)
        space.add(body, shape)
        balls.append(body)
    
    # connect balls next to each other
    for a, b in zip(balls[:-1], balls[1:]):
        spring = pymunk.DampedSpring(a, b, (0, 0), (0, 0), rest_length=10.5, stiffness=200, damping=50)
        space.add(spring)

    # connect balls with one ball in between
    for a, b in zip(balls[:-2], balls[2:]):
        spring = pymunk.DampedSpring(a, b, (0, 0), (0, 0), rest_length=2*10.5, stiffness=200, damping=50)
        space.add(spring)

    return balls

def apply_random_impulse(body_list):
    impulses_list = []
    for body in body_list:
        impulse = np.random.normal(0, 20, size=(2,))
        body.apply_impulse_at_local_point(impulse=impulse, point=(0,0))
        impulses_list.append(impulse)
    return impulses_list

def store_space_bodies(space, bodies):
    with open("space.pkl", 'wb') as f:
        pickle.dump((space, bodies), f)

def main():
    space = pymunk.Space()
    space.gravity = (0.0, 0.0)
    space.damping = 0.01

    balls = add_rope(space)

    for _ in range(100):
        apply_random_impulse(balls)
        space.step(0.02)

    # apply impulses before step
    apply_random_impulse(balls)
    store_space_bodies(space, balls)

    print_state(balls)
    space.step(0.02)
    print_state(balls)

if __name__ == "__main__":
    main()`

create a space with 15 balls connected by springs, and save the space before the step in space.pkl

`import pickle
import numpy as np

with open('space.pkl', 'rb') as f:
space, bodies = pickle.load(f)

print("print state")
for body in space.bodies:
print(body.position, body.velocity)

space.step(0.02)

print("print state")
for body in space.bodies:
print(body.position, body.velocity)`

loads the space.pkl and run one step

I expected the states after the step are the same, but the velocity have difference up to 0.1 which will accumulate to large values in about 1000 time steps.

[ShaoxiongYao]

Set np.random.seed(100), output from first script,
print state Vec2d(16.00564506776223, 12.897785950718719) Vec2d(24.688662442204407, 117.43482533622887) Vec2d(30.351513602377203, 15.823084355095963) Vec2d(39.0207300816362, 53.91865980657133) Vec2d(37.74599437735488, 6.777054107096236) Vec2d(-18.876469023389497, 11.413410812760468) Vec2d(50.99016688225109, 9.076567693053489) Vec2d(33.14868149605506, 18.7550937992691) Vec2d(58.986139354729644, 0.23981670577076403) Vec2d(5.309084272728137, 27.888768508841814) Vec2d(70.03188649891426, -0.02132835998703142) Vec2d(-25.127791515994247, -14.092461596993322) Vec2d(80.56310203377265, -0.817407133331093) Vec2d(-24.61072439332755, 6.258864629437042) Vec2d(88.64648556567055, -9.952141844527894) Vec2d(-41.81834049888057, -14.9754780533418) Vec2d(98.0731908790206, -8.018187039621528) Vec2d(21.414134678101462, 19.557754803763576) Vec2d(110.12698169701567, -15.797907502116203) Vec2d(43.670445761217465, 38.55362587268885) Vec2d(122.03715348815369, -9.152605910798696) Vec2d(20.37340094585184, 89.92230038777194) Vec2d(133.36977499209394, -11.481688975829492) Vec2d(16.73897605320624, 22.13327858587288) Vec2d(144.2901057750979, -3.809918627423708) Vec2d(20.027332138712406, 12.046147623936019) Vec2d(154.6809084304642, -9.204130787116181) Vec2d(-1.9700327427315827, -12.057040976072056) Vec2d(147.91442011193655, -20.192907960030535) Vec2d(-40.591119169977404, 40.45439258482302) print state Vec2d(16.49941831660632, 15.246482457443296) Vec2d(27.939585251866493, 102.0835875224425) Vec2d(31.131928204009927, 16.901457551227388) Vec2d(27.409600423846552, 47.49720307913304) Vec2d(37.36846499688709, 7.005322323351445) Vec2d(1.0393373889262678, 28.36109793428116) Vec2d(51.65314051217219, 9.451669569038872) Vec2d(8.388122248120046, 17.138388814956713) Vec2d(59.09232104018421, 0.7975920759476003) Vec2d(-3.3009535329172053, 15.378578990755738) Vec2d(69.52933066859438, -0.30317759192689786) Vec2d(-8.617495068318767, -15.572394050078481) Vec2d(80.0708875459061, -0.6922298407423522) Vec2d(-6.960531410659645, -5.535368629134767) Vec2d(87.79400769948022, -10.254956701935178) Vec2d(0.6078093799693675, -1.1438570805525483) Vec2d(98.51758462879535, -7.623726647205809) Vec2d(4.6394787073073624, 16.30278171648401) Vec2d(111.00039061224001, -15.026834984662425) Vec2d(14.37827585296515, 45.366641838195456) Vec2d(122.44462150707074, -7.354159903043257) Vec2d(7.504797291632363, 62.38957583447182) Vec2d(133.70455451315806, -11.039023404112035) Vec2d(5.219384222325519, 36.15007680458822) Vec2d(144.69065241787217, -3.5689956749449876) Vec2d(14.270620783715568, 21.701884582428757) Vec2d(154.64150777560957, -9.445271606637622) Vec2d(1.9931475165188135, -6.813327852954451) Vec2d(147.102597728537, -19.383820108334074) Vec2d(-29.39636803894742, 17.197126424389914)
output from the second script
print state Vec2d(16.00564506776223, 12.897785950718719) Vec2d(24.688662442204407, 117.43482533622887) Vec2d(88.64648556567055, -9.952141844527894) Vec2d(-41.81834049888057, -14.9754780533418) Vec2d(50.99016688225109, 9.076567693053489) Vec2d(33.14868149605506, 18.7550937992691) Vec2d(58.986139354729644, 0.23981670577076403) Vec2d(5.309084272728137, 27.888768508841814) Vec2d(154.6809084304642, -9.204130787116181) Vec2d(-1.9700327427315827, -12.057040976072056) Vec2d(37.74599437735488, 6.777054107096236) Vec2d(-18.876469023389497, 11.413410812760468) Vec2d(122.03715348815369, -9.152605910798696) Vec2d(20.37340094585184, 89.92230038777194) Vec2d(110.12698169701567, -15.797907502116203) Vec2d(43.670445761217465, 38.55362587268885) Vec2d(30.351513602377203, 15.823084355095963) Vec2d(39.0207300816362, 53.91865980657133) Vec2d(70.03188649891426, -0.02132835998703142) Vec2d(-25.127791515994247, -14.092461596993322) Vec2d(98.0731908790206, -8.018187039621528) Vec2d(21.414134678101462, 19.557754803763576) Vec2d(147.91442011193655, -20.192907960030535) Vec2d(-40.591119169977404, 40.45439258482302) Vec2d(144.2901057750979, -3.809918627423708) Vec2d(20.027332138712406, 12.046147623936019) Vec2d(133.36977499209394, -11.481688975829492) Vec2d(16.73897605320624, 22.13327858587288) Vec2d(80.56310203377265, -0.817407133331093) Vec2d(-24.61072439332755, 6.258864629437042) print state Vec2d(16.49941831660632, 15.246482457443296) Vec2d(27.86312246090799, 102.10085433973705) Vec2d(87.81011875569294, -10.25165140559473) Vec2d(0.8505269752209847, -1.387824928852377) Vec2d(51.65314051217219, 9.451669569038872) Vec2d(8.224578208403743, 17.147694932691934) Vec2d(59.09232104018421, 0.7975920759476003) Vec2d(-3.0812242228345, 15.54698744180916) Vec2d(154.64150777560957, -9.445271606637622) Vec2d(1.9365518833506283, -6.677213601725891) Vec2d(37.36846499688709, 7.005322323351445) Vec2d(1.0288295504373417, 28.34253349904828) Vec2d(122.44462150707074, -7.354159903043257) Vec2d(7.493675734430882, 62.38624548163456) Vec2d(111.00039061224001, -15.026834984662425) Vec2d(14.467612385306651, 45.66027611552293) Vec2d(31.131928204009927, 16.901457551227388) Vec2d(27.384245509516575, 47.51875343355925) Vec2d(69.52933066859438, -0.30317759192689786) Vec2d(-8.583767709196179, -15.4875421037615) Vec2d(98.50147357258263, -7.627031943546257) Vec2d(4.382197185564825, 16.29662579833478) Vec2d(147.102597728537, -19.383820108334074) Vec2d(-29.224107114879093, 17.231847919996145) Vec2d(144.69065241787217, -3.5689956749449876) Vec2d(14.403207838344345, 21.737432799298926) Vec2d(133.70455451315806, -11.039023404112035) Vec2d(5.022093744556288, 35.74766553599956) Vec2d(80.0708875459061, -0.6922298407423522) Vec2d(-7.052731412780005, -5.662340733885765)

[viblo]

Ok, thanks.
Did some small experiments now, and I can also see the issue. A little bit surprising for me was that even loading and running the same space.pkl file twice will give different results. Actually if I sum up all the velocity the sum is different on even the first print just after loading the pickeled file. Maybe something is going on in the pickling code.

This will need more investigation

[ShaoxiongYao]

I tried to print out impulses on springs and it shows they are 0s for space loaded from the pickle but the impulses in the original engine are not 0s. I wonder if that may cause this issue.

[viblo]

Ive made some more investigation into this.

Good news is that I found one problem which is easy to fix: Bodies and constraints in the space are kept in a set. When pickle / unpickle the set will be in different order (since sets are unordered). I committed a fix for this just now, which will make the order fixed, as long as you are using Python 3.6 or newer (since my fix is using dict instead, and they are ordered starting from Python 3.6.). I think this is a ok tradeoff, since otherwise Pymunk would have to include or depend on a 3rd party library for ordered set.

If I also set the space.collision_persistence = 0 I get completely stable result from pickle/unpickle. This shows that the remaining issue is the internal contact cache used by Pymunk/Chimpmunk. This will be more difficult to fix, I will have to do more thinking about this. Maybe you can try with space.collision_persistence = 0 and see if the simulation quality is good enough for your case without collision caching?

[viblo]

Small update: To make it fully stable the cached collision data needs to be exported out from Chipmunk. Currently Chipmunk does not provide any API for this, so it needs to be updated to allow the data to be extracted, and then Pymunk can be updated to read/write it when pickling.

Im currently in the process of modifing the binding code (Pymunk uses CFFI to talk with the c library) to better allow these kinds of modifications. But I think it will take some time to make work..

So, I will keep this issue open, and take another look at it once the binding code updating is completed.

[viblo]

After quite some time I have now implemented a fix for this issue with #221 Pickling & unpickling should work as expected with the same result each time. This will be included in the next release of Pymunk.