#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Nov 16 13:44:30 2017

@author: valerie
"""

import os
os.chdir('/users/valerie/Dropbox/ENSEIGNEMENT/PYTHON')

import numpy as np
import matplotlib.pyplot as plt
# Régression linéaire, skin cancer

# importation des données -----------------------------------------------------
tas = open('tas_1901_2015.csv', 'r')
lignes = tas.readlines()
tas.close()
tab = []
for chn in lignes:
    tab.append(chn.split(','))
    
temperature = np.zeros(len(tab)-1)
for k in range(len(tab)-1):
    temperature[k] = float(tab[k+1][0])
    
# Tracé de la série temporelle ------------------------------------------------   
tps = np.arange(1901,2016,1/12) # 12 mois par an
plt.plot(tps,temperature)  
plt.xlabel('Années')
plt.ylabel('Températures')
plt.title('Moyennes mensuelles') 
plt.show()


#  Moyennes annuelles ---------------------------------------------------------
# Solution sans utiliser de matrices, ni np.mean
def moyennes_annuelles(temperature): 
    m_yr = np.zeros(len(temperature)/12.)
    cnt = 0
    for yy in range(len(m_yr)):
        S = 0
        for k in range(12):
            S = S+temperature[cnt]
            cnt += 1
        m_yr = S/12.0
    return(m_yr)

# Solution sans utiliser de matrices
def moyennes_annuelles(temperature): 
    m_yr = np.zeros(len(temperature)/12.)
    cnt = 0
    for yy in range(len(m_yr)):
        m_yr = np.mean(temperature[(yy*12+1):(yy*12+12)])
    return(m_yr)

    
# 
def moyennes_annuelles(temperature): 
    temp = np.reshape(temperature,[12,int(len(temperature)/12)])
    return(np.mean(temp,axis=0))

yr = np.arange(1901,2016) # 12 mois par an
m_yr = moyennes_annuelles(temperature)
plt.plot(yr,m_yr)
plt.axis([1900 ,2016, 8, 15])

#  Lissage --------------------------------------------------------------------
def moyennes_mobiles(temperature,bw):
    """ bw : largeur de fenetre, entier pair"""
    temp = np.copy(temperature)
    b = int(bw/2)
    for k in range(b+1,len(temp)-b):
        fenetre = int()
        temp[k] = np.mean(temperature[k-b:k+b])
    return(temp)
m_liss = moyennes_mobiles(temperature,48)
plt.plot(tps,temperature,"k")
plt.plot(tps,m_liss,"r")
plt.show()
#plt.axis([1900 ,2016, 8, 15])


# Estimation de la tendance --------------------------------------------------- 
plt.plot(tps[25:len(tps)-25],m_liss[25:len(tps)-25])
plt.axis([1900 ,2016, 8, 15])
plt.title("Serie des moyennes mobiles, sans les bords")

class RegLin:
    def __init__(self):
        self.intercept = 0
        self.pente = 0
    
    def fit(self,x,y):
        xc = x-np.mean(x)
        yc = y-np.mean(y)
        a = np.sum(xc*yc)/np.sum(xc**2)
        b = np.mean(y)-a*np.mean(x)
        self.intercept = b
        self.pente = a
        
    def getParams(self):
        print("Pente = ", self.pente)
        print("Intercept = ", self.intercept)
        return([self.intercept,self.pente])
        
    def plotReg(self,x,y):
        plt.plot(x,y,"+")
        t = np.arange(np.min(x),np.max(x),0.1)
        plt.plot(t,self.intercept+self.pente*t,"r")
        plt.show()
    
    def predict(self,xnew):
        ynew = self.intercept+self.pente*xnew
        return(ynew)
        
        
x = tps[24:len(tps)-24] 
y = m_liss[24:len(tps)-24]
mod = RegLin()
mod.fit(x,y)
mod.getParams()
mod.plotReg(x,y)   

# Prédiction ------------------------------------------------------------------
mod.predict(np.array([2020,2050]))

# Significativité
tps_array = np.reshape(x,[int(len(x)/12),12])
temp_array = np.reshape(y,[int(len(x)/12),12])

B = 50
pente = np.zeros(50)
for b in range(B):
    ii = np.random.permutation(np.arange(temp_array.shape[0]))    
    mod1 = RegLin()
    mod1.fit(x,np.reshape(temp_array[ii,:],len(y)))
    pente[b] = mod1.getParams()[1]
plt.boxplot(pente)
mod.getParams()    # On remarque que la pente estimée 
                   # sur les données non permutées est bien plus forte
                   # que celles obtenues sur les données permutées.
                   # On conclut que la tendante au réchauffement est significative.