Open writing projects/Sage and cython a brief introduction - Revision history

Marshall Hampton at 14:13, 30 September 2010

Marshall Hampton — Thu, 30 Sep 2010 14:13:56 GMT

←Older revision		Revision as of 14:13, 30 September 2010
Line 172:		Line 172:
	return m_data		return m_data
	bdigits = RealField(16)		bdigits = RealField(16)
-	~~html(~~"~~<h3>~~ Repressilator simulator ~~</h3>~~")	+	print "Repressilator simulator"
	@interact		@interact
	def rep_run(alpha = slider(0,20,.01,10), beta = slider([10^i for i in srange(-3,1,.1)], default = .1), exp_n = slider(1,10,.1,2.1), t_end = slider([10^i for i in srange(0,3,.1)],default = 2)):		def rep_run(alpha = slider(0,20,.01,10), beta = slider([10^i for i in srange(-3,1,.1)], default = .1), exp_n = slider(1,10,.1,2.1), t_end = slider([10^i for i in srange(0,3,.1)],default = 2)):

Marshall Hampton at 14:12, 30 September 2010

Marshall Hampton — Thu, 30 Sep 2010 14:12:44 GMT

←Older revision		Revision as of 14:12, 30 September 2010
Line 172:		Line 172:
	return m_data		return m_data
	bdigits = RealField(16)		bdigits = RealField(16)
-	html('<h3> Repressilator simulator </h3>')	+	html("<h3> Repressilator simulator </h3>")
	@interact		@interact
	def rep_run(alpha = slider(0,20,.01,10), beta = slider([10^i for i in srange(-3,1,.1)], default = .1), exp_n = slider(1,10,.1,2.1), t_end = slider([10^i for i in srange(0,3,.1)],default = 2)):		def rep_run(alpha = slider(0,20,.01,10), beta = slider([10^i for i in srange(-3,1,.1)], default = .1), exp_n = slider(1,10,.1,2.1), t_end = slider([10^i for i in srange(0,3,.1)],default = 2)):

Marshall Hampton at 14:10, 30 September 2010

Marshall Hampton — Thu, 30 Sep 2010 14:10:57 GMT

←Older revision		Revision as of 14:10, 30 September 2010
Line 82:		Line 82:

	<code>		<code>
-	from Bio import SeqIO	+	from Bio import SeqIO
-	import urllib2 as U	+	import urllib2 as U
-	import scipy.stats as Stat	+	import scipy.stats as Stat
-	f = U.urlopen('http://www.d.umn.edu/~mhampton/HomoSmito.fa')	+	f = U.urlopen('http://www.d.umn.edu/~mhampton/HomoSmito.fa')
-	p = SeqIO.parse(f, 'fasta')	+	p = SeqIO.parse(f, 'fasta')
-	hsmito = p.next().seq.tostring()	+	hsmito = p.next().seq.tostring()
-	f.close()	+	f.close()
-	display_f = RealField(16)	+	display_f = RealField(16)
-	@interact	+	@interact
-	def cgplot(window_length = slider([2^i for i in range(4,12)],default = 2^9)):	+	def cgplot(window_length = slider([2^i for i in range(4,12)],default = 2^9)):
	avg = [16.0hsmito[x-window_length: x+window_length].count('CG')/(2window_length-1) for x in range(window_length, len(hsmito) - window_length)]		avg = [16.0hsmito[x-window_length: x+window_length].count('CG')/(2window_length-1) for x in range(window_length, len(hsmito) - window_length)]
	mean = display_f(Stat.mean(avg))		mean = display_f(Stat.mean(avg))
Line 112:		Line 112:
	Lets look at a toy example of Cython usage. Consider the following python function, which we can imagine is the inner loop of some more complex operation:		Lets look at a toy example of Cython usage. Consider the following python function, which we can imagine is the inner loop of some more complex operation:
	<code>		<code>
-	def slow_logistic(n,a,lamb = 3.82):	+	def slow_logistic(n,a,lamb = 3.82):
	"""		"""
	A slow logistic map.		A slow logistic map.
Line 121:		Line 121:
	return q		return q

-	time a = slow_logistic(10^6,.5)	+	time a = slow_logistic(10^6,.5)

-	Time: CPU 6.17 s, Wall: 6.17 s	+	Time: CPU 6.17 s, Wall: 6.17 s
	</code>		</code>
	So the pure python version takes about 6 seconds for 10^6 iterations. If we tell Sage we'd like to use Cython, and declare variables to C types, we get quite a speedup:		So the pure python version takes about 6 seconds for 10^6 iterations. If we tell Sage we'd like to use Cython, and declare variables to C types, we get quite a speedup:
Line 129:		Line 129:

	<code>		<code>
-	%cython	+	%cython
-	def cython_logistic(int n, float a, float lamb = 3.82):	+	def cython_logistic(int n, float a, float lamb = 3.82):
	cdef float q		cdef float q
	q = a		q = a
Line 137:		Line 137:
	return q		return q

-	time a = cython_logistic(10^6,.5)	+	time a = cython_logistic(10^6,.5)

-	CPU time: 0.04 s, Wall time: 0.04 s	+	CPU time: 0.04 s, Wall time: 0.04 s
	</code>		</code>

Line 147:		Line 147:

	<code>		<code>
-	def m_der(m,p,a,a0,n):	+	def m_der(m,p,a,a0,n):
	return -m + a/(1+p^n) - a0		return -m + a/(1+p^n) - a0
-	def p_der(m,p,b):	+	def p_der(m,p,b):
	return -b*(p-m)		return -b*(p-m)
-	def der_vec(m_list,p_list,a,a0,n,b):	+	def der_vec(m_list,p_list,a,a0,n,b):
	d_vec = []		d_vec = []
	for i in (0,1,2):		for i in (0,1,2):
Line 158:		Line 158:
	d_vec.append(p_der(m_list[i],p_list[i],b))		d_vec.append(p_der(m_list[i],p_list[i],b))
	return d_vec		return d_vec
-	def repressilator_run(m_init = [.1,.2,.3], p_init = [.1,.2,.3], a_val = 6, b_val = .1, n_val = 2.1, t_end = 100):	+	def repressilator_run(m_init = [.1,.2,.3], p_init = [.1,.2,.3], a_val = 6, b_val = .1, n_val = 2.1, t_end = 100):
	a0_val = .01		a0_val = .01
	t_0 = 0		t_0 = 0
Line 171:		Line 171:
	t = t + stepsize		t = t + stepsize
	return m_data		return m_data
-	bdigits = RealField(16)	+	bdigits = RealField(16)
-	html('<h3> Repressilator simulator </h3>')	+	html('<h3> Repressilator simulator </h3>')
-	@interact	+	@interact
-	def rep_run(alpha = slider(0,20,.01,10), beta = slider([10^i for i in srange(-3,1,.1)], default = .1), exp_n = slider(1,10,.1,2.1), t_end = slider([10^i for i in srange(0,3,.1)],default = 2)):	+	def rep_run(alpha = slider(0,20,.01,10), beta = slider([10^i for i in srange(-3,1,.1)], default = .1), exp_n = slider(1,10,.1,2.1), t_end = slider([10^i for i in srange(0,3,.1)],default = 2)):
	test = repressilator_run(a_val = alpha, n_val = exp_n)		test = repressilator_run(a_val = alpha, n_val = exp_n)
	html('beta = ' + str(bdigits(beta)) + '; alpha = ' + str(bdigits(alpha))+ '; exponent = ' + str(bdigits(exp_n)))		html('beta = ' + str(bdigits(beta)) + '; alpha = ' + str(bdigits(alpha))+ '; exponent = ' + str(bdigits(exp_n)))
Line 186:		Line 186:
	<code>		<code>
	%cython		%cython
-	cdef m_der(m,p,a,a0,n):	+	cdef m_der(m,p,a,a0,n):
	return float(-m + a/(1+p**n) - a0)		return float(-m + a/(1+p**n) - a0)
-	cdef p_der(m,p,b):	+	cdef p_der(m,p,b):
	return float(-b*(p-m))		return float(-b*(p-m))
-	cdef der_vec(m_list,p_list,a,a0,n,b):	+	cdef der_vec(m_list,p_list,a,a0,n,b):
	d_vec = []		d_vec = []
	for i in (0,1,2):		for i in (0,1,2):
Line 197:		Line 197:
	d_vec.append(p_der(m_list[i],p_list[i],b))		d_vec.append(p_der(m_list[i],p_list[i],b))
	return d_vec		return d_vec
-	def repressilator_run(m_init = [.1,.2,.3], p_init = [.1,.2,.3], a_val = 6, b_val = .1, n_val = 2.1, t_end = 100):	+	def repressilator_run(m_init = [.1,.2,.3], p_init = [.1,.2,.3], a_val = 6, b_val = .1, n_val = 2.1, t_end = 100):
	cdef float a0_val, t_0, stepsize, t		cdef float a0_val, t_0, stepsize, t
	a0_val = .01		a0_val = .01

Marshall Hampton at 14:09, 30 September 2010

Marshall Hampton — Thu, 30 Sep 2010 14:09:11 GMT

←Older revision		Revision as of 14:09, 30 September 2010
Line 47:		Line 47:

	<code>		<code>
-	def PStoRE(PrositePattern):	+	def PStoRE(PrositePattern):
	"""		"""
	Converts a PROSITE regular expression to a python r.e.		Converts a PROSITE regular expression to a python r.e.
Line 60:		Line 60:
	rePattern = rePattern.replace(')','}')		rePattern = rePattern.replace(')','}')
	return rePattern		return rePattern
-	from Bio import Fasta	+	from Bio import Fasta
-	import re	+	import re
-	import urllib2 as U	+	import urllib2 as U
-	@interact	+	@interact
-	def re_scan(fasta_file_url = 'http://www.d.umn.edu/~mhampton/PlasProtsRef.fa', pat = input_box('G - x - P - [AG] - x(2) - [LIVM] - x - [IV] ', type = str, width = 60)):	+	def re_scan(fasta_file_url = 'http://www.d.umn.edu/~mhampton/PlasProtsRef.fa', pat = input_box('G - x - P - [AG] - x(2) - [LIVM] - x - [IV] ', type = str, width = 60)):
	re_pat = re.compile(PStoRE(pat))		re_pat = re.compile(PStoRE(pat))
	parser = Fasta.RecordParser()		parser = Fasta.RecordParser()

Marshall Hampton at 14:08, 30 September 2010

Marshall Hampton — Thu, 30 Sep 2010 14:08:15 GMT

←Older revision		Revision as of 14:08, 30 September 2010
Line 45:		Line 45:
	The @interact command in the Sage notebook provides an easy way to make simple GUIs to explore data. In the example below, a user can enter the URL of a fasta-formatted protein file and a PROSITE-style regular expression. Using biopython and the "re" module of python we can search for and display matches to the pattern. For this screenshot, I used proteins from the malaria-causing Plasmodium falciparum and a fragment of the transthyretin pattern ([http://ca.expasy.org/cgi-bin/prosite-search-ac?PDOC00617 Prosite PS00768]). This is just for illustrative purposes - I do not claim any significance for these matches.		The @interact command in the Sage notebook provides an easy way to make simple GUIs to explore data. In the example below, a user can enter the URL of a fasta-formatted protein file and a PROSITE-style regular expression. Using biopython and the "re" module of python we can search for and display matches to the pattern. For this screenshot, I used proteins from the malaria-causing Plasmodium falciparum and a fragment of the transthyretin pattern ([http://ca.expasy.org/cgi-bin/prosite-search-ac?PDOC00617 Prosite PS00768]). This is just for illustrative purposes - I do not claim any significance for these matches.

-	~~<nowiki>~~	+
	<code>		<code>
	def PStoRE(PrositePattern):		def PStoRE(PrositePattern):
Line 76:		Line 76:
	print ''		print ''
	</code>		</code>
-	~~</nowiki>~~
	[[Image:rebrowse1.png]]		[[Image:rebrowse1.png]]

Marshall Hampton at 14:07, 30 September 2010

Marshall Hampton — Thu, 30 Sep 2010 14:07:24 GMT

←Older revision		Revision as of 14:07, 30 September 2010
Line 45:		Line 45:
	The @interact command in the Sage notebook provides an easy way to make simple GUIs to explore data. In the example below, a user can enter the URL of a fasta-formatted protein file and a PROSITE-style regular expression. Using biopython and the "re" module of python we can search for and display matches to the pattern. For this screenshot, I used proteins from the malaria-causing Plasmodium falciparum and a fragment of the transthyretin pattern ([http://ca.expasy.org/cgi-bin/prosite-search-ac?PDOC00617 Prosite PS00768]). This is just for illustrative purposes - I do not claim any significance for these matches.		The @interact command in the Sage notebook provides an easy way to make simple GUIs to explore data. In the example below, a user can enter the URL of a fasta-formatted protein file and a PROSITE-style regular expression. Using biopython and the "re" module of python we can search for and display matches to the pattern. For this screenshot, I used proteins from the malaria-causing Plasmodium falciparum and a fragment of the transthyretin pattern ([http://ca.expasy.org/cgi-bin/prosite-search-ac?PDOC00617 Prosite PS00768]). This is just for illustrative purposes - I do not claim any significance for these matches.

		+	<nowiki>
	<code>		<code>
	def PStoRE(PrositePattern):		def PStoRE(PrositePattern):
Line 75:		Line 76:
	print ''		print ''
	</code>		</code>
		+	</nowiki>
	[[Image:rebrowse1.png]]		[[Image:rebrowse1.png]]

Marshall Hampton at 14:06, 30 September 2010

Marshall Hampton — Thu, 30 Sep 2010 14:06:24 GMT

←Older revision		Revision as of 14:06, 30 September 2010
Line 45:		Line 45:
	The @interact command in the Sage notebook provides an easy way to make simple GUIs to explore data. In the example below, a user can enter the URL of a fasta-formatted protein file and a PROSITE-style regular expression. Using biopython and the "re" module of python we can search for and display matches to the pattern. For this screenshot, I used proteins from the malaria-causing Plasmodium falciparum and a fragment of the transthyretin pattern ([http://ca.expasy.org/cgi-bin/prosite-search-ac?PDOC00617 Prosite PS00768]). This is just for illustrative purposes - I do not claim any significance for these matches.		The @interact command in the Sage notebook provides an easy way to make simple GUIs to explore data. In the example below, a user can enter the URL of a fasta-formatted protein file and a PROSITE-style regular expression. Using biopython and the "re" module of python we can search for and display matches to the pattern. For this screenshot, I used proteins from the malaria-causing Plasmodium falciparum and a fragment of the transthyretin pattern ([http://ca.expasy.org/cgi-bin/prosite-search-ac?PDOC00617 Prosite PS00768]). This is just for illustrative purposes - I do not claim any significance for these matches.

-	<~~nowiki~~>	+	<code>
	def PStoRE(PrositePattern):		def PStoRE(PrositePattern):
	"""		"""
Line 74:		Line 74:
	html(matches)		html(matches)
	print ''		print ''
-	</~~nowiki~~>	+	</code>
	[[Image:rebrowse1.png]]		[[Image:rebrowse1.png]]

Line 80:		Line 80:
	Here is another example of the interact command, along with some of Sage's 2D plotting. We take a human mitochondrial DNA sequence and plot the fraction of CG dinucleotides in a window of variable size. This fraction is divided by 16, to normalize it against the expected value in the case of independently, uniformly distributed nucleotides (obviously a simple model in this case).		Here is another example of the interact command, along with some of Sage's 2D plotting. We take a human mitochondrial DNA sequence and plot the fraction of CG dinucleotides in a window of variable size. This fraction is divided by 16, to normalize it against the expected value in the case of independently, uniformly distributed nucleotides (obviously a simple model in this case).

-	<~~syntax type = "python"~~>	+	<code>
	from Bio import SeqIO		from Bio import SeqIO
	import urllib2 as U		import urllib2 as U
Line 96:		Line 96:
	html('Ratio of CG dinucleotides in a window of size ' + str(2*window_length) + ' to the expected fraction 1/16 <BR> in the human mitochondrion.<BR>Mean value: ' + str(mean) + '; standard deviation: ' + str(std))		html('Ratio of CG dinucleotides in a window of size ' + str(2*window_length) + ' to the expected fraction 1/16 <BR> in the human mitochondrion.<BR>Mean value: ' + str(mean) + '; standard deviation: ' + str(std))
	show(list_plot(avg, plotjoined=True), ymin = 0, ymax = max(avg))		show(list_plot(avg, plotjoined=True), ymin = 0, ymax = max(avg))
-	</~~syntax~~>	+	</code>
	[[Image:cgplot.png]]		[[Image:cgplot.png]]

Line 110:		Line 110:

	Lets look at a toy example of Cython usage. Consider the following python function, which we can imagine is the inner loop of some more complex operation:		Lets look at a toy example of Cython usage. Consider the following python function, which we can imagine is the inner loop of some more complex operation:
-	<~~syntax type="python"~~>	+	<code>
	def slow_logistic(n,a,lamb = 3.82):		def slow_logistic(n,a,lamb = 3.82):
	"""		"""
Line 123:		Line 123:

	Time: CPU 6.17 s, Wall: 6.17 s		Time: CPU 6.17 s, Wall: 6.17 s
-	</~~syntax~~>	+	</code>
	So the pure python version takes about 6 seconds for 10^6 iterations. If we tell Sage we'd like to use Cython, and declare variables to C types, we get quite a speedup:		So the pure python version takes about 6 seconds for 10^6 iterations. If we tell Sage we'd like to use Cython, and declare variables to C types, we get quite a speedup:


-	<~~syntax type="python"~~>	+	<code>
	%cython		%cython
	def cython_logistic(int n, float a, float lamb = 3.82):		def cython_logistic(int n, float a, float lamb = 3.82):
Line 139:		Line 139:

	CPU time: 0.04 s, Wall time: 0.04 s		CPU time: 0.04 s, Wall time: 0.04 s
-	</~~syntax~~>	+	</code>

	So "Cythonizing" our loop sped it up by a factor of over 100. Often by converting a few crucial functions to Cython we get all the speed we need.		So "Cythonizing" our loop sped it up by a factor of over 100. Often by converting a few crucial functions to Cython we get all the speed we need.
Line 145:		Line 145:
	For a slightly more real-world example, suppose we want to interactively explore the behavior of a system of ODEs as we change parameters. The example below uses the "repressilator" system of Elowitz and Leibler <cite>Repressilator</cite>, and a crude Euler's method (just for illustrative purposes!) to solve it. The behavior is not too bad for smaller runs but it gets a bit sluggish as we increase the simulation time.		For a slightly more real-world example, suppose we want to interactively explore the behavior of a system of ODEs as we change parameters. The example below uses the "repressilator" system of Elowitz and Leibler <cite>Repressilator</cite>, and a crude Euler's method (just for illustrative purposes!) to solve it. The behavior is not too bad for smaller runs but it gets a bit sluggish as we increase the simulation time.

-	<~~syntax type="python"~~>	+	<code>
	def m_der(m,p,a,a0,n):		def m_der(m,p,a,a0,n):
	return -m + a/(1+p^n) - a0		return -m + a/(1+p^n) - a0
Line 178:		Line 178:
	print jsmath('\\frac{d p_{i}}{dt} = -\\beta (p_i - m_i), \ \\frac{d m_{i}}{dt} = -m_i + \\frac{\\alpha}{1 + p_i^{n}} + \\alpha_0')		print jsmath('\\frac{d p_{i}}{dt} = -\\beta (p_i - m_i), \ \\frac{d m_{i}}{dt} = -m_i + \\frac{\\alpha}{1 + p_i^{n}} + \\alpha_0')
	show(sum([list_plot([[x[2],x[1][i]] for x in test], hue = i/3.0, plotjoined=True) for i in (0,1,2)]), ymin = 0)		show(sum([list_plot([[x[2],x[1][i]] for x in test], hue = i/3.0, plotjoined=True) for i in (0,1,2)]), ymin = 0)
-	</~~syntax~~>	+	</code>
	[[Image:rep1.png]]		[[Image:rep1.png]]

	If we time a typical short run, it takes about .3 seconds. The moderately Cythonized form below takes about .02 seconds - noticeably snappier (all these timings are on the same machine).		If we time a typical short run, it takes about .3 seconds. The moderately Cythonized form below takes about .02 seconds - noticeably snappier (all these timings are on the same machine).

-	<~~syntax type = "python"~~>	+	<code>
	%cython		%cython
	cdef m_der(m,p,a,a0,n):		cdef m_der(m,p,a,a0,n):
Line 210:		Line 210:
	t = t + stepsize		t = t + stepsize
	return m_data		return m_data
-	</~~syntax~~>	+	</code>

	Note that the above code is still using quite a few python objects and methods, and could be sped up further with more effort.		Note that the above code is still using quite a few python objects and methods, and could be sped up further with more effort.

Marshall Hampton at 14:04, 30 September 2010

Marshall Hampton — Thu, 30 Sep 2010 14:04:31 GMT

←Older revision		Revision as of 14:04, 30 September 2010
Line 45:		Line 45:
	The @interact command in the Sage notebook provides an easy way to make simple GUIs to explore data. In the example below, a user can enter the URL of a fasta-formatted protein file and a PROSITE-style regular expression. Using biopython and the "re" module of python we can search for and display matches to the pattern. For this screenshot, I used proteins from the malaria-causing Plasmodium falciparum and a fragment of the transthyretin pattern ([http://ca.expasy.org/cgi-bin/prosite-search-ac?PDOC00617 Prosite PS00768]). This is just for illustrative purposes - I do not claim any significance for these matches.		The @interact command in the Sage notebook provides an easy way to make simple GUIs to explore data. In the example below, a user can enter the URL of a fasta-formatted protein file and a PROSITE-style regular expression. Using biopython and the "re" module of python we can search for and display matches to the pattern. For this screenshot, I used proteins from the malaria-causing Plasmodium falciparum and a fragment of the transthyretin pattern ([http://ca.expasy.org/cgi-bin/prosite-search-ac?PDOC00617 Prosite PS00768]). This is just for illustrative purposes - I do not claim any significance for these matches.

-	<~~syntax type="python"~~>	+	<nowiki>
	def PStoRE(PrositePattern):		def PStoRE(PrositePattern):
	"""		"""
Line 74:		Line 74:
	html(matches)		html(matches)
	print ''		print ''
-	</~~syntax~~>	+	</nowiki>
	[[Image:rebrowse1.png]]		[[Image:rebrowse1.png]]

Marshall Hampton: removed some superfluous comments.

Marshall Hampton — Thu, 30 Sep 2010 14:01:08 GMT

removed some superfluous comments.

←Older revision		Revision as of 14:01, 30 September 2010
Line 22:		Line 22:

	Introduction.		Introduction.
-
-	~~Shockwave video:~~

	<html>		<html>
Line 40:		Line 38:
	</object>		</object>
	</html>		</html>
-
-	~~(Please don't ask me shockwave questions - I am merely posting this for Jose and don't know how it works!)~~

	The notebook starts by default in native Sage mode, which is almost identical to python. However, other environments are available. The screenshot below shows a Sage notebook in [http://www.r-project.org/ R] mode, with some simple R commands executed <cite>R</cite>.		The notebook starts by default in native Sage mode, which is almost identical to python. However, other environments are available. The screenshot below shows a Sage notebook in [http://www.r-project.org/ R] mode, with some simple R commands executed <cite>R</cite>.

Marshall Hampton at 13:24, 7 November 2009

Marshall Hampton — Sat, 07 Nov 2009 13:24:37 GMT

←Older revision		Revision as of 13:24, 7 November 2009
Line 132:		Line 132:

	<syntax type="python">		<syntax type="python">
-	%cython	+	%cython
	def cython_logistic(int n, float a, float lamb = 3.82):		def cython_logistic(int n, float a, float lamb = 3.82):
	cdef float q		cdef float q