
# coding: utf-8

## Introduzione al Python per il calcolo scientifico.

# 05/12/2014
# 
# docente: Francesca Mazzia
# 
# Il materiale è distribuito con licenza Creative Commons Attribuzione-Non commerciale-Non Opere Derivate 2.5 Italy License http://creativecommons.org/licenses/by-nc-nd/2.5/it/

#### Matrici e Vettori

# Scriviamo una matrice casuale con 3 righe  4 colonne 
# 

# In[2]:

import numpy as np

# np.random.rand genera numeri casuali appartenenti a una distribuzione uniforme
A = np.random.rand(3,4)*10
print(A)


# In[3]:

# questa e' una matrice rettangolare
# scriviamo una matrice quadrata con elementi casuali interi
A = np.random.randint(10,size=(3,3))
print(A)


# In[4]:

# gli elementi delle matrici si denotano
# usando l'indice di riga e di colonna
print(A[0,0])
print(A[1,1])
print(A[2,2])


# In[5]:

print(type(A))


# Trasposta

# In[6]:

print('trasposta\n',A.T)


# Prodotto riga per colonne e prodotto fra matrici

# In[7]:

print(np.dot(A,A))
print(A*A)


# Tipo matrix

# In[8]:

A = np.asmatrix(A)
print(A)    
print(type(A)) # A di tipo matrix
print(A.I)     # A.I inversa di A
print(A.I*A)    # il * e' il prodotto riga per colonne
print(np.dot(A.I,A))


# In[9]:

U=np.triu(A)   # operazioni con in input matrix in output array
print(type(U))


# In[10]:

A.shape  # numero di righe e colonne


# In[11]:

A.ndim  # dimensione della matrice, possiamo avere array a piu' dimensioni


# In[12]:

O = np.zeros( (3,3) ) # matrice nulla
print(O)
print(type(O))


# In[13]:

I = np.eye(3) # matrice identità
print(type(I))
print(I)


# In[14]:

d = np.diag(np.array([1, 2, 3, 4])) # array con elementi sulla diagonale
print(d)


# In[15]:

narray = np.arange(1,10,2) # arange array
print(narray)


# In[17]:

B = np.empty((3,4) ) # array vuoto non inizializzato
print(B)


# In[18]:

import matplotlib.pyplot as plt
get_ipython().magic('matplotlib inline')
image = np.random.rand(30, 30)
plt.imshow(image, cmap=plt.cm.hot)    
plt.colorbar()   


# In[19]:

a = np.arange(10) 
a[::-1]      # ordine inverso


# In[20]:

a = np.arange(10000)
get_ipython().magic('timeit a + 1')


# In[21]:

l = range(10000)
get_ipython().magic('timeit [i+1 for i in l]')


# In[22]:

a = np.array([1, 2, 3, 4])
b = np.array([4, 2, 2, 4])
a == b


# In[23]:

a > b


# In[24]:

np.array_equal(a,b)


# In[25]:

np.sin(a)


# In[26]:

np.sin(a[0])


# In[27]:

f = np.random.randint(10,size=(2,2,2))
print(f)
print(f.ndim)


# In[28]:

f[0,0,:]


# In[29]:

f[0,:,0]


# In[30]:

np.sum(f,axis=2)


# In[31]:

import scipy as sp


# In[33]:

import  linalgmod
import  time

dim=100
A=np.random.rand(dim,dim)
B=np.random.rand(dim,dim)
#import importlib as imp
#imp.reload(linalgmod)  #per ricaricare il modulo se si fanno modifiche


# In[34]:

t1=time.clock()
C1=linalgmod.prodmatmat(A,B)
t2=time.clock()-t1
print( "clock = ", t2 )

 
# A e B sono array il prodotto riga per colonne si esegue utilizzando 
# np.dot
t1=time.clock()
C2=np.dot(A,B)
t2d=time.clock()-t1
print( "clock = ", t2d )


print( "errore", np.max(np.abs(C1-C2)) )

# C1 e C2 contengono gli stessi elementi


# In[35]:

import scipy.linalg as la
# la function inv calcola l'inversa di A 
AI = la.inv(A)


# In[36]:

print( np.dot(AI,A) ) 

# una matrice che non ha inversa si dice
# singolare
print(la.det(A)) #  det calcola il determinate di A

A = np.random.rand(3,4)
# norma di Frobenious
print('norma frobenious= ',la.norm(A) )
# norma 1
n1 = la.norm(A,1)
print('norma 1 = ',n1)

# norma infinito
ninf = la.norm(A,np.inf)
print('norma infinito = ',ninf)
# norma 2
n2= la.norm(A,2)
print('norma 2 = ',n2 )


# In[37]:

a=np.array([-4, -3, -2, -1,  0,  1,  2,  3,  4])
b = a.reshape((3, 3))
print(b)


# Sistema lineare: trovare  x tale che A*x = b

# In[38]:

np.lookfor('solve linear system')


# In[39]:

A = np.random.rand(3,3)
b = np.array([1,2,3])
x = np.linalg.solve(A,b)
print(x)
type(x)


# In[40]:

print(type(A))
print(np.dot(A,x)-b) # residuo


# In[42]:

print(x.shape)
y = x.T
print(y.shape)
print(A*x-b)


# Risolviamo un sistema lineare con matrice dei coefficienti
# la matrice di Hilbert.La matrice viene costruita nel modulo linalgmod.
# Il termine noto lo scegliamo in modo che la soluzione x
# sia il vettore di tutti  uno.

# In[52]:

n = 7
A = linalgmod.hilb(n)
#print(A)
x = np.ones( (n,1))
b = np.dot(A,x)
#print(b)


# risolviamo il sistema, calcoliamo il residuo e l'errore assoluto

# In[53]:

xa = np.linalg.solve(A,b)
xb = la.solve(A,b)


# In[54]:

print(la.norm(x-xa,2))
print(la.norm(x-xb,2))


# Stampiamo il numero di condizione di A. Il numero di condizione nella soluzione di un sistema lineare ci dice come una perturbazione sui dati di input (A e b) si ripercuote sul risultato (xa). Il numero di condizione $$\kappa(A) = \| A \| \| A^{-1}\|$$ 

# In[55]:

np.linalg.cond(A)


# la matrice di Hilbert è mal condizionata

# In[ ]:

Possiamo anche prima fattorizzare la matrice a poi risolvere i sistemi
lineari usando la fattorizzazione. La fattorizzazione usata dagli 
algoritmi classici è la fattorizzazione LU con pivot


# In[56]:

#%%
# function lu di scipy
p,l,u=la.lu(A)
la.norm(A-np.dot(np.dot(p,l),u),1)


# In[57]:

np.linalg.cond(l)


# In[58]:

np.linalg.cond(u)


# In[59]:

# function lu della libreria lapack
x1=la.solve(A,b)
lu,piv=la.lu_factor(A)
# piv contiene il vettore delle permutazioni
# la riga i e' stata cambiata con la riga piv[i


# In[61]:

x2=la.lu_solve((lu,piv),b)
la.norm(x2-x,2)


# visualizziamo il residuo. Un residuo piccolo non implica un errore piccolo se la matrice è mal condizionata.

# In[62]:

la.norm(np.dot(A,x2)-b,2)


#### Matrici Sparse

# In[63]:

import numpy as np
import scipy.sparse as sparse
A = sparse.csr_matrix([[1, 2, 0], [0, 0, 3], [4, 0, 5]])
print(A)


# In[65]:

v = np.array([1, 0, -1])
A.dot(v)


# per eseguire il prodotto riga per colonne bisogna utilizzare metodi speciali

# In[66]:

np.dot(A, v)


# oppure trasformare la matrice in matrice densa, operazione che non conviene fare.

# In[70]:

A = sparse.coo_matrix( ([1,2,3,4],([0,1,2,3],[1,2,3,0]) ) )
print(A)
print(A.todense())
print(type(A.todense()))


# numero degli elementi diversi da zero in A

# In[71]:

A.nnz


# esistono vari formati di rappresentazione, fra cui:
# il formato per coordinate(coo), il formato compresso per colonne (csc) 
# e il formato compresso per righe (csr)

# In[72]:

B=sparse.csc_matrix(A)
print(B)


# l'algebra lineare per le matrici sparse è implementata nel pacchetto
# scipy.sparse.linalg,  risulta ottimizzata per matrici nel formato csc e csr

# In[74]:

import scipy.sparse.linalg as sla
b = np.random.rand(4,1)
x = sla.spsolve(A,b)
print(x)
x = sla.spsolve(B,b)
print(x)


# In[75]:

B.dot(x)


# In[76]:

N = 100000;
I = np.arange(0,N)
J = I+1
J[N-1]=0
A = sparse.coo_matrix( (np.arange(1,N+1),(I,J) ) )
b = I[:]
B=sparse.csc_matrix(A)
x1 = sla.spsolve(A,b)
x2 = sla.spsolve(B,b)
print(la.norm(x1-x2,np.inf))


# In[77]:

np.dot(A.todense(), v)


#### Interpolazione 

# In[ ]:

get_ipython().magic('matplotlib inline')
get_ipython().magic('run spline_dem.py')


# grafico di una funzione bidimensionale

# In[ ]:

x, y = np.mgrid[-1:1:20j, -1:1:20j]
z = (x+y) * np.exp(-6.0*(x*x+y*y))

plt.figure()
plt.pcolor(x, y, z)
plt.colorbar()
plt.title("funzione con elementi sparsi.")
plt.show()


# In[ ]:

from mpl_toolkits.mplot3d import Axes3D
from matplotlib.ticker import LinearLocator, FormatStrFormatter
import matplotlib.pyplot as plt
import numpy as np

fig = plt.figure()
ax = fig.gca(projection='3d')
surf = ax.plot_surface(x, y, z, rstride=1, cstride=1, cmap=cm.coolwarm,
        linewidth=0, antialiased=False)

ax.set_zlim(-0.3, 0.3)

ax.zaxis.set_major_locator(LinearLocator(10))
ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))

fig.colorbar(surf, shrink=0.5, aspect=5)

plt.show()
ax.zaxis.set_major_locator(LinearLocator(20))


plt.show()


# In[ ]:

from scipy import interpolate
xnew, ynew = np.mgrid[-1:1:70j, -1:1:70j]
tck = interpolate.bisplrep(x, y, z, s=0)
znew = interpolate.bisplev(xnew[:,0], ynew[0,:], tck)

plt.figure()
plt.pcolor(xnew, ynew, znew)
plt.colorbar()
plt.title("Interpolated function.")
plt.show()


# In[ ]:


fig = plt.figure()
ax = fig.gca(projection='3d')
surf = ax.plot_surface(xnew, ynew, znew, rstride=1, cstride=1, cmap=cm.coolwarm,
        linewidth=0, antialiased=False)

ax.set_zlim(-0.3, 0.3)

ax.zaxis.set_major_locator(LinearLocator(10))
ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))

fig.colorbar(surf, shrink=0.5, aspect=5)

plt.show()
ax.zaxis.set_major_locator(LinearLocator(20))



plt.show()


# In[ ]:



