[Stackless] Stackless Python and Simulation

Simon Frost sdfrost at ucsd.edu
Fri Jul 11 02:40:12 CEST 2003


Dear SimPy Users, Stackless List, and Psyco List,

(apologies for the cross-post)

I'm attaching an example of a SimPy-like model that uses Stackless Python's 
tasklets rather than generators to simulate coroutines. I've tried to keep 
the syntax as SimPy-like as possible for comparison. You'll need Stackless 
Python (available from http://www.stackless.com). For real simulations, I'd 
use an AVL tree rather than bisect to keep the list of sorted timestamps, 
and a proper random number generator.

I have a couple of problems with the code though related to Psyco:

1. The use of 'from psyco.classes import __metaclass__' results in about a 
two-fold speedup. Is there any way I can eke out more performance by 
syntactic changes? Rather strangely, I don't get this with the SimPy 
version (perhaps because it uses generators?).

2. The use of 'import psyco' and 'psyco.full()' results in an access 
violation on my Win2K machine. Can anyone else reproduce this?

Best
Simon

#!/usr/local/bin/python2.2s
# epidemic_stackless.py
# Authors: Simon Frost, Klaus Muller, Tony Vignaux
# Date: July 10th 2003
# A simple example of an epidemic process:
# 1. Generate N susceptible individuals.
# 2. For each individual, schedule a contact according to a random 
exponential variable.
# 3. Find the individual, 'A', with the next contact; if the next contact 
is after
#    'runTime', end
# 4. Choose another individual, 'B', to contact.
# 5. If individual 'B' is infected, then there is a probability 
'infectProb' that A
#    becomes infected.
# 6. If all individuals are infected, end; otherwise, go to 3.
#
# This uses Stackless Python (http://www.stackless.com) to simulate 
co-routines.

# Use Psyco if available
try:
         from psyco.classes import __metaclass__
except ImportError:
         pass

# Import libraries
import stackless
import bisect
import random
import time

class Epidemic:
         """A simple epidemic simulation."""
         def __init__(self,infectProb,contactRate,now=0,rng=random.Random(1)):
                 Epidemic.allIndividuals=[]
                 Epidemic.eventList=Epidemic.EventList()
                 Epidemic.contactRate=contactRate
                 Epidemic.infectProb=infectProb
                 Epidemic.now=now
                 Epidemic.times=[now]
                 Epidemic.numInfected=[0]
                 Epidemic.rng=rng
                 Epidemic.channel=stackless.channel()

         def activate(self):
                 for ind in Epidemic.allIndividuals:
                         stackless.tasklet(ind.live)()
                         if ind.infected==0:
                                 nextEvent=Epidemic.channel.receive()
                                 Epidemic.eventList.insert(nextEvent[0],nextEvent[1])

         def simulate(self,until):
                 while Epidemic.now <= until:
                         thisEvent=Epidemic.eventList.pop()
                         when=thisEvent[0]
                         who=thisEvent[1]
                         Epidemic.now=when
                         who.channel.send(None)
                         try:
                                 nextEvent=Epidemic.channel.receive()
                         except RuntimeError:
                                 nextEvent=Epidemic.eventList.pop()
                         Epidemic.eventList.insert(nextEvent[0],nextEvent[1])

         class EventList:
                 """A class to maintain an ordered list using a bisection 
algorithm."""
                 def __init__(self):
                         self.queue = []
                 def insert(self, timestamp,data):
                         """ Insert a new element in the queue according to 
its priority. """
                         bisect.insort(self.queue, [timestamp, data])
                 def pop(self):
                         """ Pop the highest-priority element of the queue. """
                         return self.queue.pop(0)
                 def remove(self,timestamp,data):
                         """ Remove item from queue. """
                         item_delete_point = 
bisect.bisect(self.queue,(timestamp,data))
                         del self.queue[item_delete_point-1]
                 def is_present(self,timestamp,data):
                         """ Determines whether item is present in queue. """
                         item_insert_point = bisect.bisect(self.queue, 
(timestamp,data))
                         is_present = 
self.queue[item_insert_point-1:item_insert_point] == [item]

         class Individual:
                 """An individual in a population that is either 
susceptible or infected."""
                 def __init__(self):
                         self.infected=0
                         self.channel=stackless.channel()
                         Epidemic.allIndividuals.append(self)

                 def contactInterval(self):
                         return Epidemic.rng.expovariate(Epidemic.contactRate)

                 def chooseContact(self):
                         contact = self
                         while contact == self:
                                 contact=Epidemic.rng.choice(Epidemic.allIndividuals)
                         return contact

                 def live(self):
                         while self.infected==0:
                                 ci = self.contactInterval()
                                 Epidemic.channel.send((Epidemic.now+ci,self))
                                 self.channel.receive()
                                 B=self.chooseContact()
                                 if B.infected ==1 and 
Epidemic.rng.random() <= Epidemic.infectProb:
                                         Epidemic.times.append(Epidemic.now)
                                         Epidemic.numInfected.append(Epidemic.numInfected[-1]+1)
                                         print 
str(Epidemic.now),Epidemic.numInfected[-1]
                                         self.infected = 1

def model(nrIndividuals,contactRate,infectProb,initialInfected,runTime):
         run=Epidemic(infectProb=infectProb,contactRate=contactRate)
         for i in range(nrIndividuals):
                 ind=Epidemic.Individual()
                 if i >= (nrIndividuals-initialInfected):
                         ind.infected=1
                         Epidemic.numInfected[0]+=1
         run.activate()
         run.simulate(until=runTime)

if __name__=="__main__":
         t=time.time()
         run=model(nrIndividuals=100,contactRate=3,infectProb=0.003,initialInfected=1,runTime=1000)
         print time.time()-t



_______________________________________________
Stackless mailing list
Stackless at www.tismer.com
http://www.tismer.com/mailman/listinfo/stackless




More information about the Stackless mailing list