DM559/DM545 – Linear and Integer Programming
Interactive Python¶In this course, we will use Python to carry out computation. In particular, it would be helpful to use Python interactively, trying things and learning by doing. Interactive Python is made available by the IDE Spyder and by Jupyter. Spyder is a Python Integrated Development Environment (IDE) that is installed in the Linux machines of the IMADA Computer Lab. Just type Jupyter provides a language-agnostic architecture for interactive computing. IPython is focused on Python providing a Python kernel for Jupyter. Jupyter notebook provides a python interface in the browser. The tutorial you are reading is a python notebook. To install Jupyter follow the instructions at the Jupyter page. For example, you can use Anacoda. Anaconda provides also a broad set of Python modules, among which Gurobi, that we will use in the second part of the course. Anaconda isn't the only choice in Python distributions and/or IDEs. Popular alternatives include Canopy, Eric, iep, and PyDev. However, if you are a more experienced user you can also use
Jupyter can be used from the shell by calling a console that allows the inlining of plots and text highlighting by typing: jupyter qtconsole
or as notebook from a borwser started by typing jupyter notebook
Notebook-style interfaces allow to mix executable code, text, and graphics and to create a self-documenting stream of results. Notebooks can be saved and continued later, which make them particularly well suited for prototyping and experimentation. You can create a new notebook by clicking on the You can read documentation on ipython at these links: Plotting¶Plotting in Python works with the Matplotlib plotting library. Matplotlib integrates with IPython's display system and event loop handling.
To make plots in IPython using Matplotlib, you must first enable IPython's matplotlib mode. To do this, run the In [1]:
%matplotlib inline
Functions with % are extra functions of IPython that add functionalities to the environment. They are called magic functions. Other useful magic function are %timeit to determine running time of a command and %run to run a script from a file. The command %matplotlib does not load the library in python. As usual we do: In [2]:
import matplotlib.pyplot as plt
To work with Linear Algebra we will use the In [3]:
import numpy as np
In IPython you can have access to the documentation of a function by typing '?' followed by the name of the module followed by the function: eg, Let's plot the following polynomial of degree 3: $$P_3(x)=x^3-7x+6=(x-1)(x-2)(x+3)$$ We use the a[0] * x**n + a[1] * x**(n-1) + ... + a[n-1]*x + a[n]
(Try In [4]:
a=[1,0,-7,6]
P=np.poly1d(a)
print(P)
In [5]:
x = np.linspace(-3.5, 3.5, 500)
plt.plot(x, P(x), '-')
plt.axhline(y=0)
plt.title('A polynomial of order 3');
Roots of a polynomial of degree n Let's $P_n(x)$ be a general polynomial of degree $n$: $$P_n(x)=p_n+p_{n-1}x+p_{n-2}x^2+\cdots+p_0x^n$$
In [6]:
np.roots(P)
Out[6]:
Arrays and Matrices¶We shall use In [7]:
b=np.array([1,2,3])
print(2*b)
print(b/2)
#np.norm(b)
Indexing and slices Indices start at 0. In [8]:
A=np.array([[1.,2,3],
[4,5,6],
[7,8,9]])
print("A: ", A)
print(A[0,0])
print(A[1:]) # all rows starting from the second
print(A[1]) # the second row of the matrix
print(A[:,2]) # a vector filled by the 3nd column of the matrix
print("Type: ", A.dtype)
print("Dimension:", A.shape, A.ndim)
In [9]:
A.reshape(1,9)
Out[9]:
Caution: reshape is just a view on the array, it does not copy the data. Beware: the following do shallow copy: In [10]:
N=A
N[0,0] = 0 # now A[0,0] = 1
M=A.reshape(1,9)
M[0,7]=10 # now A[2,1] is 10
A
Out[10]:
We can see the base object by: In [11]:
print(A.base)
print(M.base) # the base object of M is the matrix A
To implement a deep copy you must use: In [12]:
B=np.copy(A)
B
Out[12]:
Caution with casting: In [13]:
B=A.astype('int')
A[0,0]=1.5
B[0,1]=2.5
print(A)
print(B)
In Jupyter it is possible from command line to ask for completion via tab. This can also be used to explore which functions are available for a given module. Try for example to type
followed by a tab. You should see a list of available functions. Among
them there are two submodules that will be very useful for us:
Functions to construct arrays¶The following is a list of functions to construct matrices and vectors. In [14]:
np.zeros((3,4))
Out[14]:
In [15]:
np.ones((3,4))
np.identity(4)
np.diag([1,2,3,4])
np.random.rand(3,4)
np.random.randint(1,10,size=(3,4))
Out[15]:
In [16]:
np.arange(5)
Out[16]:
In [17]:
np.linspace(0,5,11)
Out[17]:
To concatenate arrays you have to use the concatenate function In [18]:
np.concatenate([A,A],axis=1) # same as np.column_stack() and np.hstack()
Out[18]:
In [19]:
np.concatenate([A,A],axis=0) # same as np.row_stack() and np.vstack()
Out[19]:
Note that
work also with 1-dimensional arrays. Operations with arrays¶The operations In [20]:
b=np.array([1,2,3])
v=np.array([2,0,1])
b*v
Out[20]:
In [21]:
np.dot(b,v)
Out[21]:
In [22]:
A*b
Out[22]:
In [23]:
np.dot(A,b)
Out[23]:
In [24]:
A.T # Matrix transpose
Out[24]:
Beware that like reshape, transpose does not copy the data. Vectors have one dimension hence transposing does not do anything particular. In this case it is more appropriate to use reshape: In [25]:
b.reshape(-1,1) # yelds a column vector, -1 tells to guess the size
Out[25]:
Most mathematical functions act element-wise. They are also called universal functions. In [26]:
np.cos(np.array([1,2,3]))
Out[26]:
In [27]:
2**np.arange(4) # arange includes the 0
Out[27]:
In [28]:
np.arange(4)**np.arange(4)
Out[28]:
Some functions do not act component-wise, eg, In [29]:
print(A)
np.sum(A) # on the whole matrix
np.sum(A,axis=1) # row-wise
np.sum(A,axis=0) # column-wise (this is also the result of `sum` from base python)
Out[29]:
If it is not, it is possible to make a function element-wise via In [30]:
def f(x):
return 0 if x<=5 else 1
# f(A) # error
np.vectorize(f)(A)
Out[30]:
Matrix inverse¶The inverse of a matrix can be calculated with a function from the set of linear algebra functions of numpy. In [31]:
np.linalg.det(A) # returns the determinant of A, $det(A) \neq 0$ if matrix invertible
np.linalg.inv(A)
Out[31]:
Plotting vectors¶To plot vectors we can use In [32]:
plt.quiver(0,0,3,4,scale=1,scale_units='xy',angles='xy')
plt.gca().set_xlim([0,5])
plt.gca().set_ylim([0,5]);
In [33]:
soa =np.array( [ [3,2], [1,1], [9,9]])
U,V = zip(*soa)
plt.quiver([0,0,0],[0,0,0],U,V,angles='xy',scale_units='xy',scale=1)
ax=plt.gca()
ax.set_xlim([-1,10])
ax.set_ylim([-1,10]);
Solving a Linear System¶If $A$ is a matrix and $\mathbf{b}$ is a vector, the system of linear equations $$A\mathbf{x}=\mathbf{b}$$ is solved using the In [34]:
import scipy.linalg as sl
x=sl.solve(A, b)
x
Out[34]:
In [35]:
np.allclose(np.dot(A, x), b)
Out[35]:
Easy printing for LaTex¶In [36]:
def bmatrix(a):
"""Returns a LaTeX bmatrix
:a: numpy array
:returns: LaTeX bmatrix as a string
"""
if len(a.shape) > 2:
raise ValueError('bmatrix can at most display two dimensions')
lines = str(a).replace('[', '').replace(']', '').splitlines()
rv = [r'\begin{bmatrix}']
rv += [r' ' + ' & '.join(l.split()) + '\\\\' for l in lines]
rv += [r'\end{bmatrix}']
return '\n'.join(rv)
In [37]:
A=np.random.randint(4,size=(2,2))
print(bmatrix(A))
print(bmatrix(A.T))
Fraction Mode¶It is possible to carry out operations in fraction mode using the module In [38]:
from fractions import Fraction as f
A = np.array([[f(1,1),f(2,1)],
[f(3,1),f(4,1)]])
A[0] = f(1,3)*A[0]
Some pretty printing function, which translates fractions in strings: In [39]:
def printm(a):
"""Prints the array as strings
:a: numpy array
:returns: prints the array
"""
def p(x):
return str(x)
p = np.vectorize(p,otypes=[str])
print(p(a))
In [40]:
printm(A)
|