Sage Days 14 Talk
system:sage
Sage
Creating a viable free open source alternative to
Magma, Maple, Mathematica and Matlab
William Stein, Associate Professor, University of Washington
{{{id=215|
///
}}}
{{{id=2|
///
}}}
History
- I started Sage at Harvard in January 2005.
- No existing math software good enough.
- I got very annoyed that my students and colleagues had to pay a ridiculous amount to use the code I wrote in Magma for Arithmetic Geometry (modular forms, elliptic curves, etc.).
- Sage-1.0 released February 2006 at Sage Days 1 (San Diego).
- Sage Days Workshops 1, 2, ..., 13, at UCLA, UW, Cambridge, Bristol, Austin, France, San Diego, Seattle, Georgia, etc.
- Sage won first prize in Trophees du Libre (November 2007)
- Funding from Microsoft, UW, NSF, DoD, Google, Sun, private donations, etc.
{{{id=117|
///
}}}
{{{id=214|
///
}}}
sagemath.org
Each day about 1500-2000 different real human visitors each day.
Sage Devel Headquarters: Four 24-core Sun X4450's with 128GB RAM each + 1 Sun X4540 with 24TB disk
sagenb.org
{{{id=62|
///
}}}
Some Advantages of Sage...
- Python-- a general purpose language at core of Sage; huge user base compared to Matlab, Mathematica, Maple and Magma
- Cython -- write fast compiled code in Sage
- Free and Open source
- Peer review of Code: "I do really and truly believe that the Sage refereeing model results in better code." -- John Cremona
{{{id=119|
///
}}}
{{{id=118|
///
}}}
Sage Is...
- Over 300,000 lines of new Python/Cython code
- Distribution of mathematical software; easy to build from source (over 5 million lines of code).
- About 150 developers: developer map
- Interfaces to most math software, including Magma, Macaulay 2, Singular, PHCpack, Polymake.
- Exact and numerical linear algebra, optimization, statistics (numpy, scipy, R, and gsl all included), group theory, combinatorics, cryptography.
- Number Theory: elliptic curves, modular forms, L-functions -- this is my research area
- 2d and 3d plotting
{{{id=96|
///
}}}
Examples: Calculus
Symbolic expressions:
{{{id=23|
x, y = var('x,y')
type(x)
///
}}}
{{{id=197|
reset('e')
///
}}}
{{{id=19|
a = 1 + sqrt(7)^2 + pi^e + 2/3 + x^y
///
}}}
{{{id=163|
show(a)
///
{x}^{y} + {\pi}^{e} + \frac{26}{3}
}}}
{{{id=20|
show(expand(a^2))
///
{x}^{{2 y}} + {{2 {\pi}^{e} } {x}^{y} } + \frac{{52 {x}^{y} }}{3} + {\pi}^{{2 e}} + \frac{{52 {\pi}^{e} }}{3} + \frac{676}{9}
}}}
Solve equations
{{{id=101|
var('a,b,c,x')
show(solve(x^2 + sqrt(17)*a*x + b == 0, x))
///
\begin{array}{l}[x = \frac{-\left( \sqrt{ {17 {a}^{2} } - {4 b} } \right) - {\sqrt{ 17 } a}}{2},\\
x = \frac{\sqrt{ {17 {a}^{2} } - {4 b} } - {\sqrt{ 17 } a}}{2}]\end{array}
}}}
{{{id=181|
var('a,b,c,x')
show(solve(a*x^3 + b*x + c == 0, x)[0])
///
x = {\left( \frac{{-\sqrt{ 3 } i}}{2} - \frac{1}{2} \right) {\left( \frac{\sqrt{ \frac{{{27 a} {c}^{2} } + {4 {b}^{3} }}{a} }}{{{6 \sqrt{ 3 }} a}} - \frac{c}{{2 a}} \right)}^{\frac{1}{3}} } - \frac{{\left( \frac{{\sqrt{ 3 } i}}{2} - \frac{1}{2} \right) b}}{{{3 a} {\left( \frac{\sqrt{ \frac{{{27 a} {c}^{2} } + {4 {b}^{3} }}{a} }}{{{6 \sqrt{ 3 }} a}} - \frac{c}{{2 a}} \right)}^{\frac{1}{3}} }}
}}}
Example: Some Rings
{{{id=207|
QQ[sqrt(2)]
///
Number Field in sqrt2 with defining polynomial x^2 - 2
}}}
{{{id=201|
RDF
///
Real Double Field
}}}
{{{id=202|
R = RealIntervalField(100); R
///
Real Interval Field with 100 bits of precision
}}}
{{{id=203|
R((1.01030,1.010332))
///
1.0103?
}}}
{{{id=164|
A = random_matrix(QQ, 500); v = random_matrix(QQ,500,1)
time x = A \ v
///
Time: CPU 1.56 s, Wall: 1.54 s
}}}
{{{id=3|
len(str(x[0]))
///
1484
}}}
{{{id=32|
///
}}}
Example: A Huge Integer Determinant
{{{id=210|
///
}}}
{{{id=31|
a = random_matrix(ZZ,200,200,x=-2^127,y=2^127)
time d = a.determinant()
len(str(d))
///
Time: CPU 3.78 s, Wall: 3.85 s
7786
}}}
We can also copy this matrix over to Maple and compute the same determinant there...
{{{id=167|
a[0,0]
///
-91532941219663566551217034654788002094
}}}
{{{id=27|
maple.with_package('LinearAlgebra')
B = maple(a)
t = maple.cputime()
time c = B.Determinant()
maple.cputime(t)
///
Time: CPU 0.00 s, Wall: 23.18 s
22.613
}}}
{{{id=168|
c == d
///
True
}}}
This ability to easily move objects between math software is unique to Sage.
{{{id=82|
B = magma(a)
t = magma.cputime()
time e = B.Determinant()
magma.cputime(t)
///
Time: CPU 0.00 s, Wall: 10.74 s
10.51
}}}
{{{id=208|
e == c
///
True
}}}
{{{id=112|
bernoulli
///
}}}
{{{id=25|
///
}}}
Example: A Symbolic Expression
{{{id=50|
x = var('x')
f(x) = sin(3*x)*x+log(x) + 1/(x+1)^2
show(f)
///
x \ {\mapsto}\ {x \sin \left( {3 x} \right)} + \log \left( x \right) + \frac{1}{{\left( x + 1 \right)}^{2} }
}}}
Plotting functions has similar syntax to Mathematica:
{{{id=171|
show(f.integrate(x))
///
x \ {\mapsto}\ \frac{\sin \left( {3 x} \right) - {{3 x} \cos \left( {3 x} \right)}}{9} + {x \log \left( x \right)} - \frac{1}{x + 1} - x
}}}
{{{id=109|
plot(f,(0.01,2), thickness=4) + text("Mathematica-style plotting in Sage", (1,-2), rgbcolor='black')
///
}}}
NOTE: Sage also has 2d plotting that is almost identical to MATLAB's 2d plotting:
{{{id=159|
import pylab as p
p.figure()
t = p.arange(0.01, 2.0, 0.01)
s = p.sin(2 * p.pi * t)
s = p.array([float(f(x)) for x in t])
P = p.plot(t, s, linewidth=4)
p.xlabel('time (s)'); p.ylabel('voltage (mV)')
p.title('Matlab-style plotting in Sage')
p.grid(True)
p.savefig('sage.png')
///
}}}
{{{id=110|
///
}}}
{{{id=73|
///
}}}
Example: Fast Float Evaluation
{{{id=172|
f(x,y,z) = sin(3*x)*x + log(x*y*z)^3 + 1/(1+y)^2
///
}}}
{{{id=74|
g = f._fast_float_(x,y,z)
timeit('g(4.5r,3.2r,5.7r)')
///
625 loops, best of 3: 669 ns per loop
}}}
{{{id=75|
%python
import math
def g(x,y,z): return math.sin(3*x)*x + log(x*y*z)**3 + 1/(1+y)**2
///
}}}
{{{id=77|
timeit('g(4.5r,3.2r,5.7r)')
///
625 loops, best of 3: 7.15 µs per loop
}}}
{{{id=213|
7.15/.669
///
10.6875934230194
}}}
{{{id=48|
///
}}}
{{{id=76|
///
}}}
Example: Compare Answers with Maple
{{{id=83|
var('x')
f = sin(3*x)*x+log(x) + 1/(x+1)^2
show(integrate(f))
///
\frac{\sin \left( {3 x} \right) - {{3 x} \cos \left( {3 x} \right)}}{9} + {x \log \left( x \right)} - \frac{1}{x + 1} - x
}}}
The command maple(...) fires up Maple (if you have it!), and creates a reference to a live object:
{{{id=0|
m = maple(f); m
///
sin(3*x)*x+ln(x)+1/(x+1)^2
}}}
{{{id=89|
type(m)
///
}}}
{{{id=90|
m.parent()
///
Maple
}}}
{{{id=91|
m.parent().pid()
///
67798
}}}
{{{id=92|
os.system('ps ax |grep 24038')
///
67818 s001 S+ 0:00.00 sh -c ps ax |grep 24038
67820 s001 R+ 0:00.00 grep 24038
0
}}}
Use Maple objects via a Pythonic notation:
{{{id=88|
show(m.integrate('x'))
///
1/9\,\sin \left( 3\,x \right) -1/3\,\cos \left( 3\,x \right) x+\ln
\left( x \right) x-x- \left( x+1 \right) ^{-1}
}}}
{{{id=42|
mathematica(f).Integrate(x)
///
-x - (1 + x)^(-1) - (x*Cos[3*x])/3 + x*Log[x] + Sin[3*x]/9
}}}
{{{id=43|
///
}}}
Example: Interactive Image Compression
This illustrates pylab (matplotlib + numpy), Sage plotting, html output, and @interact.
{{{id=97|
# first just play
import pylab
A = pylab.imread(DATA + 'simons.png')
graphics_array([matrix_plot(A), matrix_plot(1-A[0:,0:,1])]).show(figsize=[10,4])
///
}}}
{{{id=174|
A[0,0,]
///
array([ 0.03529412, 0.04705882, 0.01960784, 1. ], dtype=float32)
}}}
{{{id=44|
import pylab
A_image = pylab.mean(pylab.imread(DATA + 'simons.png'), 2)
@interact
def svd_image(i=(20,(1..100)), display_axes=True):
u,s,v = pylab.linalg.svd(A_image)
A = sum(s[j]*pylab.outer(u[0:,j], v[j,0:]) for j in range(i))
g = graphics_array([matrix_plot(A),matrix_plot(A_image)])
show(g, axes=display_axes, figsize=(8,3))
html('Compressed using %s eigenvalues
'%i)
///
}}}
{{{id=94|
///
}}}
{{{id=220|
///
}}}
Example: Implicit 2d Plot
{{{id=223|
var("x y")
implicit_plot(y^3 + x*y - x^2 - x, (-2,3), (-3,3), plot_points=400, cmap='hsv').show(aspect_ratio=1)
///
}}}
{{{id=55|
///
}}}
Example: 3d Plots
{{{id=84|
var('x,y')
P = plot3d(sqrt(x^2-y^2+3), (x,-3,3), (y,-2,2), color='blue', opacity=0.7)
P + sphere((0,0,0),2,color='red', aspect_ratio=[1,1,1], opacity=0.6)
///
}}}
{{{id=124|
( sphere((0,0,0), opacity=0.7) + sphere((0,1,0), color='red', opacity=0.5)
+ icosahedron((1,1,0), color='green') )
///
}}}
{{{id=224|
var('x,y,z')
T = RDF(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.77
implicit_plot3d(p, (-r, r), (-r, r), (-r, r), plot_points=40) # not released yet!!
///
}}}
{{{id=54|
var('x,y')
P = plot3d(sqrt(x^2-y^2+3), (x,-3,3), (y,-2,2), color='blue', opacity=0.7)
P + sphere((0,0,0),2,color='red', aspect_ratio=[1,1,1], opacity=0.6)
///
}}}
3d plotting (using jmol) is fast even though it does not use Java3d or OpenGL or require any special signed code or drivers.
{{{id=53|
# Yoda! -- over 50,000 triangles.
from scipy import io
X = io.loadmat(DATA + 'yodapose.mat')
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='green') + IndexFaceSet(F4,V,color='green')
Y = Y.rotateX(-1)
Y.show(aspect_ratio=[1,1,1], frame=False, figsize=4)
html('"Use the source, Luke..."')
///
"Use the source, Luke..."
}}}
{{{id=145|
///
}}}
Easily Write Very Fast Code (using Cython)
to sage-support
date Sat, Jan 31, 2009 at 11:15 AM
Hi,
I received first a MemoryError, and later on Sage reported:
uitkomst1=[]
uitkomst2=[]
eind=int((10^9+2)/(2*sqrt(3)))
print eind
for y in srange(1,eind):
test1=is_square(3*y^2+1,True)
test2=is_square(48*y^2+1,True)
if test1[0] and test1[1]%3==2: uitkomst1.append((y,(2*test1[1]-1)/3))
if test2[0] and test2[1]%3==1: uitkomst2.append((y,(2*test2[1]+1)/3))
print uitkomst1
een=sum([3*x-1 for (y,x) in uitkomst1 if 3*x-1<10^9])
print uitkomst2
twee=sum([3*x+1 for (y,x) in uitkomst2 if 3*x+1<10^9])
print een+twee
If you replace 10^9 with 10^6, the above listing works properly.
Maybe I made a mistake?
Rolandb
I rewrote Roland's code slightly so it wouldn't waste memory by constructing big lists... but the result was slow.
{{{id=150|
def f_python(n):
uitkomst1=[]
uitkomst2=[]
eind=int((n+2)/(2*sqrt(3)))
print eind
for y in (1..eind):
test1=is_square(3*y^2+1,True)
test2=is_square(48*y^2+1,True)
if test1[0] and test1[1]%3==2:
uitkomst1.append((y,(2*test1[1]-1)/3))
if test2[0] and test2[1]%3==1:
uitkomst2.append((y,(2*test2[1]+1)/3))
print uitkomst1
een=sum(3*x-1 for (y,x) in uitkomst1 if 3*x-1<10^9)
print uitkomst2
twee=sum(3*x+1 for (y,x) in uitkomst2 if 3*x+1<10^9)
print een+twee
///
}}}
{{{id=182|
time f_python(10^5)
///
28868
[(1, 1), (15, 17), (209, 241), (2911, 3361)]
[(1, 5), (14, 65), (195, 901), (2716, 12545)]
51408
Time: CPU 0.72 s, Wall: 0.77 s
}}}
{{{id=8|
time f_python(10^6)
///
288675
[(1, 1), (15, 17), (209, 241), (2911, 3361), (40545, 46817)]
[(1, 5), (14, 65), (195, 901), (2716, 12545), (37829, 174725)]
716034
Time: CPU 7.14 s, Wall: 7.65 s
}}}
Rewrite using CYTHON
{{{id=1|
%cython
from sage.all import is_square
cdef extern from "math.h":
long double sqrtl(long double)
def f(n):
uitkomst1=[]
uitkomst2=[]
cdef long long eind=int((n+2)/(2*sqrt(3)))
cdef long long y, t
print eind
for y in range(1,eind):
t = sqrtl( (3*y*y + 1))
if t * t == 3*y*y + 1:
uitkomst1.append((y, (2*t-1)/3))
t = sqrtl( (48*y*y + 1))
if t * t == 48*y*y + 1:
uitkomst2.append((y, (2*t+1)/3))
print uitkomst1
een=sum([3*x-1 for (y,x) in uitkomst1 if 3*x-1<10^9])
print uitkomst2
twee=sum([3*x+1 for (y,x) in uitkomst2 if 3*x+1<10^9])
print een+twee
///
}}}
{{{id=45|
time f(10^5)
///
28868
[(1L, 1L), (4L, 4L), (15L, 17L), (56L, 64L), (209L, 241L), (780L, 900L), (2911L, 3361L), (10864L, 12544L)]
[(1L, 5L), (14L, 65L), (195L, 901L), (2716L, 12545L)]
2
Time: CPU 0.00 s, Wall: 0.00 s
}}}
{{{id=4|
time f(10^6)
///
288675
[(1L, 1L), (4L, 4L), (15L, 17L), (56L, 64L), (209L, 241L), (780L, 900L), (2911L, 3361L), (10864L, 12544L), (40545L, 46817L), (151316L, 174724L)]
[(1L, 5L), (14L, 65L), (195L, 901L), (2716L, 12545L), (37829L, 174725L)]
2
Time: CPU 0.03 s, Wall: 0.03 s
}}}
{{{id=5|
time f(10^9)
///
288675135
[(1L, 1L), (4L, 4L), (15L, 17L), (56L, 64L), (209L, 241L), (780L, 900L), (2911L, 3361L), (10864L, 12544L), (40545L, 46817L), (151316L, 174724L), (564719L, 652081L), (2107560L, 2433600L), (7865521L, 9082321L), (29354524L, 33895684L), (109552575L, 126500417L)]
[(1L, 5L), (14L, 65L), (195L, 901L), (2716L, 12545L), (37829L, 174725L), (526890L, 2433601L), (7338631L, 33895685L), (102213944L, 472105985L)]
2
Time: CPU 25.60 s, Wall: 26.50 s
}}}
{{{id=6|
7.14/0.03
///
238.000000000000
}}}
A speedup by a factor of 238!
This is not a contrived example. This is a real world example that came up a this weekend. For C-style computations, Sage (via Cython) is as fast as C.
{{{id=151|
///
}}}
Example: Computing a Groebner Basis
{{{id=228|
P. = PolynomialRing(QQ,3, order='lex')
P
///
Multivariate Polynomial Ring in a, b, c over Rational Field
}}}
{{{id=227|
show(sage.rings.ideal.Katsura(P,3).groebner_basis())
///
\begin{array}{l}[a - 60 c^{3} + \frac{158}{7} c^{2} + \frac{8}{7} c - 1,\\
b + 30 c^{3} - \frac{79}{7} c^{2} + \frac{3}{7} c,\\
c^{4} - \frac{10}{21} c^{3} + \frac{1}{84} c^{2} + \frac{1}{84} c]\end{array}
}}}
{{{id=229|
show(sage.rings.ideal.Katsura(P,3).groebner_basis('singular:std'))
///
\begin{array}{l}[a - 60 c^{3} + \frac{158}{7} c^{2} + \frac{8}{7} c - 1,\\
b + 30 c^{3} - \frac{79}{7} c^{2} + \frac{3}{7} c,\\
c^{4} - \frac{10}{21} c^{3} + \frac{1}{84} c^{2} + \frac{1}{84} c]\end{array}
}}}
{{{id=230|
show(sage.rings.ideal.Katsura(P,3).groebner_basis('singular:slimgb'))
///
\begin{array}{l}[a - 60 c^{3} + \frac{158}{7} c^{2} + \frac{8}{7} c - 1,\\
b + 30 c^{3} - \frac{79}{7} c^{2} + \frac{3}{7} c,\\
c^{4} - \frac{10}{21} c^{3} + \frac{1}{84} c^{2} + \frac{1}{84} c]\end{array}
}}}
{{{id=231|
show(sage.rings.ideal.Katsura(P,3).groebner_basis('macaulay2'))
///
Traceback (most recent call last):
File "", line 1, in
File "/Users/wstein/talks/20090310-msri-sd14/nb/worksheets/admin/1/code/23.py", line 7, in
exec compile(ur'show(sage.rings.ideal.Katsura(P,_sage_const_3 ).groebner_basis(\u0027macaulay2\u0027))' + '\n', '', 'single')
File "/Users/wstein/build/sage-3.4.alpha0/local/lib/python2.5/site-packages/SQLAlchemy-0.4.6-py2.5.egg/", line 1, in
File "/Users/wstein/build/sage-3.4.alpha0/local/lib/python2.5/site-packages/sage/misc/cachefunc.py", line 242, in __call__
cache[key] = self.f(self._instance, *args, **kwds)
File "/Users/wstein/build/sage-3.4.alpha0/local/lib/python2.5/site-packages/sage/rings/polynomial/multi_polynomial_ideal.py", line 2458, in groebner_basis
gb = self._groebner_basis_macaulay2(*args, **kwds)
File "/Users/wstein/build/sage-3.4.alpha0/local/lib/python2.5/site-packages/sage/rings/polynomial/multi_polynomial_ideal.py", line 2058, in _groebner_basis_macaulay2
I = self._macaulay2_()
File "/Users/wstein/build/sage-3.4.alpha0/local/lib/python2.5/site-packages/sage/rings/polynomial/multi_polynomial_ideal.py", line 2020, in _macaulay2_
R._macaulay2_set_ring(macaulay2)
File "multi_polynomial_libsingular.pyx", line 880, in sage.rings.polynomial.multi_polynomial_libsingular.MPolynomialRing_libsingular._macaulay2_set_ring (sage/rings/polynomial/multi_polynomial_libsingular.cpp:8716)
File "/Users/wstein/build/sage-3.4.alpha0/local/lib/python2.5/site-packages/sage/interfaces/macaulay2.py", line 409, in ring
return self.new('%s[%s, MonomialSize=>16, MonomialOrder=>%s]'%(base_ring, varstr, order))
File "/Users/wstein/build/sage-3.4.alpha0/local/lib/python2.5/site-packages/sage/interfaces/expect.py", line 1067, in new
return self(code)
File "/Users/wstein/build/sage-3.4.alpha0/local/lib/python2.5/site-packages/sage/interfaces/expect.py", line 1002, in __call__
return cls(self, x, name=name)
File "/Users/wstein/build/sage-3.4.alpha0/local/lib/python2.5/site-packages/sage/interfaces/expect.py", line 1377, in __init__
raise TypeError, x
TypeError: Unable to start macaulay2
}}}
{{{id=232|
show(sage.rings.ideal.Katsura(P,3).groebner_basis('magma'))
///
\begin{array}{l}[a - 60 c^{3} + \frac{158}{7} c^{2} + \frac{8}{7} c - 1,\\
b + 30 c^{3} - \frac{79}{7} c^{2} + \frac{3}{7} c,\\
c^{4} - \frac{10}{21} c^{3} + \frac{1}{84} c^{2} + \frac{1}{84} c]\end{array}
}}}
{{{id=133|
I = sage.rings.ideal.Katsura(P,3)
///
}}}
{{{id=61|
F = I.groebner_fan()
F
///
Groebner fan of the ideal:
Ideal (a + 2*b + 2*c - 1, a^2 - a + 2*b^2 + 2*c^2, 2*a*b + 2*b*c - b) of Multivariate Polynomial Ring in a, b, c over Rational Field
}}}
{{{id=233|
F.weight_vectors()
///
[(4, 4, 1), (3, 2, 1), (5, 2, 3), (2, 1, 3), (2, 3, 5), (2, 3, 1), (2, 5, 3), (1, 4, 4)]
}}}
{{{id=108|
F.render()
///
}}}
{{{id=234|
///
}}} 
