CBIO (CSCI) 4835/6835: Introduction to Computational Biology
Today, we'll wrap up our module on image processing with some more in-depth examples of how to analyze biological images. By the end of this lecture, you should be able:
What does "filter" mean to you?
Instagram?
Photoshop?
Both are technically "yes", as they use the same underlying principle.
A filter applies a convolution kernel to an image.
Although this sounds fancy, this is just a generalization of a very simple idea: rather than applying a specific function to one pixel, this uses the pixel's surrounding neighborhood to apply a function.
The kernel is represented by an $n$x$n$ matrix where the target pixel is in the center (so $n$ needs to be an odd number).
The output of the filter is the sum of the products of the matrix elements with the corresponding pixels they overlap.
Confused?
Here's what it looks like. The image is the bigger $5 \times 5$ matrix; the filter is the smaller $3 \times 3$ matrix:
$$ \begin{bmatrix} 1 & 0 & 1 \\ 0 & 1 & 0 \\ 1 & 0 & 1 \end{bmatrix} $$A convolution is basically a multiplication: the smaller matrix is multiplied element-wise with the overlapped entries of the larger matrix. These products are then all summed together into a new value for the pixel at the very center of the filter.
Then the filter is moved and the process repeats. This is true for any and all filters (or kernels).
What creates the specific effect, then--edge finding, blurring, sharpening, and so on--is the specific numbers in the filter. As such, there are some common filters:
| Identity | Blur | Edge Detection |
Wikipedia has a whole page#Details) on common convolutional filters:
Of course, you don't have to design the filter and code up the convolution yourself (though you could!). Most image processing packages have default versions of these filters included.
%matplotlib inline
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
from PIL import Image
from PIL import ImageFilter
img = Image.open("Lecture22/image1.png")
plt.imshow(img)
<matplotlib.image.AxesImage at 0x10d51d4a8>
Some of the PIL built-in filters include:
Here's BLUR:
import numpy as np
blurred = img.filter(ImageFilter.BLUR)
f = plt.figure(figsize = (12, 6))
f.add_subplot(1, 2, 1)
plt.imshow(np.array(img))
f.add_subplot(1, 2, 2)
plt.imshow(np.array(blurred))
<matplotlib.image.AxesImage at 0x10ddff0f0>
Here's SHARPEN:
sharpened = img.filter(ImageFilter.SHARPEN)
f = plt.figure(figsize = (12, 6))
f.add_subplot(1, 2, 1)
plt.imshow(np.array(img))
f.add_subplot(1, 2, 2)
plt.imshow(np.array(sharpened))
<matplotlib.image.AxesImage at 0x10e2e0f28>
And FIND_EDGES (more on this later):
edges = img.filter(ImageFilter.FIND_EDGES)
f = plt.figure(figsize = (12, 6))
f.add_subplot(1, 2, 1)
plt.imshow(np.array(img))
f.add_subplot(1, 2, 2)
plt.imshow(np.array(edges))
<matplotlib.image.AxesImage at 0x10da77550>
There are other filters that don't quite abide by the "multiply every corresponding element and then sum the products" rule of convolution.
In deep learning parlance, these filters are known as "pooling" operators, because while you still have a filter that slides over the image, you choose one of the pixel values within that filter instead of computing a function over all the pixel values.
As common examples, we have
Here's max-pooling:
max_pool = img.filter(ImageFilter.MaxFilter(5)) # This means a 5x5 filter
f = plt.figure(figsize = (12, 6))
f.add_subplot(1, 2, 1)
plt.imshow(np.array(img))
f.add_subplot(1, 2, 2)
plt.imshow(np.array(max_pool))
<matplotlib.image.AxesImage at 0x10fd934e0>
(just to help things, let's convert to grayscale and try that again)
img = img.convert("L")
# Same code as before.
max_pool = img.filter(ImageFilter.MaxFilter(5)) # This means a 5x5 filter
f = plt.figure(figsize = (12, 6))
f.add_subplot(1, 2, 1)
plt.imshow(np.array(img), cmap = "gray")
f.add_subplot(1, 2, 2)
plt.imshow(np.array(max_pool), cmap = "gray")
<matplotlib.image.AxesImage at 0x110805e10>
How about median pooling?
median_pool = img.filter(ImageFilter.MedianFilter(5)) # This means a 5x5 filter
f = plt.figure(figsize = (12, 6))
f.add_subplot(1, 2, 1)
plt.imshow(np.array(img), cmap = "gray")
f.add_subplot(1, 2, 2)
plt.imshow(np.array(median_pool), cmap = "gray")
<matplotlib.image.AxesImage at 0x110917a20>
What are the benefits of filtering and pooling?
These filters clear out a lot of noise.
Protip: median filters are awesome for getting rid of tiny little specks of light you don't want, while maintaining sharp edges between objects (unlike Gaussian blur filters).
Edge-finding is an important aspect of bioimaging and image processing in general.
Edges essentially denote "image derivatives" (and are essentially calculated as such!). Edges are formed by the image pixels changing suddenly.
Edges delineate objects. The question isn't necessarily where the edges of the object are, but which edges belong to the object you're interested in.
One of the most popular algorithms for finding edges in an image is the Canny Edge Detector.
The Canny edge detector works in three distinct phases (don't worry, you won't have to implement any of these):
1: It runs a filter (!) over the image, using the Gaussian formulation, to generate a filtered image. However, this Gaussian filter is a special variant--it actually represents the first derivative of a Gaussian. Running this filter over the image essentially computes the image derivatives at each pixel.
2: The filter generates a lot of "candidate" edges that have to be pruned down; the edges are reduced until they're only 1 pixel thick.
3: Finally, a threshold is applied: a pixel is considered part of an "edge" if its derivative (as computed in step 1) exceeds a certain value. The higher the threshold, the larger the pixel derivative has to be to be considered an edge.
In action, Canny edge detectors looks something like this:
import skimage.feature as feature
img = np.array(img)
canny_small = feature.canny(img, sigma = 1)
canny_large = feature.canny(img, sigma = 3)
f = plt.figure(figsize = (12, 6))
f.add_subplot(1, 3, 1)
plt.imshow(img, cmap = "gray")
f.add_subplot(1, 3, 2)
plt.imshow(canny_small, cmap = "gray")
f.add_subplot(1, 3, 3)
plt.imshow(canny_large, cmap = "gray")
<matplotlib.image.AxesImage at 0x102e63be0>
sigma = 3
than sigma = 1
.sigma
argument essentially equates to the "width" of the filter--smaller width, smaller neighborhood; therefore, more representation by the middle pixel.Another common noise-reducing step is some combination of erosion and dilation.
These don't find edges per se, but rather they modify the edges.
Erosion will move along objects and erode them by 1 pixel (or more).
This has the benefit of utterly wiping out objects that are very small--these are usually noise anyway.
Dilation is the inverse operation: it will move along objects, padding their edges by 1 pixel (or more).
This has the effect of smoothing object edges--jagged edges are filled in.
These two effects are often used in tandem to smooth out potentially rough images.
import skimage.morphology as morphology
de1 = morphology.dilation(morphology.erosion(img))
de2 = morphology.dilation(morphology.erosion(de1))
f = plt.figure(figsize = (18, 9))
f.add_subplot(1, 3, 1)
plt.imshow(img, cmap = "gray")
f.add_subplot(1, 3, 2)
plt.imshow(de1, cmap = "gray")
f.add_subplot(1, 3, 3)
plt.imshow(de2, cmap = "gray")
<matplotlib.image.AxesImage at 0x113c55390>
Hmm, not an obvious effect in the full grayscale.
Let's try again, taking a page from yesterday's lecture--specifically, let's use a single channel and threshold the image.
import scipy.ndimage as ndimg
img = ndimg.imread("Lecture22/image1.png")
hsp = img[:, :, 1] # Green channel - heat shock protein
hsp_bin = hsp > np.mean(hsp) # Median threshold.
# Now we'll do the same code as before.
de1 = morphology.binary_dilation(morphology.binary_erosion(hsp_bin))
de2 = morphology.binary_dilation(morphology.binary_erosion(de1))
f = plt.figure(figsize = (18, 9))
f.add_subplot(1, 3, 1)
plt.imshow(hsp_bin, cmap = "gray")
f.add_subplot(1, 3, 2)
plt.imshow(de1, cmap = "gray")
f.add_subplot(1, 3, 3)
plt.imshow(de2, cmap = "gray")
<matplotlib.image.AxesImage at 0x115612b38>
You can go a step further. If you're not interested in the overall shape of the objects, but just want some measure of how many and where, you can use skeletonization.
This essentially performs the erosion operation over and over until each object is only 1 pixel wide.
skeleton = morphology.skeletonize(hsp_bin)
f = plt.figure(figsize = (12, 6))
f.add_subplot(1, 2, 1)
plt.imshow(hsp_bin, cmap = "gray")
f.add_subplot(1, 2, 2)
plt.imshow(skeleton, cmap = "gray")
<matplotlib.image.AxesImage at 0x115a09668>
What went wrong, do you think? Why didn't we just get 1-2 lines per object?
We should try to eliminate all those specks!
import skimage.filters as filters
bin_median = filters.median(hsp_bin, morphology.square(5))
bin_median[bin_median == 255] = 1
skeleton = morphology.skeletonize(bin_median)
f = plt.figure(figsize = (18, 9))
f.add_subplot(1, 3, 1)
plt.imshow(hsp_bin, cmap = "gray")
f.add_subplot(1, 3, 2)
plt.imshow(bin_median, cmap = "gray")
f.add_subplot(1, 3, 3)
plt.imshow(skeleton, cmap = "gray")
<matplotlib.image.AxesImage at 0x1145ede80>
The inverse operation of skeletonization is the to find the convex hull of objects in the image. This is basically a fancy way of saying:
This is a very good way of deriving masks for complex!
Let's give ourselves a bit of a head start by using that binarized median filtered image:
convex_hulls = morphology.convex_hull_object(bin_median)
f = plt.figure(figsize = (12, 6))
f.add_subplot(1, 2, 1)
plt.imshow(bin_median, cmap = "gray")
f.add_subplot(1, 2, 2)
plt.imshow(convex_hulls, cmap = "gray")
<matplotlib.image.AxesImage at 0x1161060b8>
We can overlay the two plots to show how the convex hulls were drawn:
import skimage
img_float = skimage.img_as_float(convex_hulls.copy())
img_float[bin_median == 1] = 2
plt.imshow(img_float, cmap = "gray")
<matplotlib.image.AxesImage at 0x1173eecc0>
Ok--we've seen so far how to identify objects (kinda) in images. It'd be great to actually get their coordinates, extract them, and learn something about their shape and size and general morphology.
This is the process of segmentation: explicitly developing a mask that pulls out the objects you're interested in.
The convex hulls are a good start, but that one blob near the middle combined a bunch of the cells.
(yes, we're purposely hobbling ourselves by not using the DAPI channel; for illustration purposes)
We can employ a specific segmentation to try and break that blob apart into multiple objects, ideally correpsonding with the blobs we see in the grayscale image.
This starts with the idea of treating your image as a contour map, where pixel intensity indicates elevation.
From the "lowest" parts of your image, you start filling your image with water.
Where water from each "basin" meets, that's a segmentation boundary.
Continue raising the water level of your image until you have fully segmented it.
This will give you an image with boundaries between all the objects.
"Wait a minute," you start, "that looks great and all, but how do you know where to start adding the water?"
Excellent question. Ideally, we'd like to place a "basin" (or "seed") inside each object. But that requires knowing where the objects are ahead of time, which kind of obviates the need for watershed in the first place.
Well, I kinda gave it away earlier--we're looking for pixels that have values at the extremes. These are known as local maxima.
We'd like to think that each object has a "center" that is bright; ideally, it's the brightest part of the object.
This isn't always the case, but we can kind of force it using a distance map. If we compute the distance of any point inside an object to its edge, the center should be the farthest, right?
Let's first see about computing the distance map for our binarized median filtered image.
distmap = ndimg.distance_transform_edt(bin_median)
f = plt.figure(figsize = (12, 6))
f.add_subplot(1, 2, 1)
plt.imshow(bin_median, cmap = "gray")
f.add_subplot(1, 2, 2)
plt.imshow(distmap, cmap = "gray")
<matplotlib.image.AxesImage at 0x1133fd0f0>
(yep, you could probably use skeletonization to achieve the same effect!)
Now, we'll find the coordinates of the local maxima and use those as seeds for watershed.
maxima = feature.peak_local_max(distmap, indices = False, min_distance = 2, labels = bin_median)
plt.imshow(maxima, cmap = "gray")
<matplotlib.image.AxesImage at 0x119ceb198>
Now we'll use these coordinates as our starting points for watershed!
markers = ndimg.label(maxima)[0]
labels = morphology.watershed(-distmap, markers, mask = bin_median)
plt.imshow(labels, cmap = "nipy_spectral_r")
<matplotlib.image.AxesImage at 0x116e918d0>
It's certainly not perfect.
There are a lot of improvements that could be made, but this was just to illustrate how to get segmented regions.
Once you have these regions, you can ask a lot of interesting questions about them:
import skimage.measure as measure
regions = measure.regionprops(labels)
print(len(regions)) # 28 regions, each with some of the following properties:
28
region1 = labels.copy()
region1[region1 != regions[1].label] = 0
plt.imshow(region1, cmap = "gray")
print(regions[1].area)
7601
print(regions[1].centroid)
(92.791080121036714, 79.788975134850688)
print(regions[1].major_axis_length)
print(regions[1].minor_axis_length)
104.25411574839332 93.36852033227652
print(regions[1].solidity)
0.933210558625
matplotlib
(finally), and Open Science


TERMS AND CONDITIONS OF USE OPEN SHELTER DATABASE
Article 1. Open Shelter Database.
International Federation of Red Cross and Red Crescent Societies- Shelter
Research Unit (“IFRC-SRU”) (“we”, “us” and “our”) welcomes
you to the Open Shelter Database (“OSD”) enabling you to upload and update
your own data (text, image and sound material) and through which you may
access without charge a selection of content and data regarding
humanitarian shelters ("OSD Content"). Set out below are
the terms and conditions for use of OSD and the OSD Content. These terms
and conditions apply whenever you access, upload and update the OSD
Content, on whatever device and by any means.
By applying to register for OSD you accept these terms and conditions. You
agree that terms and conditions may be updated by us from time to time at
our discretion without prior notice to you and you have entered into an
agreement with us on such terms.
We reserve the right to change any aspect of OSD and these terms and
conditions at any time without prior notice to you. You are therefore
responsible for regularly reviewing these terms and conditions.
No one may use OSD in a manner that is contrary to the Fundamental
Principles of the Red Cross and Red Crescent Movement (
http://www.ifrc.org/who-we-are/vision-and-mission/the-seven-fundamental-principles/
).
Article 2. Participator.
2.1 OSD is an online open-content database which is developed by a
voluntary association of individuals and groups - having a Humanitarian ID
of OCHA - (“Participator (-s)”) working to realize a common resource of
humanitarian shelters allowing to add the OSD Content, provided that any
update is only possible as to OSD Content as uploaded by yourself. Once you
upload and update the OSD Content, you are considered a Participator. By
uploading/updating the OSD Content you consent to publish the same without
any limitations and charge on OSD. OSD`s role is limited to a role of a
facilitator.
2.2 Participator is responsible and liable for the accuracy and reliability of the OSD Content provided. Participator guarantees that
the OSD Content uploaded and/or updated by you:
2.2.1 is not against the law any existing national or international law or
regulation;
2.2.2 does not create any infringement on intellectual property/copy rights
or privacy of third parties;
2.2.3 will not be otherwise unlawful.
2.3 Participator shall not (allow) use of OSD for improper, illegal and/or
criminal acts, which implies among others :
2.3.1the distribution of material/data which is contrary to any national or
international law, regulation and/or these terms and conditions;
2.3.2 the illegal/criminal distribution of non-public information;
2.3.3 the distribution of material/data which in our view is misleading, is
in violation of any national or international law and/or public
order/morality;
2.3.4 the distribution of material/data which create infringement on
intellectual property rights and/or on (privacy) rights;
2.3.5 the distribution of material/data which create infringement
on someone else's personality rights
;
2.3.6 unauthorized intrusion (‘hacking’);
2.3.7 the distribution of material/data which refer to illegal services
and/or to services which may create damages in whatever way to
persons/property;
2.3.8 damage/destruction of systems/software.
2.4 We reserve the right to change/refuse/delete your material/data which
is covered by article 2.2 respectively by article 2.3 without being liable
to pay any compensation to you and without prejudice to all of our other
rights.
Article 3. Use OSD
3.1 You shall not (allow) use of OSD for any purpose other than for which
it is intended and shall refrain from sending unsolicited/unwanted messages. If you require specific advice (for example,
medical, legal, financial or risk management), the OSD Content - subject to
article 4 - may be of use however it is your responsibility to seek a
professional assistance from a licensed practitioner in that area. Any
reliance upon OSD Content shall be at your own risk.
3.2 We are fully entitled, at our discretion to, to change, adjust, limit, update in full or in part any content in OSD and/or to discontinue OSD at any time. We shall have no obligation whatsoever to give you prior notice of such decision of the OSD.
3.3 We reserve the right to deny in our sole discretion any user access to OSD or any portion thereof without notice.
3.4 You agree to indemnify us from and against any and all actions, claims, losses, damages, liabilities and expenses (including reasonable attorneys’ fees) arising out of your breach of these terms and conditions.
Article 4. Exclusions of our liability.
4.1 Please be advised that the OSD Content has not necessarily been
reviewed by people with the expertise required to provide you with
complete, accurate or reliable information. The OSD Content is provided “as
is” basis, without warranty of any kind, either express or implied,
including, without limitation, warranties of merchantability, fitness for a
particular purpose and non-infringement. Therefore we are not liable in
whatever way for the completeness, fitness, legitimacy and correctness
neither of OSD, the OSD Content nor for the acts and omissions by
Participators. This also applies as to the OSD Content which contains any
information generated by ourselves whether or not based on the OSD Content
provided by Participators.
If you are not satisfied with OSD and/or the OSD Content your sole and
exclusive remedy is to discontinue the use of OSD.
4.2 We do not guarantee:
4.2.1 the completeness, fitness, legitimacy and correctness neither of the
OSD, the OSD Content nor for the acts and omissions by Participators;
4.2.2 that OSD is completely ‘error-free’ provided that we shall use our
reasonable endeavours to achieve that OSD shall perform correctly; and
4.2.3 that OSD is legitimately used.
4.3 TO THE MAXIMUM EXTENT PERMITTED BY THE APPLICABLE LAW, IN NO EVENT SHALL WE BE LIABLE FOR ANY DIRECT, LOST REVENUES OR PROFITS OR ANTICIPATED SAVINGS, GOODWILL, OR OTHER SPECIAL, INDIRECT, INCIDENTAL, EXEMPLARY, CONSEQUENTIAL OR PUNITIVE DAMAGES HOWEVER CAUSED AND REGARDLESS OF THEORY OF LIABILITY WHETHER ANY CLAIM FOR SUCH RECOVERY IS BASED UPON THEORIES OF CONTRACT, NEGLIGENCE OR TORT (INCLUDING STRICT LIABILITY), EVEN IF WE HAVE BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGES. YOU HEREBY WAIVE, FOR YOURSELF AND YOUR SUCCESSORS AND ASSIGNEES, ANY AND ALL CLAIMS FOR DIRECT, LOST REVENUES OR PROFITS OR ANTICIPATED SAVINGS, GOODWILL, OR OTHER SPECIAL, INDIRECT, INCIDENTAL, EXEMPLARY, CONSEQUENTIAL AND/OR PUNITIVE DAMAGES.
4.4 We shall not be held liable for any loss/damage of any data/content.
4.5 In the event that a limitation shall be held to be invalid for any reason and we become liable for loss or damage that would otherwise have been excluded, such total liability shall never exceed the amount as covered and paid out by our relevant insurance company for such specific loss or damage, provided that you agree that indirect/consequential damages/costs of whatever nature shall always be excluded.
Article 5. Privacy & Intellectual property.
5.1 We refer to the privacy notices and regulations of the OCHA which shall
apply. Furthermore:
5.2 We are committed, in the course of our activities and in compliance with the applicable law on data protection, to ensure the confidentiality, integrity and availability of your personal data, as well as to respect your privacy. We, through OSD, collect, use and process your personal data in order to provide you with a service, and for the purpose of publishing and sharing information plus documents related to humanitarian shelters.
5.3 We will not sell, rent or transfer your personal data to third parties, such as but not limited to our IT service supplier, consultants and any other third party processing personal data on our behalf, without your consent and for any other reasons consistent with the purpose for which the data were originally collected or for any other purposes authorized by law.
5.4 However, you expressly give your consent and agree to the disclosure of
your personal data by us to third parties if we:
5.4.1 need to send the information to third parties who work on behalf of
or for us to provide a product or service you have requested - unless
otherwise stated so, these third parties are prohibited to use the personal
data we provided them beyond what is necessary to provide the product or
service to you.
5.4.2 need to comply with any requests from a competent authority or court
orders.
5.4.3 find that it is necessary, as determined in our sole discretion, to
investigate, prevent or take action regarding illegal activities, suspected
fraud, emergency situations involving potential threats to the physical
safety of any persons, violations to these terms and conditions or as
otherwise required by law.
5.5 Before disclosing your personal data to third parties, we shall try to achieve that sufficient security measures and adequate data protection are guaranteed by the third party and take reasonable steps to achieve the compliance with those measures.
5.6 As part of its internal data auditing process, we may conduct regular checks on processing by third parties, especially with respect to security measures.
5.7 We do not retain personal data longer than necessary for the purposes stated and to comply with legal requirements, unless otherwise agreed to by you.
5.8 You have the right to request information on personal data related to you which are processed by us. You may request that inaccurate or incomplete personal data related to you be corrected or deleted where the purpose of data processing has lapsed or ceased to be applicable for other reasons.
5.9 All requests should be sent by regular mail along with a proof of
identity (copy of official identification document with a photo mentioning
your date & place of birth), mentioning OSD, Humanitarian ID of OCHA
and your contact details (phone number and email) to:
Tel: (+352)27558904; Mail<EMAIL_ADDRESS>5.10 We may refuse, restrict or defer your request
respectively the provision of information where we have the right to do so
under current data protection legislation.
5.11 In order to provide you safe access and use of OSD to the best of our ability, we define and implement the reasonable technical and organizational measures necessary to maintain the security of personal data, according to the nature of the personal data processed and the circumstances of the processing, with the objective of avoiding (in the realm of the possible and having regard to the state of the art) non-authorised processing or access, alteration or loss (to ensure to our best endeavours confidentiality, integrity and availability respectively).
5.12 Any of the intellectual property rights that are mentioned, used or cited in the OSD Content are and shall remain the property of their respective owners.
Article 6. Law.
Any dispute, controversy or claim arising out of or relating to the
execution of these terms and conditions shall be settled amicably through
direct negotiation. If these negotiations are unsuccessful, the matter
shall be finally submitted to the jurisdiction of the laws and courts of
Luxembourg. You hereby
consent to the exclusive personal jurisdiction of such courts and
unequivocally waive any and all defenses and/or rights to challenge the
jurisdiction or venue of such courts.


Accursed Summoner Keats:
Ribbit, ribbit, ribbit!
"I heard what happened from Sparky. I'm here to help you."


WAX, OR THE DISCOVERY OF TELEVISION AMONG THE BEES A film review by Caleb J. Howard Copyright 1992 Caleb J. Howard
Hi there. I just thought I'd share a small thing. Last night I saw a movie by the above title at the Bloor Cinema in Toronto as a part of the Festival of Grand Illusions. The film is 85 minutes, produced by Independent Electronic Cinema, and consists of a powerful mixture of computer imagery, military footage of weapons guidance systems, state of the art video post-production techniques, and a *lot* of Babylonian subterranean bees.
I'll use the terms I'm familiar with to describe the movie as it appealed to me, because that's my option.
The plot of the film traced the karmic journey of the protagonist, a designer of targeting sights on flight simulators and beekeeper, through his links backward through time to his grandfather, and forward through time to his destiny. The concept of time was illustrated as fluid and malleable at the spiritual level. The use of video processing and computer techniques to smoothly move the viewer from the spiritual perspective of the protagonist to the perspective of a huge hive of bees to the perspective of the dead souls held on reels of magnetic tape that contain the landscape data of the flight simulators. This is intended to convey the flavour of the film rather than as an in-depth review, so I will leave it as such.
Concepts of interest to me especially well illustrated in the movie include the state of consciousness of the hive of bees taken as a whole and how this state of mind reflects the global state of mind of humanity and all life taken as a whole. Also well illustrated are the interactions of mind, television, reality and virtual reality. Especially neat was the exploration of the state of mind of a semi-intelligent guidance computer in the nose of a tank-killer missile in the Gulf War.
As far as I'm concerned, the movie is absolutely inspired, and is a precursor to the future of media events. I was consumed by the imagery and enlightened by the vision of the film and the mastery with which it was executed. I also felt at the end of the film that perhaps a half a dozen others in the theatre were similarly enthusiastic, and that the majority were left completely cold. Several people left the theatre during the film (one man was really angry at having wasted his money on the film). At the end, there were only about six or seven of us applauding, and the rest were looking at us, puzzled.
I thought it was great, but would only recommend it to people with an appetite for the surreal and the symbolic.
Has anyone else seen it?
-caleb
.
The review above was posted to the
rec.arts.movies.reviews newsgroup (de.rec.film.kritiken for German reviews).
The Internet Movie Database accepts no responsibility for the contents of the
review and has no editorial control. Unless stated otherwise, the copyright
belongs to the author.
Please direct comments/criticisms of the review to relevant newsgroups.
Broken URLs inthe reviews are the responsibility of the author.
The formatting of the review is likely to differ from the original due
to ASCII to HTML conversion.
Related links: index of all rec.arts.movies.reviews reviews


| 105 | (A) | a kynge and a knawe childe consayuede at the Ere |
| 100 | (C) | A kynge & a man childe conceyued at her ere |
| 101 | (D) | A kyng and a knave child conceyvid at here ere |
| 99 | (E) | a kyng & a knaue child . conceyued be her here |
| (Ex) | ---- no text ---- | |
| 104 | (L) | a kyng & a knaue child ȝo conceyued at ere |
| 99 | (P) | A kynge and a knafe childe consayued at hir ere |
| 100 | (U) | A kyng and a knaue chiled conceiued at hir ere |
| (V) | ---- no text ---- |


Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
|
urlparse, add_params_to_uri, add_params_to_qs, ) MissingCodeException, MissingTokenException, MissingTokenTypeException, MismatchingStateException, )
scope=None, state=None, **kwargs): """Prepare the authorization grant request URI.
The client constructs the request URI by adding the following parameters to the query component of the authorization endpoint URI using the ``application/x-www-form-urlencoded`` format:
:param uri: The authorize endpoint to fetch "code" or "token". :param client_id: The client identifier as described in `Section 2.2`_. :param response_type: To indicate which OAuth 2 grant/flow is required, "code" and "token". :param redirect_uri: The client provided URI to redirect back to after authorization as described in `Section 3.1.2`_. :param scope: The scope of the access request as described by `Section 3.3`_. :param state: An opaque value used by the client to maintain state between the request and callback. The authorization server includes this value when redirecting the user-agent back to the client. The parameter SHOULD be used for preventing cross-site request forgery as described in `Section 10.12`_. :param kwargs: Extra arguments to embed in the grant/authorization URL.
An example of an authorization code grant authorization URL::
/authorize?response_type=code&client_id=s6BhdRkqt3&state=xyz &redirect_uri=https%3A%2F%2Fclient%2Eexample%2Ecom%2Fcb
.. _`Section 2.2`: https://tools.ietf.org/html/rfc6749#section-2.2 .. _`Section 3.1.2`: https://tools.ietf.org/html/rfc6749#section-3.1.2 .. _`Section 3.3`: https://tools.ietf.org/html/rfc6749#section-3.3 .. _`section 10.12`: https://tools.ietf.org/html/rfc6749#section-10.12 """ params = [ ('response_type', response_type), ('client_id', client_id) ]
if redirect_uri: params.append(('redirect_uri', redirect_uri)) if scope: params.append(('scope', list_to_scope(scope))) if state: params.append(('state', state))
for k in kwargs: if kwargs[k]: params.append((to_unicode(k), kwargs[k]))
return add_params_to_uri(uri, params)
"""Prepare the access token request. Per `Section 4.1.3`_.
The client makes a request to the token endpoint by adding the following parameters using the ``application/x-www-form-urlencoded`` format in the HTTP request entity-body:
:param grant_type: To indicate grant type being used, i.e. "password", "authorization_code" or "client_credentials". :param body: Existing request body to embed parameters in. :param redirect_uri: If the "redirect_uri" parameter was included in the authorization request as described in `Section 4.1.1`_, and their values MUST be identical. :param kwargs: Extra arguments to embed in the request body.
An example of an authorization code token request body::
grant_type=authorization_code&code=SplxlOBeZQQYbYS6WxSbIA &redirect_uri=https%3A%2F%2Fclient%2Eexample%2Ecom%2Fcb
.. _`Section 4.1.1`: https://tools.ietf.org/html/rfc6749#section-4.1.1 .. _`Section 4.1.3`: https://tools.ietf.org/html/rfc6749#section-4.1.3 """ params = [('grant_type', grant_type)]
if redirect_uri: params.append(('redirect_uri', redirect_uri))
if 'scope' in kwargs: kwargs['scope'] = list_to_scope(kwargs['scope'])
if grant_type == 'authorization_code' and 'code' not in kwargs: raise MissingCodeException()
for k in kwargs: if kwargs[k]: params.append((to_unicode(k), kwargs[k]))
return add_params_to_qs(body, params)
"""Parse authorization grant response URI into a dict.
If the resource owner grants the access request, the authorization server issues an authorization code and delivers it to the client by adding the following parameters to the query component of the redirection URI using the ``application/x-www-form-urlencoded`` format:
**code** REQUIRED. The authorization code generated by the authorization server. The authorization code MUST expire shortly after it is issued to mitigate the risk of leaks. A maximum authorization code lifetime of 10 minutes is RECOMMENDED. The client MUST NOT use the authorization code more than once. If an authorization code is used more than once, the authorization server MUST deny the request and SHOULD revoke (when possible) all tokens previously issued based on that authorization code. The authorization code is bound to the client identifier and redirection URI.
**state** REQUIRED if the "state" parameter was present in the client authorization request. The exact value received from the client.
:param uri: The full redirect URL back to the client. :param state: The state parameter from the authorization request.
For example, the authorization server redirects the user-agent by sending the following HTTP response:
.. code-block:: http
HTTP/1.1 302 Found Location: https://client.example.com/cb?code=SplxlOBeZQQYbYS6WxSbIA &state=xyz
""" query = urlparse.urlparse(uri).query params = dict(urlparse.parse_qsl(query))
if 'code' not in params: raise MissingCodeException()
if state and params.get('state', None) != state: raise MismatchingStateException()
return params
"""Parse the implicit token response URI into a dict.
If the resource owner grants the access request, the authorization server issues an access token and delivers it to the client by adding the following parameters to the fragment component of the redirection URI using the ``application/x-www-form-urlencoded`` format:
**access_token** REQUIRED. The access token issued by the authorization server.
**token_type** REQUIRED. The type of the token issued as described in Section 7.1. Value is case insensitive.
**expires_in** RECOMMENDED. The lifetime in seconds of the access token. For example, the value "3600" denotes that the access token will expire in one hour from the time the response was generated. If omitted, the authorization server SHOULD provide the expiration time via other means or document the default value.
**scope** OPTIONAL, if identical to the scope requested by the client, otherwise REQUIRED. The scope of the access token as described by Section 3.3.
**state** REQUIRED if the "state" parameter was present in the client authorization request. The exact value received from the client.
Similar to the authorization code response, but with a full token provided in the URL fragment:
.. code-block:: http
HTTP/1.1 302 Found Location: http://example.com/cb#access_token=2YotnFZFEjr1zCsicMWpAA &state=xyz&token_type=example&expires_in=3600 """ fragment = urlparse.urlparse(uri).fragment params = dict(urlparse.parse_qsl(fragment, keep_blank_values=True))
if 'access_token' not in params: raise MissingTokenException()
if 'token_type' not in params: raise MissingTokenTypeException()
if state and params.get('state', None) != state: raise MismatchingStateException()
return params |


Simulation
The system
object can generate one or multiple simulation
objects, which determine how to simulate the system
.
This includes among other things the diffusion integrator, the reaction handler, observables.
The initial positions of particles are also set on the simulation
object.
Given a system
one can generate a simulation
by invoking
simulation = system.simulation()
The function takes a number of arguments that influence the way the simulation is executed:
simulation = system.simulation(
kernel="SingleCPU",
output_file="",
integrator="EulerBDIntegrator",
reaction_handler="Gillespie",
evaluate_topology_reactions=True,
evaluate_forces=True,
evaluate_observables=True,
skin=0
)
Except for the kernel
argument, all of these arguments can also be modified by setting properties on the
simulation object. The configuration of the reaction diffusion system
is copied into the simulation object,
so subsequent changes to the reaction diffusion system will not propagate into the simulation.
Selecting a kernel
Currently there are two different kernels that are supported: the SingleCPU
and the CPU
kernel. As the name
suggests, the SingleCPU
kernel is single-threaded whereas the CPU
kernel attempts to parallelize as much as
possible, thus making use of more cores.
It can be expected that the SingleCPU
implementation is roughly as fast as the CPU
implementation with two threads,
as it applies Newton’s third law for calculating pairwise forces and evaluates reactions per particle pair, which are
the two major performance bottlenecks.
Configuration
In the following it will be explained how to add particles, add topologies, configure specifics of the selected kernel, how to record a trajectory, and how to perform checkpointing.
Adding particles
For adding particles to a system there are two separate methods. One is can be to place a single particle, one
is for bulk insertion.
Adding a single particle of type A
to a simulation box amounts to
simulation.add_particle(type="A", position=pos)
where pos
can be a list [x, y, z]
, tuple (x, y, z)
, or a numpy.ndarray: np.array([x, y, z])
with three entries
representing the x,y,z components.
When one wants several particles of a certain type to the simulation, one can can exchange multiple calls to
simulation.add_particle
by the better performing variant
X = np.random.random((100, 3))
simulation.add_particles(type="A", positions=X)
taking a (N, 3)
-shaped numpy array as position argument, resulting in N
particles with their respective positions
being added to the simulation. In this example, 100 particles of type A
would be placed uniformly at random
in $[0,1)^3$.
Adding topologies
A topology can be added by invoking
my_topology = simulation.add_topology(topology_type="My topology type", particle_types="T",
positions=np.random.random((100, 3)))
which requires a “My topology type” topology type and a topology species “T”
to be registered in the ReactionDiffusionSystem
. In this example the
topology will contain 100 randomly placed topology particles of type “T” that are for now disconnected.
If the topology should contain several different particle types one can pass a list of particle types to the particle_types
argument
that contains types for all the positions:
my_topology = simulation.add_topology(
topology_type="My topology type",
particle_types=["T1", "T2", "T3", "T2", "T1"],
positions=np.random.random((5, 3))
)
Unless the topology consists out of only one particle, one still needs to set up the connectivity graph before running
the simulation. The returned object my_topology
is a topology object as the ones described in
topology reactions. Edges in the graph can be introduced like
my_graph = my_topology.get_graph()
for i in range(len(graph.get_vertices())-1):
my_graph.add_edge(i, i+1)
where the indices that go into add_edge
correspond to the particle positions that were entered in add_topology
.
The simulation can only be executed if the graph of each topology is connected,
- i.e., there are no independent components (between each pair of vertices there is at least one path in the graph that connects them),
and for each edge there is a bond,
- i.e., all particle type pairs that are contained in the graph have at least one entry in the topology potential configuration.
Should one of these two conditions be not fulfilled, starting the simulation will raise an error.
Kernel configuration
In case of selecting the CPU kernel with a parallelized implementation of the ReaDDy model, one can change certain aspects of its behavior:
- The number of threads to be used can be selected by
simulation.kernel_configuration.n_threads = 4
-
Mainly due to pairwise interactions and bimolecular reactions there is a neighbor list to reduce the time needed for evaluating these. The neighbor list imposes a spatial discretization on the simulation box into cuboids. In the default case, each of these cuboids has an edge length of at least the maximal cutoff radius / reaction radius. This means that instead of naively looping over all particle pairs ($\mathcal{O}(N^2)$), one can assign each particle to its cuboid and then loop over all particles in a cuboid and its 26 neighboring cuboids to find particle pairs.
When collecting particle pairs in this fashion one effectively approximates a sphere with cuboids. The number of potential interaction or reaction partners can be further reduced by using only a fraction of the edge length but increasing the search radius of the neighboring boxes so that one still covers at least the cutoff radius in each spatial dimension.
Reducing the edge length usually comes with a price, at some point the bookkeeping of neighboring boxes dominates the runtime.
The edge length and therefore search radius can be controlled by
simulation.kernel_configuration.cell_linked_list_radius = 4
which would yield cuboids with $\frac{1}{4}$ the edge lengths of the default case.
Recording a trajectory
ReaDDy records trajectories and observable data in HDF5 files. For doing so one needs to set an output file
simulation.output_file = "my_trajectory.h5"
and instruct the simulation to record a trajectory:
simulation.record_trajectory(stride=1, name="", chunk_size=1000)
The stride
arguments causes the trajectory to be recorded every stride
time steps. If a name
(other than
the default one) is given, the trajectory data will be stored in a different data set. The chunk_size
is mainly
a performance parameter that has an effect on how large every chunk of data in the binary file is going to be,
influencing the time needed for IO during the simulation and the resulting file size.
For reading back the trajectory data, please refer to post-processing.
NOTE: When running long simulations on a cluster it can happen that the process runs into a timeout causing the already recorded data to be corrupted. This can possibly be mitigated by configuring the job manager to send SIGINT before KILLing the process. That way the file can still be properly closed (see issue #220, thanks @jansteinkuehler). Please make sure this works in your environment before running long simulations.
Checkpointing
Checkpoints in ReaDDy consist out of the particles’ and topologies’ configurations at specific points in simulation time. They can be enabled by calling
simulation.make_checkpoints(stride=1000, output_directory="checkpoints/", max_n_saves=10)
which causes checkpoints to be made every 1000 steps. Each checkpoint is a separate file and all checkpoint files will be
saved to the directory specified by output_directory
. The option max_n_saves
decides how many checkpoint files
are allowed to be saved to the directory, e.g. if max_n_saves=10
then only the last 10 most recent checkpoints
are kept.
Once the simulation has run its course and checkpoint files have been written, they can be listed by
simulation.list_checkpoint_files('checkpoints/')
A particular checkpoint file can in principle also contain multiple checkpoints. They can be inspected by
simulation.list_checkpoints('checkpoints/checkpoint_10000.h5')
and a system’s state can be restored by a call to
simulation.load_particles_from_checkpoint('checkpoints/checkpoint_10000.h5')
amounting to restoring the latest checkpoint of that particular file. If the file contains multiple checkpoints, let’s say 5, you can
select the 5th checkpoint by supplying the optional argument n=4
(enumeration starts at n=0
per file).
Oftentimes you just need the last checkpoint of all checkpoint files in a certain directory. This can be achieved by
simulation.load_particles_from_latest_checkpoint('checkpoints/')
It should be noted that if the simulation should be continued and the output_directory
for the new checkpoints is the
same as of the original simulation, the old checkpoint files will be overwritten. If you want to keep the checkpoints
of the original simulation, specify another output_directory
.
Observables
The currently available observables are:
- Radial distribution function
- Particles
- Particle positions
- Number of particles
- Energy
- Forces
- Reactions
- Reaction counts
- Virial
- Pressure
There are three things that all observables have in common: The evaluation can be strided, they can have a callback
function and they can be saved to the simulation.output_file
.
The callback
is usually a function that takes one argument being the current value of the observable. During the
course of the simulation this callback function will be evaluated whenever the particular observable is evaluated.
Per default, whenever an output_file
is given, the registered observables’ outputs are saved to that file. Each
observable has a certain place in the group hierarchy of the HDF5 file, however this place can be modified so that,
e.g., multiple observables of the same type can be recorded into the same file.
To this end, the save
argument of the respective observable can be modified. By providing
None
orFalse
writing the results to file can be switched off,- providing a dictionary with keys
'name'
and'chunk_size'
can modify the name under which the observable data is stored in the group hierarchy and the hdf5 chunk size. Thechunk_size
is always to be considered into the “temporal direction”, i.e., if an observable yields data in the form of3x3
matrices each time it is evaluated, a chunk would be of shape(3, 3, chunk_size)
.
Radial distribution function
The radial distribution function for certain particle types can be observed by
def rdf_callback(current_value):
print(current_value)
simulation.observe.rdf(
stride=1,
bin_borders=np.linspace(0., 2., 10),
types_count_from=["A"],
types_count_to=["B"],
particle_to_density=1./system.box_volume,
callback=rdf_callback)
which causes the observable to be evaluated in each time step (stride=1
) and print the value (callback=rdf_callback
).
The RDF is determined by calculating the distance from all particles of a type contained in types_count_from
to
all particles of a type contained in types_count_to
and then binning the distance into a histogram as given by bin_borders
.
The histogram is normalized with respect to $g(r) = 4\pi r^2\rho dr$, where $\rho$ is the number density of particles
with types contained in types_count_to
, reflected by particle_to_density
.
Particles
This observable records all particles in the system, as in: Each particle’s type, (unique) id, and position. It can be registered by
def particles_callback(particles):
types, ids, positions = particles
print("Particle 5 has type {}, id {}, and position {}."
.format(types[5], ids[5], positions[5])
simulation.observe.particles(
stride=5,
callback=particles_callback,
save=False
)
where the argument of the callback function is a 3-tuple containing a list of types, unique ids, and positions
corresponding to each particle in the system. In this example the callback function prints these properties of the
fifth particle every fifth time step, the output of the observable is not saved into the trajectory file (save=False
).
Particle positions
The particles’ positions can be recorded by
simulation.observe.particle_positions(
stride=200,
types=None,
callback=lambda x: print(x)
)
which makes this observable very similar to the particles
one, however one can select specific types of particles
that are recorded. In case of types=None
, all particle positions will be recorded, in case of types=["A", "B"]
only positions of particles with type A
or B
are returned.
In this case the callback will simply print x
every 200 steps, where x
is a list of three-dimensional vectors.
Since save
is not explicitly set to False
or None
the observed data will be recorded into the trajectory file
if n simulation.output_file
was configured.
Number of particles
When one is only interested in the sheer number of particles then one can use this observable. Depending on the input, it will either observe the total number of particles or the number of particles per selected type:
simulation.observe.number_of_particles(
stride=1,
types=["A", "B", "C"],
callback=lambda x: print(x),
save=False
)
This example records the numbers of particles with types A
, B
, C
in each time step. The callback takes
a list with three elements as argument where each element corresponds to a particle type as given in types
and
contains the respective counts. If types=None
was given, the observable would record the total number of particles,
regardless of their types.
Energy
The system’s current potential energy can be observed by
simulation.observe.energy(
stride=123,
callback=lambda x: print("Potential energy is {}".format(x)),
save=False
)
where stride=123
indicates that the observable will be evaluated every 123rd time step. The argument of the callback
function is a scalar value and the observable’s output is not saved to a potentially configured trajectory file.
Forces
The forces acting on particles can be observed by
simulation.observe.forces(
stride=1,
types=["A"],
callback=lambda x: print(sum(x))
)
yielding an observable that is evaluated every time step (stride=1
) and collects the forces acting on all particles
of type A
. If types=None
, all types are considered. The callback function takes a list of vectors as argument.
Since save
is not further specified, this observable would be recorded into the trajectory file.
Reactions
This observable records all occurred reactions in the system for a particular time step. It can be registered by invoking
def reactions_callback(reactions):
for r in reactions:
print("{} reaction {} occurred: {} -> {}, position {}"
.format(r.type, r.reaction_label, r.educts, r.products, r.position))
print("----")
simulation.observe.reactions(
stride=5,
callback=reactions_callback
)
where stride=5
indicates that the observable is evaluated every fifth time step. The callback takes a list of
reaction records as argument, where each reaction record stores information about the
- type of reaction (
r.type
), i.e., one of conversion, fission, fusion, enzymatic, decay, - reaction name (
r.reaction_label
), i.e., the name under which the reaction was registered in thesystem
, - educt unique particle ids (
r.educts
) as in the particles observable, - product unique particle ids (
r.products
), - and position (
r.position
) of the reaction event which is evaluated to the midpoint between educts in case of a bimolecular reaction.
Since the save
argument was left out, it is defaulted and given a simulation.output_file
, the observed reactions are
recorded.
Reaction counts
Instead of recording all reaction events one can also record the number of occurred reactions per registered reaction per time step. This can be achieved by
simulation.observe.reaction_counts(
stride=2,
callback=lambda x: print(x),
save=False
)
where stride=2
causes the observable to be evaluated every second time step. The callback function takes
a dictionary as argument where the keys are the reaction names as given in the
system configuration and the values
are the occurrences of the corresponding reaction in that time step.
Virial
This observable evaluates the pressure virial according to pair-wise forces by evaluating
where $\mathbf{r}_{ij}$ is the vector difference between the positions of particle i and j, $\mathbf{f}$ the pair-wise force.
It can be registered by
simulation.observe.virial(
stride=5,
callback=lambda x: print(x)
)
which causes it to be evaluated every fifth time step (stride=5
). The callback function takes a numpy array of shape (3,3)
as argument.
Pressure
The pressure of a system can be understood in terms of the trace of a pressure tensor
which is defined via the virial tensor
where $V$ is the system’s volume, $N_t$ the number of (physical) particles that somehow partake in the dynamics of the system, and $k_BT$ the thermal energy.
Observing the pressure in a simulation amounts to
simulation.observe.pressure(
stride=1,
physical_particles=["A", "B", "C"],
callback=lambda p: print("current pressure: {}".format(p))
)
where stride=1
causes the pressure to be evaluated in every time step, physical_particles
are set to be particles of
type A
, B
, or C
(physical_particles=None
causes all particles to be considered physical), and the callback
function takes a scalar value as argument.
Internally, the pressure observable builds up on the virial observable
and the number_of_particles observable and when writing it to file
not the actual pressure is recorded but the output of these other two observables.
In the HDF5 group hierarchy the observable’s group is per default postfixed by _pressure
, causing the virial
to be stored under virial_pressure
and the number of particles under n_particles_pressure
.
Changing the postfix amounts to providing a dictionary to the save
argument with
pressure_save_options = {
'name': "empty_or_other_postfix",
'chunk_size': 500
}
simulation.observe.pressure(1, ["A", "B", "C"], save=pressure_save_options)
where the chunk_size
refers to the HDF5 chunk size of data sets.
Running the simulation
Per default, the simulation loop looks like below
This means that all observables are evaluated at the initial state regardless of their stride
, after which
the actual loop is performed.
A simulation is started by invoking
simulation.run(n_steps=1000, timestep=1e-5)
where n_steps
is the number of time steps to perform and the timestep
is the time step. Per default an overview
of the system configuration is printed upon simulation start, this can be disabled by providing the
argument show_system=False
.
One can influence portions of the loop through the simulation
object:
- Per default a progress bar is shown when the simulation is executed, however it can be hidden
print(simulation.show_progress) simulation.show_progress = False
- If one does not want to evaluate topology reactions at all, one can invoke
simulation.evaluate_topology_reactions = False
- Evaluation of forces can be deactivated by invoking
simulation.evaluate_forces=False
- Evaluation of observables can be switched off altogether by
simulation.evaluate_observables=False
Note that setting this to
False
also causes the trajectory to not be recorded. - In case of a large simulation box but small interaction radii one can sometimes boost performance by artifically
increasing the cuboid size of the neighbor list’s spatial discretization by setting
simulation.skin = s
where $s\geq 0$ is a scalar that will be added to the maximal cutoff.
Furthermore, one can select the
Integrator
Currently the only available integrator is the EulerBDIntegrator
which is selected by default and can be selected by
a call to
simulation.integrator = "EulerBDIntegrator"
It integrates the isotropic brownian dynamics
with an Euler-Maruyama discretization
where
Reaction handler
Reactions in ReaDDy are evaluated per time step, meaning that each particle might have multiple possible reaction partners. In order to solve this, one can chose between two different reaction handlers:
- The
UncontrolledApproximation
reaction handler is the less rigorous of the two. It performs as follows:- A list of all possible reaction events is gathered.
- For each element of this list a decision is made whether it is to be evaluated based on the time step and its rate as described in the section about reactions.
- The filtered list is shuffled.
- For each event in the list evaluate it unless its educts have already been consumed by another reaction event.
A problem of this reaction handler is that it does not give preference to reaction events based on their rate in case of a conflict, i.e., two events that share educts. However in the limit of small time steps this problem disappears.
This reaction handler can be selected by invoking
simulation.reaction_handler = "UncontrolledApproximation"
-
The
Gillespie
reaction handler is the default choice and is statistically exact in the sense that it imposes a Gillespie reaction order on the events of a particular time step:A list of all possible reaction events is gathered. Then
- Each reaction event is assigned its normalized reaction rate $\lambda_i/\sum_i\lambda_i$
- A random event is picked from the list based on its normalized reaction rate, i.e., a uniform random number between 0 and 1 is drawn and then used together with the cumulative normalized reaction rates to determine the event
- Depending on its rate the reaction described by the event might be performed:
- if not the event is simply removed from the list
- if it was performed it is also removed and any other event that shared educts with this particular event
- if there are still events in the list go back to 1., otherwise stop
An example of conflicting reaction events with expected outcome might be
$$ \left\{ \begin{array}{rcc} A + B &\xrightarrow{\lambda_1}& C\\ A &\xrightarrow{\lambda_2}& D \end{array} \right. \xrightarrow{\lambda_1 \gg \lambda_2} \left\{ \begin{array}{rcc} A + B &\xrightarrow{\lambda_1}& C\\ &\mathrm{ignored} \end{array} \right. $$This reaction handler can be selected by invoking
simulation.reaction_handler = "Gillespie"