*guest post by David Tanzer*

Petri nets are a simple model of computation with a range of modelling applications that include chemical reaction networks, manufacturing processes, and population biology dynamics. They are graphs through which entities flow and have their types transformed. They also have far-reaching mathematical properties which are the subject of an extensive literature. See the network theory series here on Azimuth for a panoramic and well-diagrammed introduction to the theory and applications of Petri nets.

In this article, I will introduce Petri nets and give an expository presentation of a program to simulate them. I am hoping that this will be of interest to anyone who wants to see how scientific concepts can be formulated and actualized in computer software. Here the simulator program will be applied to a “toy model” of a chemical reaction network for the synthesis and dissociation of H_{2}O molecules. The source code for this program should give the reader some “clay” to work with.

This material will prepare you for the further topic of *stochastic* Petri nets, where the sequencing of the Petri net is probabilistically controlled.

### Definition of Petri nets

A **Petri net** is a diagram with two kinds of nodes: container nodes, called **species** (or “places”, “states”), which can hold zero or more **tokens**, and **transition** nodes, which are “wired” to some containers called its **inputs** and some called its **outputs**. A transition can have multiple input or output connections to a single container.

When a transition **fires**, it removes one token from each input container, and adds one token to each output container. When there are multiple inputs or outputs to the same container, then that many tokens get removed or added.

The total state of the Petri net is described by a **labelling** function which maps each container to the number of tokens it holds. A transition is **enabled** to fire in a labelling if there are enough tokens at its input containers. If no transitions are enabled, the it is **halted**.

The sequencing is non-deterministic, because in a given labelling there may be several transitions that are enabled. Dataflow arises whenever one transition sends tokens to a container that is read by another transition.

Petri nets represent *entity conversion networks*, which consist of entities of different species, along with reactions that convert input entities to output entities. Each entity is symbolized by a token in the net, and all the tokens for a species are grouped into an associated container. The conversion of entities is represented by the transitions that transform tokens.

### Example 1: Disease processes

Here’s an example discussed earlier on Azimuth. It describes the virus that causes AIDS. The species are **healthy** cell, **infected** cell, and **virion** (the technical term for an individual virus). The transitions are for infection, production of healthy cells, reproduction of virions within an infected cell, death of healthy cells, death of infected cells, and death of virions.

Here the species are yellow circles and the transitions are aqua squares. Note that there are three transitions called “death” and two called “production.” They are disambiguated by a Greek letter suffix.

• **production** () describes the birth of one healthy cell, so it has no input and one **healthy** as output.

• **death** () has one **healthy** as input and no output.

• **death** () has one **infected** as input and no output.

• **death** () has one **virion** as input and no output.

• **infection** () takes one **healthy** and one **virion** as input, and has one **infected** cell as output.

• **production** () describes the reproduction of the virus in an infected cell, so it has one **infected** as input, and one **infected** and one **virion** as output.

### Example 2: Population dynamics

Here the tokens represent organisms, and the species are biological species. This simple example involving two species, **wolf** and **rabbit**:

There are three transitions: **birth**, which inputs one **rabbit** and outputs **rabbit** + **rabbit** like asexual reproduction), **predation**, which converts **rabbit** plus **wolf** to **wolf** plus **wolf**, and **death**, which inputs one **wolf** and outputs nothing.

### Example 3: Chemical reaction networks

Here the entities are chemical units: molecules, isolated atoms, or ions. Each chemical unit is represented by a token in the net, and a container holds all of the tokens for a chemical species. Chemical reactions are then modeled as transitions that consume input tokens and generate output tokens.

We will be using the following simplified model for *water formation and dissociation*:

• The species are H, O, and H_{2}O.

• Transition **combine** inputs two H atoms and one O atom, and outputs an H_{2}O molecule.

• Transition **split** is the reverse process, which inputs one H_{2} and outputs two H and one O.

Note that the model is intended to show the *idea* of chemical reaction Petri nets, so it is not physically realistic. The actual reactions involve H_{2} and O_{2}, and there are intermediate transitions as well. For more details, see Part 3 of the network theory series.

### A Petri net simulator

The following Python script will simulate a Petri net, using parameters that describe the species, the transitions, and the initial labelling. It will run the net for a specified number of steps. At each step it chooses randomly among the enabled transitions, fires it, and prints the new labelling on the console.

#### Download

Here is the first Petri net simulator.

#### Running the script

The model parameters are already coded into the script. So let’s give it a whirl:

python petri1.py

This produced the output:

H, O, H2O, Transition 5, 3, 4, split 7, 4, 3, split 9, 5, 2, combine 7, 4, 3, combine 5, 3, 4, split 7, 4, 3, split 9, 5, 2, split 11, 6, 1, combine 9, 5, 2, split 11, 6, 1, split 13, 7, 0, ...

We started out in a state with 5 H’s, 3 O’s, and 4 H_{2}O’s, then a split took place, which increased H by 2, increased O by 1, and decreased H_{2}O by one, then…

Running it again gives a different execution sequence.

#### Software structures used in the program

Before performing a full dissection of the code, I’d like to make a digression to discuss some of the basic software constructs that are used in this program. This is directed in particular to those of you who are coming from the science side, and are interested in learning more about programming.

This program exercises a few of the basic mechanisms of object-oriented and functional programming. The Petri net logic is bundled into the definition of a single *class* called PetriNet, and each PetriNet *object* is an *instance* of the class. This logic is grouped into *methods*, which are functions whose first parameter is bound to an instance of the class.

For example, here is the method IsHalted of the class PetriNet:

def IsHalted(this): return len(this.EnabledTransitions()) == 0

The first parameter, conventionally called “this” or “self,” refers to the PetriNet class instance. len returns the length of a list, and EnabledTransitions is a method that returns the list of transitions enabled in the current labelling.

Here is the syntax for using the method:

if petriNetInstance.IsHalted(): ...

If it were the case that IsHalted took additional parameters, the definition would look like this:

def IsHalted(this, arg1, ..., argn):

and the call would look like this:

if petriNetInstance.IsHalted(arg1, ..., argn)

Here is the definition of EnabledTransitions, which shows a basic element of functional programming:

def EnabledTransitions(this): return filter(lambda transition: transition.IsEnabled(this.labelling), this.transitions)

A PetriNet instance holds this list of Transition objects, called this.transitions. The expression “lambda: transition …” is an **anonymous function** that maps a Transition object to the boolean result returned by calling the IsEnabled method on that transition, given the current labelling this.labelling. The filter function takes an anonymous boolean-valued function and a list of objects, and returns the sublist consisting of those objects that satisfy this predicate.

The program also gives an example of *class inheritance*, because the class PetriNet inherits all fields and methods from a base class called PetriNetBase.

#### Python as a language for pedagogical and scientific programming

We continue with our digression, to talk now about the choice of language. With something as fundamental to our thinking as a language, it’s easy to see how the topic of language choice tends to become partisan. Imagine, for example, a debate between different nationalities, about which language was more beautiful, logical, etc. Here the issue is put into perspective by David Tweed from the Azimuth project:

Programming languages are a lot like “motor vehicles”: a family car has different trade-offs to a sports car to a small van to a big truck to a motorbike to …. Each of these has their own niche where they make the most sense to use.

The Petri net simulator which I wrote in Python, and which will soon to be dissected, could have been equivalently written in any modern object-oriented programming language, such as C++, Java or C#. I chose Python for the following reasons. First, it is a scripting language. This means that the casual user need not get involved with a compilation step that is preparatory to running the program. Second, it does provide the abstraction capabilities needed to express the Petri net model. Third, I like the way that the syntax rolls out. This syntax has been described as executable pseudo-code. Finally, the language has a medium-sized user-base and support community.

Python is good for proof of concept programs, and it works well as a pedagogical programming language. Yet it is also practical. It has been widely adopted in the field of scientific programming. David Tweed explains it in these terms:

I think a big part of Python’s use comes from an the way that a lot of scientific programming is as much about “management” — parsing files of prexisting structured data, iterating over data, controlling network connections, calling C code, launching subprograms, etc — as much as “numeric computation”. This is where Python is particularly good, and it’s now acquired the NumPy & SciPy extensions to speed up the numerical programming elements, but it’s primarily the higher level elements that make it attractive in science.

Because variables in Python are neither declared in the program text, nor inferred by the byte-code compiler, the types of the variables are known only at run time. This has a negative impact on performance for data intensive calculations: rather than having the compiler generate the right code for the data type that is being processed, the types need to be checked at run time. The NumPy and SciPy libraries address this by providing bulk array operations, e.g., addition of two matrices, in a native code library that is integrated with the Python runtime environment. If this framework does not suffice for your high-performance numerical application, then you will have to turn to other languages, notably, C++.

All this being said, it is now time to return to the main topic of this article, which is the *content* of Petri net programming.

#### Top-level structure of the Petri net simulator script

At a high level, the script constructs a Petri net, constructs the initial labelling, and then runs the simulation for a given number of steps.

First construct the Petri net:

# combine: 2H + 1O -> 1H2O combineSpec = ("combine", [["H",2],["O",1]], [["H2O",1]]) # split: 1H2O -> 2H + 1O splitSpec = ("split", [["H2O",1]], [["H",2],["O",1]]) petriNet = PetriNet( ["H","O","H2O"], # species [combineSpec,splitSpec] # transitions )

Then establish the initial labelling:

initialLabelling = {"H":5, "O":3, "H2O":4}

Then run it for twenty steps:

steps = 20 petriNet.RunSimulation(steps, initialLabelling)

#### Program code

Species will be represented simply by their names, i.e., as strings. PetriNet and Transition will be defined as object classes. Each PetriNet instance holds a list of species names, a list of transition names, a list of Transition objects, and the current labelling. The labelling is a dictionary from species name to token count. Each Transition object contains an input dictionary and an input dictionary. The input dictionary map the name of a species to the number of times that the transition takes it as input, and similarly for the output dictionary.

##### Class PetriNet

The PetriNet class has a top-level method called RunSimulation, which makes repeated calls to FireOneRule. FireOneRule obtains the list of enabled transitions, chooses one randomly, and fires it. This is accomplished by the method SelectRandom, which uses a random integer between 1 and N to choose a transition from the list of enabled transitions.

class PetriNet(PetriNetBase): def RunSimulation(this, iterations, initialLabelling): this.PrintHeader() # prints e.g. "H, O, H2O" this.labelling = initialLabelling this.PrintLabelling() # prints e.g. "3, 5, 2" for i in range(iterations): if this.IsHalted(): print "halted" return else: this.FireOneRule() this.PrintLabelling(); print "iterations completed" def EnabledTransitions(this): return filter(lambda transition: transition.IsEnabled(this.labelling), this.transitions) def IsHalted(this): return len(this.EnabledTransitions()) == 0 def FireOneRule(this): this.SelectRandom(this.EnabledTransitions()).Fire (this.labelling) def SelectRandom(this, items): randomIndex = randrange(len(items)) return items[randomIndex]

##### Class Transition

The Transition class exposes two key methods. IsEnabled takes a labeling as parameter, and returns a boolean saying whether the transition can fire. This is determined by comparing the input map for the transition with the token counts in the labeling, to see if there is sufficient tokens for it to fire. The Fire method takes a labelling in, and updates the counts in it to reflect the action of removing input tokens and creating output tokens.

class Transition: # Fields: # transitionName # inputMap: speciesName -> inputCount # outputMap: speciesName -> outputCount # constructor def __init__(this, transitionName): this.transitionName = transitionName this.inputMap = {} this.outputMap = {} def IsEnabled(this, labelling): for inputSpecies in this.inputMap.keys(): if labelling[inputSpecies] < this.inputMap[inputSpecies]: return False # not enough tokens return True # good to go def Fire(this, labelling): print this.transitionName for inputName in this.inputMap.keys(): labelling[inputName] = labelling[inputName] - this.inputMap[inputName] for outputName in this.outputMap.keys(): labelling[outputName] = labelling[outputName] + this.outputMap[outputName]

##### Class PetriNetBase

Notice that the class line for PetriNet declares that it inherits from a base class PetriNetBase. The base class contains utility methods that support PetriNet: PrintHeader, PrintLabelling, SelectRandom, and the constructor, which converts the transition specifications into Transition objects.

class PetriNetBase: # Fields: # speciesNames # Transition list # labelling: speciesName -> token count # constructor def __init__(this, speciesNames, transitionSpecs): this.speciesNames = speciesNames this.transitions = this.BuildTransitions(transitionSpecs) def BuildTransitions(this, transitionSpecs): transitions = [] for (transitionName, inputSpecs, outputSpecs) in transitionSpecs: transition = Transition(transitionName) for degreeSpec in inputSpecs: this.SetDegree(transition.inputMap, degreeSpec) for degreeSpec in outputSpecs: this.SetDegree(transition.outputMap, degreeSpec) transitions.append(transition) return transitions def SetDegree(this, dictionary, degreeSpec): speciesName = degreeSpec[0] if len(degreeSpec) == 2: degree = degreeSpec[1] else: degree = 1 dictionary[speciesName] = degree def PrintHeader(this): print string.join(this.speciesNames, ", ") + ", Transition" def PrintLabelling(this): for speciesName in this.speciesNames: print str(this.labelling[speciesName]) + ",",

### Summary

We’ve learned about an important model for computation, called Petri nets, and seen how it can be used to model a general class of entity conversion networks, which include chemical reaction networks as a major case. Then we wrote a program to simulate them, and applied it to a simple model for formation and dissociation of water molecules.

This is a good beginning, but observe the following limitation of our current program: it just randomly picks a rule. When we run the program, the simulation makes a kind of random walk, back and forth, between the states of full synthesis and full dissociation. But in a real chemical system, the rates at which the transitions fires are _probabilistically determined_, and depend, among other things, on the temperature. With a high probability for formation, and a low probability for dissociation, we would expect the system to reach an equilibrium state in which H_{2}O is the predominant “token” in the system. The relative concentration of the various chemical species would be determined by the relative firing rates of the various transitions.

This gives motivation for the next article that I am writing, on *stochastic* Petri nets.

### Programming exercises

1. Extend the constructor for PetriNet to accept a dictionary from transitionName to a number, which will give the relative probability of that transition firing. Modify the firing rule logic to use these values. This is a step in the direction of the stochastic Petri nets covered in the next article.

2. Make a prediction about the overall evolution of the system, given fixed probabilities for synthesis and dissociation, and then run the program to see if your prediction is confirmed.

3. If you like, improve the usability of the script, by passing the model parameters from the command line. Use sys.argv and eval. You can use single quotes, and pass a string like “{‘H':5, ‘O':3, ‘H2O':4}” from the command line.

By the way: if you do this exercises, please post a comment including your code, or a link to your code! To make Python code look pretty on this blog, see this.

### Acknowledgements

Thanks to David Tweed and Ken Webb for reviewing the code and coming up with improvements to the specification syntax for the Petri nets. Thanks to Jacob Biamonte for the diagrams, and for energetically reviewing the article. John Baez contributed the three splendid examples of Petri nets. He also encouraged me along the way, and provided good support as an editor.

### Appendix: Notes on the language interpreter

The sample program is written in Python, which is a low-fuss scripting language with abstraction capabilities. The language is well-suited for proof-of-concept programming, and it has a medium-sized user base and support community. The main website is www.python.org.

You have a few distributions to choose from:

• If you are on a Linux or Mac type of system, it is likely to already be installed. Just open a shell and type “python.” Otherwise, use the package manager to install it.

• In Windows, you can use the version from the python.org web site. Alternatively, install cygwin, and in the installer, select Python, along with any of your other favorite Linux-style packages. Then you can partially pretend that you are working on a Linux type of system.

• The Enthought Python distribution comes packaged with a suite of Python-based open-source scientific and numeric computing packages. This includes ipython, scipy, numpy and matplotlib.

The program is distributed as a self-contained script, so in Linux and cygwin you can just execute ./petri1.py. When running this way under cygwin, you can either refer to the cygwin version of the Python executable, or to any of the native Windows versions of interpreter. You just need to adjust the first line of the script to point to the python executable.

Great post, David! I hope we pick up enough steam to write some stochastic Petri net programs that tackle really interesting real-world issues, like:

• how likely is it that a random stochastic Petri net displays bistability, like an ‘on-off switch’?

or

• how likely is it that a random stochastic Petri net displays periodic behavior, like a ‘biological clock’?

or

• what rules for

evolutionof stochastic Petri nets push them toward interesting behavior?or

• take some stochastic Petri nets that are studied in the biology / chemistry literature, and answer some interesting open questions about them.

But for now, here’s a naive programming question: what’s an ‘anonymous function’?

Anonymous functions are what most functions always were before we forgot they needn’t all have names.

Jim, below suggests lambda-abstracts as one way to get anonymous functions. Strictly speaking, any composite term that can be used as a “function” in the program-writing sense might as well be considered an anonymous function.

Languages like C tend to discourage function anonymity, and they get around the attendent difficulties (if they are difficulties) with things like “callbacks” and clever manipulation of machine state. C particularly is designed with a byte-chunked RAM machine in mind, and doing things directly to that machine; the things that are used as functions in C are thought of as addresses where a particular bunch of machine instructions start. These addresses can be used as parameters in calling other functions, or returned as values, but returned values won’t be “new” functions or combinates.

Other languages, especially those with a “functional programming” style are more encouraging of anonymous functions; lisp comes to mind particularly. Of course, one can compile lisp into C (e.g.), so it ought to be possible to write in a functional style in C, directly; though the things that would be functions in lisp probably wouldn’t be C functions.

Interesting questions about random stochastic Petri nets that you posed. They lead to further questions.

1. By what method would we construct a random Petri net? I see definitions on the web of a random graph, where you fix the set of nodes, and then assign a probability to each possible edge — then use these probabilities to build a graph. With Petri nets we have two types of nodes, plus the possibility that an edge from A to B occurs multiple times.

The following questions will betray my lack of knowledge of the further reaches of Petri net theory. I’m guessing my way through the terrain, so please correct and/oor elaborate.

2. Bistable and periodic refer to the solutions to the rate equation, right? We could think of the state space of a Petri net with k species either as the non-negative region of or , depending on whether you use a continuous approximation or not. In the continuous case we use a differential equation, the rate equation. ( Is there a discrete analog?)

The solutions to the equations give the trajectories in state space, given an initial state. So we can look at attractor points and basins of attraction.

Does bistable mean that there are two attractors and the whole (non-negative) region of is divided into the basins of attraction for these points?

Does periodic means that there are points that cycle back to themselves, or is it a broader notion, meaning that there are points whose trajectories are bounded, but do not converge to an attractor?

Can a system be bistable and periodic, in the sense that there are two attractors, but there are also points that are periodic?

Basically I’m asking about the character of the solutions to the rate equation, which I have not yet investigated. Can you summarize the main theorems regarding the solutions.

3. Given a stochastic Petri net, is there a decision procedure to say whether it is “bistable,” or “periodic”? Or does one have to resort to a probabilistic test, involving selections of an initial starting vector, and testing to see whether (1) it appears to be headed for an attractor (or infinity), or (2) whether it appears to be a periodic point. Clearly “appears” needs to be defined.

I realize these are somewhat broad and sketchy questions, but any expert hints you can provide would be great!

Thanks

David Tanzer asked:

Right.

Right. My book with Jacob Biamonte wound up spending a lot of time on the rate equation. We explained that when a Petri net has ‘deficiency zero’ and is ‘weakly reversible’, people have a good understanding of solutions of the rate equation. There are just as many equilibrium solutions as you’d expect, no more and no less. They’re all stable, and there are no periodic solutions.

Even better, in this case, the ‘Anderson–Craciun–Kurtz theorem’ lets you quickly obtain equilibrium solutions of the

masterequation from equilibrium solutions of the rate equation.But this case is, in a sense, the

boringcase!Yes, but people haven’t studied this as much. In fact, I don’t know if I’ve seen it anywhere, but I could easily describe it. It might be fun to generalize the zero deficiency theorem to this case, if nobody has yet. I think I could do it.

What I’m calling a ‘stable equilibrium’ is exactly the same as what you’re calling an attractor. When the ‘zero deficiency’ condition holds, we have a good handle on the attractors, and there’s only one attractor in each ‘stoichiometric compatibility class’.

No, but close. It means there are two attractors in some stoichiometric compability class.

I’m reluctant to explain ‘stoichiometric compatibility class’ in general, having already done so here, but if you think about the water formation and dissociation example you described, you’ll quickly get the idea.

There’s one stoichiometric compatibility class for each choice of these two quantities:

1) number of H

_{2}O’s plus the number of O’s2) twice the number of H

_{2}O’s plus the number of H’sThese quantities can never change, since hydrogen and oxygen atoms can’t be created or destroyed by the reactions in this Petri net. So, a solution of the rate equation

ormaster equation can only wander around in one stoichiometric compatibility class.Thanks to the deficiency zero theorem, the rate equation has one attractor in each stoichometric compatibility class. Where this attractor is depends on the rate constants for the formation and dissociation of water.

This is a nice example of the zero deficiency theorem.

Yes. The time it takes to come back is called the ‘period’.

No, that broader notion would includes not just periodic solutions but quasiperiodic and chaotic ones.

What’s periodic is not the system but a particular solution. A system can have any number of attractors and any number of periodic solutions.

The word ‘bistable’ focuses undo attention on the number

two, which is interesting just because it’s the first number bigger than one. A light switch is bistable, and people are interested in ‘switches’ in biochemistry, but a general switch could have many stable settings.I summarized the deficiency zero theorem and Anderson-Craciun-Kurtz theorem rather vaguely above. I did it quite precisely in the book. There are also lots of other theorems, like the deficiency one theorem.

There are theorems—most notably the deficiency zero theorem and deficiency one theorem–that give necessary or sufficient theorems for these properties, but I don’t believe a full-fledged decision procedure is known.

People can already do a lot better than this brute-force method, thanks in part to all the theorems people have shown, but I think there’s a huge unexplored territory here. Mathematical chemists are interested in this problem, and it seems to be tricky.

A naive programming answer :

Creating an anonymous function is like cooking a cake in a pastry shop. The customer is the executor of the function. He will refer to it by name (like myPrecious, theYummy) and he uses it as he wishes (eating it the same day, eating half tomorrow and half today, eating it now, giving it as a present).

The person who cooks it will not use it so there is no reason to give it a name.

Some more :

A pure anonymous function has no side effects. If the cooker has forgotten his keys inside he will not be able to enter his house after selling it. This cake (anonymous function) is not pure. But if it has no side effects (after giving it away nothing happens) then it is pure.

Another way of putting it is the difference between someone you introduce yourself to, say a new acquaintance and the person behind the till in a shop you’ll visit once. There’s no reason you can’t avoid asking the name of an acquaintance,but it makes things harder,and likewise you could find the name of the server,but its not clear it’s worth it even for them as they’re probably being tracked on throughout. Sometimes naming isn’t worth doing,even in programs.

The lambda in lambda calculus is an anonymous function except for its name “lambda”

An example is:

\x -> x**2 where x is a variable and \ is a lambda..

producer input (fine) – process – output

producer input (virus) – process – output

reproducer input (fine) – process – output

reproducer input (virus) – process – output

repeating…

basically input (fine)- process wont be infected- output(fine)

input (virus) – – process will be infected- output(virus)

Can you elaborate a bit in words the thought process that you are going though here?

Okay, thanks everyone: I think I get the idea of an anonymous function. If I understand right,

is ‘nonymous’ (my silly term), since we’ve given it the name , but

is an ‘anonymous’ description of the same function in lambda-calculus terminology, and lots of mathematicians would use the anonymous description

I think you may be missing something here. Your discussion gets quite deep into the problems of building a petri-net simulator precisely in Python, but reading it, it is not clear to me that the concepts you find in Python are particularly relevant for understanding Petri-nets (this does not mean that you can’t build a good petri-net simulator in Python, only that what is interesting is the result, not the details of the implementation).

When I think about how to implement a Petri-net simulator (speaking as someone who first encountered them at least 30 years ago, when I was about 17) I don’t think in terms of object/class structures at all, but in terms of possible paths through a relational algebra. (I had the same feeling reading James Sethna’s book on statistical physics a while ago, where he provides a lot of class-based (python) code for working with graphs, when sparse matrixes and an algebraic approach would have been much more effective).

The result would be more concise and map much better onto the problem, so that intutions on one side could be reinforced on the other. I’m sceptical that this is the case with something built around the limitations of Python.

Object-orientation is a bit like the gibbs vector calculus: it works and it is better than some of the alternatives, but algebraic geometry can do it better!

[you know, blog comment boxes are too short for this sort of point, it suddenly occurs to me]

Hi Sean,

Thanks for your reflective comments. I’d be very interested to hear more about your relational-algebra interpretation of Petri net structures.

In terms of the article itself, I do not believe that I have missed the mark in terms of the aim that I set out, which was to address a three-fold audience of programming people who want to learn about science applications, science people who want to learn about programming, and any other curious folks. I wrote it as, for lack of a better term, a “praxis paper,” in the sense of practical application of a field of study. So first I introduced the concepts, and then moved to a pragmatic construction stage. For the language, I choose one that is both capable of abstraction, and which has a substantial application base. To give the idea of this genre in simpler terms, consider an article that introduces regular polyhedra, and the shows how to build them with common household materials.

I would further maintain that this “household materials” programming language can faithfully represent the mathematical concepts of many application domains, which are presented in terms of objects, functions and relations. Take for example a definition of a Petri net. A Petri net P consists of a pair (S,T) where S is a set of species names, and T is a set of transitions. Each T is defined by a pair (I,O), where I: S -> N is the input degree map, and O: S -> N is the output degree map. A labelling L: S -> N is a map that counts the number of tokens for each species. A transition T is enabled in a labelling L if for each s in S, L(s) >= InputMap(T)(s).

The definition itself introduces nouns, which means, in the first place, objects. Functions are involved in the definition. This one happens not to involve any relations. So any language that allows us to work with objects, functions, and relations is a Useful Thing. If you look at the definitions in the program I gave, sure there is some syntactic baggage, but if you look through that, you will see an isomorphic copy of the tuples and functions used in the above mathematical definition. Consider, for example, the method IsEnabled of a Transition. It is precisely a boolean valued function of a Transition (with a parameter that is conventionally called “this” or “self”) and a labelling, and the body of the function says exactly what I said above in natural language.

The beauty of this kind of formal language is that it is mathematics that runs.

Now it is true, as you pointed out, that I gave some discussion of Python per se, which is incidental to the topic of Petri net programming. I was aware of this, and that is why I introduced those sections as part of a “digression.” I put this in to get as many readers on board as possible. Since the program text was to be presented in Python, if only because that is a native language of the author, then why not prep the uninitiated. Finally, in the spirit of hands-on and “literate” coding, we _do_ want to present the code and dissect it. This is not merely an “implementation.” We’ve learned something about a theory, and now we want to learn the method of translating this theory into something that can do cool things. We look _into_ the code, to see the isomorphic copy of theory.

* * *

Your view of Petri nets in terms of “possible paths through a relational algebra,” is intriguing, but I’d need to hear more about it to really understand what you mean. Let’s suppose that your conceptualization — which undoubtedly is based on far more experience with Petri nets than I have — gives some kind of structural insight in the semantics of Petri nets, which can be put to good use in simulator. I have no reason to doubt this.

Here is an analogy to explain the relation of your approach to mine. Suppose we gave an introductory “praxis paper” on Abelian group theory. A group is first defined, an then an “Abelian” group is introduced as one which has the constraint that the group operation is commutative. After giving some examples, we could then introduce some theory-isomorphic code to test if the group operator is commutative, perform multiplications, using the simplifications that come along with commutativity, etc. This could be some useful pedagogical code, that readers could tinker with. Now, based on further analysis of the domain, along comes the Structure Theorem for finite Abelian groups, which gives a powerful, simplifying representation for such groups: they are isomorphic to direct products of cyclic groups. This leads to a more incisive representation for the groups, and new an better data structures to be used in a computational system. But this latter advance doesn’t negate the pedagogical or experimental value of the original code, which is an executable expression of group theory, using the general terms of group theory.

Furthermore, once the new theory based on cyclic groups is put forward, it would be very natural then to translate these definitions into isomorphic forms in the household computing language, and use that for an expository presentation of the more advanced simulator.

So if you share with us the definitions behind your relational-algebraic approach to Petri net computation, we might ask for your permission to express it in a household abstract programming language and blog about it — with due credit to you of course!

Seriously though, in addition to giving us some more color on your idea, check out the discussions taking place on the Forum of the Azimuth Project. There is a lot of ongoing, active interest in Petri nets, and their potential application to areas such as climate modelling. (Theme question: do we live in a bistable climate?)

It’s clear that you have an experienced eye for this problem domain, and if you are interested, you can make a real contribution to our ongoing research. To become a posting member of the group only takes a few steps.

Regards,

Dave

P.S. – This invitation is generally extended to anyone who wants to discuss and work on the computational, mathematical and scientific challenges associated with environmental problems. This is a focal point, but is not taken in a restrictive sense. For example, we have recently been discussing Qinglan Xia’s gorgeous paper “The Formation of a Tree Leaf,” at http://www.ma.utexas.edu/users/qlxia/Research/leaf.pdf, which presents a simple physical model that accounts for the shape of tree leaves, in terms of an optimization process that seeks to minimize transport costs. from the root to the cells.

Hmmm… My point wasn’t really about petri nets, and more about the

relevance of the specific programming language. However, now that I

am confronted with the problem, I suppose I should put my money where

my mouth is. How would I implement petri-nets?

Call it Graham’s principle that the structure visible in a program

should correspond to the problem, and not to the programming

language. The ideal for this is usually Lisp (for which see Paul

Graham passim), or, say, for large scale linear algebra, APL (alas

only once upon a time) or (today) Matlab. For working with relations

on tuples I think Prolog is likely to provide a clean map between

concrete and abstract.

A petri-net is a bipartite graph G plus a set of counters on one of

the types of nodes, T, so that the vector of these counters defines a

state. G defines a relation R, of transitions from state to state;

thus a particular run of a simulation of a petri-net has the form

T1 -R-> T2 -R-> T3 -R-> T4 -R-> …

Prolog allows me to model tuple transitions very easily. If I model

numbers ‘peano-style’, as s(s(…(0))), and T has the form

t(n1, n2, n3, ….)

Then I can model a rule directly in Prolog as:

r(rulename,

t(n1.old, n2.old, n3.old, ….),

t(n1.new, n2.new, n3.new, ….)).

which is simply a declared relation between states (I’ve also included

the name – which gives me a bit more control and documentation) and I

am done.

So how does this work in practice?

Well I can write down the the ‘Chemical reaction’ net as:

% t(O, H, H2O).

r(split,

t( O , H , s(H2O)),

t(s(O), s(s(H)), H2O ).

r(combine,

t(s(O), s(s(H)), H2O ),

t( O , H , s(H2O)).

Note that rules are directly declarative and I have completely

separated them from any other interpretive framework – I can modify my

non-deterministic model, or implement breadth first or bounded

depth-first search on the space easily as extensions. I can even run

my model backwards, should I so want.

A basic simulator (with a bound on the number of transitions) then

looks like:

simulate(T, T, [], 0).

simulate(Tinit, Tfinal, [(Tinit, R) | L], N) :-

choose(R),

r(R, Tinit, Tnext),

PrecN is N – 1,

simulate(Tnext, Tfinal, L, PrecN).

simulate(T, T, [], N) :-

not r(_, T, _),

N > 0.

choose is a bit more messy, since it should enumerate the possible

rules in random order. A state specific version (we can do much better

– more general – with a bit more space) would be

choose(R) :-

(0 =:= rand(2) ->

(X = split ; X = combine)

| (X = combine; X = split ).

A call to the simulator, with a bound of 20 transitions, is then:

simulate(t(s(s(s(s(s(0))))), s(s(s(0))), s(s(s(s(0))))), Tfinal, L, 20).

Not sure there is enough code here to count as a program (<= 19

lines). I don't have a prolog interpreter to hand – my employers are

very unenthusiastic about my installing my own software, so there may

be typos above, but I hope the idea is clear.

Note the argument here is not in favor of Prolog (a similar-flavored

– though not quite identical and not quite so concise – implementation

would be easy in a language like Scheme), but in favor of choosing

your tools so as to move you close to the problem rather than trying

to move the problem closer to your tools, and thus, in this case, no

objects/classes; if you were designing a user-interface, then an

object/class model would be close to your problem.

How does this compare to something like the Stochastic Pi Machine (http://research.microsoft.com/en-us/um/cambridge/groups/science/tools/spim/spim.htm)?

Keith,

Thank you for the link to this magnificent program!

For the comparison, right off the bat we have that:

My program is a proof-of-concept implementation of a simulator for Petri nets. It is meant to teach the concepts, and to serve as a kind of “clay” that actually works.

That program is a full-fledged, incredibly well presented implementation of a simulation system for reactions that use the formalism of pi-calculus, not Petri nets.

I know next to nothing about the pi-calculus, but now that I see what it can do, I will learn more about it!!

Can anyone recommend a good tutorial introduction to the pi-calculus? I find the Wikipedia presentation to be formalistic and unintuitive. Thanks.

There’s surely better actual tutorials, but here’s a lightning introduction by Mike Stay.

Wow. I sort of get the idea, but not well.

If the stochastic pi-machine is pi-calculus in action, then that is proof that the calculus is cool and useful.

So I’m trying to make sense of the code that have for the sample reactions. Here a document that explains the some of the examples in the Sliverlight reaction simulator:

http://research.microsoft.com/en-us/projects/spim/gillespie.pdf

It includes radioactive decay, and Lotka reactions. Here is the Lotka reaction logic:

directive sample 5.0 1000

directive plot Y()

new c1@5.0:chan

new c2@0.0025:chan

let X() = ?c1; X()

let Y() =

do !c1; (Y() | Y())

or !c2

or ?c2

run (X() | 10 of Y())

The following questions are directed generally to anyone here. How much of this is pure pi-calculus, and how much is in the stochastic extension? The channels have rate constants associated with them, but is that the extent of the extension?

Can anyone give a detailed analysis and explanation of this code, to explain how it is acheiving the marvelous simulation result? Remember, I don’t really know what pi-calculus is, so you’ll have to start from a point of few assumptions.

Here is the the language definition document for stochastic pi-calculus:

http://research.microsoft.com/en-us/projects/spim/language.pdf

Thanks

The pi calculus is meant as a distilled framework for describing ‘concurrency’—roughly, computer networks where nodes can create channels to other nodes and send messages along these nodes—much as the lambda calculus is a distilled framework for describing computations done sequentially by a single processor.

Personally I find these distilled frameworks, that seek to get everything done with as few primitives as possible, less clear than a category-theoretic description of

large classesof frameworks. When you try to distill everything down to a few primitives you wind up making a lot of arbitrary choices, and the resulting framework tends to look cryptic.Of course, category theory also looks cryptic to those who haven’t studied it… but it has the advantage of being so generally useful that once you learn it and formulate some idea in this way, a bunch of connections to other ideas are instantly visible—ideas in math, logic, physics, and computation, for starters.

Right now I have just one grad student: Mike Stay, who works at Google. He’s working with me on categories and computation. In that post David Tweed mentions, he was endeavoring to explain the pi calculus to me, in response to my plea for help. As a result, he got interested in describing

By now he’s gone a lot further. See:

• Mike Stay, Higher categories for concurrency,

n-Category Café.for the state of his thinking a year and a half ago. Since then he’s been working on this subject with Jamie Vicary and making even more progress.

By now it’s clear to me that bigraphs, the pi calculus and other frameworks for concurrency would really profit from a category-theoretic description, and that when this is done it’ll look a bit like a

categorificationof Petri net theory, with symmetric monoidal bicategories replacing symmetric monoidal categories.In other words: instead of describing a bunch of ‘particles’ interacting in time, these fancier frameworks describe a bunch of

labelled graphsinteracting and changing topology in time. The vertices of the graph are what I called ‘nodes’, and the edges are what I called ‘channels’. The nodes can do things (e.g. compute) and send messages to each other along the edges.This sort of framework might be useful in biology and ecology, too. It would certainly be relevant to the ‘smart grid’.

While I don’t feel this way of running Petri nets is very useful, I used the opportunity of the proposed exercise to extend the engine for the case of transitions with relative rates. This, to better show the flexibility of Python as a base to create DSL – domain specific languages. In this case a DSL for petri nets, which for the three examples of the post, allows to encode them entirely as

The optional number and vertical bar at the beginning of a rule expresses a weight or unnormalized rate different from unity.

The code is written for Python 3.3 but except for the use of greek letters in the “AIDS” petri net, should also run on Python 2.7. The implementation of the syntax in 65 LOCs is minimalistic, as no step is taken ensure graceful errors when the DSL code strays from the syntax illustrated.

The way to invoke the code is

And here is the code. I just hope the sourcecode quoting works as advertised :)

Cool. Can you say a bit about how these DSL language constructs are working together to make this program run?

I.e., How does it work? I am not familiar with these 2.7 language constructs.

Looking it over some more, I’m getting some hazy pciture of it.. The syntax “@petrinet ___” is a way to pass the defined set of transitions to the function “petrinet.”

What is the type of the defined objects like “chemical_reaction.” Are they something like rulesets? And what is the general role of the symbols |, _ and >>. Yet you call the object AIDS like a function.

Above and beyond these language specifics, can you post here a small amount of dissection of the code? Saying essentially, here is what I wanted to achieve, and these are the mechanisms I used, which work together correctly because of XYZ. Thanks.

The only thing that’s 2.7 is the Counter library class (itself a subclass of dict) which saves me some definitions, such as that of the & operator by a special method of the user-defined subclass Bag. That subclass itself (through special __methods__ that are the Python way of overloading operators) serves to implement the syntax for transitions. This is comparable in power and limitations to making up user syntax by defining operators in Prolog. “Practicality beats purity” is a motto of Python, and this is reflected in using that single general-purpose subclass Bag to represent at once species, transition inputs, transitions outputs, and transitions (given this, the name of the subclass is debatable – as it doesn’t properly fit the last specialization).

The implementation is a hack, and can be largely characterized as using reflection to re-target Python’s own compilation engine. The syntax @petrinet is what’s called a decorator (iirc it was borrowed from Java around python 2.4), and the way it works is that the result of compiling the function definition that follows, gets passed to the function petrinet, and the return value bound instead to the name of the defined function in the outer scope.

The function petrinet() takes apart the compiled function object which allows it to create bindings for all unbound variables meant to name species, to values that are singleton Bags with 1 exemplary of the name or species token counted in. The code of the function is then evaluated while the Bag class of these bound values provides appropriate definitions for the algebraic syntax of transitions. The most flaky is the cheap way transition names get associated with transitions; this is done by simply mapping evaluations of the >> operator in the sandboxed function evaluation, in order, to the names of local variables listed by the compiled function object. This boils down to the unenforced requirement that the decorated petri net function is formed of simple assignment statements to distinct variables, of the form

transition_name = (expression involving a single call to >>)

that is, with the RHS an expression involving exactly one evaluation of the >> operator. Also, the function must be declared with an empty parameter list or things will break.

The function petrirun is there to demonstrate that the intended semantics is captured (in a more realistic setting the procedural interpretation of the petri net would be better decoupled). petrirun() gets instantiated by petrinet() for each petri net as the executable object. It is actually a generator because it returns values with a yield statement, what allows to pull out of it, subsequent values on demand, for instance as the generator in a for-loop as demonstrated above. The list of transitions and the list of species are attached to it as user function attributes for allowing later user access. The petrirun instance should also be renamed to the decorated function name (to appear under that name when coming out in the REPL or tracebacks) but that was left out. petrirun() is also designed to exploit named parameters function call syntax to serve the stipulation of the initial labelling of the run.

In the previous article, I explored a simple computational model called Petri nets [...]