{{{id=302| 1+2+3 # testing /// }}} {{{id=301| /// }}}


Sage:  Creating a Viable Open Source Alternative to
Magma, Maple, Matlab, and Mathematica

SACNAS 2012, October 10, 2012

William Stein (Professor of Mathematics, University of Washington)






 

{{{id=194| /// }}} {{{id=174| /// }}} {{{id=177| /// }}} {{{id=171| /// }}} {{{id=170| /// }}} {{{id=101| /// }}}

 

Sage: Mission Statement

"Create A Viable Free Open Source Alternative to Magma, Maple, Mathematica, and Matlab"

 

Sage is not a drop-in replacement: does not run programs written in the custom languages of the Ma's.

Sage is not like Octave (versus Matlab).  

Sage's culture, architecture, programming language, and feel are different than the Ma's. 

{{{id=272| /// }}} {{{id=271| /// }}} {{{id=158| /// }}} {{{id=133| /// }}} {{{id=132| /// }}} {{{id=94| /// }}}

Why not Magma, Maple, Matlab, Mathematica?

  1. Commercial: Expensive for my collaborators and students.  Not community owned.

  2. Closed: Implementation of algorithms often secret

  3. Frustrating: Tight control of development

  4. Copy protection: Hard to run on supercomputer or my new laptop or after my 1-year license expires.

  5. Programming language: All use a special math-only language

  6. Bugs: Bug tracking is secret

  7. Compiler: Lack of compilers for their math-only languages
{{{id=242| /// }}} {{{id=241| /// }}} {{{id=247| /// }}} {{{id=246| /// }}} {{{id=153| /// }}} {{{id=105| /// }}} {{{id=93| /// }}} {{{id=106| /// }}}

  1. Python: Sage uses a mainstream general purpose programming language (with a compiler: Cython)

  2. Distribution: about 100 open source packages (written by you and your colleagues!)

  3. Interfaces: smoothly tie together all these libraries and packages

  4. New library: implements novel algorithms; over a half million lines; worldwide community of several hundred people.

{{{id=180| # Sage uses Python code def foo(n, m): if n > m: print("%s is bigger!"%n) foo(10, 5) /// }}} {{{id=178| foo('william', 'jon') /// }}} {{{id=253| /// }}} {{{id=252| /// }}} {{{id=154| /// }}} {{{id=92| /// }}}

Distribution

... each with their own rich history!
{{{id=134| /// }}} {{{id=91| /// }}}

Hundreds of Sage Developers

(There are at least 247 contributors in 167 different places from all around the world.)

 

{{{id=8| /// }}} {{{id=135| /// }}} {{{id=83| /// }}}

History

See this article for more details about the (pre-)history of Sage.

{{{id=98| /// }}} {{{id=157| /// }}} {{{id=156| /// }}} {{{id=181| /// }}} {{{id=185| /// }}} {{{id=184| /// }}}

Random Question Break

????

Some random examples involving Sage:

{{{id=275| factorial(3) /// }}} {{{id=274| A = random_matrix(QQ, 4); b = A^(-1); b /// }}} {{{id=305| show(b) /// }}} {{{id=273| latex(b) /// }}} {{{id=304| show(integrate(sin(x)*cos(x), x)) /// }}} {{{id=190| G = plot(sin(x)*cos(x), (0, 4), figsize=3) G /// }}} {{{id=260| G.save('a.pdf') /// }}} {{{id=259| a = octave('rand(5)'); a /// }}} {{{id=258| a.eig() /// }}} {{{id=262| /// }}} {{{id=276| /// }}} {{{id=261| /// }}} {{{id=256| /// }}} {{{id=320| /// }}}

How to Use Sage

{{{id=319| /// }}} {{{id=318| /// }}} {{{id=317| /// }}} {{{id=316| /// }}} {{{id=255| /// }}} {{{id=254| /// }}} {{{id=202| /// }}} {{{id=204| /// }}}

Interactive Image Compression

(using numpy) 

{{{id=29| import pylab, numpy X = pylab.imread(DATA + 'sacnas.png') A_image = numpy.mean(X, 2) u,s,v = numpy.linalg.svd(A_image) S = numpy.zeros(A_image.shape) S[:len(s),:len(s)] = numpy.diag(s) n = A_image.shape[0] @interact def svd_image(i = ("Eigenvalues (quality)", (20,(1..A_image.shape[0]//2)))): A_approx = numpy.dot(numpy.dot(u[:,:i], S[:i,:i]), v[:i,:]) g = graphics_array([matrix_plot(A_approx), matrix_plot(A_image)]) show(g, axes=False, figsize=8) html("%sx%s image compressed to %.1f%% of size using %s eigenvalues."%( A_image.shape[0], A_image.shape[1],100*(2.0*i*n+i)/(n*n), i)) /// }}} {{{id=278| /// }}} {{{id=277| /// }}} {{{id=82| /// }}}

Number Theory

{{{id=115| factor(2009201020112012) # makes use of PARI /// }}} {{{id=11| # Jon Bober - Rademacher's formula time number_of_partitions(10^6) /// }}} {{{id=102| @interact def _(n=(25..10000)): plot(prime_pi, 0, n, gridlines='minor').show(figsize=[8,3]) /// }}} {{{id=117| /// }}} {{{id=279| /// }}} {{{id=6| /// }}}

Graph Theory

"Sage's graph theory crushes anything I have met from the point of view of methods implemented.  I would say: 'if you found a proprietary graph library and you are convinced that it is better for your needs than Sage, your license is on me.' But those licenses are really expensive :-D"

{{{id=306| graphs. /// }}} {{{id=216| G = graphs.BuckyBall() G.plot().show(figsize=8) /// }}} {{{id=76| set_random_seed(1) G = graphs.RandomLobster(8, .6, .3) show(G, figsize=7) /// }}} {{{id=77| G.automorphism_group() /// }}} {{{id=159| G.chromatic_number() /// }}} {{{id=196| G.shortest_path(13,20) /// }}} {{{id=74| /// }}} {{{id=113| /// }}} {{{id=288| /// }}} {{{id=287| /// }}} {{{id=299| /// }}} {{{id=297| /// }}}

Cython

{{{id=111| def python_sum(n): s = int(0) for i in xrange(1, n+1): s += i*i return s /// }}} {{{id=123| python_sum(3) /// }}} {{{id=142| time python_sum(2*10^6) /// }}} {{{id=228| timeit('python_sum(2*10^6)') /// }}} {{{id=229| def python_sum2(n): return sum(i*i for i in xrange(1,n+1)) /// }}} {{{id=230| time python_sum2(2*10^6) /// }}} {{{id=120| %cython def cython_sum(long n): cdef long long i, s = 0 for i in range(1, n+1): s += i*i return s /// }}} {{{id=124| cython_sum(3) /// }}} {{{id=141| time cython_sum(2*10^6) /// }}} {{{id=109| timeit('cython_sum(2*10^6)') /// }}} {{{id=122| 165/.663 /// }}} {{{id=140| %cython def cython_sum2(long n): cdef long long i return sum(i*i for i in range(1,n+1)) /// }}} {{{id=233| time cython_sum2(2*10^6) /// }}} {{{id=232| /// }}} Of course, it is better to choose a different algorithm: {{{id=234| var('k, n') factor(sum(k^2, k, 1, n)) /// }}} {{{id=121| def sum2(n): return n*(2*n+1)*(n+1)/6 /// }}} {{{id=146| sum2(2*10^6) /// }}}

Even then, Cython provides a speedup:

{{{id=143| %cython def c_sum2(long long n): return n*(2*n+1)*(n+1)/6 /// }}} {{{id=197| c_sum2(3) /// }}} {{{id=236| c_sum2(2*10^6) /// }}} {{{id=237| n = 2*10^6 timeit('sum2(n)') /// }}} {{{id=238| timeit('c_sum2(n)') /// }}} {{{id=235| 2.01/.218 /// }}} {{{id=240| /// }}} {{{id=239| /// }}}

But at a cost!

{{{id=145| c_sum2(2*10^6) # WARNING: overflow -- it's just like C... /// }}} {{{id=147| n=2*10^6; n*(2*n+1)*(n+1) > 2^63 /// }}} {{{id=284| /// }}} {{{id=283| /// }}} {{{id=137| /// }}}

Solving Equations

Solve a cubic equation:

{{{id=15| x = var('x'); show(solve(x^3 + x - 1==0, x)[0]) /// }}}

Solve a system of two linear equations with one unknown coefficient $\alpha$:

{{{id=17| var('alpha, y') show(solve([3*x + 7*y == 2, alpha*x + 3*y == 8], x,y)[0]) /// }}}

Solve a system of 500 linear equations exactly over the rational numbers:

{{{id=127| n = 500 A = random_matrix(QQ,n,n,num_bound=10, den_bound=10) v = random_matrix(QQ,n,1,num_bound=10, den_bound=10) A.visualize_structure() /// }}} {{{id=148| # IML is used -- http://www.cs.uwaterloo.ca/~astorjoh/iml.html time w = A \ v /// }}} {{{id=150| print w.str() /// }}} {{{id=149| /// }}} {{{id=310| /// }}} {{{id=151| /// }}}

Symbolic Calculus

Symbolic Calculus makes use of Maxima and Ginac under the hood.

{{{id=26| var('x') f = sqrt(1 - x^2) plot(f, thickness=3, color='red', aspect_ratio=1, fill=True, figsize=4) /// }}} {{{id=163| var('t') assume(t+1 > 0) g = integrate(f, (x, -1, t)); show(g) /// }}} {{{id=162| show(g(t=1) - g(t=-1)) /// }}} {{{id=311| /// }}} {{{id=160| /// }}} {{{id=161| /// }}}

Plotting in 2D

{{{id=71| plot(sin(1/x^2), (x,.1,.5), figsize=4) /// }}} {{{id=166| G = sum(circle((random(), random()), random()/3, color=hue(random())) for i in range(25)) G.show(aspect_ratio=1, frame=True, figsize=4) /// }}} {{{id=165| /// }}} {{{id=32| /// }}}

Plotting in 3D

{{{id=30| f(x,y) = sin(x - y) * y * cos(x) plot3d(f, (x,-3,3), (y,-3,3), opacity=.9, color='red', figsize=3) /// }}} {{{id=322| # user viewer='tachyon' on the ipad... f(x,y) = sin(x - y) * y * cos(x) plot3d(f, (x,-3,3), (y,-3,3), opacity=.9, color='red', figsize=3, viewer='tachyon') /// }}} {{{id=168| G = sum(sphere((random(), random(), random()), random()/4, color=hue(random()), opacity=.6) for i in range(20)) G.show(aspect_ratio=1, frame=True, figsize=3) /// }}} {{{id=167| /// }}} {{{id=65| /// }}} {{{id=64| /// }}}

Plotting a 3D Model

See http://www.davidson.edu/math/chartier/Starwars/

{{{id=49| # Yoda! (53,756 vertices) from scipy import io yoda = io.loadmat(DATA + 'yodapose.mat') from sage.plot.plot3d.index_face_set import IndexFaceSet V = yoda['V']; F3=yoda['F3']-1; F4=yoda['F4']-1 Y = IndexFaceSet(F3,V,color=Color('#444444')) + IndexFaceSet(F4,V,color=Color('#007700')) Y = Y.rotateX(-1) Y.show(aspect_ratio=1, frame=False, zoom=1.2) /// }}} {{{id=96| /// }}}

Questions

????

{{{id=314| /// }}} {{{id=313| /// }}}

If you made it this far, next try the Sage Tutorial. 

{{{id=312| /// }}} {{{id=131| /// }}}