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=219|
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=220|
///
}}}
Example: Implicit 2d Plot
{{{id=222|
var("x y")
implicit_plot(x^2+y^2-2, (-3,3), (-3,3)).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=116|
L = []
@interact
def random_list(number_of_points=(10..50), dots=True):
n = normalvariate
global L
if len(L) != number_of_points:
L = [(n(0,1), n(0,1), n(0,1)) for i in range(number_of_points)]
G = list_plot3d(L,interpolation_type='nn', texture=Color('#ff7500'),num_points=120)
if dots: G += point3d(L)
G.show()
///
}}}
{{{id=54|
///
}}}
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|
///
}}}
Real World Example: Super 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|
///
}}}
Numerical Matrix Algebra
{{{id=192|
a = random_matrix(RDF, 3); show(a)
///
\left(\begin{array}{rrr}
-0.898388350554 & -0.245114654439 & 0.630541996552 \\
0.654081247919 & 0.696119569249 & -0.957945963168 \\
0.00662634420862 & -0.704817828713 & -0.517608714138
\end{array}\right)
}}}
{{{id=193|
show(a.eigenvalues())
///
\left[0.961673429437,
-0.4808854176,
-1.20066550728\right]
}}}
{{{id=195|
a = random_matrix(RDF, 1000)
time v = a.eigenvalues()
///
Time: CPU 5.81 s, Wall: 6.21 s
}}}
{{{id=196|
show(points(v), axes=False, aspect_ratio=1, figsize=4)
///
}}}
{{{id=184|
a = random_matrix(RDF, 200) # read doubles in [-1,1]
G = graphics_array([matrix_plot(a^i) for i in [-3,-2,-1,1,2,3]], 2, 3)
show(G, axes=False)
///
}}}
{{{id=187|
@interact
def f(n=(10..400), field=[RDF,CDF]):
G = points(random_matrix(field,n).eigenvalues())
show(G, axes=False, frame=True)
///
}}}
{{{id=186|
///
}}}
{{{id=144|
///
}}}
{{{id=133|
///
}}}
{{{id=61|
///
}}}
{{{id=108|
///
}}} 
