
# coding: utf-8

## Introduzione al Python per il calcolo scientifico.

# 28/11/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/

## Python

# Linguaggio di programmazione  ad alte prestazioni che può essere utilizzato in modo efficiente anche  per il calcolo scientifico.
#  
# Le istruzione Python possono essere eseguito utilizzando alcune IDE come spyder  o IPython e associandolo a pacchetti quali numpy, scipy, simpy,  diventa un potenete strumento di calcolo scientifico, confrontabile con Matlab, Octave o R per quanto riguarda la facilità di utilizzo ma con prestazioni che possono essere  molto più efficienti.  Fra i linguaggi di alto livello l'unico che lo supera per prestazioni
# è il nuovissimo linguaggio di nome Julia. Per il calcolo scientifico il Fortran rimane comunque il linguaggio più efficienti al pari del c, ma il livello di difficoltà  per  scrivere un codice in Fortran o c è sicuramente molto più elevato di quello richiesto per Python.

# Possiamo eseguire alcune istruzioni Python usando una interfaccia come Ipython o Spyder,
# entrati in questo ambiente compare il prompt, il che significa che il calcolatore \`e pronto a ricevere le
# istruzioni Python ed ad eseguirle.  
# Se vogliamo eseguire la somma 3+2 basta
# scrivere l'istruzione e premere il tasto invio. 
# Otteniamo immediatamente la risposta.

# In[167]:

3+2


# e il calcolatore è di nuovo pronto a ricevere un nuovo comando. Se vogliamo assegnare il risultato dell'operazione alla variabile di nome var
# basta scrivere:

# In[168]:

var = 3+2
var


# Questa istruzione si chiama **istruzione di assegnazione**. Il risultato dell'operazione non viene visualizzato. Per visualizzare il risultato dobbiamo digitare al prompt il nome della variabile oppure scrivere una istruzione di stampa:
# 

# In[169]:

print(var)


# Possiamo assegnare a una variabile il valore di una espressione più complicata.  Per
# far questo dobbiamo sapere a quali simboli corrispondono le principali
# operazioni aritmetico-logiche:
# 
# 
# * $+$    addizione
# * $-$    sottrazione 
# * $*$    moltiplicazione 
# * $/$     divisione 
# * $//$   divisione intera
# * $**$    elevamento a potenza 
# * &     and logico 
# * |     or logico 
# * ~     not logico 
# * ==    uguale  
# * !     diverso 
# * $<>$    diverso
# * $<$     minore 
# * $<=$    minore uguale 
# * $>$     maggiore
# * $>=$    maggiore uguale
# 
# Le espressioni aritmetiche e logiche seguono le classiche regole di precedenza,
# che è possibile modificare utilizzando le parentesi tonde.
# 
# Per avere informazioni su tutte le operazioni che possono essere eseguite in
# Python si può consultare The Python Language Reference, che in spyder si ottiene digitando 
# dal tasto $?$  il link alla Python documentation. La funzione help() permette di aprire un help interattivo sulle console di IPython e 
# di spyder. In IPython si può digitare anche il nome di una funzione seguito dal ?.
# 
# 

# In[170]:

get_ipython().magic('pinfo help')


# In[172]:

help('==')


# Possiamo ora eseguire, per calcolare $$\frac{5 \cdot 2}{(3+2)^2}$$
# 

# In[173]:

var=5*3/(3+2)**2 
print(var)
type(var)


# se state usando 2.7 il risultato e' zero la divisione viene eseguita fra interi perchè tutte le costanti sono intere. Questa e' una caratteristica di Python 2.7.  Per avere un risultato reale e' importante usare costanti reali. Calcoliamo $$\frac{5 \cdot 2}{(3+2)^2}$$

# In[174]:

var=5.*3./(3.+2.)**2. 
print(var)
type(var)


# Per poter lavorare in python 2.7 come in Python 3 si  utilizza l'istruzione divisione analoga a quella di Python 3.

# In[175]:

from __future__ import division
var=5*3/(3+2)**2 
print(var)
type(var)


# Un'altra sostanziale differenza fra Python 2.7 e Python 3 e' nella funzione print. Se usate Python 2.7 importate  la funzione print di Python 3:

# In[179]:

from __future__ import print_function
get_ipython().magic('pinfo print')


# In[180]:

var=5*3/(3+2)**2 
print('var = ',var)
type(var)


# Si possono  calcolare i valori delle principali funzioni. Ad
# esempio sono predefinite le funzioni in una variabile:
# 
# 
# * abs(x)  valore assoluto di x     
# * log(x)  logaritmo naturale          
# * exp(x)  esponenziale  $e^x$
# * sin(x)  seno,                        
# * cos(x)  coseno,                      
# * tan(x)  tangente.
# 
# La maggior parte delle funzioni aritmetiche fanno parte del modulo numpy.
# Per poter utilizzare il modulo numpy bisogna, dopo averlo installato importarlo,
# se si usa python in ambiente spyder numpy è automaticamente  importato all'avvio della console.
# Le istuzioni per importare il modulo numpy sono, per esempio:
# 

# In[182]:

import numpy as np


# 
# Tutti i metodi e le funzioni  di numpy devono essere precedute da np. per poter essere eseguite.

# In[183]:

var = (5 * 3/(3+np.cos(5)))**np.log(2) * 100 
print('var = ',var)
type(var)


# In[184]:

var = (5 * 3/(3+cos(5)))**log(2) * 100 
var
type(var)


# 
# Sono di particolare utilità anche alcune  *costanti*. Le
# costanti sono dei contenitori di valori che non cambiano nel corso
# dell'esecuzione di un algoritmo. In numpy vi sono delle costanti predefinite,
# come, ad esempio, la costante np.pi  che contiene il valore di $\pi$:
# 

# In[185]:

np.pi


# Esiste anche il modulo math che permette di lavorare con le funzioni standard. Anche questo modulo va importato prima di essere utilizzato:
# 

# In[186]:

import math as mt


# In[187]:

var = (5 * 3/(3+mt.cos(5)))**mt.log(2) * 100 
print('var = ', var)
type(var)


# Le istruzioni possone essere inserite in un file di tipo script

# In[190]:

get_ipython().magic('run stampavar')


# Le variabili possono essere di tipi diversi. Le operazioni automaticamente cambiano i tipi di dati. La function type permette di controllare il tipo del dato. 

#### Tipo Intero

# In[191]:

n=3
print(n)
type(n)


# In[192]:

n=10**50
print(n)
type(n)


# In[193]:

n = np.int(3)
type(n)


# In[194]:

n = np.int(10**30)
type(n)


#### Reali 

# In[195]:

type(3.9)


# In[196]:

v=np.float32(1.0)
v


# In[197]:

type(v)


# In[198]:

x=np.float64(1.0)
type(x)


# Possiamo decidere il comportamento da assumere quando si verifica una operazione invalida, di default si ignora l'errore

# In[201]:

x=np.float32(0.0)
y=np.float32(1.0)
np.seterr(all='ignore')
z= y/x
print('z = ',z )  


# In[202]:

s= x/x
print( 's = %(s)10.00f' % {"s":s})


# La precisione di macchina e' il più piccolo numero positivo  che aggiunto ad uno ci genera il successivo numero di macchina, usando 32 bit   è:

# In[203]:

np.finfo(x).eps


# Il più piccolo numero floating point positivo  usando la rappresentazione a 32 bit è

# In[204]:

np.finfo(x).tiny


# Il massimo numero floating point positivo usando la rappresentazione a 32 bit è:

# In[205]:

np.finfo(x).max


# In[206]:

type(np.finfo(x).max)


# In questo caso  è possibile usare dei numeri più grandi, infatti si passa automaticamente alla rappresentazione floating point a 64 bit

# In[207]:

print(np.finfo(x).max*2)
type(np.finfo(x).max*2)


# Per la rappresentazione a 64 bit la precisione di macchina è:
# 

# In[208]:

x = np.float64(0.0)
np.finfo(x).eps


# In questo caso usare numeri più grandi del massimo numero rappresentabile genera INF oppure un errore di overflow:

# In[209]:

np.finfo(x).max


# non possiamo superare il massimo numero rappresentabile, questa operazione genera inf

# In[211]:

np.finfo(x).max*2


# possiamo tranquillamente usare numeri più piccoli del più piccolo numero floating point positivo. Questi numeri vengono detti numeri denormali e consentono l'underflow graduale.

# In[212]:

np.finfo(x).tiny


# In[217]:

np.finfo(x).tiny/2**53


#### Tipo Complesso

# In[218]:

print(complex(-1.0))


# In[219]:

a =-1+0j
b = 0+1j
print(a+b)


# In[220]:

print(np.sqrt(-1.0))


# In[221]:

print(np.sqrt(complex(-1.0)))


# In[222]:

import scipy as sp
print(sp.sqrt(-1.0))


# In[223]:

import cmath
print(cmath.sqrt(-1.0))


#### Stringhe

# In[224]:

c = "la mia prima stringa"
print(c)


# In[225]:

c="""la mia seconda stringa 
su più righe
in tutto tre righe"""
print(c)


#### Liste

# In[227]:

v = [1, 2, 3]
v
type(v)


# In[228]:

v[1]


# In[229]:

v[0]


# In[230]:

v.append(4)
print(v)


# In[231]:

g=v.pop()
print(g)
print(v)


# L'istruzione whos permette di visualizzare le variabili in uso nel workspace:

# In[232]:

get_ipython().magic('whos')


# IPython supporta le  "magic functions" facendo preceder il comando dal carattere %. Per esempio run e whos sono magic functions. 
# Essendo per default settato il comando automagic si può anche non mettere il simbolo % prima delle magic functions.
# 

# In[233]:

del var


# In[234]:

whos


# In spyder per visualizzare le variabili si usa  la finestra variable explorer. 

#### Istruzioni di assegnazione

# In[235]:

a = 3
b = 2*a
print(type(b))
print(b)
d=a*b
print(d)


# In[236]:

b = 'Ciao'
print(type(b))

c = b + b
print(c)

d=2*b
print(d)


# Per le liste l'istruzione di assegnazione fa condividere alle due variabili le stesse locazioni di memoria:

# In[237]:

a =[1,2,'ciao']
b=a
b[2]=10
print(a) 


# Il modo più semplice per fare la copia di una lista e copiare non l'oggetto ma tutti gli elementi

# In[238]:

a =[1,2,'ciao']
b=a[:]
b[2]=10
print(a)
print(b)


# La lunghezza di una lista si trova con la funzione len

# In[239]:

len(a)


#### Tuple

# In[240]:

a=(1,2,'ciao')
b=a
a[0]


# Non si possono aggiornare gli elementi di una tupla come quelli di una lista

# In[242]:

a[0]=10


# la  tupla si può spacchettare

# In[243]:

(g1,g2,g3)=a
print(g1)
print(g2)
print(g3)


#### Dictionary

# In[244]:

voti = {"Maria":10,"Antonio":7, "Giovanni":9}
print(voti)


# In[245]:

voti.values()


# In[246]:

voti.keys()


# In[247]:

voti["Giovanni"]


# In[248]:

len(voti)


#### Funzioni

# In[249]:

def somma(a,b=1):
    s=a+b
    return(s)


# In[250]:

somma(2,3)


# In[251]:

def somma(a,b=1):
    s=a+b
    a=4
    return(s)


# In[252]:

a=3
print(somma(a))
print(a)


# In[253]:

somma2 = somma  # le funzioni sono oggetti posso usarli in operazioni di assegnazione 
print(somma2(2,3))


# In[254]:

somma2=2*somma # non posso però moltplicare una funzione con uno scalare


# In[255]:

somma2=somma*somma # e neanche due funzioni fra loro


# Funzioni semplici si possono implementare usando le cosidette lambda function

# In[256]:

print("lambda function utilizzo = ",  (lambda a,b:a+b)(2,3)  )


# In[257]:

sommal =  (lambda a,b:a+b)  # assegnazione 
print(sommal(2,3))


# Esempio semplice ciclo di selezione

# In[258]:

def somma(vl):
    """
    somma dei primi due elementi di una lista
    input  lista
    output  tupla (somma , messaggio di errore)
    """
    n=len(vl)
    if n == 2:
       s=vl[0]+vl[1]
       errormessage = "nessun errore"
    elif n > 2:
       s=vl[0]+vl[1]
       errormessage = "solo i primi due elementi sono stati sommati"
    else:
       s=vl[0]
       errormessage = "s contiene un singolo elemento"         
    vl[0]=5
    return(s, errormessage)


# In[259]:

vl = [1,2,3]
print(somma(vl))
print(vl)


# In[260]:

help(somma)


# Esempio ciclo di ripetizione while.
# 
# Costruiamo una funzione che abbia come dati di input la base usata dal calcolatore e in output la precisione di macchina. Ricordo che la precisione di macchina e'  il numero, potenza della base, che aggiunto ad 1 genera il successivo numero di macchina. Un numero più piccolo della precisione di macchina aggiunto ad uno genera zero.

# In[261]:

import numpy as np
def epsmach(base):
    epsmach = 1
    k = 0
    while ((1+epsmach) > 1):
        epsmach = epsmach/base
        k = k+1
    return (epsmach*2, k-1)

print(epsmach(2))
print(np.finfo(1.0).eps)
epsmio, kmio = epsmach(2)
print('eps =',epsmio,',  k =',kmio)


# Esempio ciclo di ripetizione for.
# Funzione che calcola per n volte la radice quadrata di un numero e per n volte il quadrato del numero risultante. Esempio di ciclo for, esempio di errori di arrotondamento

# In[262]:

def  esempio_radici(xt,n):
    x=xt
    for i in range(0,n,1):
      x = np.sqrt(x)
     
    for i in range(0,n):
      x = x**2

    errore_relativo = abs(xt-x)/abs(xt)
    return (x, errore_relativo)


# In[266]:

n = 55
x, errore = esempio_radici(100,n)
print('x =',x,',n = ', n,', errore= ',errore)


#### Vettori e Matrici

# In[267]:

A = np.array([3., 5., 6.]) # per utilizzare vettori e matrici facciamo ricorso a numpy
print(A)


# In[268]:

A.size


# In[269]:

A[0] + A[1] + A[2]


# In[270]:

np.sum(A)


# In[271]:

np.min(A)


# In[272]:

np.argmin(A)


# In[273]:

A[np.argmin(A)]


# In[274]:

np.argmax(A)


# In[275]:

A[np.argmax(A)]


# In[276]:

A=np.array([[  3.,   5.,   7.],
           [  8.,   9.,  10.],
           [  2.,   3.,   5.]])
A


# In[277]:

A[0,0]


# In[278]:

A[:,0]


# In[279]:

A[0,:]


# In[280]:

A[1:3,1:3]


# In[281]:

np.linspace(0,5,20)


# In[282]:

np.logspace(0,5,10)


# In[283]:

get_ipython().magic('matplotlib inline')
import matplotlib.pyplot as plt

x=np.logspace(0,9,500)
r1=x*(np.sqrt(x**2+1)-x)
r2=x*np.sqrt(x**2+1)-x**2
r3=x/(np.sqrt(x**2+1)+x)
v=plt.semilogx(x,r1,'r.',x,r2,'b.',x,r3,'g-',lw=2)


# In[284]:

import matplotlib.pyplot as plt
get_ipython().magic('matplotlib qt  # funzione di Ipython')
x=np.logspace(0,9,500)

r1=x*(np.sqrt(x**2+1)-x)
r2=x*np.sqrt(x**2+1)-x**2
r3=x/(np.sqrt(x**2+1)+x)

v=plt.semilogx(x,r1,'r.',x,r2,'b.',x,r3,'g-',lw=2)
print(v)
print(v[0])


# Approssimazione della derivata usando le differenze in avanti

# In[285]:

#%%
h = np.array([0.1, 0.01, 0.001, 1e-4, 1e-7, 1e-12])
# approssimazione della derivata di exp(1) differenze in avanti
dfh = (np.exp(1+h)-np.exp(1))/h
# valore esatto
df = np.exp(1)
# errore assoluto
errass = np.abs(df-dfh)
errrel = errass/np.abs(df)
# grafico in scala logaritmica
plt.loglog(h,errrel,'o-b')
plt.show()


# In[286]:

#%%
# costruiamo un vettore con elementi
# che variano da 10**(-16) a 1
h = np.logspace(-19,0,500)
# approssimazione della derivata
dfh = (np.exp(1+h)-np.exp(1))/h
# valore esatto
df = np.exp(1)
# errore assoluto
errass = np.abs(df-dfh)
errrel = errass/np.abs(df)
# grafico in scala logaritmica
plt.loglog(h,errrel,'*-r')
plt.show()


# In[287]:


# esempio di problema mal condizionato
# trovare le radici multiple di polinomi
# (x-1)^4

p4 = np.array([1, -4, 6, -4, 1])
r4 = np.roots(p4)
# roots calcola le radici
# errore relativo
err4 = abs(r4-1)
print ("errore assoluto polinomio di grado 4" )
for rad in err4:
    print (rad)
# calcoliamo le radici di (x-1)^6

p6 = np.array([ 1, -6, 15, -20, 15, -6, 1])
r6 = np.roots(p6)
# roots calcola le radici 
# errore relativo
err6 = abs(r6-1)

print( "errore assoluto polinomio di grado 6" )
for rad in err6:
    print (rad)


# Grafico di un polinomio in un intervallo specificato

# In[288]:

import matplotlib.pyplot as plt
get_ipython().magic('matplotlib qt')
p = np.array([1,-6,15,-20,15,-6,1]) # coefficienti del polinomio  (x-1)^6
x = np.linspace(0.995,1.005,500)  # intervallo [0.995,1.005] 
y = (x-1)**6
ye = np.polyval(p,x)#funzione per valutare un polinomio in un punto
plt.plot(x,y,'-b',lw=4)
plt.plot(x,ye,'og')


#### Moduli

# I moduli sono dei file che contengono oggetti Python come funzioni o classi o variabili che possono essere importati e utilizzati da altri moduli e possono far parte di pacchetti più complessi.
# Scriviamo il modulo nel file primo_modulo.py

# In[ ]:

import primo_modulo as pm


# In[ ]:

print(pm.a)
print(pm.b)
print(pm.somma(2,3))
print(pm.prodotto(2,2))


# Se modifichiamo il codice del modulo pm  aggiungendo, per esempio, una funzione dobbiamo ricaricarlo con l'istruzione

# In[ ]:

import importlib as imp
imp.reload(pm)


#### Files 

# In[ ]:

pf =open("provafile.txt","w")
pf.write("CIAO, provo a scrivere nel file\n")
pf.write("CIAO, provo ancora a scrivere nel file\n")
pf.close()


# In[ ]:

pf =open("provafile.txt","r")
print(pf.read())
pf.close()
pf =open("provafile.txt","r")
print(pf.readline())
print(pf.readline())


# In[ ]:

pf =open("provafile.txt","w")
pf.write("CIAO; provo a scrivere nel file; 1\n")
pf.write("CIAO; provo a scrivere nel file; 2\n")
pf.writelines(["CIAO; provo a scrivere nel file; 3\n","CIAO; provo a scrivere nel file; 4\n" ])
pf.close()


# In[ ]:

pf =open("provafile.txt","r")
file = pf.read()
print(file)
pf.close()


# Alcuni esempi di operazioni sulle stringhe

# In[ ]:

for k in file:
    print (k)


# In[ ]:

sfile = set(file)
print(sfile)


# In[ ]:

print(file.upper())


# In[ ]:

print(file.lower())


# In[ ]:

print(file.split())


# In[ ]:

print(file.split(';'))


#### Classi

# In[ ]:

class Punto: 
  def __init__(self, x=0, y=0): 
    self.x = x 
    self.y = y 
  def __str__(self): 
    return '(' + str(self.x) + ', ' + str(self.y) + ')' 


# In[ ]:

P1= Punto(2,3)
print(P1)


# In[ ]:

class Punto: 
    """Classe punto"""
    def __init__(self, x=0, y=0): 
        self.x = x 
        self.y = y     
    def __add__(self,AltroPunto):
        return Punto(self.x + AltroPunto.x, self.y + AltroPunto.y) 
    def __str__(self): 
        return '(' + str(self.x) + ', ' + str(self.y) + ')' 


# In[ ]:

P1= Punto(2,3)
P2=Punto(3,4)
P3= P1+P2
print(P3)


# In[ ]:

print(Punto.__doc__)

