Size: 18252
Comment:
|
Size: 2906
Comment:
|
Deletions are marked like this. | Additions are marked like this. |
Line 3: | Line 3: |
Post code that demonstrates the use of the interact command in Sage here. It should be easy for people to just scroll through and paste examples out of here into their own sage notebooks. | Post code that demonstrates the use of the interact command in Sage here. It should be easy to just scroll through and paste examples out of here into their own sage notebooks.If you have suggestions on how to improve interact, add them [:interactSuggestions: here] or email [email protected]. |
Line 5: | Line 5: |
We'll likely restructure and reorganize this once we have some nontrivial content and get a sense of how it is laid out. [[TableOfContents]] |
* [:interact/graph_theory:Graph Theory] * [:interact/calculus:Calculus] * [:interact/diffeq:Differential Equations] * [:interact/linear_algebra:Linear Algebra] * [:interact/algebra:Algebra] * [:interact/number_theory:Number Theory] * [:interact/web:Web Applications] * [:interact/bio:Bioinformatics] * [:interact/graphics:Drawing Graphics] |
Line 10: | Line 16: |
== Profile a snippet of code == {{{ html('<h2>Profile the given input</h2>') import cProfile; import profile @interact def _(cmd = ("Statement", '2 + 2'), do_preparse=("Preparse?", True), cprof =("cProfile?", False)): if do_preparse: cmd = preparse(cmd) print "<html>" # trick to avoid word wrap if cprof: cProfile.run(cmd) else: profile.run(cmd) print "</html>" }}} attachment:profile.png |
|
Line 14: | Line 38: |
Line 23: | Line 46: |
=== A Random Walk === | |
Line 24: | Line 48: |
== Graph Theory == === Automorphism Groups of some Graphs === by William Stein (I spent less than five minutes on this): {{{ @interact def _(graph=['CycleGraph', 'CubeGraph', 'RandomGNP'], n=selector([1..10],nrows=1), p=selector([10,20,..,100],nrows=1)): print graph if graph == 'CycleGraph': print "n (=%s): number of vertices"%n G = graphs.CycleGraph(n) elif graph == 'CubeGraph': if n > 8: print "n reduced to 8" n = 8 print "n (=%s): dimension"%n G = graphs.CubeGraph(n) elif graph == 'RandomGNP': print "n (=%s) vertices"%n print "p (=%s%%) probability"%p G = graphs.RandomGNP(n, p/100.0) print G.automorphism_group() show(plot(G)) }}} attachment:autograph.png == Calculus == === A contour map and 3d plot of two inverse distance functions === by William Stein {{{ @interact def _(q1=(-1,(-3,3)), q2=(-2,(-3,3)), cmap=['autumn', 'bone', 'cool', 'copper', 'gray', 'hot', 'hsv', 'jet', 'pink', 'prism', 'spring', 'summer', 'winter']): x,y = var('x,y') f = q1/sqrt((x+1)^2 + y^2) + q2/sqrt((x-1)^2+(y+0.5)^2) C = contour_plot(f, (-2,2), (-2,2), plot_points=30, contours=15, cmap=cmap) show(C, figsize=3, aspect_ratio=1) show(plot3d(f, (x,-2,2), (y,-2,2)), figsize=5, viewer='tachyon') }}} attachment:mountains.png === A simple tangent line grapher === by Marshall Hampton {{{ html('<h2>Tangent line grapher</h2>') @interact def tangent_line(f = input_box(default=sin(x)), xbegin = slider(0,10,1/10,0), xend = slider(0,10,1/10,10), x0 = slider(0, 1, 1/100, 1/2)): prange = [xbegin, xend] x0i = xbegin + x0*(xend-xbegin) var('x') df = diff(f) tanf = f(x0i) + df(x0i)*(x-x0i) fplot = plot(f, prange[0], prange[1]) print 'Tangent line is y = ' + tanf._repr_() tanplot = plot(tanf, prange[0], prange[1], rgbcolor = (1,0,0)) fmax = f.find_maximum_on_interval(prange[0], prange[1])[0] fmin = f.find_minimum_on_interval(prange[0], prange[1])[0] show(fplot + tanplot, xmin = prange[0], xmax = prange[1], ymax = fmax, ymin = fmin) }}} attachment:tangents.png == Differential Equations == === Euler's Method in one variable === by Marshall Hampton. This needs some polishing but its usable as is. {{{ def tab_list(y, headers = None): ''' Converts a list into an html table with borders. ''' s = '<table border = 1>' if headers: for q in headers: s = s + '<th>' + str(q) + '</th>' for x in y: s = s + '<tr>' for q in x: s = s + '<td>' + str(q) + '</td>' s = s + '</tr>' s = s + '</table>' return s var('x y') @interact def euler_method(y_exact_in = input_box('-cos(x)+1.0', type = str, label = 'Exact solution = '), y_prime_in = input_box('sin(x)', type = str, label = "y' = "), start = input_box(0.0, label = 'x starting value: '), stop = input_box(6.0, label = 'x stopping value: '), startval = input_box(0.0, label = 'y starting value: '), nsteps = slider([2^m for m in range(0,10)], default = 10, label = 'Number of steps: '), show_steps = slider([2^m for m in range(0,10)], default = 8, label = 'Number of steps shown in table: ')): y_exact = lambda x: eval(y_exact_in) y_prime = lambda x,y: eval(y_prime_in) stepsize = float((stop-start)/nsteps) steps_shown = max(nsteps,show_steps) sol = [startval] xvals = [start] for step in range(nsteps): sol.append(sol[-1] + stepsize*y_prime(xvals[-1],sol[-1])) xvals.append(xvals[-1] + stepsize) sol_max = max(sol + [find_maximum_on_interval(y_exact,start,stop)[0]]) sol_min = min(sol + [find_minimum_on_interval(y_exact,start,stop)[0]]) show(plot(y_exact(x),start,stop,rgbcolor=(1,0,0))+line([[xvals[index],sol[index]] for index in range(len(sol))]),xmin=start,xmax = stop, ymax = sol_max, ymin = sol_min) if nsteps < steps_shown: table_range = range(len(sol)) else: table_range = range(0,floor(steps_shown/2)) + range(len(sol)-floor(steps_shown/2),len(sol)) html(tab_list([[i,xvals[i],sol[i]] for i in table_range], headers = ['step','x','y'])) }}} attachment:eulermethod.png == Linear Algebra == === Numerical instability of the classical Gram-Schmidt algorithm === by Marshall Hampton (tested by William Stein, who thinks this is really nice!) {{{ def GS_classic(a_list): ''' Given a list of vectors or a matrix, returns the QR factorization using the classical (and numerically unstable) Gram-Schmidt algorithm. ''' if type(a_list) != list: cols = a_list.cols() a_list = [x for x in cols] indices = range(len(a_list)) q = [] r = [[0 for i in indices] for j in indices] v = [a_list[i].copy() for i in indices] for i in indices: for j in range(0,i): r[j][i] = q[j].inner_product(a_list[i]) v[i] = v[i] - r[j][i]*q[j] r[i][i] = (v[i]*v[i])^(1/2) q.append(v[i]/r[i][i]) q = matrix([q[i] for i in indices]).transpose() return q, matrix(r) def GS_modern(a_list): ''' Given a list of vectors or a matrix, returns the QR factorization using the 'modern' Gram-Schmidt algorithm. ''' if type(a_list) != list: cols = a_list.cols() a_list = [x for x in cols] indices = range(len(a_list)) q = [] r = [[0 for i in indices] for j in indices] v = [a_list[i].copy() for i in indices] for i in indices: r[i][i] = v[i].norm(2) q.append(v[i]/r[i][i]) for j in range(i+1, len(indices)): r[i][j] = q[i].inner_product(v[j]) v[j] = v[j] - r[i][j]*q[i] q = matrix([q[i] for i in indices]).transpose() return q, matrix(r) html('<h2>Numerical instability of the classical Gram-Schmidt algorithm</h2>') @interact def gstest(precision = slider(range(3,53), default = 10), a1 = input_box([1,1/1000,1/1000]), a2 = input_box([1,1/1000,0]), a3 = input_box([1,0,1/1000])): myR = RealField(precision) displayR = RealField(5) html('precision in bits: ' + str(precision) + '<br>') A = matrix([a1,a2,a3]) A = [vector(myR,x) for x in A] qn, rn = GS_classic(A) qb, rb = GS_modern(A) html('Classical Gram-Schmidt:') show(matrix(displayR,qn)) html('Stable Gram-Schmidt:') show(matrix(displayR,qb)) }}} attachment:GramSchmidt.png === Linear transformations === by Jason Grout A square matrix defines a linear transformation which rotates and/or scales vectors. In the interact command below, the red vector represents the original vector (v) and the blue vector represents the image w under the linear transformation. You can change the angle and length of v by changing theta and r. {{{ @interact def linear_transformation(theta=slider(0, 2*pi, .1), r=slider(0.1, 2, .1, default=1)): A=matrix([[1,-1],[-1,1/2]]) v=vector([r*cos(theta), r*sin(theta)]) w = A*v circles = sum([circle((0,0), radius=i, rgbcolor=(0,0,0)) for i in [1..2]]) print jsmath("v = %s,\; %s v=%s"%(v.n(4),latex(A),w.n(4))) show(v.plot(rgbcolor=(1,0,0))+w.plot(rgbcolor=(0,0,1))+circles,aspect_ratio=1) }}} attachment:Linear-Transformations.png === Singular value decomposition === by Marshall Hampton {{{ import scipy.linalg as lin var('t') def rotell(sig,umat,t,offset=0): temp = matrix(umat)*matrix(2,1,[sig[0]*cos(t),sig[1]*sin(t)]) return [offset+temp[0][0],temp[1][0]] @interact def svd_vis(a11=slider(-1,1,.05,1),a12=slider(-1,1,.05,1),a21=slider(-1,1,.05,0),a22=slider(-1,1,.05,1),ofs= selector(['On','Off'],label='offset image from domain')): rf_low = RealField(12) my_mat = matrix(rf_low,2,2,[a11,a12,a21,a22]) u,s,vh = lin.svd(my_mat.numpy()) if ofs == 'On': offset = 3 fsize = 6 colors = [(1,0,0),(0,0,1),(1,0,0),(0,0,1)] else: offset = 0 fsize = 5 colors = [(1,0,0),(0,0,1),(.7,.2,0),(0,.3,.7)] vvects = sum([arrow([0,0],matrix(vh).row(i),rgbcolor = colors[i]) for i in (0,1)]) uvects = Graphics() for i in (0,1): if s[i] != 0: uvects += arrow([offset,0],vector([offset,0])+matrix(s*u).column(i),rgbcolor = colors[i+2]) show(my_mat) html('<h3>Singular value decomposition: image of the unit circle and the singular vectors</h3>') image_ell = parametric_plot(rotell(s,u,t, offset),0,2*pi) graph_stuff=circle((0,0),1)+image_ell+vvects+uvects graph_stuff.set_aspect_ratio(1) show(graph_stuff,frame = False,axes=False,figsize=[fsize,fsize]) }}} attachment:svd1.png == Number Theory == === Continued Fraction Plotter === by William Stein {{{ @interact def _(number=e, ymax=selector([None,5,20,..,400],nrows=2), clr=Color('purple'), prec=[500,1000,..,5000]): c = list(continued_fraction(RealField(prec)(number))); print c show(line([(i,z) for i, z in enumerate(c)],rgbcolor=clr),ymax=ymax,figsize=[10,2]) }}} attachment:contfracplot.png === Illustrating the prime number thoerem === by William Stein {{{ @interact def _(N=(100,(2..2000))): html("<font color='red'>$\pi(x)$</font> and <font color='blue'>$x/(\log(x)-1)$</font> for $x < %s$"%N) show(plot(prime_pi, 0, N, rgbcolor='red') + plot(x/(log(x)-1), 5, N, rgbcolor='blue')) }}} attachment:primes.png === Computing Generalized Bernoulli Numbers === by William Stein (Sage-2.10.3) {{{ @interact def _(m=selector([1..15],nrows=2), n=(7,(3..10))): G = DirichletGroup(m) s = "<h3>First n=%s Bernoulli numbers attached to characters with modulus m=%s</h3>"%(n,m) s += '<table border=1>' s += '<tr bgcolor="#edcc9c"><td align=center>$\\chi$</td><td>Conductor</td>' + \ ''.join('<td>$B_{%s,\chi}$</td>'%k for k in [1..n]) + '</tr>' for eps in G.list(): v = ''.join(['<td align=center bgcolor="#efe5cd">$%s$</td>'%latex(eps.bernoulli(k)) for k in [1..n]]) s += '<tr><td bgcolor="#edcc9c">%s</td><td bgcolor="#efe5cd" align=center>%s</td>%s</tr>\n'%( eps, eps.conductor(), v) s += '</table>' html(s) }}} attachment:bernoulli.png === Computing modular forms === by William Stein {{{ j = 0 @interact def _(N=[1..100], k=selector([2,4,..,12],nrows=1), prec=(3..40), group=[(Gamma0, 'Gamma0'), (Gamma1, 'Gamma1')]): M = CuspForms(group(N),k) print j; global j; j += 1 print M; print '\n'*3 print "Computing basis...\n\n" if M.dimension() == 0: print "Space has dimension 0" else: prec = max(prec, M.dimension()+1) for f in M.basis(): view(f.q_expansion(prec)) print "\n\n\nDone computing basis." }}} attachment:modformbasis.png === Computing the cuspidal subgroup === by William Stein {{{ html('<h1>Cuspidal Subgroups of Modular Jacobians J0(N)</h1>') @interact def _(N=selector([1..8*13], ncols=8, width=10, default=10)): A = J0(N) print A.cuspidal_subgroup() }}} attachment:cuspgroup.png === A Charpoly and Hecke Operator Graph === |
|
Line 326: | Line 51: |
# Note -- in Sage-2.10.3; multiedges are missing in plots; loops are missing in 3d plots | html('<h1>A Random Walk</h1>') vv = []; nn = 0 |
Line 328: | Line 54: |
def f(N = prime_range(11,400), p = selector(prime_range(2,12),nrows=1), three_d = ("Three Dimensional", False)): S = SupersingularModule(N) T = S.hecke_matrix(p) G = Graph(T, multiedges=True, loops=not three_d) html("<h1>Charpoly and Hecke Graph: Level %s, T_%s</h1>"%(N,p)) show(T.charpoly().factor()) if three_d: show(G.plot3d(), aspect_ratio=[1,1,1]) |
def foo(pts = checkbox(True, "Show points"), refresh = checkbox(False, "New random walk every time"), steps = (50,(10..500))): # We cache the walk in the global variable vv, so that # checking or unchecking the points checkbox doesn't change # the random walk. html("<h2>%s steps</h2>"%steps) global vv if refresh or len(vv) == 0: s = 0; v = [(0,0)] for i in range(steps): s += random() - 0.5 v.append((i, s)) vv = v elif len(vv) != steps: # Add or subtract some points s = vv[-1][1]; j = len(vv) for i in range(steps - len(vv)): s += random() - 0.5 vv.append((i+j,s)) v = vv[:steps] |
Line 339: | Line 76: |
show(G.plot(),figsize=7) | v = vv L = line(v, rgbcolor='#4a8de2') if pts: L += points(v, pointsize=10, rgbcolor='red') show(L, xmin=0, figsize=[8,3]) |
Line 341: | Line 81: |
attachment:randomwalk.png | |
Line 342: | Line 83: |
attachment:heckegraph.png === Demonstrating the Diffie-Hellman Key Exchange Protocol === by Timothy Clemans (refereed by William Stein) |
=== 3D Random Walk === |
Line 348: | Line 86: |
def diffie_hellman(button=selector(["New example"],label='',buttons=True), bits=("Number of bits of prime", (8,12,..512))): maxp = 2^bits p = random_prime(maxp) k = GF(p) g = k.multiplicative_generator() a = ZZ.random_element(10, maxp) b = ZZ.random_element(10, maxp) print """ <html> <style> .gamodp { background:yellow } .gbmodp { background:orange } .dhsame { color:green; font-weight:bold } </style> <h2>%s-Bit Diffie-Hellman Key Exchange</h2> <ol style="color:#000;font:12px Arial, Helvetica, sans-serif"> <li>Alice and Bob agree to use the prime number p=%s and base g=%s.</li> <li>Alice chooses the secret integer a=%s, then sends Bob (<span class="gamodp">g<sup>a</sup> mod p</span>):<br/>%s<sup>%s</sup> mod %s = <span class="gamodp">%s</span>.</li> <li>Bob chooses the secret integer b=%s, then sends Alice (<span class="gbmodp">g<sup>b</sup> mod p</span>):<br/>%s<sup>%s</sup> mod %s = <span class="gbmodp">%s</span>.</li> <li>Alice computes (<span class="gbmodp">g<sup>b</sup> mod p</span>)<sup>a</sup> mod p:<br/>%s<sup>%s</sup> mod %s = <span class="dhsame">%s</span>.</li> <li>Bob computes (<span class="gamodp">g<sup>a</sup> mod p</span>)<sup>b</sup> mod p:<br/>%s<sup>%s</sup> mod %s = <span class="dhsame">%s</span>.</li> </ol></html> """ % (bits, p, g, a, g, a, p, (g^a), b, g, b, p, (g^b), (g^b), a, p, (g^ b)^a, g^a, b, p, (g^a)^b) |
def rwalk3d(n=(50,1000), frame=True): pnt = [0,0,0] v = [copy(pnt)] for i in range(n): pnt[0] += random()-0.5 pnt[1] += random()-0.5 pnt[2] += random()-0.5 v.append(copy(pnt)) show(line3d(v,color='black'),aspect_ratio=[1,1,1],frame=frame) |
Line 382: | Line 96: |
attachment:dh.png === Plotting an elliptic curve over a finite field === {{{ E = EllipticCurve('37a') @interact def _(p=slider(prime_range(1000), default=389)): show(E) print "p = %s"%p show(E.change_ring(GF(p)).plot(),xmin=0,ymin=0) }}} attachment:ellffplot.png == Bioinformatics == === Web app: protein browser === by Marshall Hampton (tested by William Stein) {{{ import urllib2 as U @interact def protein_browser(GenBank_ID = input_box('165940577', type = str), file_type = selector([(1,'fasta'),(2,'GenPept')])): if file_type == 2: gen_str = 'http://www.ncbi.nlm.nih.gov/entrez/viewer.fcgi?db=protein&sendto=t&id=' else: gen_str = 'http://www.ncbi.nlm.nih.gov/entrez/viewer.fcgi?db=protein&sendto=t&dopt=fasta&id=' f = U.urlopen(gen_str + GenBank_ID) g = f.read() f.close() html(g) }}} attachment:biobrowse.png === Coalescent simulator === by Marshall Hampton {{{ def next_gen(x, selection=1.0): '''Creates the next generation from the previous; also returns parent-child indexing list''' next_x = [] for ind in range(len(x)): if random() < (1 + selection)/len(x): rind = 0 else: rind = int(round(random()*(len(x)-1)+1/2)) next_x.append((x[rind],rind)) next_x.sort() return [[x[0] for x in next_x],[x[1] for x in next_x]] def coal_plot(some_data): '''Creates a graphics object from coalescent data''' gens = some_data[0] inds = some_data[1] gen_lines = line([[0,0]]) pts = Graphics() ngens = len(gens) gen_size = len(gens[0]) for x in range(gen_size): pts += point((x,ngens-1), hue = gens[0][x]/float(gen_size*1.1)) p_frame = line([[-.5,-.5],[-.5,ngens-.5], [gen_size-.5,ngens-.5], [gen_size-.5,-.5], [-.5,-.5]]) for g in range(1,ngens): for x in range(gen_size): old_x = inds[g-1][x] gen_lines += line([[x,ngens-g-1],[old_x,ngens-g]], hue = gens[g-1][old_x]/float(gen_size*1.1)) pts += point((x,ngens-g-1), hue = gens[g][x]/float(gen_size*1.1)) return pts+gen_lines+p_frame d_field = RealField(10) @interact def coalescents(pop_size = slider(2,100,1,15,'Population size'), selection = slider(-1,1,.1,0, 'Selection for first taxon'), s = selector(['Again!'], label='Refresh', buttons=True)): print 'Population size: ' + str(pop_size) print 'Selection coefficient for first taxon: ' + str(d_field(selection)) start = [i for i in range(pop_size)] gens = [start] inds = [] while gens[-1][0] != gens[-1][-1]: g_index = len(gens) - 1 n_gen = next_gen(gens[g_index], selection = selection) gens.append(n_gen[0]) inds.append(n_gen[1]) coal_data1 = [gens,inds] print 'Generations until coalescence: ' + str(len(gens)) show(coal_plot(coal_data1), axes = False, figsize = [8,4.0*len(gens)/pop_size], ymax = len(gens)-1) }}} attachment:coalescent.png == Miscellaneous Graphics == === Interactive rotatable raytracing with Tachyon3d === {{{ C = cube(color=['red', 'green', 'blue'], aspect_ratio=[1,1,1], viewer='tachyon') + sphere((1,0,0),0.2) @interact def example(theta=(0,2*pi), phi=(0,2*pi), zoom=(1,(1,4))): show(C.rotate((0,0,1), theta).rotate((0,1,0),phi), zoom=zoom) }}} attachment:tachyonrotate.png [[Anchor(eggpaint)]] === Somewhat Silly Egg Painter === by Marshall Hampton (refereed by William Stein) {{{ var('s,t') g(s) = ((0.57496*sqrt(121 - 16.0*s^2))/sqrt(10.+ s)) def P(color, rng): return parametric_plot3d((cos(t)*g(s), sin(t)*g(s), s), (s,rng[0],rng[1]), (t,0,2*pi), plot_points = [150,150], rgbcolor=color, frame = False, opacity = 1) colorlist = ['red','blue','red','blue'] @interact def _(band_number = selector(range(1,5)), current_color = Color('red')): html('<h1 align=center>Egg Painter</h1>') colorlist[band_number-1] = current_color egg = sum([P(colorlist[i],[-2.75+5.5*(i/4),-2.75+5.5*(i+1)/4]) for i in range(4)]) show(egg) }}} attachment:eggpaint.png |
attachment:randomwalk3d.png |
Sage Interactions
Post code that demonstrates the use of the interact command in Sage here. It should be easy to just scroll through and paste examples out of here into their own sage notebooks.If you have suggestions on how to improve interact, add them [:interactSuggestions: here] or email [email protected].
- [:interact/graph_theory:Graph Theory]
- [:interact/calculus:Calculus]
- [:interact/diffeq:Differential Equations]
- [:interact/linear_algebra:Linear Algebra]
- [:interact/algebra:Algebra]
- [:interact/number_theory:Number Theory]
- [:interact/web:Web Applications]
- [:interact/bio:Bioinformatics]
- [:interact/graphics:Drawing Graphics]
Miscellaneous
Profile a snippet of code
html('<h2>Profile the given input</h2>') import cProfile; import profile @interact def _(cmd = ("Statement", '2 + 2'), do_preparse=("Preparse?", True), cprof =("cProfile?", False)): if do_preparse: cmd = preparse(cmd) print "<html>" # trick to avoid word wrap if cprof: cProfile.run(cmd) else: profile.run(cmd) print "</html>"
attachment:profile.png
Evaluate a bit of code in a given system
by William Stein (there is no way yet to make the text box big):
@interact def _(system=selector([('sage0', 'Sage'), ('gp', 'PARI'), ('magma', 'Magma')]), code='2+2'): print globals()[system].eval(code)
attachment:evalsys.png
A Random Walk
by William Stein
html('<h1>A Random Walk</h1>') vv = []; nn = 0 @interact def foo(pts = checkbox(True, "Show points"), refresh = checkbox(False, "New random walk every time"), steps = (50,(10..500))): # We cache the walk in the global variable vv, so that # checking or unchecking the points checkbox doesn't change # the random walk. html("<h2>%s steps</h2>"%steps) global vv if refresh or len(vv) == 0: s = 0; v = [(0,0)] for i in range(steps): s += random() - 0.5 v.append((i, s)) vv = v elif len(vv) != steps: # Add or subtract some points s = vv[-1][1]; j = len(vv) for i in range(steps - len(vv)): s += random() - 0.5 vv.append((i+j,s)) v = vv[:steps] else: v = vv L = line(v, rgbcolor='#4a8de2') if pts: L += points(v, pointsize=10, rgbcolor='red') show(L, xmin=0, figsize=[8,3])
attachment:randomwalk.png
3D Random Walk
@interact def rwalk3d(n=(50,1000), frame=True): pnt = [0,0,0] v = [copy(pnt)] for i in range(n): pnt[0] += random()-0.5 pnt[1] += random()-0.5 pnt[2] += random()-0.5 v.append(copy(pnt)) show(line3d(v,color='black'),aspect_ratio=[1,1,1],frame=frame)
attachment:randomwalk3d.png