Post a video for visualizing the neighbors of English word (http://en.wikipedia.org/wiki/Word_ladder). It shows what can be done now with minimum change to IPython 0.13-dev source + some simple monkey patches.
I will write down what I think where we can go from here later. The IPython notebook can be download from here. The monkey patches that make this working can be downloaded from my fork of the IPython source code from the GitHub site too.
I gave a lightening talk in a recent Bay Area Python meet-up. I went over some of my recent hacks on combining ipython notebook and d3.js. What I wanted to show was how to mix python code and javascript code to create a dynamic programming/data analysis notebook. I created yet another example to demonstrate the great potential on combining the powerful tools.
If you are interested, you can try the this ipython notebook. You will need to download the development branch of the ipython v 0.13 to see the notebook. The notebook itself includes some of the explanation on how to run it and how it is done. I did not spend too much polishing the code and the motion chart, but it got the basic ingredients. If you want to peek it, here is a short screen recoding to show it looks like.
Thanks for Brian Granger pointing out how to make bi-directional communication from javascript in a ipython-notebook front-end to back-end ipython kernel using the existing websocket/zmq channel architecture in ipython (see the thread ). I have been hacking around to see how to do it. I need to modified a few lines of the ipython-notebook javascript to make it work ( see my github commit ). I wrote some example to show how it works ( the ipython notebook and a screen shot ). Pretty cool that it works. It seems that one can develop a widget library to avoid hand-crafting both the javascript and python code for such communication. All right, one more small step toward to building some cool interactive visualization / analysis tools with ipython.
While I have been using ipython for a long time, I never really it more than just checking whether some code snippets working as expected. (Well, I tried to play with the parallel computing framework with ipython, but I never put it into production.) Just recently, I start to look into the ipython web-based notebook feature more carefully. It is great and make me think the ipython will make a python programmer or someone uses python for data analysis much more productive. (I used to envy the “RStudio” in the R-lang land, now, we python programmer finally have something more competitive.)
The cool thing using a web page as front-end is there are a lot potential using web interface for some cool visualization. I played with protovis.js a while ago. Recently, I went to a visualization meets-up, d3.js was mentioned a numbers of time. Then the idea comes to my mind “is it possible to combine the best of two world, python and d3.js?” After consulting some more experience users in the ipython-dev mailing list to see what is possible, I decided to spend some of my weekend time to hack it around. In the meantime, I get the chance to play with tornado, zero-mq and websocket, all the fun stuff these days. At the end, I am able to pass some javascript code written within the ipython notebook to get the browser to execute it and show some animation with d3.js. This will enable to create more fancier visualization in an interactive way all in a browser.
My weekend hacking results are hosted at github . I think there is a great potential to make thing like this working better. (For example, can we have a pythonic backend of d3.js? 
) It definitely worth to mess it around to see more use like this.
Jun 26
I read Richard Stallman’s story about how he was able to fix a printing problem in the MIT AI lab long time ago. The main point is that if you know how to modify a device or an instrument then you have the ability to extend the functionality well and you can make things more useful. In order to do so, you will have to understand the underlying mechanism about how a thing works. In the old days before any computer code is popular, one might learn how a machine works by just looking how it works if one has good amount domain knowledge. However, the universal applications using computer program to control machines makes modification harder. Computer is a wonderful tool to process complicated information and procedures but it also allows concealing about how a device or an instrument works. When there is no source code or some simple API document allowing you to get some understand how the computer controlling a device or an instrument, then some reverse engineering is un-avoidable. This is my story about how I learn to “hack” a machine in order to do some interesting science.
I got a Ph.D in theoretical physics. My thesis was about stochastic fluctuations in non-equilibrium system. My major interests were to use some statistical mechanics framework to explain interesting natural phenomena. What kind of non-equilibrium system can be more interesting than life itself? I started working as a post-doc in UCSF after graduate school. Initially, I mostly focused on some biological network analysis of protein-protein interactions and gene regulations in the beginning. Trained as a theorist with some computational skills, my initial research plan followed my advisor’s paradigm, which was to gather data in public domain and did analysis on the data to publish new analysis results. However, while there are a few successful cases, such approach is in general not favored in the biological research field in my opinion. First of all, unless you know all the details of the initial experiment that generates the data, it is very likely one might not have enough information to interpret the new analysis results or the interpretation might be simply wrong. Secondly, the researchers who generate the data should have done some “first-order” analysis that reveals most significant information with straightforward analysis with the data set. Most low hanging fruits should be presented in the initial study. Secondly analysis on such data set is just harder and may not have the proportional impact as the initial analysis. Moreover, the biological research field generally values original experiments more than fancy theoretical analysis, although there is a new trend after more and more quantitative emphasis on biological study recently. These lines of thought made me to decide I would need to do my own experiments to gather my own data. But, how?
My advisor in UCSF was also a theoretical physicist who became a bioinformatist and built a career in UCSF providing his quantitative skill to the research community in UCSF. However, as far as I could tell, he had very limited experience in experimental sciences. That has been said, the setup in the new UCSF mission campus was indeed great for me who only had mostly theoretical background to learn new tricks in experimental biological sciences. Our lab was surrounded by some of the best molecule biology labs in the country. I started to learn some basic molecule techniques from the post-docs and graduate students in those labs. I still remember some generous graduate students and post-docs spent their time to show me how to do the basic like how to do precision pipetting. My first few attempts to do simple PCR failed miserably because I was not really able to pipet the 1 ul enzyme precisely for the reaction. Sometimes, I felt so embarrassing since I was not even able to calculate the amount regents for PCR and carried the reaction out correctly even though I though I had no problem to attack some difficult and highly mathematical problems in physics if I wanted to. The learning curve was indeed quite steep, however, I had some other advantage: trained as a physicist, I knew some of the fundamental principles on how certain instruments worked. Moreover, it was more or less easy for me to understand how software could be written and how the software controlling the instruments worked. These skills allowed me to do something interesting and went beyond to do some experiments that the vendor could not support.
When I was analyzing some of the time course microarray data, one question we asked was about how to explain certain pattern in the time course. Were there some fundamental and new principles in governing the dynamics of gene expression? Although microarray allowed monitoring a large number of gene expression at the mRNA level at once, the resolution in time and expression level was quite rough. Could I get a better measuring system for understanding gene regulation dynamics with much better time resolution and expression level quantitation? At that time, GFP-fused proteins were used to measure gene expression at the protein level. Two labs in UCSF worked hard to tag several thousands of proteins in S. cerevisiae with GFP and study their static expression under several conditions and the localization of the proteins in a cell. The high throughput measurements of protein expression were done by using a flow cytometer. A flow cytometer uses a flow stream to separate individual cell such that high throughput optical measurement for each individual cells can be done. Although you can do tons of single cell measurement at once, the typically experiment setup is that you still need to collect cells, put them into a tube, put the tube into the flow cytometer and click a number of buttons on the control software to get the measurement done. I did try to take time course this way manually. For a simple low time resolution experiment, like measurement every 20 minute for several hours, I will have run between the incubator of my cell culture and the instrument room of the flow cytometer all the time. I will have to check the cell density and wash the cell and the flow cytometer intake every 20 minutes. Well, to do large-scale experiments with this way was just probably not my style. I needed to find a better way to do it.
The previous study using the same flow cytometer that I was using was indeed done with a robotic setup. The interesting part was that the robotic system was synchronized to the flow cytometer software by simulating UI events to trigger the measurement software to start and stop recording. The flow cytometer vendor did not support such feature yet. I think it was Joe deRisi figured out how to do that and wrote some program that controlled a sample-delivering robot and trigger the measurement software automatically. Instead of the robotic sample delivering system, I was also setting up a different sampling delivering system that used a bioreactor to deliver live culture continuously to the flow cytometer. Even though I could using same UI simulation hack that Joe did, I was also trying to find some other to collect the data more efficiently since pulling the data out of the instrument software was also quite painful to do. It would be nice I could have controlling system and data collection system in one single program that I had more or less full control on the control and data flow. Joe did notice that the vendor’s software running on a PC talked to the embedding system on the instrument through a standard network interface. Basically, the instrument was controlling by an embedding system and data was collected by the same system. The software running on a PC was merely a frontend and a client. If we could figure out the data stream between the PC and the embedding system, we could do much more than the vendor’s software could do. By the way, we did try to talk with the vendor to see whether we could get some documentation about it. Well, as usual, the vendor was quite reluctant to reveal more. (Now, I am on the vendor side. I probably understand their concern better. But, we did not ask too much. What we needed was some simple description about the communication protocol. It had to be already somewhere in their software organization. Otherwise, software engineer could not write the code.)
I did not remember why I knew how to sniff the network communication. I probably learnt that because I installed some on-the-edge linux software that did not work after installation and I needed to do some tcpdump or so to see what was going on to fix the problem as a hobby. So, I installed “ethereal” (it is now called Wireshark) and start to dump the TCP/IP data stream between the vendor’s software on the PC and the embedding system. Of course, I was hopping they simply used some clear text protocol. Unfortunately, they were using some sort of binary protocol. However, when I was staring at those hex numbers of the output from an “ethereal” session. I did notice some regularity. This means, at least, they did not encrypt the data. Great, at least, I could start to guess what those hex numbers represented. Some portion of the data stream had some regularity. The numbers seems repeating similar pattern every 4 bytes. What could they be? Well, there were a lot of things in computer represented with 32 bits. On the other hand, the most common one for such instrument was probably some floating numbers. Yes, I took those 4 byte chucks and figured out the right byte orders, treated them as IEEE float number format, then I wrote some short code to translate them. The number I got totally make sense!! I could read those numbers and understand what they were. Moreover, I could start to tell which parts are control commands, and which parts are data. The rest of that week, I was running between the lab doing some experiments, collected more TCP/IP data streams under different conditions and went to a hospital to stay with my wife who was expecting a new baby in the antepartum department of hospital. I took my laptop with me and worked on deciphering the data and code while my wife was resting.
Eventually, I got to a point that I could decipher and read the data stream to get everything I needed without using the vendor’s controlling software. The next thing is to combine such knowledge to write my own software (GUI application for control and data processing) to bypass the control application from the vendor. I used the same old version Borland C/C++ to build a GUI interface under Windows such that other people can use it too. (If I were the only user, some scripting through a command line interface was probably good enough.) The results were much useful than the original software that the vendor provides. We could combine the data acquisition and an external robot seamlessly to make the measurement more reliable. Actually, there were a few glitches here and there initially since I was almost the sole developer for the project and there were no comprehensive testing and user case studies. Most of them got fixed and could be fixed easily. However, I was always afraid something might be broken after a vendor software/firmware upgrade. It is possible that the vendor can totally do a firmware upgrade that changing the protocol that I have reversed engineered. If so, I will need to restart the whole process.
There were also some side-effect after people can really see the “raw data”. The flowcytometer integrated the electric signals from photo multiplier tubes to get total fluorescence. Anyone who had some experience in such numerical integration method would not be surprised with negative values. My code faithfully reported some of the negative values reported by the firmware of the machine. However, it caused me some trouble. I could not remember how many times I needed to explain the meaning of such negative values to my colleagues who did not know how such machine and algorithm worked.
To be able to control the acquisition for the flow cytometer is an important step stone for what I wanted to do: monitoring GFP expression of a few proteins in real time. I wrote a few simple program with few simple parts (a chemostats as an incubator and a syringe pump to deliver the sample), I was able to see the time course of gene induction by staving the yeast cells for some essential ingredient. I though it was a great system to really measure how gene circuits response to environmental perturbation to get a better understanding and modeling about gene regulations. We wrote a paper using the system to study a particular amino acid synthesis pathway. Unfortunately, due to the hacky nature of the code and lack of general instrumentation mind set in biological research field. We might be the only group that was doing such thing. Although we can share our code, it will be much easier if the vendor opens their protocols such that we can make sure the code will work across different machines. If so, there might be much more other researchers who will like to use such powerful tools to study gene expression in real time.
Since then, I proudly quitted my postdoc job, and got “real job” in few bio-tech companies working on mostly on data analysis, software and algorithm development. I am hoping that the company I am working for now or in the further will embrace openness and I like to see people hacking a machine beyond its original purpose to generate interesting scientific research results.
Apr 10
William McVey left a message about my code for the “DiggsInAbox” and this. I never really intended to release the code. I don’t mind to share but it is just I have never really motivated to polish the code so it can be released as a professional written code. Anyway, I am pasting my code here so I hope it can be useful for some one. I wrote this code for fun and to learn the algorithm to generate the Treemap. There are a lot of Treemap implementation these days. You might be able to find better code.
DiggsInABox.py
#!/usr/bin/env python
"""
Copyright 2007-2011 Jason Chin, All right reserved
Redistribution and use in source and binary forms, with or without modification,
are permitted provided that the following conditions are met:
* Redistributions of source code must retain the above copyright notice,
this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright notice,
this list of conditions and the following disclaimer in the documentation
and/or other materials provided with the distribution.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 'AS IS'
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE."""
from Treemap import Treemap, Node
class diggTreemap(Treemap):
def __init__(self, rootNode):
self.rootNode = rootNode;
self.setWidthHeight(300,200)
self.setPosition(0,0)
self._cMap = self._colorMap()
def _colorMap(self):
cMap = ["#fc0","#cc0","#f63","#39c","#696","#f93",
"#6f0","#c69","#cf9","#36c","#393","#c0f"]
i = 0
while True:
yield cMap[i]
i = i+1
if i > len(cMap) - 1: i = i % len(cMap)
def writeCSS(self):
return u"""\
<style type="text/css" media="screen">
body {
background: black;
}
div.group2 {
position:absolute;
overflow:hidden;
margin: 3px;
z-index: 1;
}
div.leave {
position:absolute;
text-align:center;
overflow:hidden;
vertical-align:middle;
opacity: 1;
border: 1px outset #000;
margin: 3px;
z-index:2;
}
div.leave > a > div {
display: table-cell;
position: static;
vertical-align: middle;
padding:3px;
color:#000000;
}
div.leave > a {
text-decoration:none;
}
div.leave > a:link > div {
color: black;
}
div.leave > a:visited > div {
color: #bcc;
}
div.index {
position:absolute;
left: 810px;
width: 90px;
text-align:center;
overflow:hidden;
vertical-align:middle;
opacity: 1;
border: 2px #fff outset;
}
div.index > a {
text-decoration:none;
}
div.index > a:link > div {
color: black;
}
div.index > a:visited > div {
color: black;
}
div.index > a > div {
display: table-cell;
position: static;
vertical-align: middle;
padding:3px;
color:#000000;
width: 90px;
-moz-user-select: none;
}
</style>\n"""
def writeJS(self):
return """\
<script language='javascript'>
function hiliteBlock(blockId) {
elm = document.getElementById(blockId);
elm.style.width = parseInt(elm.style.width) - 6 + 'px';
elm.style.height = parseInt(elm.style.height) - 6 + 'px';
elm.style.borderWidth = '3px';
elm.style.borderStyle = 'dashed';
elm.style.borderColor = '#fff';
return true;
}
function removeHilite(blockId) {
elm = document.getElementById(blockId);
elm.style.width = parseInt(elm.style.width) + 6 + 'px';
elm.style.height = parseInt(elm.style.height) + 6 + 'px';
elm.style.borderWidth = '0px';
elm.style.borderStyle = '';
elm.style.borderColor = '';
return true;
}
function showSummary(event,leaveId) {
if (document.getElementById('summaryPopup')) {
popup = document.getElementById('summaryPopup');
document.body.removeChild(popup);
}
//elm = document.getElementById(leaveId);
popup = document.createElement('div');
popup.id = 'summaryPopup';
x = event.pageX;
y = event.pageY;
document.body.appendChild(popup);
popup.style.position = 'absolute';
popup.style.left = x + 'px';
popup.style.top = y + 'px';
popup.style.width = '250px';
popup.style.height = '150px';
popup.style.zIndex = '100';
popup.style.background = '#cec';
//still working on this function, ajax maybe needed
return true;
}
function adjustFont() {
tmElm = document.getElementById('treemap');
for (idx=0; idx < tmElm.childNodes.length; idx++) {
elm = tmElm.childNodes[idx];
if (elm.className != 'leave' && elm.className != 'index') {
continue;
}
while ( elm.scrollWidth > elm.clientWidth ||
elm.scrollHeight > elm.clientHeight) {
curFontsize = parseInt(elm.style.fontSize);
newFontSizeInPx = parseInt(elm.style.fontSize) - 1;
if (newFontSizeInPx <= 2) {
newFontSizeInPx = 2;
elm.style.fontSize = newFontSizeInPx + 'px';
break;
}
elm.style.fontSize = newFontSizeInPx + 'px';
}
}
}
</script>
"""
def writeAll(self):
outStr = u""
outStr += u'<html><head>'
outStr += self.writeCSS()
outStr += self.writeJS()
outStr += u"</head><body onload='adjustFont();'><div style='font-size:30px;color:white;margin-left:40px'>Diggs in a Box<br> <span style='font-size:0.3em;'>(v. 0.02 by Chen-Shan Chin)</span></div>"
outStr += u"<div id='treemap' style='position:relative;left:40px; top:5px;width:%dpx;height:%dpx;'>" % (self.width+5, self.height+5)
#write feed navigator
rssMapLabel = ['All','Technology','Science','Business','Sports',
'Entertainment','Gamming']
outStr += "<div id='nav1' style='position:absolute;right:0px;top:-30px;height:30px;width:600px'>"
outStr += "<table align=right><tr>"
for label in rssMapLabel:
col = '#fff'
if self.rootNode.name == label:
col = '#fff'
else:
col = '#888'
outStr += "<td><a href='http://infoecho.net/Sandbox/DiggsInABox.py?feed=%s' style='text-decoration:none;'><div style='color:%s;border:1px #888 solid;align:right;padding:3px;text-decoration:none;'>%s</div></a></td>" % (label,col,label)
outStr += "</table></div>"
#write index navigator
y = 0;
cMap = self._colorMap()
totalWeight = 0
for node in self.rootNode.children:
if node.weight > 10:
totalWeight += node.weight
else:
totalWeight += 10
dhMap = {}
dhSum = 0
for node in self.rootNode.children:
dh = 1.0 * node.weight / totalWeight * (self.height+3)
if dh > 20:
dhMap[node] = dh
else:
dhMap[node] = 20
dhSum += dhMap[node]
for node in self.rootNode.children:
dh = 1.0 * dhMap[node] / dhSum * (self.height+3)
fs = 12
if dh > 24:
fs *= 1.8
if max([len(w) for w in node.name.split(' ')]) > 10:
fs = (1.0*max([len(w) for w in node.name.split(' ')]) / 10)
if len(node.name.split(' ')) > 1 and dh < 48:
fs *= 0.7
if fs < 12: fs = 12
if fs > 24: fs = 24
outStr += u"<div id='%s_index' class='index' \
style='top:%dpx; height:%dpx; background:%s;font-size:%dpx' onmouseover='hiliteBlock(\"%s\");' onmouseout='removeHilite(\"%s\");'><a target='_blank' href='%s'><div style='height:%dpx;'>%s</div></a></div>"\
% (node.name, y, dh-2, cMap.next(), fs, node.name, node.name, node.properties['link'], dh-2, node.name)
y += dh
#write all nodes
outStr += self.writeNodes(self.rootNode)
outStr += u"</div></body></html>"
return outStr
def writeNodes(self, node):
outStr = self.writeNode(node)
for n in node.children:
outStr += self.writeNodes(n)
return outStr
def writeNode(self, node):
outStr = u""
if "group2" in node.properties:
x,y,dw,dh = node.rect
x = int(round(x))
y = int(round(y))
dw = int(round(dw))
dh = int(round(dh))
color = self._cMap.next()
outStr += u"<div id='%s' class='group2' \
style='left:%dpx; top:%dpx; width:%dpx; height:%dpx; background:%s;'></div>"\
% (node.name, x, y, dw, dh, color)
return outStr
elif node.properties['is_leave'] == False:
return outStr
x,y,dw,dh = node.rect
x = int(round(x+3))
y = int(round(y+3))
dw = int(round(dw-9))
dh = int(round(dh-9))
label = node.properties['data']['title'].strip()
if len(label) > 60:
label = label[:60]+" ..."
if dw > 20 and dh >20:
fs = node.area**0.5 / 7;
if len(label) > 40:
fs *= 0.75
if len(label) < 20:
fs *= 1.25
if max([len(w) for w in label.split(' ')]) * fs * 0.7 > dw:
fs = 2 * dw / (max([len(w) for w in label.split(' ')]))
fs = int(fs)
outStr += u"<div id='%s' class='leave' \
style='left:%dpx; top:%dpx; width:%dpx; height:%dpx;\
font-size:%fpx;'>" % (node.name, x, y, dw, dh, fs)
outStr += u"<a target='_blank' href='%s'><div style='width:%d;height:%d;'>%s</div></a></div>" % (node.properties['data']['link'], dw, dh, label)
else:
outStr += u"<a target='_blank' href='%s'><div id='%s' class='leave' \
style='left:%dpx; top:%dpx; width:%dpx; height:%dpx;'></div></a>" % (node.properties['data']['link'], label, x, y, dw, dh)
return outStr
###########################################################
import cgitb; cgitb.enable()
import cgi
import feedparser
form = cgi.FieldStorage()
print "Content-Type: text/html"
print
feed = ""
if form.has_key('feed'):
feed = form['feed'].value.strip()
rssMap = {'All':'http://digg.com/rss/index.xml',
'Technology':'http://digg.com/rss/containertechnology.xml',
'Science':'http://digg.com/rss/containerscience.xml',
'Business':'http://digg.com/rss/containerworld_business.xml',
'Sports':'http://digg.com/rss/containersports.xml',
'Entertainment':'http://digg.com/rss/containerentertainment.xml',
'Gamming':'http://digg.com/rss/containergaming.xml'}
if feed not in rssMap:
feed = 'All'
data = feedparser.parse(rssMap[feed])
entries = data['entries']
term2entries = {}
term2link = {}
for e in entries:
term = e['digg_category']
if term not in term2entries:
term2entries[term] = []
term2entries[term].append( {'id':e['id'],
'diggcount':int(e['digg_diggcount']),
'link':e['link'],'title':e['title'],
'commentcount':int(e['digg_commentcount']),
'summary':e['summary']} )
if term not in term2link:
term2link[term] = "/".join(e['link'].split('/')[:-1])
rootNode = Node(feed)
for term in term2entries:
n2 = Node(term)
n2.properties['is_leave'] = False
n2.properties['group2'] = True
n2.properties['link'] = term2link[term]
n2.weight = 0
for entry in term2entries[term]:
n3 = Node(entry['id'].split("/")[-1])
n3.properties['is_leave'] = True
n3.properties['data'] = entry
n3.weight = entry['diggcount']
n2.addAChild(n3)
n2.weight = n2.weight + n3.weight
rootNode.addAChild(n2)
rootNode.weight = rootNode.weight + n2.weight
rootNode.properties['is_leave'] = False
#for n in rootNode.children:
# print n.name, [n2.name for n2 in n.children]
rootNode.sortChildrenByWeight()
TM = diggTreemap(rootNode)
TM.setWidthHeight(800,540)
TM.layout()
outStr = TM.writeAll()
print outStr.encode('utf-8')
Treemap.py
#!/usr/bin/env python
"""
Copyright 2007-2011 Jason Chin, All right reserved
Redistribution and use in source and binary forms, with or without modification,
are permitted provided that the following conditions are met:
* Redistributions of source code must retain the above copyright notice,
this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright notice,
this list of conditions and the following disclaimer in the documentation
and/or other materials provided with the distribution.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 'AS IS'
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE."""
import sys
class Node:
def __init__(self, name, weight=1):
self.name = name
self.properties = {}
self.children = []
self.parentNode = []
self.area = 1.0 #need to normalized such that sum(children.area) = the area assigned
self.rect = []
self.weight = weight
def addAChild(self, aNode):
self.children.append(aNode)
aNode.parentNode.append(self)
def addChildren(self, Nodes):
for n in Nodes:
self.children.append(n)
n.parentNode.append(self)
def numOfChildren(self):
return len(children)
def addAProperty(self, pk, pv):
self.properties[pk] = pv
def sortChildrenByWeight(self):
if len(self.children) == 0:
return
tmpNodes = [ [-c.weight, c] for c in self.children]
tmpNodes.sort()
self.children = 1 for c in tmpNodes]
for c in self.children:
c.sortChildrenByWeight()
def normalizeChildrenArea(self, totalArea):
sw = 1.0 * sum([n.weight for n in self.children])
for n in self.children:
n.area = n.weight / sw * totalArea;
class Treemap:
def __init__(self, rootNode):
self.rootNode = rootNode
self.setHeightWidth(300,200)
self.setPosition(0,0)
pass
def setWidthHeight(self, width,height):
self.height = height
self.width = width
def setPosition(self,left,top):
self.left = left
self.top = top
def _worst(self, rowOfNodes, mw):
rowBlockAreas = [n.area for n in rowOfNodes];
s = sum(rowBlockAreas)
rmin = min(rowBlockAreas)
rmax = max(rowBlockAreas)
return max([ (1.0*mw*mw * rmax)/(s*s), (1.0*s*s)/(mw*mw * rmin)] )
def _squarified(self, nodes, rowOfNodes, top, left, mw, mh, layoutDir=1):
if mw > mh:
mw, mh = mh, mw
layoutDir = -layoutDir
nodesToPlot = nodes;
nodesInRow = []
while nodesToPlot:
c = nodesToPlot[0];
if not nodesInRow:
nodesInRow.append(c)
nodesToPlot = nodesToPlot[1:]
continue
if self._worst(nodesInRow, min([mw,mh])) >= self._worst(nodesInRow+1, min([mw,mh])):
nodesInRow = nodesInRow+1
nodesToPlot = nodesToPlot[1:]
continue
else:
dh = 1.0*sum([n.area for n in nodesInRow])/mw;
self._layoutRowOfNodes(nodesInRow, left, top, dh, layoutDir)
if layoutDir == 1:
top = top + dh
else:
left = left + dh
mh = mh - dh;
if mh < mw:
mw, mh = mh, mw
layoutDir = -layoutDir
nodesInRow = []
dh = 1.0*sum([n.area for n in nodesInRow])/mw;
self._layoutRowOfNodes(nodesInRow, left, top, dh, layoutDir)
def _layoutRowOfNodes(self,rowOfNodes,left,top,mh,ld):
x = left
y = top
for n in rowOfNodes:
r = n.area
if ld == 1:
dw = 1.0 * r / mh
dh = 0
n.rect = [x,y,dw,mh]
#print "Rect(%f,%f,%f,%f);" % (x,y,dw,mh)
else:
dw = 0
dh = 1.0 * r / mh
n.rect = [x,y,mh,dh]
#print "Rect(%f,%f,%f,%f);" % (x,y,mh,dh)
x = x + dw
y = y + dh
def _layoutANode(self, aNode, left, top, w, h):
if len(aNode.children)==0: return
aNode.normalizeChildrenArea(w*h)
#print aNode,[n for n in aNode.children]
self._squarified(aNode.children, [], top, left, w, h)
for n in aNode.children:
x,y,w,h = n.rect
self._layoutANode(n, x, y, w, h)
def layout(self):
w = self.width
h = self.height
self.rootNode.area = w * h;
self.rootNode.rect = [0, 0, w, h];
#self.rootNode.normalizeChildrenArea(w*h)
self._layoutANode(self.rootNode, 0, 0, w, h)
def writeAll(self, outputStream=sys.stdout):
self.outputStream = outputStream
self.printNodes(self.rootNode)
pass
def writeNodes(self, node):
self.writeNode(node)
for n in node.children:
self.writeNodes(n)
def writeNode(self, node):
outputStream.write(node.name)
outputStream.write("\n")
outputStream.write(node.rect)
outputStream.write("\n")
pass
class CanvasTreemap(Treemap):
import random
def __init__(self, rootNode):
self.rootNode = rootNode;
self.setWidthHeight(300,200)
self.setPosition(0,0)
def printAll(self):
print "function plotCanvas(cId){"
print "ctx = document.getElementById(cId).getContext('2d');"
print "ctx.globalAlpha=0.2;"
self.printNodes(self.rootNode);
print "}"
def printNodes(self, node):
self.printNode(node)
for n in node.children:
self.printNodes(n)
def printNode(self,node):
x,y,dw,dh = node.rect
x = x+1
y = y+1
dw = dw-1
dh = dh-1
level = node.properties['level']
style = {0:'rgb(255,0,0)', 1:'rgb(255,255,0)', 2: 'rgb(0,255,0)', 3: 'rgb(0,255,255)', 4: 'rgb(255,0,255)'}
print "ctx.strokeStyle = '%s';" % style[int(random.uniform(0,4.9))]
print "ctx.fillStyle = '%s';" % style[int(random.uniform(0,4.9))]
print "ctx.lineCap='round';"
print "ctx.lineWidth = %d;" % (10-2*int(random.uniform(0,4)))
print "ctx.beginPath();"
print "ctx.moveTo(%.1f,%.1f);" % (x, y);
print "ctx.lineTo(%.1f,%.1f);" % (x+dw, y);
print "ctx.lineTo(%.1f,%.1f);" % (x+dw, y+dh);
print "ctx.lineTo(%.1f,%.1f);" % (x, y+dh);
print "ctx.lineTo(%.1f,%.1f);" % (x, y);
if random.uniform(0,1) < 0.5:
print "ctx.stroke();"
else:
print "ctx.fill();"
#print "ctx.stroke();"
class DivTreemap(Treemap):
def __init__(self, rootNode):
self.rootNode = rootNode;
self.setWidthHeight(300,200)
self.setPosition(0,0)
def printAll(self):
print "<html>"
print '<head><style type="text/css" media="screen">'
print """div.node {
position:absolute;
text-align:center;border:2px solid #bbaaaa;
overflow:hidden;
vertical-align:middle;
background:#f0ffff
}"""
print '</style><head>'
print "<body>This is a test.<br/>"
print "<div id='treemap' style='position:relative;left:30px; top:60px;width:%dpx;height:%dpx;'>" % (self.width, self.height)
self.printNodes(self.rootNode);
print "</div></body></html>"
def printNodes(self, node):
self.printNode(node)
for n in node.children:
self.printNodes(n)
def printNode(self,node):
if node.name == "root": return
x,y,dw,dh = node.rect
level = node.properties['level']
fs = node.area**0.5 / 100;
if fs < 0.75: fs = 0.75
print "<div id='%s' class='node' \
style='left:%dpx; top:%dpx; width:%dpx; height:%dpx;\
line-height:%dpx; font-size:%fem;' onclick='alert(\"%s\")'>%s</div>"\
% (node.name, x, y, dw-5, dh-5, dh-5, fs, \
"Do you like to eat "+node.name+"?", node.name)
import random
def testCanvas():
nodes = []
for i in range(0,5):
n1 = Node('node:%d' % i)
n1.properties['level'] = 1
n1.weight = random.uniform(1,20)
for j in range(0,5):
n2 = Node('node:%d-%d' %(i,j))
n2.properties['level'] = 2
n2.weight = random.uniform(1,20)
for k in range(0,10):
n3 = Node('node:%d-%d-%d' % (i,j,k))
n3.properties['level'] = 3
n3.weight = random.uniform(1,20)
for l in range(0,50):
n4 = Node('node:(%d-%d-%d-%d)' % (i,j,k,l))
n4.properties['level'] = 4
n4.weight = random.uniform(1,20)
n3.addAChild(n4)
n2.addAChild(n3)
n1.addAChild(n2)
nodes.append(n1)
root = Node('root')
root.addChildren(nodes)
root.sortChildrenByWeight()
root.properties['level']=0
TM = CanvasTreemap(root)
TM.setWidthHeight(800,800)
TM.layout()
print "<html><head>"
print "<script>"
TM.printAll()
print "</script><head>"
print """<body onload="plotCanvas('c')">
<canvas id="c" width=810 height=810></canvas></body></html>"""
def testDiv():
tagArray = {"apples": 12,
"oranges": 38,
"pears" : 10,
"mangos" : 24,
"grapes" : 18,
"bananas" : 56,
"watermelons" : 80,
"lemons" : 12,
"limes" : 12,
"pineapples" : 15,
"strawberries" : 20,
"coconuts" : 43,
"cherries" : 20,
"raspberries" : 8,
"peaches" : 25
}
nodes = []
for tag in tagArray:
n = Node(tag)
n.weight = tagArray[tag]
n.properties['level']=1
nodes.append(n)
root = Node('root')
root.addChildren(nodes)
root.properties['level']=0
root.sortChildrenByWeight()
TM = DivTreemap(root)
TM.setWidthHeight(800,250)
TM.layout()
TM.printAll()
if __name__ == '__main__':
#testCanvas();
testDiv();
Apr 10
I am still fascinated about the programming style using co-routine. Actually, it is possible to implement the Needleman–Wunsch alignment algorithm by purely message passing fashion. The following code shows how to implement the algorithm using co-routines again. I modify the code from my previous post such that the alignment array itself is also generated dynamically. We can completely remove those setting up loops. This code is also annotated to show how it is done. If any reader is interested and have any comment, I do like to hear.
# @author Jason Chin
#
# Copyright (C) 2011 by Jason Chin
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
"""
This is an example to implement Needleman-Wunsch sequence algorithm using python's
co-routine. One of the most interest aspect of such implementation is that there
is no explicitly loop. You can not find either the "for" nor "while" keywords
in this code. Each alignment cell is a co-routine and the calculation of alignement
score and backtracking that generates the alignment string are done with a message
passing fashion. The alignment cells are also generated in a dynamic way. A
banded alignment can be done by limiting not generate the whole alignment array but
only the banded part of the array.
Is it useful? I am not sure, but it is definitely fun to show it is possible.
--Jason Chin, Apr. 10, 2011
"""
### Set up the alignment score scheme
matchScore, mismatchScore, gapScore = 4, -3, -4
### Two testing string for alingment
seq1 = "TTAAGTGTAGCCTTGTGTGACATGTATTTTTAT"
seq2 = "TTTCTAGGTAGTTGTGGTGAGTTTAGTTGATAT"
### cellMap is a dictionary that maps integer pairs to the co-routines
cellMap = {}
### For tracking the global best alignment cell
globalBestCellScore = [None, -100000]
def getAnAlignCell(x, y, seq1, seq2):
"""
This function returns a co-routine the represents an alignment cell at position
x and y. The alignment strings are passed explicilty for simplicity.
"""
def alnCell():
"""
This is the co-routine for an alignment cell. A alignment cell co-routine is
excuted in roughly two stage. The first stage it collects the alignment score
from the cells at (x-1,y-1), (x-1,y), and (x, y-1) and calculate the best
alignment score. Depending the alignment path through the alignment cell, a new
alignment score is generated and passed to the cells at (x+1, y+1),
(x+1,y), and (x, y+1). If any of those cell has not be generated, it will
generate the co-routine and regisiter them with the cellMap dictionary. After
this it waits for the backtracking caculation. If a cell is in the best alignment
path, it will pass the best alignment pair to next cell in the best alignment
path.
"""
global globalBestCellScore
global cellMap
b1, b2 = seq1[x], seq2[y]
mx, my = len(seq1), len(seq2)
cellData = []
# if the cell is on the top or the left side of the alignment, they only have
# to wait for one other cell to pass in the alignment score. Otherwise, they
# need to collect three messages from those (x-1,y), (x,y-1), and (x-1, y-1)
# before they can do any calculation.
if x == 0 or y == 0:
cellId, s = yield
cellData.append( (cellId, s) )
else:
cellId, s = yield
cellData.append( (cellId, s) )
cellId, s = yield
cellData.append( (cellId, s) )
cellId, s = yield
cellData.append( (cellId, s) )
# find the best cell that gives the best alignment score
cellData.sort( key=lambda x: -x[1] )
bestCell, bestScore = cellData[0]
if bestScore > globalBestCellScore[1]:
globalBestCellScore = [ (x,y), bestScore ]
# pass the new alignment score to (x+1, y+1)
if x+1 < mx and y+1 < my:
# generate the cell at (x+1, y+1) if necessary
if (x+1, y+1) not in cellMap:
cellMap[ (x+1, y+1) ] = getAnAlignCell( x+1, y+1, seq1, seq2 )()
cellMap[ (x+1, y+1) ].next()
if b1 == b2: # a match, seq1[x] == seq[2], new_score = bestScore + matchScore
cellMap[ (x+1, y+1) ].send( ((x,y), bestScore + matchScore) ) # pass the new score to cell (x+1, y+1)
else: # a mismatch, seq1[x] != seq[2], new_score = bestScore + mismatchScore
cellMap[ (x+1, y+1) ].send( ((x,y), bestScore + mismatchScore) ) # pass the new score to cell (x+1, y+1)
# pass the new alignment score to (x+1, y), namely, the base seq1[x] is aligned to a gap
if x+1 < mx:
# generate the cell at (x+1, y) if necessary
if (x+1, y) not in cellMap:
cellMap[ (x+1, y) ] = getAnAlignCell( x+1, y, seq1, seq2 )()
cellMap[ (x+1, y) ].next()
cellMap[ (x+1, y) ].send( ((x,y), bestScore + gapScore) )
# pass the new alignment score to (x, y+1), namely, the base seq2[y] is aligned to a gap
if y+1 < my:
# generate the cell at (x, y+1) if necessary
if (x, y+1) not in cellMap:
cellMap[ (x, y+1) ] = getAnAlignCell( x, y+1, seq1, seq2 )()
cellMap[ (x, y+1) ].next()
cellMap[ (x, y+1) ].send( ((x,y), bestScore + gapScore) )
path = yield # wait, if the cell is on the best path, the co-routine will resume
# generate the alignment pair according the best alinged cells
if bestCell[0] >= 0 and bestCell[1] >=0 :
if path == None:
path = []
if bestCell[0] - x == 0:
c1 = "-"
else:
c1 = seq1[x-1]
if bestCell[1] - y == 0:
c2 = "-"
else:
c2 = seq2[y-1]
path.extend( [ (c1, c2) ] )
# send calculated partial path to the best alingment cell to this cell
cellMap[ bestCell ].send( path )
# return the best path if bestCell[0] = -1 or bestCell[1] = -1
yield path
return alnCell
# initialize the cell at (0,0)
cellMap[ (0,0) ] = getAnAlignCell( 0, 0, seq1, seq2 )()
# prime it
cellMap[(0,0)].next()
# start the whole execution by sending in the initial score to cell at (0,0)
cellMap[(0,0)].send( ( (-1, -1), 0 ) )
# get the best global cell
bestCell = globalBestCellScore[0]
# continue to excute the best cell co-routine to get the alignment path
bestPath = cellMap[bestCell].next()
bestPath.reverse()
# some simple mechinary to print out the alignment path
alnRes = zip(*bestPath)
print "".join(alnRes[0])
print "".join(alnRes[1])
The result:
$ python coAlign_v2.py -TT-AAGTGTAGCCTTGT-GTGACATGTA-TTTTTA TTTCTAG-GTAG--TTGTGGTGA-GTTTAGTTGATA
It might be a totally useless python hack. Yes, it is possible to implement dynamic programming using message passing style python co-routine with the enhanced python generator. Here is the code. I will write some details about how this piece code works. However, the main idea is simple (although you might need some background knowledge about sequence alignment algorithm.) We create a co-routine for each alignment cell. The alignment score is generated by passing the best score around the neighboring cells. The backtracking is also implemented as message passing backward.
matchScore, mismatchScore, gapScore = 4, -5, -3
seq1 = "AGTGTAGTTGTGTGAATGTATTTTTAT"
seq2 = "AGGTAGTTGTGGTGATTTAGTTGATAT"
cellMap = {}
globalBestCellScore = [None, -100]
def getAnAlignCell(x, y, p):
def f():
global globalBestCellScore
global cellMap
b1, b2 = p
cell1Id, s1 = yield
cell2Id, s2 = yield
cell3Id, s3 = yield
cellData = [ (cell1Id, s1), (cell2Id, s2), (cell3Id, s3) ]
cellData.sort( key=lambda x: -x[1] )
bestCell, bestScore = cellData[0]
if bestScore > globalBestCellScore[1]:
globalBestCellScore = [ (x,y), bestScore ]
if x+1 < len(seq1) and y+1 < len(seq2):
if b1 == b2:
cellMap[ (x+1, y+1) ].send( ((x,y), bestScore + matchScore) )
else:
cellMap[ (x+1, y+1) ].send( ((x,y), bestScore + mismatchScore) )
if x+1 < len(seq1):
cellMap[ (x+1, y) ].send( ((x,y), bestScore + gapScore) )
if y+1 < len(seq2):
cellMap[ (x, y+1) ].send( ((x,y), bestScore + gapScore) )
path = yield
if bestCell[0] >= 0 and bestCell[1] >=0 :
if path == None:
path = []
path.extend( [ (x,y) ] )
cellMap[ bestCell ].send( path )
yield path
return f
for x in range(len(seq1)):
for y in range(len(seq2)):
cellMap[ (x,y) ] = getAnAlignCell( x, y, (seq1[x], seq2[y]) )()
cellMap[ (x,y) ].next()
for x in range(len(seq1)):
cellMap[ (x,0) ].send( ( (x, -1), 0 ) )
cellMap[ (x,0) ].send( ( (x-1, -1), 0 ) )
for y in range(len(seq2)):
if y != 0:
cellMap[ (0,y) ].send( ( (-1, y), 0 ) )
cellMap[ (0,y) ].send( ( (-1, y-1), 0 ) )
cellMap[(0,0)].send( ( (-1, -1), 0 ) )
bestCell = globalBestCellScore[0]
bestPath = cellMap[bestCell].next()
bestPath.reverse()
s1 = []
s2 = []
px, py = bestPath[0]
for x,y in bestPath[1:]:
if x - px != 0:
s1.append(seq1[px])
else:
s1.append("-")
if y - py != 0:
s2.append(seq21)
else:
s2.append("-")
px, py = x, y
print "".join(s1)
print "".join(s2)
The result seems to be correct
$ python coAlign.py GTGTAGTTGTGTGAATGTATTT--TT-A G-GTAGTTGTG-G--TG-ATTTAGTTGA
Some thought on the interview puzzles from a dot-com/DNA sequencing data processing company
Oct 22
Posted by Jason Chin in comment, programming, Uncategorized | No Comments
Run into this earlier today. If you are interested in solving computational puzzles, I think they are good ones. It is tempting to write some real code to solve them but it is more interesting and important for me to spend my time analyzing some real data to solve real scientific puzzles this weekend. Anyway, to think and find the clues to the answers is straightforward, especially during and after my morning shower time.
Regardless the verbose description of the questions, the route to solve question is quite straightforward more or less. A real implementation might be slightly more complicated. The following are my tips on “how to solve” these puzzles. (I could be giving wrong answers/tips. If you want a job from DNAnexus, you are on your own.)
(1) Insane in the Membrane
Obviously, the question has nothing really to do with any biological membrane. It is actually more related to “maze solving algorithm” or “finding shortest path” in a graph. One way to solve the problem is to convert the lattice to a graph where each node in the graph represents each empty space (“o”). The edges connect all nodes satisfying the constrain “Each successive position is only 1 nm away from the one before it”. You can start the search with a seed node that only has an edge. Then, using the standard BFS algorithm to find the longest path that connect to the seed node. If there is no node with single edge attached, one can pick any node that is not visited during the graph search as the new seed node. If you find one path that is longer the “Danny Dendrite’s genome”, output the solution. If not, try other seed node until you test all seed nodes. If you can not find any path from any seed node that is longer than the “Danny Dendrite’s genome”, output “impossible”.
(2) Hungry Hungry Coders
Think the “enjoyment values” as a matrix of M by N. What you try to do to find maximum sum when picking one single element per row but non of the element can have a share column. One obvious initial state is to pick the maximum value for each row. If there is no overlapping column, Done. If such assignment is not possible, one need to find potential other assignment which minimize the reduction of the sum. I have not to encounter such problem in my work or research yet, although I can see such algorithm is super-useful on solving practical resource allocation. A quick search shows the “correct” solution is probably the “Hungarian algorithm“. There are several different variants to solve such problem. It would be interested to know which one is mostly efficient especially in the case that there might be a lot of degenerated solutions. Also, it might be interesting to see whether it really works sociologically asking engineers to writing down 10 numbers on 10 menu items for every lunch. It seems quite a way to waste of time.
(3) Genome Search
The major constrain of the problem is “You must not load the entire genome into memory; furthermore, you may read through the genome sequence only once.” Namely, you have to stream the reference genome and doing sequence comparison at the same time. Like all good exact sequencing match problem, using hash values is the way to go. For streaming approach, I think the answer is in some sort of rolling hash. One can calculate the hash values of the “K sequences, each of length M”, into K hash values with one of the rolling hash algorithms. Then, calculate the rolling hash values of the 670G bps sequence with the various lengths of the K sequences while streaming the 670G bps sequence file through the memory. If there is a matched hash value, double check the strings do match and the matched hash value is not due to collision. By the way, good luck on sequencing and doing assembly on the 670 Gbps potentially highly repetitive genome.
Oh, well, while it could be fun to implement these algorithms to see that I indeed get the “correct” answer, I would be more interested to see how they are used in solving real genomics problems. While doing exact match string in a smart way is cool, to build a computing infrastructure to be able to do not-exact matches over and over again is way more useful, e.g., NCBI’s blast server. I assume that is one of what DNAnexus’ directions in the future. I do hope there is indeed a good HPC computation platform to help the scientific communities to solve large data / large analysis problems. In the meantime, I do also believe that the innovation on the DNA detection/sequencing technologies remains one of the most important parts driving biological/medical research to solve important problems. Maybe large data analysis with well-known analytics is only part of the equation. Using data and good analytics to solve basic technology, chemistry, signal processing problems, and to understand the nature of different kind of DNA sequence data is really fun and important.