#
# Python introduction, Oct 19th-20th, 2015
#


# write a first script
print('\nThis is my first Python script.\n')


# write an interactive script
print('Use Python to sum numbers')
a = input('a = ')
b = input('b = ')
c = a+b
print('c = a + b = ' + str(c))


# booleans, logical expressions
a==b
a!=b
a<b
a<=b
a>b
a>=b


# if construct (conditional statement)
if a<0:
  print('a is negative') # the body of the if/elif and else statements are indented
elif a>0:                # the elif and else statements have the same indent like the if statement
  print('a is positive')
else:
  print('a is zero')


# while loop
i = 0
while i<=10:
  print('i = ' + str(i)) # the body of the loop is indented
  i+= 1                  # increment i


# for loop
for i in range(5):
  print('i = ' + str(i))


# sequences
print('range(5)   = ' + str(range(5)))
print('range(0,5) = ' + str(range(0,5)))
print('range(1,6) = ' + str(range(1,6)))


# import numpy module
import numpy as np


# data format of numpy arrays
A = np.zeros((2,2),int)
print(A)
print(type(A[0,0]))
A = np.zeros((2,2),float)
print(A)
print(type(A[0,0]))
A = np.zeros((2,2))
print(A)
print(type(A[0,0]))


# multiplication of numpy arrays 
A = np.array([[1,2],[3,4]])
B = np.identity(2)
A*B         # componentwise multiplication
np.dot(A,B) # normal matrix matrix multiplication


# matrix vector products
a = np.array([5,6])
B*a         # componentwise multiplication
np.dot(B,a) # normal matrix vector multiplication
np.dot(a,B) # also works (-> a does not have dimension 2x1 or 1x2, it has dimension 2)

b = np.array([[5,6]])
np.dot(B,b) # does not work, dimensions mismatch
np.dot(b,B) # works

c = np.array([[5],[6]]) # = b.transpose()


# solve linear systems
x = np.linalg.solve(A,c)     # uses LAPACK implementation of LU decomposition
np.allclose(np.dot(A, x), c)


# eigenvalues and -vectors
e = np.linalg.eig(A)
eigenvalues  = e[0]
eigenvectors = e[1]
print('eigenvalues of A:  ' + str(eigenvalues))
print('eigenvectors of A: ' + str(eigenvectors))


# sparse matrices
print('sparse matrices')
import scipy.sparse as sp
S = sp.eye(2)             # identity matrix 


# multiplication of sparse matrices
np.dot(S,c)           # does not work as numpy's dot function on;ly works with (dense) numpy arrays 
np.dot(S.todense(),c) # works but sparse matrix is coverted to a dense matrix, which is usually too expensive
S.dot(c)              # works
S*c                   # works as well (no componentwise multiplication!)

# indexing of sparse matrices
S[0,0]        # does not work
S = S.tolil() # convert to, e.g., list of lists format
S[0,0]        # works

# solve sparse systems
import scipy.sparse.linalg as spla  
x = spla.spsolve(S,c) # works but with warning
S = S.tocsr()         # convert to compressed sparse row format (or to csc)
x = spla.spsolve(S,c) # works


# graphical output using matplotlib
import matplotlib
matplotlib.__version__ # print version number of matplotlib on screen (should be larger than 1.2 to offer all functions needed in this course)


# import matplotlib.pyplot
import matplotlib.pyplot as plt


# plot structure of sparse matrix
A = sp.rand(10,10,0.7)
plt.spy(A)
plt.show()


# plot sin
x = np.linspace(0,2*np.pi,1000)
plt.plot(x,np.sin(x))
plt.xlabel('x')
plt.ylabel('y')
plt.show()


# plot sin and cos with legend
plt.plot(x,np.sin(x),label='sin')
plt.plot(x,np.cos(x),label='cos')
plt.legend()
plt.xlabel('x')
plt.ylabel('y')
plt.show()


# write and import an own module
import mymodule as mymod
x = 2.0
y = mymod.squareroot(x)
print(str(y) + ' is the square root of ' + str(x))
# mymodule.py:
# from numpy import sqrt
#
# def squareroot(x):
#   return sqrt(x)


# reload a module in Python2
reload(mymod)


# reload a module in Python3
import imp 
imp.reload(mymod)