DM559/DM545 – Linear and Integer Programming

Jupyter and IPython

Jupyter provides a language-agnostic architecture for interactive computing. IPython is focused on Python providing a Python kernel for Jupyter.

To install Jupyter follow the instructions at the Jupyter page. For example, you can use Anacoda, a widely-used Python platform that includes an Integrated Development Environment (Spyder), a notebook-style interface (Jupyter), and 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 pip. In the course we will use Python 2.7 because Gurobi does not support Python 3 yet. However, the behaviour of Python 3 can be approximated in Python by importing the module __future__. For example:

from __future__ import division, print_function, unicode_literals

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 New icon (in the upper right) and choosing one of the Notebook options. Once your new notebook starts, you can type standard Python commands or Gurobi Interactive Shell commands directly into the In window.

Alternatively, you can use the Spyder Python IDE. Simply type spyder in a terminal window.

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 (there is also a feature to load Matplotlib examples %load URL_to_py eg: %load http://matplotlib.org/mpl_examples/style_sheets/plot_ggplot.py)

To make plots using Matplotlib, you must first enable IPython's matplotlib mode. To do this, run the %matplotlib command to enable plotting in the current Notebook. This takes an optional argument that specifies which Matplotlib backend should be used. Most of the time, in the Notebook, you will want to use the inline backend, which will embed plots inside the Notebook:

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
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, ?np.polyd

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 numpy.poly1d function to represent the polynomial. It takes an array of coefficients a of length $n+1$ such that the polynomial is described by:

a[0] * x**n + a[1] * x**(n-1) + ... + a[n-1]*x + a[n]

(Try numpy.polyfit? to know more).

In [3]:
a=[1,0,-7,6]
P=np.poly1d(a)
print P

   3
1 x - 7 x + 6
In [4]:
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$$

numpy.roots(p) Return the roots of a polynomial with coefficients given in P. For the polynomial of 3rd degree seen above:

In [5]:
np.roots(P)
Out[5]:
array([-3.,  2.,  1.])

Arrays and Matrices

In python we shall use np.arrays to represent matrices and vectors. Althoguh very similar to lists they differ in that they store only the same numeric type and there is no append method.

In [6]:
b=np.array([1,2,3])
print 2*b
print b/2
#np.norm(b)

[2 4 6]
[0 1 1]

Indexing and slices

Indices start at 0.

In [7]:
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

A =  [[ 1.  2.  3.]
 [ 4.  5.  6.]
 [ 7.  8.  9.]]
1.0
[[ 4.  5.  6.]
 [ 7.  8.  9.]]
[ 4.  5.  6.]
[ 3.  6.  9.]
Type:  float64
Dimension: (3, 3) 2
In [8]:
A.reshape(1,9)
Out[8]:
array([[ 1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9.]])

Caution: reshape is just a view on the array, it does not copy the data. Hence:

In [9]:
M=A.reshape(1,9)
M[0,7]=10 # now A[2,1] is 10
A
Out[9]:
array([[  1.,   2.,   3.],
       [  4.,   5.,   6.],
       [  7.,  10.,   9.]])
In [10]:
print A.base
print M.base

None
[[  1.   2.   3.]
 [  4.   5.   6.]
 [  7.  10.   9.]]

Caution with casting:

In [11]:
B=A.astype('int')
A[0,0]=1.5
B[0,1]=2.5
A, B
Out[11]:
(array([[  1.5,   2. ,   3. ],
        [  4. ,   5. ,   6. ],
        [  7. ,  10. ,   9. ]]), array([[ 1,  2,  3],
        [ 4,  5,  6],
        [ 7, 10,  9]]))

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

import numpy as np
np.

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: random and linalg. The first implements function to generate random numbers and matrices. The second implements functions from linear algebra.

Functions to construct arrays

The following is a list of functions to construct matrices and vectors.

In [12]:
np.zeros((3,4))
Out[12]:
array([[ 0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.]])
In [13]:
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[13]:
array([[3, 5, 4, 9],
       [8, 9, 9, 5],
       [8, 6, 8, 2]])
In [14]:
np.arange(5)
Out[14]:
array([0, 1, 2, 3, 4])
In [15]:
np.linspace(0,5,11)
Out[15]:
array([ 0. ,  0.5,  1. ,  1.5,  2. ,  2.5,  3. ,  3.5,  4. ,  4.5,  5. ])

To concatenate arrays you have to use the concatenate function

In [16]:
np.concatenate([A,A],axis=1) # same as np.column_stack() and np.hstack()
Out[16]:
array([[  1.5,   2. ,   3. ,   1.5,   2. ,   3. ],
       [  4. ,   5. ,   6. ,   4. ,   5. ,   6. ],
       [  7. ,  10. ,   9. ,   7. ,  10. ,   9. ]])
In [17]:
np.concatenate([A,A],axis=0) # same as np.row_stack() and np.vstack()
Out[17]:
array([[  1.5,   2. ,   3. ],
       [  4. ,   5. ,   6. ],
       [  7. ,  10. ,   9. ],
       [  1.5,   2. ,   3. ],
       [  4. ,   5. ,   6. ],
       [  7. ,  10. ,   9. ]])

Note that concatenate works with arrays that are at least 2-dimensional. The methods:

  • column_stack and hstack to stack 1-D arrays as columns into a 2-D array and
  • row_stack and vstack to stack arrays in sequence vertically (row wise)

work also with 1-dimensional arrays.

Operations with arrays

The operations +,*,/,- are elementwise. For the operations of scalar multiplication and matrix multiplication the dot method must be used. Compare the following two cells:

In [18]:
b=np.array([1,2,3])
v=np.array([2,0,1])
b*v
Out[18]:
array([2, 0, 3])
In [19]:
np.dot(b,v)
Out[19]:
5
In [20]:
A*b
Out[20]:
array([[  1.5,   4. ,   9. ],
       [  4. ,  10. ,  18. ],
       [  7. ,  20. ,  27. ]])
In [21]:
np.dot(A,b)
Out[21]:
array([ 14.5,  32. ,  54. ])
In [22]:
A.T # Matrix transpose
Out[22]:
array([[  1.5,   4. ,   7. ],
       [  2. ,   5. ,  10. ],
       [  3. ,   6. ,   9. ]])

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 [23]:
b.reshape(-1,1) # yelds a column vector, -1 tells to guess the size
Out[23]:
array([[1],
       [2],
       [3]])

Most mathematical functions act element-wise. They are also called universal functions.

In [24]:
np.cos(np.array([1,2,3]))
Out[24]:
array([ 0.54030231, -0.41614684, -0.9899925 ])
In [25]:
2**np.arange(4) # arange includes the 0
Out[25]:
array([1, 2, 4, 8])
In [26]:
np.arange(4)**np.arange(4)
Out[26]:
array([ 1,  1,  4, 27])

Some functions do not act component-wise, eg, min, max, sum. These functions may operate on the whole matrix or column-wise or row-wise.

In [27]:
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)

[[  1.5   2.    3. ]
 [  4.    5.    6. ]
 [  7.   10.    9. ]]
Out[27]:
array([ 12.5,  17. ,  18. ])

If it is not, it is possible to make a function element-wise via vectorize. This yields often an improvement in performance with respect to for loops. Compare for example:

In [28]:
def f(x):
    return 0 if x<=5 else 1
# f(A) # error
np.vectorize(f)(A)
Out[28]:
array([[0, 0, 0],
       [0, 0, 1],
       [1, 1, 1]])

Matrix inverse

The inverse of a matrix can be calculated with a function from the set of linear algebra functions of numpy.

In [29]:
np.linalg.det(A) # returns the determinant of A, $det(A) \neq$ if matrix invertible 
np.linalg.inv(A)
Out[29]:
array([[-3.33333333,  2.66666667, -0.66666667],
       [ 1.33333333, -1.66666667,  0.66666667],
       [ 1.11111111, -0.22222222, -0.11111111]])

Plotting vectors

To plot vectors we can use matplotlib.pyplot.quiver. In the basic form quiver(X,Y,U,V), X and Y are the coordinates of the tail of the vector and U and V the components of the vector. The plot needs some tweak for the scales.

In [30]:
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 [31]:
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 mtrix and $b$ is a vector, the system of linear equations $$Ax=b$$ is solved using the solve method from numpy.linalg or scipy.linalg (the second being more efficient):

In [32]:
import scipy.linalg as sl
x=sl.solve(A, b)
x
Out[32]:
array([ -1.11022302e-16,   7.77156117e-17,   3.33333333e-01])
In [33]:
np.allclose(np.dot(A, x), b)
Out[33]:
True

Easy printing for LaTex

In [34]:
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 [35]:
A=np.random.randint(4,size=(2,2))
print(bmatrix(A)) 
print(bmatrix(A.T))

\begin{bmatrix}
  1 & 2\\
  1 & 1\\
\end{bmatrix}
\begin{bmatrix}
  1 & 1\\
  2 & 1\\
\end{bmatrix}