{"id":52,"date":"2015-10-30T10:09:10","date_gmt":"2015-10-30T17:09:10","guid":{"rendered":"http:\/\/www.howardklein.net\/?p=52"},"modified":"2015-11-20T17:20:42","modified_gmt":"2015-11-21T01:20:42","slug":"simulating-an-mm1-queue-using-generators-and-greenlets","status":"publish","type":"post","link":"https:\/\/www.howardklein.net\/?p=52","title":{"rendered":"Simulating an M\/M\/1 Queue Using Generators and Greenlets"},"content":{"rendered":"<p>My last two posts (<a href=\"https:\/\/www.howardklein.net\/?p=41\">here <\/a>and <a href=\"https:\/\/www.howardklein.net\/?p=46\">here<\/a>) introduced three Python-based alternatives for expressing simulation process or agent-based behavior that takes place over (or blocks for) simulated time:<\/p>\n<ul>\n<li>Generators, a core feature of the Python language<\/li>\n<li>Greenlets, provided by the third-party <a href=\"http:\/\/greenlet.readthedocs.org\/en\/latest\/\" target=\"_blank\">greenlet <\/a>Python extension package<\/li>\n<li>Tasklets, provided by two standalone Python distributions, Stackless Python and <a href=\"http:\/\/doc.pypy.org\/en\/latest\/\" target=\"_blank\">PyPy<\/a><\/li>\n<\/ul>\n<p>While those posts included a few brief code snippets, perhaps now is a good time to add meat to the bones and present the full code for one (albeit simple) simulation model.\u00a0 We\u2019ll model the classic <a href=\"http:\/\/www.win.tue.nl\/~iadan\/que\/h4.pdf\" target=\"_blank\">M\/M\/1 queue<\/a>, a system with a single queue and a single server.\u00a0 Both interarrival and service times are exponentially distributed.\u00a0 There is no upper bound on the queue length.\u00a0 We\u2019ll implement a process-oriented simulation in two ways; one using generators, the other, greenlets.\u00a0 (With small modifications, we could convert the greenlet-based code to tasklets.)<\/p>\n<p><!--more-->Most of the code is common to both approaches. \u00a0First a few basics, starting with a clock:<\/p>\n<pre class=\"brush: python; title: ; notranslate\" title=\"\">\r\nclass SimClock(object):\r\n    current_time = 0\r\n    \r\n    @staticmethod\r\n    def now():\r\n        return SimClock.current_time\r\n    \r\n    @staticmethod\r\n    def advance_to(tm):\r\n        SimClock.current_time = tm\r\n<\/pre>\n<p>Every simulation needs a clock \ud83d\ude42 .\u00a0 For now, we\u2019ll keep it simple and unit-less. We could keep it even simpler and just define <code><span style=\"font-family: 'courier new', courier, monospace;\">current_time<\/span><\/code> as a global variable, but we\u2019ll wrap it inside a class instead, where <code><span style=\"font-family: 'courier new', courier, monospace;\">now()<\/span><\/code> returns the current simulated time and <code><span style=\"font-family: 'courier new', courier, monospace;\">advance_to()<\/span><\/code> advances the simulated time.<\/p>\n<pre class=\"brush: python; title: ; notranslate\" title=\"\">\r\neventlist = &#x5B;]\r\ncounter = itertools.count()\r\ndef initialize():\r\n    eventlist.clear()\r\n    SimClock.current_time = 0\r\n<\/pre>\n<p>The simulation must maintain a list of the events to process.\u00a0 That list functions as a priority queue, where priority is the simulated time the event is scheduled to be fired\/handled; if two events are scheduled for the same time, they should be handled first-in, first-out.<\/p>\n<p>A single simulated M\/M\/1 queue is not going to terribly stress the simulator\u2019s event manager, so I generally wouldn\u2019t put a great deal of effort into optimizing this code right now.\u00a0 Fortunately, the Python standard library provides a <a href=\"https:\/\/docs.python.org\/3.5\/library\/heapq.html#\" target=\"_blank\">heap<\/a>, which makes for a pretty decent <a href=\"https:\/\/docs.python.org\/3.5\/library\/heapq.html#priority-queue-implementation-notes\" target=\"_blank\">priority queue implementation<\/a> with O(log n) performance. \u00a0With almost no additional effort, we\u00a0can keep events on a \u201cheapified\u201d global list, and maintain a counter to maintain the secondary FIFO priority.<\/p>\n<p>The <code><span style=\"font-family: 'courier new', courier, monospace;\">initialize()<\/span><\/code> function simply sets (or resets) both the event list and clock prior to starting a simulation run.<\/p>\n<p>We now have an event list, so we need events:<\/p>\n<pre class=\"brush: python; title: ; notranslate\" title=\"\">\r\nclass SimEvent(object):\r\n    def __init__(self, simtime):\r\n        self.time = simtime\r\n        \r\n    def register(self):\r\n        heappush(eventlist, (self.time, next(counter), self))\r\n        \r\n    @abstractmethod\r\n    def handle(self):\r\n        &quot;&quot;&quot;\r\n        Process the event.  Implemented by subclasses\r\n        &quot;&quot;&quot;\r\n        pass\r\n<\/pre>\n<p><code><span style=\"font-family: 'courier new', courier, monospace;\">SimEvent<\/span><\/code> is an abstract base class for the concrete event classes described below (<code><span style=\"font-family: 'courier new', courier, monospace;\">ArrivalEvent<\/span> <\/code>and <code><span style=\"font-family: 'courier new', courier, monospace;\">ResumeEvent<\/span><\/code>).\u00a0 The base class includes one piece of data: the time at which the event should be handled\/processed.\u00a0 The actual handling\/processing is performed by the <code><span style=\"font-family: 'courier new', courier, monospace;\">handle()<\/span><\/code> method, which must be implemented by the subclasses.\u00a0 The <code><span style=\"font-family: 'courier new', courier, monospace;\">register()<\/span><\/code> method simply pushes the event on to the heapified list (priority queue) as part of a time\/count\/event tuple; the event time and counter value are be used to prioritize this event within the queue.<\/p>\n<p>And to finish up with preliminaries, our simulator needs an event processor:<\/p>\n<pre class=\"brush: python; title: ; notranslate\" title=\"\">\r\ndef processEvents(until_time):\r\n    while eventlist:\r\n        next_event = heappop(eventlist)&#x5B;2]\r\n        SimClock.advance_to(next_event.time)\r\n        if SimClock.now() &amp;gt; until_time:\r\n            break\r\n        else:\r\n            next_event.handle()\r\n<\/pre>\n<p>The <code><span style=\"font-family: 'courier new', courier, monospace;\">processEvents()<\/span><\/code> function executes the simulation run for a specified length of simulated time. A seed <code><span style=\"font-family: 'courier new', courier, monospace;\">ArrivalEvent<\/span> <\/code>(see below) must be registered prior to calling <code>processEvents()<\/code>; otherwise, the initial event list will be empty and the simulation stops.\u00a0 Otherwise, <code>processEvents<\/code> runs a loop until either we run out of events or the clock is advanced passed the specified <code>until_time<\/code>.\u00a0 The loop body simply pops\u00a0the next scheduled event, advances the clock to that event\u2019s scheduled execution time, and then processes\/executes the event via it\u2019s <code>handle()<\/code> method.<\/p>\n<p>On to the aforementioned <code>ArrivalEvent<\/code>:<\/p>\n<pre class=\"brush: python; title: ; notranslate\" title=\"\">\r\nclass ArrivalEvent(SimEvent):\r\n    def __init__(self, interarrival_mean, service_mean, server,                   \r\n                 process_cls):\r\n        super().__init__(0)\r\n        self.interarrival_mean = interarrival_mean\r\n        self.service_mean = service_mean\r\n        self.server = server\r\n        self.process_cls = process_cls\r\n        self.rng = random.Random()\r\n        self.rng.seed(1000420099314159)\r\n        self.time = self.rng.expovariate(1\/self.interarrival_mean)\r\n        \r\n    def handle(self):\r\n        service_time = self.rng.expovariate(1\/self.service_mean)\r\n        process = self.process_cls(self.server, service_time)\r\n        start_process_event = ResumeEvent(process, SimClock.now())\r\n        start_process_event.register()\r\n        next_ia_time = self.rng.expovariate(1\/self.interarrival_mean)\r\n        self.time = SimClock.now() + next_ia_time\r\n        self.register()\r\n<\/pre>\n<p>The <code>ArrivalEvent<\/code> looks a bit more complicated. \u00a0Its initializer is passed the basic M\/M\/1 model parameters (mean interarrival and service time) along with a reference to the single server and a process class.\u00a0 A Generator-based simulation uses the <code>GeneratorProcess<\/code> class, while a greenlet-based simulation uses <code>GreenletProcess<\/code> (as discussed below).\u00a0 The event\u2019s \u00a0<code>handle()<\/code> method does the following:<\/p>\n<ol>\n<li>It instantiates a new process object using the process class passed to the initializer<\/li>\n<li>The service time is sampled from an exponential distribution (with a mean specified in the initializer) via a pseudo-random number generator (attribute <code>rng<\/code>, created in the initializer with a hard-coded seed.) The service time is also passed to the new process object, along with a reference to the server object.<\/li>\n<li>The new process is started by instantiating and scheduling a <code>ResumeEvent<\/code> for right now (in simulated time)<\/li>\n<li>Finally, the <code>ArrivalEvent<\/code> reschedules itself by sampling from another exponential distribution (with mean equal to the initialized interarrival time), setting a new event time, and reregistering itself. Rinse, lather, repeat.<\/li>\n<\/ol>\n<p>As noted above, a single <code>ArrivalEvent<\/code> must be created and registered prior to running <code>processEvents()<\/code>.\u00a0 The <code>ArrivalEvent<\/code> essentially seeds the simulation, providing an initial event to execute.<\/p>\n<pre class=\"brush: python; title: ; notranslate\" title=\"\">\r\nclass ResumeEvent(SimEvent):\r\n    def __init__(self, process, resumetime):\r\n        super().__init__(resumetime)\r\n        self.process = process\r\n        \r\n    def handle(self):\r\n        self.process.resume()\r\n<\/pre>\n<p>The <code>ResumeEvent<\/code> starts execution of a new process, or restarts a process that is currently blocked\/not executing while either acquiring the server or waiting for the service time to pass.\u00a0 This is done through the process object\u2019s <code>resume()<\/code> method.\u00a0 (The <code>GeneratorProcess<\/code> and <code>GreenletProcess<\/code> classes each have their own implementation of <code>resume()<\/code>.)<\/p>\n<p>Now that we\u2019re through with event classes, it\u2019s on to the server (or resource, in alternative jargon):<\/p>\n<pre class=\"brush: python; title: ; notranslate\" title=\"\">\r\nclass Server(object):\r\n    def __init__(self):\r\n        self.requests = deque()\r\n        self.total_queue_time = 0\r\n        self.process_count = 0\r\n        self.current_process = None\r\n        \r\n    def request(self, process):\r\n        self.process_count += 1\r\n        if self.current_process:\r\n            self.requests.append(process)\r\n            return False\r\n        else:\r\n            self.current_process = process\r\n            return True\r\n    \r\n    def release(self, process):\r\n        if self.requests:\r\n            self.current_process = self.requests.popleft()\r\n            event = ResumeEvent(self.current_process, SimClock.now())\r\n            event.register()\r\n            qtime = SimClock.now() - self.current_process.start_time\r\n            self.total_queue_time += qtime\r\n        else:\r\n            self.current_process = None  \r\n                         \r\n    @property\r\n    def mean_queue_time(self):\r\n        if self.process_count:\r\n            return self.total_queue_time \/ self.process_count\r\n        else:\r\n            return None \r\n\r\n@property\r\n    def queueSize(self):\r\n        return len(self.requests)\r\n<\/pre>\n<p>The <code>Server<\/code> class implements a server resource through <code>request()<\/code> and <code>release()<\/code> methods.\u00a0 Internally, it maintains a reference to the process object that currently holds the server, and a FIFO queue of processes that are waiting to acquire the server. \u00a0In another moment of gratuitous optimization, we\u2019ll implement that queue as a <code>deque<\/code> rather than a standard Python list.\u00a0 When the server is released, it pulls the next requesting process of that queue, makes it the currently-being-served process, and tells the prcoess to resume (via a <code>ResumeEvent<\/code>).<\/p>\n<p>The server instance also tracks queue time data (time in queue, waiting for a request to be fulfilled), and returns the mean queue time via a property. If this were a real model we would certainly be tracking any number of other outputs, but this is enough to start.<\/p>\n<p>Since we\u2019re building a process-oriented model, we need a process class:<\/p>\n<pre class=\"brush: python; title: ; notranslate\" title=\"\">\r\nclass SimProcess(object):\r\n    def __init__(self, server, servicetime):\r\n        self.server = server\r\n        self.servicetime = servicetime\r\n        self.start_time = SimClock.now()\r\n        \r\n    @abstractmethod\r\n    def resume(self):\r\n        pass\r\n        \r\n    @abstractmethod\r\n    def acquire(self, server):\r\n        pass\r\n        \r\n    @abstractmethod\r\n    def wait_for(self, delaytime):\r\n        pass\r\n        \r\n    @abstractmethod\r\n    def execute(self):\r\n        pass\r\n<\/pre>\n<p><code>SimProcess<\/code> is the abstract base class defining the interface for the generator and greenlet-based process classes mentioned above:<\/p>\n<ul>\n<li><code>resume()<\/code>, as noted above, is invoked by the ResumeEvent<\/li>\n<li><code>execute()<\/code> is invoked immediately after process creation; as its name implies, this method encapsulates the execution of the simulation process (for one simulation entity) over simulated time<\/li>\n<li><code>acquire()<\/code> and <code>wait_for()<\/code> are typically invoked by subclass implementations of <code>execute()<\/code><\/li>\n<\/ul>\n<p>Finally, we have the generator and greenlet-specific subclasses of <code>SimProcess<\/code>, which each also include M\/M\/1-specific processing logic:<\/p>\n<pre class=\"brush: python; title: ; notranslate\" title=\"\">\r\nclass GeneratorProcess(SimProcess):    \r\n    def __init__(self, server, servicetime):\r\n        super().__init__(server, servicetime)\r\n        self.generator = self.execute()\r\n        \r\n    def resume(self):\r\n        try:\r\n            next(self.generator)\r\n        except StopIteration as e:\r\n            pass\r\n                \r\n    def acquire(self, server):\r\n        acquired = server.request(self)\r\n        if not acquired:\r\n            yield\r\n        \r\n    def wait_for(self, delaytime):\r\n        event = ResumeEvent(self, SimClock.now() + delaytime)\r\n        event.register()\r\n        yield\r\n        \r\n    def execute(self):\r\n        yield from self.acquire(self.server)\r\n        yield from self.wait_for(self.servicetime)\r\n        self.server.release(self)\r\n\r\n\r\nclass GreenletProcess(SimProcess):\r\n    maingreenlet = greenlet.getcurrent() # the root greenlet\r\n    \r\n    def __init__(self, server, servicetime):\r\n        super().__init__(server, servicetime)\r\n        self.greenlet = greenlet(self.execute)\r\n        \r\n    def suspend(self):\r\n        self.maingreenlet.switch()\r\n        \r\n    def resume(self):\r\n        self.greenlet.switch()\r\n        \r\n    def acquire(self, server):\r\n        acquired = server.request(self)\r\n        if not acquired:\r\n            self.suspend()\r\n        \r\n    def wait_for(self, delaytime):\r\n        event = ResumeEvent(self, SimClock.now() + delaytime)\r\n        event.register()\r\n        self.suspend()\r\n        \r\n    def execute(self):\r\n        self.acquire(self.server)\r\n        self.wait_for(self.servicetime)\r\n        self.server.release(self)\r\n<\/pre>\n<p>A few notes:<\/p>\n<ul>\n<li>As noted in my earlier posts, generator-based process methods (<code>execute()<\/code>) are started and resumed through a <code>next()<\/code> call on the generator created in the process initializer. Similarly, greenlet-based process execution is started\/resumed through a <code>switch()<\/code> call on the process\u2019s <code>greenlet<\/code> (also created in the initializer)<\/li>\n<li>For a generator-based process, execution is suspended (returning control to the event processor) through a <code>yield<\/code> statement. For greenlet-based processes, this is done through a switch back to the main (parent) greenlet.<\/li>\n<li>A greenlet-based process simply calls <code>acquire()<\/code> and <code>wait_for()<\/code> like any other method; a generator-based implementation must remember to use a <code>yield from<\/code> statement.<\/li>\n<\/ul>\n<p>Last but not least, a bit of code to actually initialize and execute two simulation runs (one for each <code>SimProcess<\/code> implementation):<\/p>\n<pre class=\"brush: python; title: ; notranslate\" title=\"\">\r\nmean_interarrival_time = 10\r\nmean_service_time = 8\r\nclasses = (GreenletProcess, GeneratorProcess)\r\nfor cls in classes:    \r\n    initialize()\r\n    server = Server()\r\n    arrival_event = ArrivalEvent(mean_interarrival_time, \r\n                                 mean_service_time, server, cls)\r\n    arrival_event.register()\r\n    processEvents(1000000)\r\n        \r\n    print(&quot;Process Class:&quot;, cls.__name__)\r\n    print(&quot;total entries processed:&quot;, server.process_count)\r\n    print(&quot;mean queue time:&quot;, server.mean_queue_time)\r\n<\/pre>\n<p>This code simply defines the parameters of the M\/M\/1 model (mean interarrival and service time) and runs the simulation twice, with each simulation being run for one million simulated time units.\u00a0 The resulting output is:<\/p>\n<pre class=\"brush: plain; title: ; notranslate\" title=\"\">\r\nProcess Class: GreenletProcess\r\ntotal entries processed: 99803\r\nmean queue time: 31.48626676016153\r\nProcess Class: GeneratorProcess\r\ntotal entries processed: 99803\r\nmean queue time: 31.48626676016153\r\n<\/pre>\n<p>It\u2019s always reassuring to see that the two implementations generate identical results.<\/p>\n<p>Queuing theory allows us to independently calculate the expected mean queue time analytically:<\/p>\n<p style=\"padding-left: 30px;\">\u03bb = 1 \/ mean interarrival time1 = 0.125<\/p>\n<p style=\"padding-left: 30px;\">\u00b5 = 1 \/ mean service time = 0.1<\/p>\n<p style=\"padding-left: 30px;\">\u03c1 (utilization) = \u03bb \/ \u00b5 = 0.8<\/p>\n<p style=\"padding-left: 30px;\">mean queue time = \u03c1 \/ (\u00b5 &#8211; \u03bb) = 32<\/p>\n<p>I\u2019ll declare our simulation models valid based on the simulated mean queue time of just under 31.5 \u00a0\ud83d\ude42 .<\/p>\n<p>&nbsp;<\/p>\n","protected":false},"excerpt":{"rendered":"<p>My last two posts (here and here) introduced three Python-based alternatives for expressing simulation process or agent-based behavior that takes place over (or blocks for) simulated time: Generators, a core feature of the Python language Greenlets, provided by the third-party greenlet Python extension package Tasklets, provided by two standalone Python distributions, Stackless Python and PyPy&#8230;&nbsp;(<a href=\"https:\/\/www.howardklein.net\/?p=52\">read more<\/a>)<\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"nf_dc_page":"","_jetpack_newsletter_access":"","_jetpack_dont_email_post_to_subs":false,"_jetpack_newsletter_tier_id":0,"_jetpack_memberships_contains_paywalled_content":false,"_jetpack_feature_clip_id":0,"_jetpack_memberships_contains_paid_content":false,"footnotes":"","jetpack_post_was_ever_published":false},"categories":[3,5,2],"tags":[],"class_list":["post-52","post","type-post","status-publish","format-standard","hentry","category-modeling","category-python","category-simulation"],"jetpack_sharing_enabled":true,"jetpack_featured_media_url":"","_links":{"self":[{"href":"https:\/\/www.howardklein.net\/index.php?rest_route=\/wp\/v2\/posts\/52","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/www.howardklein.net\/index.php?rest_route=\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/www.howardklein.net\/index.php?rest_route=\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/www.howardklein.net\/index.php?rest_route=\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/www.howardklein.net\/index.php?rest_route=%2Fwp%2Fv2%2Fcomments&post=52"}],"version-history":[{"count":25,"href":"https:\/\/www.howardklein.net\/index.php?rest_route=\/wp\/v2\/posts\/52\/revisions"}],"predecessor-version":[{"id":102,"href":"https:\/\/www.howardklein.net\/index.php?rest_route=\/wp\/v2\/posts\/52\/revisions\/102"}],"wp:attachment":[{"href":"https:\/\/www.howardklein.net\/index.php?rest_route=%2Fwp%2Fv2%2Fmedia&parent=52"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/www.howardklein.net\/index.php?rest_route=%2Fwp%2Fv2%2Fcategories&post=52"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/www.howardklein.net\/index.php?rest_route=%2Fwp%2Fv2%2Ftags&post=52"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}