DM559/DM545 – Linear and Integer Programming

Tutorial4Exam

Using IPython for interactive work

You can work interactively with Python using:

  • Spyder3
  • Jupyter (IPython)

See this tutorial Tutorial for an introduction to interactive Python. This is a syntetic tutorial only about the tools that can be useful for the digital tests during the course.

Preamble for working with Arrays and Fractions

Copy this preamble in your script to load the needed modules and implement a function printm and verbatim for pretty printing. (You can also put the code in a file, for example, utils.py and import that file from your script from utils import * or, if you are using Jupyter, using the magic function: %run utils.py.

In [1]:
from numpy import *
from fractions import Fraction as f
set_printoptions(precision=3,suppress=True)

def printm(a):
    """Prints the array as strings 
    :a: numpy array
    :returns: prints the array
    """
    def p(x):
        return str(x)
    p = vectorize(p,otypes=[str])
    print(p(a))
    
def tableau(a,W=7):
    """Returns a string for verbatim printing
    :a: numpy array
    :returns: a string
    """
    if len(a.shape) != 2:
        raise ValueError('verbatim displays two dimensions')
    rv = []
    rv+=[r'|'+'+'.join('{:-^{width}}'.format('',width=W) for i in range(a.shape[1]))+"+"]
    rv+=[r'|'+'|'.join(map(lambda i: '{0:>{width}}'.format("x"+str(i+1)+" ",width=W), range(a.shape[1]-2)) )+"|"+
         '{0:>{width}}'.format("-z ",width=W)+"|"
         '{0:>{width}}'.format("b ",width=W)+"|"]
    rv+=[r'|'+'+'.join('{:-^{width}}'.format('',width=W) for i in range(a.shape[1]))+"+"]
    for i in range(a.shape[0]-1):
        rv += [r'| '+' | '.join(['{0:>{width}}'.format(str(a[i,j]),width=W-2) for j in range(a.shape[1])])+" |"]
    rv+=[r'|'+'+'.join('{:-^{width}}'.format('',width=W) for i in range(a.shape[1]))+"+"]
    i = a.shape[0]-1
    rv += [r'| '+' | '.join(['{0:>{width}}'.format(str(a[i,j]),width=W-2) for j in range(a.shape[1])])+" |"]
    rv+=[r'|'+'+'.join('{:-^{width}}'.format('',width=W) for i in range(a.shape[1]))+"+"]
    print('\n'.join(rv))

Elementary row operations

To do elementary row operations on a matrix one can do as follows:

In [2]:
A = array( [[1, f(9,4), 2, f(4,5)],
             [f(3,4), 1, f(5,4), 0]])
A[0,:] = f(4,9) * A[0,:]
printm(A)

[['4/9' '1' '8/9' '16/45']
 ['3/4' '1' '5/4' '0']]
In [3]:
tableau(A)

|-------+-------+-------+-------+
|    x1 |    x2 |    -z |     b |
|-------+-------+-------+-------+
|   4/9 |     1 |   8/9 | 16/45 |
|-------+-------+-------+-------+
|   3/4 |     1 |   5/4 |     0 |
|-------+-------+-------+-------+

Note that indices start from 0, hence a vector of size $n$ has indices from $0$ to $n-1$.

Matrix Operations

Construct an array of zeros

In [4]:
zeros((3,1))
Out[4]:
array([[0.],
       [0.],
       [0.]])

Concatenate an identity matrix to a matrix

In [5]:
A_1 = concatenate([A,identity(2)],axis=1)
printm(A_1)

[['4/9' '1' '8/9' '16/45' '1.0' '0.0']
 ['3/4' '1' '5/4' '0' '0.0' '1.0']]

Insert a row in between other two:

In [6]:
A_2 = concatenate([ A[0:1,:], array( [[0,2,1,0]] ), A[1:2,:] ],axis=0)
printm(A_2)

[['4/9' '1' '8/9' '16/45']
 ['0' '2' '1' '0']
 ['3/4' '1' '5/4' '0']]

Ask for the dimensions of an array:

In [7]:
A.shape
Out[7]:
(2, 4)

Multiply matrices

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 [8]:
A = ones((4,2)) # a matrix of 1s of size 4x2
b = array([1,2])
A*b
Out[8]:
array([[1., 2.],
       [1., 2.],
       [1., 2.],
       [1., 2.]])

which is element-wise and hence not matrix-vector multiplication. The correct way to do matrix-vector and matrix-matrix multiplication is:

In [9]:
dot(A,b)
Out[9]:
array([3., 3., 3., 3.])

Inverse

In [10]:
A = random.randint(1,10,size=(4,4))
A_i = linalg.inv(A)
around( dot(A,A_i), decimals=2) # let's verify
Out[10]:
array([[ 1., -0.,  0.,  0.],
       [ 0.,  1.,  0.,  0.],
       [ 0., -0.,  1.,  0.],
       [-0., -0.,  0.,  1.]])

Transpose

In [11]:
A.T
Out[11]:
array([[5, 2, 7, 6],
       [6, 1, 8, 1],
       [6, 3, 7, 5],
       [7, 7, 3, 5]])

Rank

In [12]:
linalg.matrix_rank(A)
Out[12]:
4

Reduced row echelon form

You should calculate the reduced row echelon form by applying elementary row operations

In [13]:
A=array([[f(3,1),1,1,1],
         [2,4,1,1],
         [2,1,2,2]])
b=array([2,2,1])

AA=column_stack([A,b])
AA[0,:] = f(1,3) * AA[0,:]  # remember indices start from 0
AA[1,:] = -2 * AA[0,:] + AA[1,:]
AA[2,:] = -2 * AA[0,:] + AA[2,:]
printm(AA)
AA[1,:] = f(3,10)* AA[1,:]
AA[2,:] += -f(1,3) * AA[1,:]
printm(AA)
AA[2,:] = f(10,13) * AA[2,:]
printm(AA)
AA[1,:] += -f(1,10) * AA[2,:]
AA[0,:] += -f(1,3) * AA[2,:]
printm(AA)
AA[0,:] += -f(1,3) * AA[1,:]
printm(AA)

[['1' '1/3' '1/3' '1/3' '2/3']
 ['0' '10/3' '1/3' '1/3' '2/3']
 ['0' '1/3' '4/3' '4/3' '-1/3']]
[['1' '1/3' '1/3' '1/3' '2/3']
 ['0' '1' '1/10' '1/10' '1/5']
 ['0' '0' '13/10' '13/10' '-2/5']]
[['1' '1/3' '1/3' '1/3' '2/3']
 ['0' '1' '1/10' '1/10' '1/5']
 ['0' '0' '1' '1' '-4/13']]
[['1' '1/3' '0' '0' '10/13']
 ['0' '1' '0' '0' '3/13']
 ['0' '0' '1' '1' '-4/13']]
[['1' '0' '0' '0' '9/13']
 ['0' '1' '0' '0' '3/13']
 ['0' '0' '1' '1' '-4/13']]

You can check your result by importing the module for symbolic computation sympy. You are however discouraged to use this for your calculations:

In [14]:
import sympy as sy
# np.linalg.solve(A,b)

sy.Matrix(AA).rref()
Out[14]:
(Matrix([
 [1, 0, 0, 0,  9/13],
 [0, 1, 0, 0,  3/13],
 [0, 0, 1, 1, -4/13]]), (0, 1, 2))