4403.sagews

Tutorial on the scientific use of Python

My name is Fabio Grazioso, and here is my picture:

One of the key features of SAGE is the ability to build "notebooks", where computing code is shown along with its output, and with other "material" assembled by the author.

This is the reason why I start with my picture: to show how to insert styled text and "external" image files in a notebook.

To insert a "Text cell" shift-click the blue inserting line

To insert an image file (from local hard drive) in a "text cell":

1) upload the file on SAGE server using the "Data" menu

2) use the "edit/insert image" button in the toolbar of the "text cell"

(this can be interesting to import a plot from other software)

more details can be found here: http://wiki.sagemath.org/quickref?action=AttachFile&do=get&target=quickref.pdf

 

Basic Python functions

N(log(12))  # N() asks for the numerical value

\newcommand{\Bold}[1]{\mathbf{#1}}2.48490664978800

N(log(e)) #let's check the base

\newcommand{\Bold}[1]{\mathbf{#1}}1.00000000000000

print N(log(12) , 1000)  # have a lot of decimal figures

2.48490664978800031022970947983887884079849082654325995997605435262428153715799839808534240880656946387197299110709923720973904469744678093511814337403742933446461791084012423755998777609214430684598927054752490743162383286809603449691907047782504850125493615385291175134283980681385134645366622146078

to have the logarithm in base 2 we have to convert

N(log(12) / log(2) )

\newcommand{\Bold}[1]{\mathbf{#1}}3.58496250072116

loading SciPy

if we want something more sophisticated, we have to load "SciPy", the scientific library:

import scipy as sp

now we can use dirctly log2(), which is defined in SciPy and returns logarithm in base 2:

sp.log2(12)

\newcommand{\Bold}[1]{\mathbf{#1}}3.58496250072

in general, to find help on SciPy visit the official website:

www.scipy.org

I also use google to search within that: typing "log2 site:www.scipy.org" into gooooogle

Arrays in SciPy

Python has a "basic" object which is the list, denoted by square brackets:

myList = [1,2,3,4,5]

and it is possible to do nice things with lists:

mySecondList = [3+x for x in myList]
print(mySecondList)

[4, 5, 6, 7, 8]

but to do more sophisticated things we use scipy arrays:

myArray = sp.array([1,2,3,4,5])
print(myArray)

[1 2 3 4 5]

my2DArray = sp.array([[0,1,2,3,4,5,6,7,8,9,10],[0,10,20,30,40,50,60,70,80,90,100],[0,11,22,33,44,55,66,77,88,99,111]])
print(my2DArray)

[[ 0 1 2 3 4 5 6 7 8 9 10] [ 0 10 20 30 40 50 60 70 80 90 100] [ 0 11 22 33 44 55 66 77 88 99 111]]

print(sp.transpose(my2DArray))

[[ 0 0 0] [ 1 10 11] [ 2 20 22] [ 3 30 33] [ 4 40 44] [ 5 50 55] [ 6 60 66] [ 7 70 77] [ 8 80 88] [ 9 90 99] [ 10 100 111]]

mySlice = my2DArray[0:2, :]  # NB first included, last excluded!
print(mySlice)

[[ 0 1 2 3 4 5 6 7 8 9 10] [ 0 10 20 30 40 50 60 70 80 90 100]]

at this link there is a nice page with a summary about SciPy arrays:

http://pages.physics.cornell.edu/~myers/teaching/ComputationalMethods/python/arrays.html

control structures: for cycle syntax

for i in mySecondList:
    print i

4 5 6 7 8

for i in range(2,10):    # remember: first included, last excluded
    print i

2 3 4 5 6 7 8 9

for i in range(20, 0, -2):
    print i

20 18 16 14 12 10 8 6 4 2

Linear algebra, solving system of linear equations

to do some linear algebra, another library is needed, so we load it:

coefficients = sp.array([[2,3,6], [3,4,0], [2,4,6]])
knowns = sp.array([-1,7,0])
print(coefficients)
print(knowns)

[[2 3 6] [3 4 0] [2 4 6]] [-1 7 0]

from scipy import linalg #for some reason linalg has to be imported explicitly
x = sp.linalg.solve(coefficients, knowns)
print(x)

[ 1. 1. -1.]

coefficients = sp.matrix(coefficients)
x = sp.matrix(x)
coefficients * sp.transpose(x)

\newcommand{\Bold}[1]{\mathbf{#1}}$$\begin{array}{l} \verb|[[-1.]|\\ \phantom{x}\verb|[|\phantom{x}\verb|7.]|\\ \phantom{x}\verb|[|\phantom{x}\verb|0.]]| \end{array}$$

here we are solving the system      «coefficients * x = knowns»

plotting

Although SAGE has its own tools for plotting,  advanced scientific plotting in python is done using the library "matplotlib" (see http://matplotlib.sourceforge.net/)

Here we show how to use matplotlib in SAGE.

import scipy as sp
import matplotlib.pyplot as plt

a = sp.arange(0,3,.02)    
b = sp.arange(0,3,.02)
c = sp.exp(a)
d = c[::-1]

fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(a,c,'k--',a,d,'k:',a,c+d,'k')    #this is the actual plotting command
leg = ax.legend(('Model length', 'Data length', 'Total message length'),
           'upper center', shadow=True)
ax.set_ylim([-1,20])
ax.grid(False)
ax.set_xlabel('Model complexity --->')
ax.set_ylabel('Message length --->')
ax.set_title('Minimum Message Length')

ax.set_yticklabels([])
ax.set_xticklabels([])

# set some legend properties.  All the code below is optional.  The
# defaults are usually sensible but if you need more control, this
# shows you how

# the matplotlib.patches.Rectangle instance surrounding the legend
frame  = leg.get_frame()
frame.set_facecolor('0.80')    # set the frame face color to light gray

# matplotlib.text.Text instances
for t in leg.get_texts():
    t.set_fontsize('small')    # the legend text fontsize

# matplotlib.lines.Line2D instances
for l in leg.get_lines():
    l.set_linewidth(1.5)  # the legend line width
plt.savefig('fig')        #here is the only difference with the use of matplotlib in python: 
                          #in python this last line has to be replaced by "plt.show()" (with no quotes)

more examples:

import scipy as sp
import matplotlib.pyplot as plt
from matplotlib.font_manager import FontProperties

Xrange=sp.linspace(0, 5, 501)
# 
def f1(x): return sp.log2(x)
def f2(x): return sp.sin(x)

 # legend font size
fontP = FontProperties() #legend font
fontP.set_size('small') #legend font
  
fig=plt.figure()
 
# # shrink the "plot box" to make room for the legend
plotbox = plt.subplot(111)
box = plotbox.get_position()
plotbox.set_position([box.x0, box.y0, box.width, box.height  * 0.8])
 
var=3
plt.plot(Xrange, f1(Xrange), linewidth=2.0, label='$f(x)$')  #
plt.plot(Xrange, f2(Xrange), label='n'+ str(var))        #
 
plt.axes(plotbox).tick_params(labelsize=16) # this is to change the font size of the axes tick labels

plt.legend(loc='best', ncol=2, labelspacing=0.01, columnspacing=2, prop=fontP, bbox_to_anchor=(1, 1.2)) #

plt.savefig('fig')

Warning: divide by zero encountered in log2

Interactive plotting

 

here we have the very useful "interact" feature, where we change a parameter using a slider

import scipy as sp
import matplotlib.pyplot as plt


const12 = 0.5 * sp.ones(501)    # definition of a constant in the plot

@interact
def _(xi=(0,1), xMax=(1,20), yMax=(1,20)):      #here the "sliders" parameters are defined. xi=(1..2) is a discrete parameter
    Xrange=sp.linspace(0, xMax, 501)    # definition of the abscissa range and resolution
    def eZ(d): return d/(0.25 + 2*d)
    def eZQ(d,xi): return d/(0.25 +(1-xi)*d)
    fig=plt.figure()
    ax = fig.add_subplot(111)
    ax.set_ylim([0, yMax])            #here limits in y are defined
    plt.plot(Xrange, eZ(Xrange))
    plt.plot(Xrange, eZQ(Xrange, xi))
    plt.plot(Xrange,Xrange, color='red')
    plt.plot(Xrange,const12, color='yellow')
    plt.savefig('fig')

xi  xMax  yMax  [removed] [removed] [removed] [removed]

3D plotting

from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
import scipy as sp

@interact
def _(pitch=(0, 180), yaw=(0, 360), xMin=(-10,0), xMax=(0, 10), xMeshStep=(0.1, 0.5), 
                        yMin=(-10,0), yMax=(0, 10), yMeshStep=(0.1, 0.5)):
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')

    X, Y = sp.mgrid[xMin:xMax:xMeshStep, yMin:yMax:yMeshStep]		# this creates the "mesh" 
						# (two 2D arrays with the (x,y) values)

    Z = sp.cos(X) + sp.sin(Y)		#this is the function


    ax.plot_wireframe(X, Y, Z, rstride=5, cstride=5)    # cstride and rstride are the size of the "wireframe"
						    # cell in "mesh units"

    ax.view_init(pitch, yaw)				# those are the "pitch" and "yaw" wiew angles
    plt.savefig('fig')

numerical integration

from scipy import integrate   #again, first import the object explicitly from
def f(x, a, b):
    return a * x**2 + b

args = (2, 1.5)
    
result = sp.integrate.quad(f, 0, 5, args)

print 'integral=',result[0],'    with error', result[1]

integral= 90.8333333333 with error 1.0084525807e-12

minimization: finding local minima for functions of multiple variables

here is the documentation for the minimization function:

http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fmin_slsqp.html

import scipy.optimize
def rosen(x):
	"""The Rosenbrock function"""
	return sum(100.0*(x[1:]-x[:-1]**2.0)**2.0 + (1-x[:-1])**2.0)
	
def constrains(x):
	return sp.sum(x) < 1000
	

x0 = sp.array([1.3, 0.7, 0.8, 1.9, 1.2, 1.1])

bounds_list = [(0,5)]* sp.alen(x0)

results = sp.optimize.fmin_slsqp(rosen, x0, iprint=2, bounds=bounds_list, iter=200, acc=1e-9, full_output=1)

NIT FC OBJFUN GNORM 1 8 8.598200E+02 2.217809E+03 2 17 2.938209E+02 1.065307E+03 3 27 2.104964E+02 8.078765E+02 4 36 1.312137E+02 6.556575E+02 5 45 2.387009E+01 2.815240E+02 6 54 1.373233E+01 2.355823E+02 7 64 1.331185E+01 2.307380E+02 8 73 9.042670E+00 1.739320E+02 9 82 6.330557E+00 1.358729E+02 10 91 1.160474E+00 3.343308E+01 11 99 6.411282E-01 1.102618E+01 12 107 5.573775E-01 9.066741E+00 13 115 4.907074E-01 9.020361E+00 14 123 4.179102E-01 2.413114E+01 15 131 3.175876E-01 7.551915E+00 16 139 2.491147E-01 5.125283E+00 17 148 2.084760E-01 7.818404E+00 18 156 1.720741E-01 1.237367E+01 19 164 1.323767E-01 6.865034E+00 20 172 8.443046E-02 4.479579E+00 21 180 5.572290E-02 8.431220E+00 22 188 2.972724E-02 2.506456E+00 23 196 1.586080E-02 2.265520E+00 24 204 1.111980E-02 3.891739E+00 25 212 4.659902E-03 1.093651E+00 26 220 2.362658E-03 5.751417E-01 27 228 3.438655E-04 6.736684E-01 28 236 3.356882E-05 6.450326E-02 29 244 8.832935E-07 3.732706E-02 30 252 9.941508E-08 1.666719E-02 31 260 1.119800E-08 4.975348E-03 32 268 4.561139E-10 4.653546E-04 Optimization terminated successfully. (Exit mode 0) Current function value: 4.56113906466e-10 Iterations: 32 Function evaluations: 268 Gradient evaluations: 32