| Download

All published worksheets from http://sagenb.org

Project: sagenb.org published worksheets

Path: pub / 2201-2301 / 2233-MSRI: compressed using %s eigenvalues.sagews

Views: ¹⁶⁸⁷⁵⁴
Image: ubuntu2004

Creating an Open Source Alternative to
Magma, Maple, Mathematica, and MATLAB

William Stein, Professor of Mathematics, University of Washington

Mission Statement

Create an open source alternative to Magma, Maple, Mathematica, and Matlab.

Firefox <--> Internet Explorer, Opera

Open Office, Latex <--> Microsoft Office

Linux <--> Microsoft Windows

PostgreSQL, MySQL <--> Oracle, Microsoft SQLserver

GIMP <--> Photoshop

Sage <--> Magma, Maple, Mathematica, Matlab

History of Sage

I started Sage at Harvard in January 2005.
Sage-1.0 released February 2006 at Sage Days 1 (UC San Diego).
Nearly 30 Sage Days Workshops (!) at UCLA, UW, Cambridge, Bristol, Austin, France, San Diego, Seattle, MSRI, ..., Barcelona, ...
Sage won first prize in the Trophees du Libre (November 2007)
Funding from Microsoft, Univ of Washington, UC San Diego, NSF, DoD, Google, Sun, private donations, etc.
Hundreds of other people subsequently got involved in Sage's development, and the scope of the project has widened to cover all mathematical computation. There is interest from all areas of mathematics, physics, engineering, etc. (Developer mailing list has 1184 subscribers and about 50 messages/day... and sage-support has 1718 subscribers.)
About 8,000 downloads per month. (Probably most users do not download sage themselves.)

What is Sage?

Sage = Python + Math
A unified self-contained distribution of open source mathematical software.
Nearly a half million lines of new Python (and Cython) code/documentation that implements new capabilities and algorithms.
Over 100,000 lines of automated tests.
A "cloud" application like GMail or Google Docs: http://sagenb.org (over 33,000 accounts); of course, Sage also runs on your desktop.

Tour of the http://sagemath.org website

Google page search
Documentation
Live Chat: http://sagemath.org/help-irc.html
Downloading and Installing Sage
Developer map

Sage is about building the car instead of reinventing the wheel

Sage uses Python, a mainstream programming language, instead of inventing a custom mathematics language
Use straightforward method to link programs together -- C library and pseudotty's, instead of XML servers/OpenMath. We implement all conversion routines, instead of expecting upstream to do it: we make them communicate with Sage, whether they want to or not.
Give copious credit to contributors and be very developer friendly (easily build from source).
Reuse, improve, and contribute to existing libraries and projects (e.g., Scipy, Numpy, R, ATLAS, CVXopt, GSL), instead of starting over and competing with them.
Make the GUI using a web browser: the world of Java and Javascript is available and Sage integrates with the web.

And now a demo...

Interfaces

2 + 3

5

The following opens a single persistent running copy of Mathematica:

%mathematica
f := Sin[x^2]

%mathematica
Integrate[ f, x ]

Pi                2
         Sqrt[--] FresnelS[Sqrt[--] x]
              2                 Pi

Lisp, Gap, etc.

%lisp
(+ 2 2)

4

%gap
2 + 2

4

Moving from Mathematica to Maple (via Sage):

for i in range(10):
    mathematica.eval('Expand[ (x+%s)^2 ]'%i)

2
         x
                    2
         1 + 2 x + x
                    2
         4 + 4 x + x
                    2
         9 + 6 x + x
                     2
         16 + 8 x + x
                      2
         25 + 10 x + x
                      2
         36 + 12 x + x
                      2
         49 + 14 x + x
                      2
         64 + 16 x + x
                      2
         81 + 18 x + x

mathematica.eval('f := x^2 + 1')

mathematica('f')

1 + x^2

g = mathematica('f')

type(g)

<class 'sage.interfaces.mathematica.MathematicaElement'>

g = mathematica('f').sage(); g

x^2 + 1

type(g)

<type 'sage.symbolic.expression.Expression'>

maple(g)

x^2+1

You can also open a connection to a remote copy of Mathematica (or Maple, etc.):

mathematica2 = Mathematica(server='sage.math.washington.edu')

No remote temporary directory (option server_tmpdir) specified, using /tmp/ on sage.math.washington.edu

mathematica2('2+3')

5

The above just happened in Seattle. But otherwise it is nearly indistinguishable from working locally.

Very Big Numbers

Very big numbers: In less than a second, Sage exactly computes $(10^6)!$ , which has over 5 million digits.

time n = factorial(10^6)

Time: CPU 0.68 s, Wall: 0.70 s

n.ndigits()

5565709

Multiplying huge numbers is also fast:

time m = n * (n+1)

Time: CPU 0.38 s, Wall: 0.39 s

and a matrix with a big determinant:

ZZ

Integer Ring

a = random_matrix(ZZ,200, x=2^128); a

200 x 200 dense matrix over Integer Ring (type 'print a.str()' to see all of the entries)

matrix_plot(a)

a[0,0]

78356817574460480744974597920662893126

time a.determinant()

-4233443366271166716754348275917633475423721960507493691709596122653827719580313904897658084879068634101493921864950543403842524209732714287233079280388377315897132340402688046633477042704953523834864346680068972301433505458966931902263581499719026625512795012963048168566970049121913635581632782072617842614181479990815688380733844006508509666987549120493128423270162870832419106637625589095550906469075902067640271467668502157621091646689938913265360039374714745888376194946027250439151644150609714002698788287825066110365670496039739153224087469601080808408819678292824620793794034643627850999330892145760899879428654646896687791561038321507616746511314031451739368427965338548872410701324281865929758073443960916483246671771498418183634916034242158323347240197054999635631364443684086301687896541545428625927450120241992308542677627358363365362220295243060658884402430290920293692926441666064271022985025648492760576232853933428499711485920843885715379416880189646942995432715971564671668086038695163047034139714916552115600507157858084323450005368652734108361957309697776117039244000234706961624045796299044270419543029813446117045435799708295317122272561841406044201863010039294830427991056086113361518508365839666234316620048044530872671901147649952766130331416332294146729187814898999447041143947424807866581250719876119770462117654063368870286812709650770946564180979638110808678385611600247314377083467091256358946092536806877100789896899724232566220607750445728573756335250141126316338213563874054327091366095615728335744629153241840706809146770644873115878664816415815068392324113608342343922819848307137102615089566712384301030728396410484191054099791260900352361137480866026845166657235964798745728683981150263071851547941151628174892396522361852495912237541907808414596641192444714570199635458250848404160579084240120984557749349177014680741231056873678205842717693962365620188600093149915520844339087659239298180655774548579457196440402637036418560070039032671453538322122154556594325868034644945989482171592361675910032303294298747270721689881267768850915584972982765648828947304417297898442413147160858742052929095470887599221576606340543823144435126194206411523084135072741626833093447519450129285831933489266877187772996441141701073047301966692929275571895413902397908577818758316286869305731064844201673629789712759490172873151382684249890395479948345232446753915292077000071699878375796368241480761504711398570425586415529276933025168584248292106607315400328646253652606543713463012676318027288216954681293481313619467034212191661676772675967478398583441077969491512889444154106241660955431994008808123275678188173462604083654678689236439684777695690956542393829582242044458695869371573021834559467946248801283571879522426942522337888973320228240991277804633018072369505475839446592772753438041429415707183203466921367497788350317057173559189252810095439971338711028003206591430495532838012852710398572581924552648495776801966065209690377568843526673161114373727525218201161576184720971431150965974456296758622667778641503671894364986884507685260258923841641390372934880474886459544635248534923737471231400273282667279142889829823086182253138793813636015889712988403124652745992119484099001007174466801357188101031007924970878712938608470913313496581746320211422563184406969637319613259184311899656614912694864640866430832913476966021147079595924417300712082146224045826371046587300053556397440561552658588654032528972462037480677948053956496235884290759719277816326879554737432195461933378592501600807692622835183829057370073217134133120063816571822051527337808413868155746791495855510754035137868395160960000352020496182074638619166275249163703843935786878543773096463013849811198092744437953857199740004632679831006474038252538307624351252476083255340168517850795850338841113997731568089841452570019644659667945681095462042268184395773784635883494820521385861701577323883055436076313293889572967897223321849220386522402681694813027455190644644491047348604219139612476672146101509547491741821742385666320479616297662276517336058601562328033191600637328393157041734830417388169998365403701323468603814968770931949007304715462505720242751094352556287224507457106570202800494793217855501633872707969410884521518973745764962135999322979675984579140668552083312165505994084811214119910222202955734255562913625817567229045154433606397214533144873617369231944339882834487124665878207306911086205506885392643003955213229427405235640780444015386702764777751208238200230283916559697890941452905846872158480287961561049840696008326831288701465300343280320195728616211160541356027321056447451427991505633700072748697642373707729725658525447288357141430292752438191542873989306280358796514505003880122674548614432919305671678195515235978019685020969138785238192222245733419407285734316570344638411818766283130074411173184656019987067537933482887777275381291534907773326957276867410898732401180542861385270951262922812436684879953979646161530319619151041014916082297038295012494986473324596579880181290910303648049703545256557256449693532144714564428730217612713970981545417082542158417969346198728273616945310531804733276366338360394581004729394341609395895759115939767326508414351463158260492220593932361573099478819211776702751253124816747221853326308702101181369073448138379705404416621474418419088799275276405466916156359420533480115552967876054816898375271809263685056067724151423976640745344796372199522415988983523859955471692427723394031183078143979285672650488573031342399448558685844879373709990977867208754118985570039101140632630069524850600494672558101253937028407463620599495334230959730099603202293501500047051885890367576201855232416870993089602642990135125622413836657214173160952922630712346070733145106748024654895744286552775225664696882135863705322400140875301849830561899036333484358662055461542253160423026595455937904749584319566783655539880623857655563702413587199677828542190143854907582895282029300576770029408259137221541500505212374968003030256521927450617915671398175054068903829235011108094051427407375760954569379645467835260532751684927987691837538502040481973110018890714209220107715146575204282913909889792698541980921351295535827120868476398792407648457466833858374790606644503738516741557872013888854386282250585054538862870286810906701084031015688559560262693586031975033822509131298072257493757563980656551896966557117752488543143963974168686381560929315252509208153193778662815682030325138715854308778977768553315147297803857263306065804012024190936523101052739035774455779863347121826801531511995802400942818640427061849361602313398360868481178301049784634179035305402367346360448861349729923029633825651107699573507289218517291054066945314247556396293643952291573542403798180810677333184197152923727578768691581874771137498356303241817118822542485742023909157694299089692931344013481874792364018022509235637387385788077865251515739898426890464236149932609825446722742148315523082759992202714369208465040995266556159349161086775129198686407521415114480491475507147597535644226279761160356236219101552987256164615375722059630772680652483032614273962235466449522247355591830704011051290465234525523424173226398560086291133373915967439737836048538215514793139865402013701022456582951482474682706614416137998807182825582669537709237176729243503270860478115337011708191925974711746517759275594528657601045097768266046522303621402606054676071626378511347336321382182967592845951875688270883364092140566392852660888532759201334629643197435355227810774889367211260012495619396355190238621882822256308228060237826471104972884663427338962561481432858800280991514216476300419951304116080914507416148173508550891295041781622325314200696809729955500989672077177976355633587149318196518237058388223724954444291506283439112041318630812402375510857306
Time: CPU 2.37 s, Wall: 2.50 s

Symbolic Calculus

Sage does Calculus:

reset()

f = sin(x)^2*cos(x)*exp(x)
show(f)

\newcommand{\Bold}[1]{\mathbf{#1}}e^{x} \sin\left(x\right)^{2} \cos\left(x\right)

g = f._fast_float_(x)

g.op_list()

['load 0', 'call exp(1)', 'load 0', 'call sin(1)', 'dup', 'mul', 'mul', 'load 0', 'call cos(1)', 'mul']

plot(f, 0, 3, thickness=5, color='brown')

show(f.taylor(x, 0, 5))

\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{2}{3} \, x^{5} - \frac{1}{3} \, x^{4} + x^{3} + x^{2}

show(integrate(f, x))

\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{3}{40} \, e^{x} \sin\left(3 \, x\right) + \frac{1}{8} \, e^{x} \sin\left(x\right) - \frac{1}{40} \, e^{x} \cos\left(3 \, x\right) + \frac{1}{8} \, e^{x} \cos\left(x\right)

Graph Theory

Sage can do graph theory:

g = graphs.CubeGraph(5); g

5-Cube: Graph on 32 vertices

g.plot()

g = graphs.FlowerSnark(); g

Flower Snark: Graph on 20 vertices

graph_editor(g)

g.plot(graph_border=True)

Sage contains many unique and deep algorithms:

G = g.automorphism_group(); G

Permutation Group with generators [(2,11)(3,12)(4,13)(5,8)(6,9)(7,10)(16,19)(17,18), (1,4,7,10,13)(2,5,8,11,14)(3,6,9,12,20)(15,16,17,18,19), (1,14)(2,4)(5,7)(8,10)(11,13)]

G.order()

20

g.cliques_maximal()

[[0, 1], [0, 14], [0, 15], [2, 1], [2, 3], [2, 7], [3, 4], [3, 16], [4, 5], [4, 14], [5, 6], [5, 10], [6, 7], [6, 17], [7, 8], [8, 9], [8, 13], [9, 10], [9, 18], [10, 11], [11, 1], [11, 12], [12, 13], [12, 19], [13, 14], [16, 15], [16, 17], [17, 18], [18, 19], [19, 15]]

Statistics

Sage includes R, scipy.stats, and GSL (=Gnu Scientific Library).

reset('r')

%r

D <- data.frame(x=c(1,2,3,1), y=c(7,19,2,2))

# Sort on x
indexes <- order(D$x)
D[indexes,]

%r
# Print out sorted dataset, sorted in reverse by y
D[rev(order(D$y)),]

There is also a C library interface to R:

import rpy2.robjects
str(rpy2.robjects.r('mean( c(1,2,5,19) )'))

'[1] 6.75'

Scipy is a large Python library included in Sage with statistics functions:

import scipy.stats

scipy.stats.

scipy.stats.norm

<scipy.stats.distributions.norm_gen object at 0x1102cba90>

GSL is a C library in Sage, which also has a lot of highly optimized statistical functions. It can only be used from Cython, though.

Sage also has its own native statistics functionality:

t = stats.TimeSeries(10^6).randomize('normal')
plot(t.sums(), frame=True, figsize=[8,3])

t.standard_deviation()

1.0001354251864949

t.plot_histogram()

mean(t)

-0.00019249918463875795

std(t)

0.99978723260878721

show(std([1,pi,e]))

\newcommand{\Bold}[1]{\mathbf{#1}}\sqrt{\frac{1}{18} \, {\left(\pi - 2 \, e + 1\right)}^{2} + \frac{1}{18} \, {\left(\pi + e - 2\right)}^{2} + \frac{1}{18} \, {\left(2 \, \pi - e - 1\right)}^{2}}

Sage has a new Generalized Hidden Markov Models Library.

This is a complete new implementation in Sage -- useful for math biology, computer learning, etc.:

hmm.

m = hmm.DiscreteHiddenMarkovModel([[0.4,0.5,0.1],[0.1,0.8,0.1],[.2,.3,.5]], [[0.1,0.9],[0.5,0.5],[1,0]],
       [.2,.5,.3], emission_symbols=['a','b']); m

Discrete Hidden Markov Model with 3 States and 2 Emissions
Transition matrix:
[0.4 0.5 0.1]
[0.1 0.8 0.1]
[0.2 0.3 0.5]
Emission matrix:
[0.1 0.9]
[0.5 0.5]
[1.0 0.0]
Initial probabilities: [0.2000, 0.5000, 0.3000]
Emission symbols: ['a', 'b']

m.baum_welch(['a','b','b','a','b','b','a','b','b'])

(-4.6196788513299064e-05, 20)

Discrete Hidden Markov Model with 3 States and 2 Emissions
Transition matrix:
[ 2.30981274908e-05     0.999976901873  1.89806678835e-23]
[7.14454696757e-221 1.98424978494e-232                1.0]
[               1.0                0.0                0.0]
Emission matrix:
[6.75306548194e-69               1.0]
[              0.0               1.0]
[              1.0               0.0]
Initial probabilities: [0.0000, 0.0000, 1.0000]
Emission symbols: ['a', 'b']

m.sample(10)

['a', 'b', 'b', 'a', 'b', 'b', 'a', 'b', 'b', 'a']

g = m.graph()

show(plot(m.graph(),graph_border=1), figsize=3)

A  = [[0.5,0.5],[0.5,0.5]]
B  = [[(0.9,(0.0,1.0)), (0.1,(1,10000))],[(1,(1,1)), (0,(0,0.1))]]
m = hmm.GaussianMixtureHiddenMarkovModel(A, B, [1,0]); m

Gaussian Mixture Hidden Markov Model with 2 States
Transition matrix:
[0.5 0.5]
[0.5 0.5]
Emission parameters:
[0.9*N(0.0,1.0) + 0.1*N(1.0,10000.0), 1.0*N(1.0,1.0) + 0.0*N(0.0,0.1)]
Initial probabilities: [1.0000, 0.0000]

plot(m.sample(1000).sums())

Fortran!

Example here: http://sagenb.org/home/pub/1708/

Matrices

Exact linear algebra: great support for very sophisticated algorithms over $\mathbf{Z}$ , finite fields, $\mathbf{F}_2$ , etc.

a = random_matrix(GF(2),10000,20000)
#a.visualize_structure()

time b = a.echelon_form()
#b.visualize_structure()

Time: CPU 4.50 s, Wall: 4.70 s

And numerical linear algebra:

reset()

import pylab
import numpy
MSRI = pylab.imread(DATA + 'msri.png')
    
A_image = numpy.mean(MSRI, 2)
u,s,v = numpy.linalg.svd(A_image)
S = numpy.zeros( A_image.shape )
S[:len(s),:len(s)] = numpy.diag(s)

@interact
def svd_image(i = ("Eigenvalues (quality)",(20,(1..A_image.shape[0])))):
    # Thus the following expression involves only the first i eigenvalues,
    # and the first i rows of v and the first i columns of u.
    # Thus we use i*m + i*n + i = i*(m+n+1) storage, compared to the
    # original image that uses m*n storage.
        
    A_approx = numpy.dot(numpy.dot(u[:,:i], S[:i,:i]), v[:i,:])
    
    # Now plot both matrices using Sage's matrix_plot command.
    g = graphics_array([matrix_plot(A_approx), matrix_plot(A_image)])
    
    # And show the plot, with no axes:
    show(g, axes=False, figsize=(8,3))
    
    # And add a caption
    html("<h2>MSRI: compressed using %s eigenvalues</h2>"%i)
    html("Compressed to %.1f%% of original size."%(100*(2.0*i*n+i)/(n*n)))

Very Deep, Cutting Edge Mathematics

Number theory
Combinatorics
Algebraic topology

Number theory - Verify that the 5-part of the Shafarevich-Tate group of a rank 2 curve is trivial:

S = EllipticCurve('389a').sha(); S

Shafarevich-Tate group for the Elliptic Curve defined by y^2 + y = x^3 + x^2 - 2*x over Rational Field

S.p_primary_bound(5)

0

Algebraic topology -- Examples of computing homology groups:

S = simplicial_complexes.Sphere(2) # the 2-sphere
S

Simplicial complex with vertex set (0, 1, 2, 3) and facets {(0, 2, 3), (0, 1, 2), (1, 2, 3), (0, 1, 3)}

S.homology()

{0: 0, 1: 0, 2: Z}

X = simplicial_complexes.SurfaceOfGenus(3); X

Simplicial complex with 15 vertices and 38 facets

X.homology()

{0: 0, 1: Z^6, 2: Z}

X.graph().plot()

FemHub -- a customized Sage

See http://femhub.org/

3-Dimensional Plots

Sage can draw 3d plots:

var('x,y')
b = 2.2
show((plot3d(sin(x^2-y^2),(x,-b,b),(y,-b,b), opacity=.7) + 
    plot3d(0, (x,-b,b), (y,-b,b), color='red')))

var('x,y,z')
T = golden_ratio
p = 2 - (cos(x + T*y) + cos(x - T*y) + cos(y + T*z) + cos(y - T*z) + cos(z - T*x) + cos(z + T*x))
r = 4.78
implicit_plot3d(p, (x, -r, r), (y, -r, r), (z, -r, r), plot_points=50)

Sage can plot Yoda:

from scipy import io
x = io.loadmat(DATA + 'yodapose.mat', struct_as_record=True)
from sage.plot.plot3d.index_face_set import IndexFaceSet
V = x['V']; F3 = x['F3']-1; F4 = x['F4']-1
Y = (IndexFaceSet(F3, V, color = Color('#00aa00')) + 
    IndexFaceSet(F4, V, color = Color('#00aa00')))
Y = Y.rotateX(-1)
Y.show(aspect_ratio = [1,1,1], frame = False, figsize = 4)

1/0

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "_sage_input_43.py", line 10, in <module>
    exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("MS8w"),globals())+"\\n"); execfile(os.path.abspath("___code___.py"))' + '\n', '', 'single')
  File "", line 1, in <module>
    
  File "/private/var/folders/FE/FEo498bGEOeewp0B4AIIP++++TM/-Tmp-/tmpacRLtw/___code___.py", line 3, in <module>
    exec compile(u'_sage_const_1 /_sage_const_0 ' + '\n', '', 'single')
  File "", line 1, in <module>
    
  File "element.pyx", line 1527, in sage.structure.element.RingElement.__div__ (sage/structure/element.c:11973)
  File "integer.pyx", line 1560, in sage.rings.integer.Integer._div_ (sage/rings/integer.c:11163)
  File "integer_ring.pyx", line 279, in sage.rings.integer_ring.IntegerRing_class._div (sage/rings/integer_ring.c:5022)
ZeroDivisionError: Rational division by zero

reset()

@interact
def f(z=sin(x), n=(1..13), m=range(1,13)):
    1/0
    plot(z(x=n*m*x)).show()

And if you're really serious, the language of Sage is Python, so you can also use:

Mayavi: Python interface to VTK
Chaco: interact 2d graphics system developed at Enthought (local Austin company)
VisIt: https://wci.llnl.gov/codes/visit/manuals.html (python bindings)
ParaView: http://www.cmake.org/Wiki/ParaView/Python_Scripting
VPython: http://vpython.org/
PyGame: http://www.pygame.org/news.html

Key take-away points:

FOSS: 100% free open source GPL software: good for clusters, sharing, research, collaboration
Python: mainstream, scientific-computing friendly programming language
Cython: optimized Python to C/C++ compiler we develop with support for C/C++ datatypes
Interfaces: to control Matlab, Octave, Mathematica, etc., from Sage
Self contained: can have many copies of Sage at once, so no "dependency hell"
Peer reviewed: get your code refereed
Web based: notebook interface

Questions?

Creating an Open Source Alternative to Magma, Maple, Mathematica, and MATLAB