#!/usr/bin/env python
# coding: utf-8

# # Programming with Python
# #### Stefan Güttel, [guettel.com](http://guettel.com)
# 
# ## Lab classes: Modules

# **Problem 1.** Write a module `symtoep` for creating various [symmetric](http://en.wikipedia.org/wiki/Symmetric_matrix) [Toeplitz matrices](http://en.wikipedia.org/wiki/Toeplitz_matrix) with the following functions:
# 
# * `general(L)` that returns a symmetric Toeplitz matrix with the first row equal to the elements of the list `L`.
# 
# * `tridiagonal(n, d, sd)` that returns a [tridiagonal](http://en.wikipedia.org/wiki/Tridiagonal_matrix) symmetric Toeplitz matrix of order `n`, with the float number `d` on the diagonal and the float number `sd` right below and above it (the rest of the elements are zero).
# 
# * `filled(n, d, nd)` that returns a symmetric Toeplitz matrix of order `n`, with the float number `d` on the diagonal and all the non-diagonal elements equal to `nd`.
# 
# * `menu()` that displays a choice
#    ```
#    Load a symmetric Toeplitz matrix:
#      1. General symmetric Toeplitz matrix
#      2. Tridiagonal symmetric Toeplitz matrix
#      3. Filled symmetric Toeplitz matrix
#    Your choice (1-3):
#    ```
#    and asks the user to choose. Then,
#    * if the user chooses `"1"`, the function asks for a list `L` of numbers (preferably as a string of comma-separated floats, but you can use some other method as well), and then returns `general(L)`,
#    * if the user chooses `"2"`, the function asks for an integer `n` and floats `d` and `sd`, and returns `tridiagonal(n, d, sd)`.
#    * if the user chooses `"3"`, the function asks for an integer `n` and floats `d` and `nd`, and returns `filled(n, d, sd)`.
#    * if the user chooses anything else, the choice is redisplayed and user is asked to choose again.
# 
# When run like a program, this module should do nothing.
# 
# Toeplitz matrices are those that have constant elements on all diagonals.
# * An example of a symmetric Toeplitz matrix of order $3$:
# $$\begin{bmatrix} 1 & 2 & 3 & 4 & 5 \\ 2 & 1 & 2 & 3 & 4 \\ 3 & 2 & 1 & 2 & 3 \\ 4 & 3 & 2 & 1 & 2 \\ 5 & 4 & 3 & 2 & 1 \end{bmatrix}$$
# * An example of a tridiagonal (all elements other than those on the main and its neighbour diagonals are zero) symmetric Toeplitz matrix of order $3$:
# $$\begin{bmatrix} 1 & 2 & 0 & 0 & 0 \\ 2 & 1 & 2 & 0 & 0 \\ 0 & 2 & 1 & 2 & 0 \\ 0 & 0 & 2 & 1 & 2 \\ 0 & 0 & 0 & 2 & 1 \end{bmatrix}$$
# * An example of a "filled" symmetric Toeplitz matrix of order $3$:
# $$\begin{bmatrix} 1 & 2 & 2 & 2 & 2 \\ 2 & 1 & 2 & 2 & 2 \\ 2 & 2 & 1 & 2 & 2 \\ 2 & 2 & 2 & 1 & 2 \\ 2 & 2 & 2 & 2 & 1 \end{bmatrix}$$
# 
# All the functions have to return the matrices as `NumPy`'s `matrix` type, which is easily created from a list or any other two-dimensional array-like type, using the function [`numpy.matrix`](http://docs.scipy.org/doc/numpy/reference/generated/numpy.matrix.html). For example, the first of the three matrices shown above can be created as follows:

# In[1]:


import numpy as np
L = list(range(1,6))
res = list()
for i in range(5):
    row = list()
    for j in range(5):
        row.append(L[abs(i-j)])
    res.append(row)
res = np.matrix(res)
print(res)


# It can be done more easily using a generator expression:

# In[2]:


import numpy as np
L = list(range(1,6))
res = np.matrix([[L[abs(i-j)] for i in range(5)] for j in range(5)])
print(res)


# You may implement the functions `tridiagonal(n, d, sd)` and `filled(n, d, nd)` either on their own or simply by creating `L` and calling `general(L)`. The former approach is prefered, in which case [`numpy.zeros`](http://docs.scipy.org/doc/numpy/reference/generated/numpy.zeros.html), [`numpy.ones`](http://docs.scipy.org/doc/numpy/reference/generated/numpy.ones.html), [`numpy.fill_diagonal`](http://docs.scipy.org/doc/numpy/reference/routines.indexing.html), and [`numpy.tri`](http://docs.scipy.org/doc/numpy/reference/generated/numpy.tri.html#numpy.tri) may help you, along with the matrix addition (done with the usual plus `+` operator) and multiplication with a constant (done with the usual multiplication `*` operator). These functions are much faster than setting the elements index by index, but you are free to work without them.

# **Problem 2.** Using the module `symtoep`, write a program that loads a matrix (one of the 3 supported types), and then prints it and its eigenvalues.
# 
# To compute the eigenvalues of a symmetric matrix, use [`scipy.linalg.eigvalsh`](http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.linalg.eigvalsh.html). If, for some reason, SciPy doesn't work, use [`numpy.linalg.eigvalsh`](http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.eigvalsh.html).
# 
# Ideally, make your program run properly with either (but this is not a must). For example,
# ```python
# try:
#     import scipy.linalg as la
# except ImportError:
#     import numpy.linalg as la
# ...
# print(la.eigvalsh(A))
# ```

# **Problem 3.** Using the module `symtoep`, write a program that loads a matrix (one of the 3 supported types), and then prints it and a message explaining if it is positive semidefinite or not.
# 
# To test positive semidefiniteness, try to compute the [Cholesky factorization](http://en.wikipedia.org/wiki/Cholesky_decomposition), using either [`scipy.linalg.cholesky`](http://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.cholesky.html) or [`numpy.linalg.cholesky`](http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.cholesky.html). If the computation fails (i.e., a `numpy.linalg.linalg.LinAlgError` exception is raised), the matrix is not positive semidefinite. Otherwise, it is.
# 
# The same remarks from Problem 2 regarding the import of SciPy/NumPy modules apply.