How to parallelize this situation with robots

Question

Hello,

I'm working on a robotic problem. The situation is something like this:

1. There are N number of robots (generally N>100) initially all at rest.
2. Each robot attracts all other robots which are with in its radius r.
3. I've set of equations with which I can compute acceleration, velocity & hence the position of the robot after time deltat. Simply put, I can find the position of each robot after deltat time.
4. All I need to do is for a given deltat. I need to display position of each robot for every deltat.

Problem is actually very simple. Algo will be something like:

del_t = ;its given
initialPositions = ;its given
num_robots = ;its given

The following code executes for every del_t

robots = range(1,no_robots)
for R in robots:
    for r in robots:
        if  distanceBetween(r,R) <= radius and r is not R:
            acceleration_along_X[R] += xAcceleration( position(r), position(R) )
            acceleration_along_Y[R] += yAcceleration( position(r), position(R) )
    currVelocity_along_X[R] = prevVelocity_along_X[R] + acceleration_along_X[R] * del_t
    currVelocity_along_Y[R] = prevVelocity_along_Y[R] + acceleration_along_Y[R] * del_t
    curr_X_coordinate[R] = prev_X_coordinate[R] + currVelocity_along_X[R] * del_t
    curr_Y_coordinate[R] = prev_Y_coordinate[R] + currVelocity_along_Y[R] * del_t
    print 'Position of robot ' + str(R) + ' is (' + curr_X_coordinate[R] + ', ' + curr_Y_coordinate[R] +' ) \n'
    prev_X_coordinate[R] = curr_X_coordinate[R]
    prev_Y_coordinate[R] = curr_Y_coordinate[R]
    prevVelocity_along_X[R] = currVelocity_along_X[R]
    prevVelocity_along_Y[R] = currVelocity_along_Y[R]

Now I need to parallelize the algorithm and set up the Cartesian grid of MPI processes.

1. Because computation for each robot is an independent task. computation for each Robot can be done by an independent thread. Right?
2. I don't know anything about MPI. What does this "Cartesian grid of MPI processes" mean? How can I setup this grid? I've no clue about this.

EDIT:

Now the problem turned interesting. Actually, it isn't as simple as I thought. After reading Unode's answer. I went on to apply his method two by parallelizing using multiprocessing.

This is the code. printPositionOfRobot is my serial algo. Basically, it is supposed to print the position of robot (with id *robot_id*) t=1,2,3,4,5,6,7,8,9,10. (Here *del_t* is taken as 1. *num_iterations = 10*. Each of the robot prints message like this: Robot8 : Position at t = 9 is (21.1051065245, -53.8757356694 )

There is bug in this code. t=0 locations of bots are given by position() for determining xAcceleration & yAcceleration. We need use the positions of previous iterations of all other particles.

from multiprocessing import Pool
import math


def printPositionOfRobot(robot_id):
    radius = 3
    del_t = 1
    num_iterations = 10
    no_robots = 10

    prevVelocity_along_X = 0
    prevVelocity_along_Y = 0
    acceleration_along_X = 0
    acceleration_along_Y = 0
    (prev_X_coordinate,prev_Y_coordinate) = position(robot_id)#!!it should call initialPosition()
    for i in range(1,num_iterations+1):
        for r in range(no_robots):
            if  distanceBetween(r,robot_id) <= radius and r is not robot_id:
                acceleration_along_X += xAcceleration( position(r), position(robot_id) ) #!! Problem !!
                acceleration_along_Y += yAcceleration( position(r), position(robot_id) )#!! Problem !!
        currVelocity_along_X = prevVelocity_along_X + acceleration_along_X * del_t
        currVelocity_along_Y = prevVelocity_along_Y + acceleration_along_Y * del_t
        curr_X_coordinate = prev_X_coordinate + currVelocity_along_X * del_t
        curr_Y_coordinate = prev_Y_coordinate + currVelocity_along_Y * del_t
        print 'Robot' + str(robot_id) + ' : Position at t = '+ str(i*del_t) +' is (' + str(curr_X_coordinate) + ', ' + str(curr_Y_coordinate) +' ) \n'
        prev_X_coordinate = curr_X_coordinate
        prev_Y_coordinate = curr_Y_coordinate
        prevVelocity_along_X = currVelocity_along_X
        prevVelocity_along_Y = currVelocity_along_Y

def xAcceleration((x1,y1),(x2,y2)):
    s = distance((x1,y1),(x2,y2))
    return 12*(x2-x1)*( pow(s,-15) - pow(s,-7) + 0.00548*s )

def yAcceleration((x1,y1),(x2,y2)):
    s = distance((x1,y1),(x2,y2))
    return 12*(y2-y1)*( pow(s,-15) - pow(s,-7) + 0.00548*s )

def distanceBetween(r,robot_id):
    return distance(position(r), position(robot_id))

def distance((x1,y1),(x2,y2)):
    return math.sqrt( (x2-x1)**2 + (y2-y1)**2 )

def Position(r): #!!name of this function should be initialPosition
    k = [(-8.750000,6.495191) , (-7.500000,8.660254) , (-10.000000,0.000000) , (-8.750000,2.165064) , (-7.500000,4.330127) , (-6.250000,6.495191) , (-5.000000,8.660254) , (-10.000000,-4.330127) , (-8.750000,-2.165064) , (-7.500000,0.000000) ]
    return k[r]

if __name__ == "__main__":
    no_robots = 10  # Number of robots you need
    p = Pool(no_robots)  # Spawn a pool of processes (one per robot in this case)
    p.map(printPositionOfRobot, range(no_robots))

the position function in acceleration_along_X & acceleration_along_Y should return the latest position of the robot.By latest I mean the position at the end of that previous iteration. So, each processes must inform other processes about its latest position. Until the latest position of the bot is know the process must wait.

Other way can be that all processes edit a global location.(I wonder if its possible, because each process have its own Virtual address space). If a process has not yet reached that iteration all other processes must wait.

Any ideas about how to go about it? I guess this is why MPI was suggested in the problem.

Answer 1

Note: Python's threads still run on the same processor. If you want to use the full range of processors of your machine you should use multiprocessing (python2.6+).

Using MPI will only bring you clear benefits if the computation is going to be spread over multiple computers.

There are two approaches to your problem. Since you have completely independent processes, you could simply launch the algorithm (passing a unique identifier for each robot) as many times as needed and let the operating system handle the concurrency.

1 - A short Linux shell script (or something equivalent in Windows BATCH language):

#!/bin/sh
for i in {0..99}; do
    echo "Running $i"
    python launch.py $i &
done

Note: the & after the launch.py this ensures that you actually launch all processes in consecutive way, rather than waiting for one to finish and then launch the next one.

2 - If instead you want to do it all in python, you can use the following simple parallelization approach:

from multiprocessing import Pool

def your_algorithm(robot_id):
    print(robot_id)

if __name__ == "__main__":
    robots = 100  # Number of robots you need
    p = Pool(robots)  # Spawn a pool of processes (one per robot in this case)
    p.map(your_algorithm, range(robots))

The map function takes care of dispatching one independent operation per robot.

If you do require the use of MPI I suggest mpi4py.

As for information on what Cartesian grid stands for, try this

Answer 2

My solution too I similar to Unode but I prefer using apply_async method in multiprocessing as it's asynchronous.

from multiprocessing import Pool

    def main():
        po = Pool(100) #subprocesses
        po.apply_async(function_to_execute, (function_params,), callback=after_processing)
        po.close() #close all processes
        po.join() #join the output of all processes
        return    

Answer 3

So the trick here is that at each step, all of the robots have to see the data at some point from all the other robots. This makes efficient parallelizations hard!

One simple approach is to have each process chuging away on its own robots, calculating the self-interactions first, then getting one by one the data from the other processors and calculating those forces. Note that this is not the only approach! Also, real-world solvers for this sort of thing (molecular dynamics, or most astrophsical N-body simulations) take a different tack entirely, treating distant objects only approximately since far away objects don't matter as much as near ones.

Below is a simple python implementation of that approach using two mpi processes (using mpi4py and matplotlib/pylab for plotting). The generalization of this would be a pipeline; each processor sends its chunk of data to the next neighbour, does the force calcs, and this process repeats until everyone has seen everyone's data. Ideally you'd update plot as the pipeline progressed, so that no one has to have all of the data in memory at once.

You'd run this program with mpirun -np 2 ./robots.py ; note that you need the MPI libraries installed, and then the mpi4py needs to know where to find these libraries.

Note too that I'm being very wasteful in sending all of the robot data along to the neighbour; all the neighbour cares about is the positions.

#!/usr/bin/env python
import numpy
import pylab
from mpi4py import MPI

class Robot(object):
    def __init__(self, id, x, y, vx, vy, mass):
        self.id = id
        self.x = x
        self.y = y
        self.vx = vx
        self.vy = vy
        self.ax = 0.
        self.ay = 0.
        self.mass = mass
    def rPrint(self):
        print "Robot ",self.id," at (",self.x,",",self.y,")"
    def interact(self, robot2):
        dx = (self.x-robot2.x)
        dy = (self.y-robot2.y)
        eps = 0.25
        idist3 = numpy.power(numpy.sqrt(dx*dx +dy*dy + eps*eps),-3)
        numerator = -self.mass*robot2.mass
        self.ax += numerator*dx*idist3
        self.ay += numerator*dy*idist3
        robot2.ax -= numerator*dx*idist3
        robot2.ay -= numerator*dy*idist3
    def updatePos(self, dt):
        self.x += 0.5*self.vx*dt
        self.y += 0.5*self.vy*dt
        self.vx += self.ax*dt
        self.vy += self.ay*dt
        self.x += 0.5*self.vx*dt
        self.y += 0.5*self.vy*dt
        self.ax = 0.
        self.ay = 0.



def init(nRobots):
    myRobotList = []
    vx = 0.
    vy = 0.
    mass = 1.
    for i in range(nRobots):
        randpos = numpy.random.uniform(-3,+3,2)
        rx = randpos[0]
        ry = randpos[1]
        myRobotList.append(Robot(i, rx, ry, vx, vy, mass))
    return myRobotList

def selfForces(robotList):
    nRobots = len(robotList)
    for i in range(nRobots-1): 
       for j in range (i+1, nRobots):
            robotList[i].interact(robotList[j])

def otherRobotForces(myRobotList, otherRobotList):
    for i in myRobotList:
        for j in otherRobotList:
            i.interact(j)

def plotRobots(robotList):
    xl = []
    yl = []
    vxl = []
    vyl = [] 
    for i in robotList:
       xl.append(i.x)
       yl.append(i.y)
       vxl.append(i.vx)
       vyl.append(i.vy)
    pylab.subplot(1,1,1)
    pylab.plot(xl,yl,'o')
    pylab.quiver(xl,yl,vxl,vyl)
    pylab.show()

if __name__ == "__main__":
    comm = MPI.COMM_WORLD
    nprocs = comm.Get_size()
    rank   = comm.Get_rank()

    if (nprocs != 2):
        print "Only doing this for 2 for now.."
        sys.exit(-1)
    neigh = (rank + 1) %  nprocs

    robotList = init(50)

    for i in range (10):
        print "[",rank,"] Doing step ", i
        selfForces(robotList)

        request = comm.isend(robotList, dest=neigh, tag=11)
        otherRobotList = comm.recv(source=neigh, tag=11)

        otherRobotForces(robotList,otherRobotList)

        request.Wait()

        for r in robotList:
            r.updatePos(0.05)



    if (rank == 0):
        print "plotting Robots"
        plotRobots(robotList + otherRobotList)