
# coding: utf-8

## Introduzione al Python per il calcolo scientifico.

# 12/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/
# 

#### Interpolazione

# In[2]:

import numpy as np
import matplotlib.pyplot as plt
# estremi dell'intervallo
a = -1
b =  1
#%%

xx = np.linspace(a,b,500)
f = np.exp
yy = f(xx)

n=8
xnodi = np.linspace(a,b,n+1)
fnodi = f(xnodi)


# In[4]:

import scipy.interpolate as sinterp
import matplotlib.pyplot as plt
get_ipython().magic('matplotlib inline')
s1 = sinterp.interp1d(xnodi, fnodi, 1) # 'linear' , interpolazione lineare a tratti
plt.plot(xx, s1(xx), xnodi, fnodi, 'o', xx, yy, '--')


# In[5]:

get_ipython().magic('matplotlib inline')
get_ipython().magic('run spline_dem.py # demo interpolazione lineare a tratti e spline cubica')
# esempi di concatenazione di array e di creazione di structured array


# In[7]:

# convergenza spline lineare  
np.log2(risultati["errore"][0:5]/risultati["errore"][1:6])



# In[8]:

# convergenza spline cubica
np.log2(risultati["errore"][6:11]/risultati["errore"][7:12])


# risultati
# 

# In[10]:

findcub=risultati["spline"]==b'cubica'
print( risultati[findcub])


##### Grafico di una funzione bidimensionale

# In[11]:

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[13]:

from matplotlib import cm
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[14]:

from scipy import interpolate
xnew, ynew = np.mgrid[-1:1:70j, -1:1:70j]
#rappresentazione B-spline della superficie 
tck = interpolate.bisplrep(x, y, z, s=0)
# valutazione della B-spline nella nuova griglia
znew = interpolate.bisplev(xnew[:,0], ynew[0,:], tck)

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


# In[15]:

# grafico della B-spline
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()


#### Zeri di funzione

# In[17]:

# calcolo degli zeri di una funzione
import numpy as np
import scipy.optimize as opt
import matplotlib.pylab as plt


a1=1.5
b1=2.5
# f1 funzione test
def f1(x):
    fx = x**3-2*x-5
    return fx
    
def df1(x):
    dfx = 3*x**2-2
    return dfx
  
a2 = 0.0
b2 = 1.0
# f2 funzione test
def f2(x):
    fx = 2-np.exp(x)
    return fx
    
def df2(x):
    dfx = -np.exp(x)
    return dfx
  
a3=0
b3=1
# f3 funzione test
def f3(x):
    fx = 2*x*np.exp(-15)-2*np.exp(-15*x)+1 
    return fx
    
# df derivata prima di f
def df3(x):
    dfx = 2*np.exp(-15)+30*x*np.exp(-15*x) 
    return dfx
    
# utilizzo la funzione test 1   
a = a1
b = b1
f = f1
df = df1    
    
x = np.linspace(a,b,100)
plt.plot(x,f(x))

# metodo delle bisezioni
bsci = opt.bisect(f,a,b,full_output=True)

print ('soluzione calcolata da bisect')
print (bsci [0])
np.seterr(all='ignore')
print ('bisect iterations',bsci[1].iterations)
print ('bisect converged', bsci[1].converged)
print ('bisect root', bsci[1].root)
print ('bisect function calls',bsci[1].function_calls)
print ('bisect flag',bsci[1].flag)
print ('%%%%%%%%%%%%%%%%%%%%%%')

# metodo di Newton
nsci = opt.newton(f,a,df,tol=1e-10,maxiter=100)

print( 'soluzione calcolata da Newton partendo da a')
print( nsci)
print( '%%%%%%%%%%%%%%%%%%%%%%')

nsci = opt.newton(f,b,df,tol=1e-10,maxiter=100)

print( 'soluzione calcolata da Newton partendo da b')
print( nsci)
print( '%%%%%%%%%%%%%%%%%%%%%%')
plt.show()


#### Approssimazione

# In[18]:

import numpy as np
import matplotlib.pylab as plt


# In[19]:

# Dati da approssimare
xnodi = np.array([ 0, 1, 2, 3, 4, 5])
fnodi = np.array([10, 25, 51, 66, 97, 118])


# In[21]:


# costruiamoci la retta di approssimazione nel senso dei minimi 
# quadrati (retta di best-fitting)

p = np.polyfit(xnodi,fnodi,1)
xx = np.linspace(0,5,100)
yp = np.polyval(p, xx)


# In[22]:

get_ipython().magic('matplotlib inline')
plt.plot(xx,yp,xnodi,fnodi,'o')
plt.show()


# In[23]:

# dati relativi alla conducibilita' termica del ferro
xnodi = np.array( [100, 200, 300, 400, 500, 600, 700, 800, 900, 1000])
fnodi = np.array([1.32, 0.94 ,0.835, 0.803, 0.694, 0.613 ,0.547, 0.487 ,0.433, 0.38])


# In[25]:

p1 = np.polyfit(xnodi,fnodi,1)
xx = np.linspace(0,1000,100)
yp1 = np.polyval(p1,xx)
#
# grafico della retta #grafico dei dati
plt.plot(xx,yp1,xnodi,fnodi,'o')
plt.show()


# In[26]:

# costruiamo il polinomio di 2 grado di approssimazione nel senso dei minimi 
# quadrati 
p2 = np.polyfit(xnodi,fnodi,2);
yp2 = np.polyval(p2,xx)
#
# grafico 
plt.plot(xx,yp2,xnodi,fnodi,'o')
plt.show()


# In[27]:

# costruiamo il polinomio di 3 grado di approssimazione nel senso dei minimi 
# quadrati 
p3 = np.polyfit(xnodi,fnodi,3)
yp3 = np.polyval(p3,xx)
#
# grafico 
plt.plot(xx,yp3,xnodi,fnodi,'o')
plt.show()


#### Integrazione

# In[28]:

import scipy.integrate as ints

fun = lambda x : 1-(1-x)**(0.5)
a=0.0
b=1.0

xx = np.linspace(a,b,500)
yy = fun(xx)    
plt.plot(xx,yy)


# In[29]:

# calcolo dell'intergale di una funzione, interfacci alla libreria Fortran QUADPACK
Iq,abserr,info=ints.quad(fun,a,b,full_output=1)
print(Iq)
print(abserr)


# In[31]:

get_ipython().magic('matplotlib inline')
import numpy as np
import matplotlib.pyplot as plt
def f(x):
    return   1-(1-x)**(0.5) #(x-3)*(x-5)*(x-7)+85

x = np.linspace(0, 1, 200)
y = f(x)
a, b = 0.1, 0.9
xint = x[np.logical_and(x>=a, x<=b)][::30]
yint = y[np.logical_and(x>=a, x<=b)][::30]
plt.plot(x, y, lw=2)
plt.axis([0, 1, 0, 1])
plt.fill_between(xint, 0, yint, facecolor='gray', alpha=0.8)
plt.text(0.5 * (a + b), 0.5,r"$\int_a^b f(x)dx$", horizontalalignment='center', fontsize=20);


#### Equazioni Differenziali Ordinarie

# In[32]:

# soluzione numerica di una equazione differenziale ordinaria
# problema ai valori iniziali
yin = np.array([2.0,0.0])
x = np.arange(0, 2.0, 0.01)
t = x
def vdpol(y,t):
    return [ y[1],((1-y[0]**2)*y[1]-y[0])/1.0e-1]

# funzione odeint interfacciao alla libreria Fortran lsoda di ODEPACK 
y = ints.odeint(vdpol, yin, t)
plt.plot(t,y)     


# In[38]:

yin = [2.0,0.0]
x = np.arange(0, 2.0, 0.01)
t = x
def vdpolode(t,y):
    return [y[1],((1-y[0]**2)*y[1]-y[0])/1.0e-6]


# In[40]:

# function ode interfaccia per utilizzare le librerie FORTRAN
# vode, zvode, lsoda, dop853, dopri5
r=ints.ode(vdpolode).set_integrator('lsoda')# vode,method='bdf') #dop853 #dopri5
r.set_initial_value(yin, 0)
tf = 2
dt = 0.01
ts =[0]
ys =[yin]
while r.successful() and r.t < tf:
     r.integrate(r.t+dt)
     ys.append(r.y)
     ts.append(r.t)       
plt.plot(ts,ys)
        


#### Pandas per l'analisi dei dati

# In[42]:

import pandas 
data = pandas.read_csv("data/weather_year.csv")
data


# http://synesthesiam.com/posts/an-introduction-to-pandas.html
# This gives us a lot of information. First, we can see that there are 366 rows (entries)
#  -- a year and a day's worth of weather. Each column is printed along with however many 
# "non-null" values are present. We'll talk more about null (or missing) values in pandas 
# later, but for now we can note that only the "Max Gust SpeedMPH" and "Events" columns 
# have fewer than 366 non-null values. Lastly, the data types (dtypes) of the columns are 
# printed at the very bottom. We can see that there are 4 float64, 16 int64, and 3 object
#  columns.

# In[44]:

print(len(data))

#nomi delle colonne
data.columns


# In[45]:


print(data["EDT"])

print(data.EDT)



# In[46]:


#prime cinque righe 
data.EDT.head()


# In[47]:

#ultime cinque righe
data.EDT.tail()


# In[49]:


data["Mean TemperatureF"].head()


# In[50]:


data.columns = ["date", "max_temp", "mean_temp", "min_temp", "max_dew",
                "mean_dew", "min_dew", "max_humidity", "mean_humidity",
                "min_humidity", "max_pressure", "mean_pressure",
                "min_pressure", "max_visibilty", "mean_visibility",
                "min_visibility", "max_wind", "mean_wind", "min_wind",
                "precipitation", "cloud_cover", "events", "wind_dir"]

data


# In[51]:

# i valori omessi sono tralasciati           
data.mean_temp.std()


# In[52]:


# istogramma
data.mean_temp.hist()


# In[53]:

data.mean_temp.plot()


# In[54]:



#deviazione standard su tutte le colonne 
data.std()


# In[56]:

# http://www.ilmeteo.it/portale/archivio-meteo/Bari/2013/Giugno             
datibarigiugno = pandas.read_csv("data/Bari-2013-Giugno.csv",sep=';')
datibarigiugno.columns


# In[57]:

mesi=["Gennaio","Febbraio","Marzo","Aprile","Maggio","Giugno","Luglio",
       "Agosto","Settembre","Ottobre", "Novembre","Dicembre"]
   

datibarif =  pandas.DataFrame()     
for nomemese in mesi:      
    file = "data/Bari-2013-"+nomemese+".csv"
    print(file)
    datibarimese = pandas.read_csv(file,sep=';')
    datibarif = (datibarif.append(datibarimese,ignore_index=True))



# In[58]:

datibarif


# In[59]:


datibarif.columns =["localita", "data", "mean_temp", "min_temp", "max_temp",
                "puntorugiada", "umidita", "visibilita_km", "ventomedia_km_h",
                "ventomax_km_h", "raffica_km_h", "pressioneslm_mb",
                "pressionemedia_mb", "pioggia_mm", "fenomeni",
                 ]
         


# In[60]:

datibarif.index=range(0,365)
datibarif.mean_temp.plot()


# In[ ]:



