text:
# Zonal analysis with SciPy

import os
import numpy as np
import scipy.stats
from osgeo import gdal

def get_bins(data):
    """Return bin edges for all unique values in data."""
    bins = np.unique(data)
    return np.append(bins, max(bins) + 1)

os.chdir(r'D:\osgeopy-data\Utah')
landcover_fn = 'landcover60.tif'
ecoregion_fn = 'utah_ecoIII60.tif'
out_fn = r'D:\Temp\histogram.csv'

# Read in the ecoregion data and get appropriate bins.
eco_ds = gdal.Open(ecoregion_fn)
eco_band = eco_ds.GetRasterBand(1)
eco_data = eco_band.ReadAsArray().flatten()
eco_bins = get_bins(eco_data)

# Read in the landcover data and get appropriate bins.
lc_ds = gdal.Open(landcover_fn)
lc_band = lc_ds.GetRasterBand(1)
lc_data = lc_band.ReadAsArray().flatten()
lc_bins = get_bins(lc_data)

# Calculate the histogram.
hist, eco_bins2, lc_bins2, bn = \
    scipy.stats.binned_statistic_2d(
        eco_data, lc_data, lc_data, 'count',
        [eco_bins, lc_bins])

# Add bin information to the histogram so the output file
# is more useful.
hist = np.insert(hist, 0, lc_bins[:-1], 0)
row_labels = np.insert(eco_bins[:-1], 0, 0)
hist = np.insert(hist, 0, row_labels, 1)

# Save the output
np.savetxt(out_fn, hist, fmt='%1.0f', delimiter=',')

max_stars_repo_name: hohsieh/osgeopy-code

max_issues_repo_name: sthagen/osgeopy-code

max_forks_repo_name: sthagen/osgeopy-code

reason: import numpy,import scipy


text:
from PIL import Image
import numpy as np
import copy, math

def pil_loader(path):
  # open path as file to avoid ResourceWarning (https://github.com/python-pillow/Pillow/issues/835)
  with open(path, 'rb') as f:
    with Image.open(f) as img:
      return img.convert('RGB')

ext: py

reason: import numpy


text:
# -*- coding: utf-8 -*-
"""Progetto_Data_Intensive.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1d5-w-eV5cf4S4Qf0MmpTzamtBUvwaTeK

<font size="5">Programmazione di Applicazioni Data Intensive 2020/2021 - Progetto

- Romagnoli Giacomo

- Parrinello Angelo

# Introduzione al problema

Ci troviamo sempre più spesso nella condizione di affidarci al giudizio di esperti per decretare la qualità di un vino. Fortunatamente il consumatore, scarsamente dotato in gusto e olfatto, può affidarsi a metri di giudizio scientifici per rendersi autonomo nella scelta.

Il progetto ha come **obbiettivo** quello di definire la qualità (buona o scadente) di un vino in base alle sue caratteristiche fisiche.

In materia il problema è definito di **classificazione**.

#Caricamento Librerie

Per prima cosa carichiamo le librerie per effettuare operazioni sui dati

*   **Numpy** per creare e operare su array a N dimensioni
*   **Pandas** per caricare e manipolare dati tabulari
*   **Matplotlib** per creare grafici

Importiamo le librerie usando i loro alias convenzionali e abilitando l'inserimento dei grafici inline.
"""

# Commented out IPython magic to ensure Python compatibility.
import csv
import numpy as np
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt

# %matplotlib inline

"""#Descrizione del problema

- Carichiamo dalla repo di Github il dataset contenente la matrice dei vini con le loro caratteristiche.
"""

import os

wine_zip = "wine.zip"

if not os.path.exists(wine_zip):
    from urllib.request import urlretrieve
    urlretrieve("https://github.com/AngeloParrinello/Wine-Quality-Classification/blob/main/wine.zip?raw=true", wine_zip)

from zipfile import ZipFile

!unzip "wine.zip"

"""- Usiamo la funzione `read_csv` per importare i dati contenuti nel file "wine.csv", specificando il separatore."""

data = pd.read_csv("wine.csv", delimiter=";")

"""- Analizziamo ora i dati, partendo dalla dimensione della matrice."""

data.shape

"""- Le colonne, ovvero le features, da cui è composta."""

data.columns

"""- Ed una piccola sezione della matrice appena importata: in questo caso visualizziamo 5 dati partendo dal fondo."""

data.tail()

"""Come si nota, il dataset contiene 6440 righe e 12 colonne.
Ogni riga costituisce un *vino registrato*. Questi dati sono stati raccolti da un'azienda vinicola portoghese.
Di ogni vino abbiamo le seguenti *feature fisiche*:

- **fixed acidity** : *acidità fissa*. Misurato in g/L.
- **volatile acidity** : *acidità volatile*. Misurato in g/L.
- **citric acid** : *acido citrico*. Misurato in g/L.
- **residual sugar** : *zucchero residuo*. Misurato in g/L.
- **chlorides** : *cloruri*. Misurato in g/L.
- **free sulfur dioxide** : *anidride solforosa libera*. Misurata mg/L.
- **total sulfur dioxide** : *anidride solforosa totale*. Misurata mg/L.
- **density** : *densità*. Misurata g/mL.
- **pH** : *pH*.
- **sulphates** : *solfati*. Misurata in g/L.
- **alcohol** : *alcool*. Misurata in percentuale, sul totale della soluzione.

Utilizziamo il metodo *info* per vedere la memoria occupata dal DataFrame e i tipi delle colonne assegnati da Pandas, al fine di conoscere il tipo delle variabili con cui lavoreremo.
"""

data.info()

"""Per le finalità del progetto la colonna **Quality** fornisce informazioni eccessive. E' quindi necessario operare una trasformazione su quest'ultima. Nello specifico sostituiremo il valore intero attuale con un booleano: se la qualità è maggiore o uguale a 6 allora verrà etichettato come "Good" altrimenti come "Bad"; assumendo due valori, dunque, possiamo dire che è una variabile categorica e possiamo convertirla nel tipo appropriato risparmiando memoria, anche se è un dataset di piccole dimensioni. Anche la variabile *alcohol* è di tipo object e anch'essa va trasformata, questa volta in float.

Andiamo ora a mappare la variabile da predire ovvero *quality*.

- Creiamo quindi una banale funzione che adempi questo scopo.
"""

def convert_value(i):
  if i >= 6:
    return 'Good'
  else:
    return 'Bad'

"""- Attraverso la funzione *apply* di Pandas, applichiamo la funzione appena creata alla colonna **Quality**."""

data.quality = data.quality.apply(convert_value)

data['quality'].unique()

data.tail()

"""- Una volta mappata *quality*, rendiamola del corretto tipo. Così come *alcohol*."""

data = data.astype({'quality': 'category', 'alcohol': 'float64'})

"""- Tutti gli altri attributi sono variabili continue di tipo float e notiamo anche che non abbiamo valori nulli."""

data.info()

"""#Analisi Esplorativa dei Dati

Iniziamo ora *l'analisi esplorativa dei dati* ove questi vengono analizzati, esplorati e si apprendono le prime informazioni generiche a riguardo, che guideranno i passi successivi. In questa fase cercheremo anche eventuali correlazioni tra essi e risolveremo imprecisioni del dataset.

- Utilizziamo la funzione *describe*() per iniziare a comprendere come sono distribuiti i dati.

- Per evitare disambiguità al lettore precisiamo che la riga con scritto 25% (25-esimo), 50% (50-esimo) e 75% (75-esimo) stanno ad indicare il *percentile*. In breve, considerando la prima colonna e la riga del 25-esimo percentile, 6.4 è il valore sotto il quale si trovano il 25% di osservazioni.
"""

data.describe()

"""- Per analizzare la distribuzione dei valori nella colonna nominale *quality*, possiamo utilizzare *value_counts*, specificando anche *normalize=True*, se vogliamo visualizzare la percentuale."""

data["quality"].value_counts()

data["quality"].value_counts(normalize=True)

data["quality"].value_counts().plot.pie();

"""- Dai grafici possiamo notare che le classi sono **sbilanciate**, cioè sono presenti molti più vini buoni rispetto a quelli cattivi. Questa peculiarità del dataset potrebbe creare problemi in fase di learning del modello: andranno bilanciate le classi.

- Visualizziamo la distribuzione dei valori delle feature del nostro dataset in degli **istogramma**.
"""

import warnings
warnings.filterwarnings("ignore")

plt.figure(figsize=(30, 20));
for n, value in enumerate(["fixed acidity", "volatile acidity", "citric acid", "residual sugar",
                           "chlorides", "free sulfur dioxide", "total sulfur dioxide", "density", "pH", "sulphates", "alcohol"], start=1):
  data[value].plot.hist(ax=plt.subplot(4,3,n), bins=30);
  plt.subplot(4, 3, n).set_title(value);

"""- Per maggiore dettaglio gli istogrammi, ma anche altri grafici in questo notebook, saranno ripetuti per ogni feature del dataset.

- Possiamo notare, ad esempio, che la maggior parte di vini ha una densità compresa tra 0 e 20 g/mL mentre molto pochi sono estremamente densi.
"""

#percentuale dei vini estremamente densi
(len(data[data.density > 20]) / len(data.density)) * 100

"""- Possiamo anche farci un'idea più numerica della distribuzione di queste feature. Quindi eseguiamo il **binning o discretizzazione** delle colonne numeriche, utilizzando la funzione *cut*, suddividendo così i valori in fasce di uguale ampiezza per contarne poi i valori in ciascuna. """

datino = data.columns.drop("quality")
for value in datino:
  print(pd.cut(data[value],4).value_counts()) #4 fasce uguali
  print("\n")

"""- E' facile notare come molte proprietà dei nostri vini hanno divisione dei valori simili tra loro: la distribuzione non è molto omogenea.

- Visualizziamo ora le statistiche di base sui valori delle feature mediante *boxplot*.
"""

data.info()

np.warnings.filterwarnings('ignore', category=np.VisibleDeprecationWarning)
plt.figure(figsize=(16, 35))
data.boxplot(ax=plt.subplot(6,2,1), column="fixed acidity", by="quality", showmeans=True);
data.boxplot(ax=plt.subplot(6,2,2), column="volatile acidity", by="quality", showmeans=True);
data.boxplot(ax=plt.subplot(6,2,3), column="citric acid", by="quality", showmeans=True);
data.boxplot(ax=plt.subplot(6,2,4), column="residual sugar", by="quality", showmeans=True);
data.boxplot(ax=plt.subplot(6,2,5), column="chlorides", by="quality", showmeans=True);
data.boxplot(ax=plt.subplot(6,2,6), column="free sulfur dioxide", by="quality", showmeans=True);
data.boxplot(ax=plt.subplot(6,2,7), column="total sulfur dioxide", by="quality", showmeans=True);
data.boxplot(ax=plt.subplot(6,2,8), column="density", by="quality", showmeans=True);
data.boxplot(ax=plt.subplot(6,2,9), column="pH", by="quality", showmeans=True);
data.boxplot(ax=plt.subplot(6,2,10), column="sulphates", by="quality", showmeans=True);
data.boxplot(ax=plt.subplot(6,2,11), column="alcohol", by="quality", showmeans=True);

"""- **ALCOHOL**: nei vini di scarsa qualità la percentuale media di alcool è più bassa rispetto ai vini di alta qualità, che però hanno una maggiore distribuzione.
- **SULPHATES**: fra buoni e cattivi vini non notiamo una grosse differenze riguardo la media, unica discriminante la maggior presenza di outliners nei vini cattivi spostati più verso il massimo.
- **pH**: il pH nei vini di maggior pregio è leggermente più alto, quindi un vino neutro ma tendente al basico, rispetto ai vini peggiori.
- **DENSITY**: il boxplot relativo alla densità non ci è molto di aiuto: infatti la stragrande maggioranza dei vini ha una bassissima densità, come già avevamo notato.
- **TOTAL SULFUR DIOXIDE**: l'anidride solforosa totale è tendezialmente più bassa nei vini "good". Come si vede, infatti, nei vini di scarsa qualità arriviamo a toccare i 400 mg/L.
- **FREE SULFUR DIOXIDE**: l'anidride solforosa libera è simile in entrambe le categorie di vini. Notare il picco di un dato ben oltre il baffo superiore del boxplot.
- **CHLORIDES**: osservando il grafico relativo ai cloridi notiamo che ambo le tipologie di vini che stiamo analizzando hanno diversi outliers, ma in media i cloridi sono più bassi nei vini buoni.
- **RESIDUAL SUGAR**: lo zucchero residuo è simile in entrambi i tipi di vino.
- **CITRIC ACID**: l'acido citrico è simile in entrambi i tipi di vino.
- **VOLATILE ACIDITY**: l'acidità volatile è sensibilmente inferiore nel massimo per i vini buoni.
- **FIXED ACIDITY**: l'acidità fissa è simile in entrambi i tipi di vino.

- Tra le features presentate quelle che ipotiziamo siano più rilevanti sono: l'alcool, l'anidride solforosa totale e l'acidità volatile. Questo perchè notiamo sostanziali differenze nelle due categorie di vino.

- Con l'operazione di *pivoting* possiamo sdoppiare ciscuna colonna del dataframe, suddividendo i valori relativi alle classi good e bad, al fine di visualizzare meglio quanto le variabili predittive siano correlate con la classe da predire.
"""

data.pivot(columns="quality")

"""- Mostriamo ora i grafici relativi ai valori più incisivi del dataframe pivotato."""

data.pivot(columns="quality")["alcohol"].plot.hist(bins=20, stacked=True);

data.pivot(columns="quality")["total sulfur dioxide"].plot.hist(bins=20, stacked=True);

"""- In enologia, l'anidride solforosa è utilizzata sin dalle primissime fasi della produzione del vino, a partire dal mosto fino all'imbottigliamento. Nell'usare l'anidride solforosa, è opportuno sapere che una parte di questo gas si combina con alcuni componenti del mosto o del vino, mentre la restante parte resta libera, cioè non combinata. Sarà proprio la parte libera a svolgere gli importanti effetti antiossidanti e antisettici: per questo motivo è indispensabile che l'anidride solforosa si combini il meno possibile. L'anidride solforosa combinata è comunque utile, poiché nel caso in cui la frazione libera si disperda - durante le operazioni di travaso, per esempio - una piccola parte di quella combinata si libera sostituendola. (fonte: https://www.agraria.org/viticoltura-enologia/anidride-solforosa.htm)."""

data.pivot(columns="quality")["volatile acidity"].plot.hist(bins=20, stacked=True);

"""- Si definisce acidità volatile in un vino, la quantità di acido acetico presente in un vino e viene espressa in g/l di acido acetico. E’ importante determinare l’acidità volatile del vino perché la quantità di acido acetico presente è indice dello stato di sanità dell’uva, di come procede la fermentazione e poi, in generale, dello stato di conservazione del vino. (fonte: https://www.rivistadiagraria.org/articoli/anno-2009/determinazione-dellacidita-volatile-di-un-vino/)

- Possiamo visualizzare i valori medi di tutte le feature, distinti per qualità, raggruppandoli per quest'ultimo.
"""

data_by_quality = data.groupby("quality")
data_by_quality.mean()

"""- Per valutare visivamente la correlazione tra due variabili, possiamo selezionarne alcune ed evidenziare anche le due diverse classi nel *diagramma a dispersione*, differenziando i punti ad esempio per colore mediante un dizionario che associ green ai vini buoni e rosso ai vini cattivi."""

quality_color_map = {"Good": "green", "Bad": "red"}
quality_colors = data["quality"].map(quality_color_map)

plt.figure(figsize=(16, 20))
data.plot.scatter(ax=plt.subplot(4,2,1), x="residual sugar", y="chlorides", c=quality_colors);
data.plot.scatter(ax=plt.subplot(4,2,2), x="chlorides", y="free sulfur dioxide", c=quality_colors);
data.plot.scatter(ax=plt.subplot(4,2,3), x="free sulfur dioxide", y="total sulfur dioxide", c=quality_colors);
data.plot.scatter(ax=plt.subplot(4,2,4), x="total sulfur dioxide", y="pH", c=quality_colors);
data.plot.scatter(ax=plt.subplot(4,2,5), x="pH", y="alcohol", c=quality_colors);
data.plot.scatter(ax=plt.subplot(4,2,6), x="citric acid", y="residual sugar", c=quality_colors);

chlorides = data["chlorides"]
residual_sugar = data["residual sugar"]
free_sulfur_dioxide = data["free sulfur dioxide"]
total_sulfur_dioxide = data["total sulfur dioxide"]
citric_acid = data["citric acid"]

"""- Negli terzo grafico notiamo che c'è una dipendenza fra le variabili. A conferma, calcoliamo il *coefficiente di correlazione di Pearson*."""

np.mean((free_sulfur_dioxide-free_sulfur_dioxide.mean()) * (total_sulfur_dioxide-total_sulfur_dioxide.mean())) / (free_sulfur_dioxide.std() * total_sulfur_dioxide.std())

"""- Il valore del coefficiente di Pearson è vicino a 1 (correlazione diretta): al raddoppiare di una, raddoppia anche l'altra.

- Mentre le altre variabili sono scorrelate tra loro, questo è deducibile visivamente guardando i grafi. Per dimostrarlo, calcoliamo alcuni coefficenti di Pearson di queste variabili.
"""

np.mean((chlorides-chlorides.mean()) * (residual_sugar-residual_sugar.mean())) / (chlorides.std() * residual_sugar.std())

np.mean((free_sulfur_dioxide-free_sulfur_dioxide.mean()) * (chlorides-chlorides.mean())) / (free_sulfur_dioxide.std() * free_sulfur_dioxide.std())

np.mean((citric_acid-citric_acid.mean()) * (residual_sugar-residual_sugar.mean())) / (citric_acid.std() * residual_sugar.std())

"""#Normalizzazione dei dati

Procediamo ora con la normalizzazione/standardizzazione delle variabili.

In generale le variabili coinvolte in un modello di regressione
possono utilizzare scale di valori molto diverse. Questo rende difficile l'interpretazione dei modelli che andremo ad addestrare. Una soluzione più semplice è però normalizzare i dati in modo che tutte le variabili abbiano valori in un medesimo intervallo.

- Nel nostro caso abbiamo necessità di normalizzare i dati come si nota dall'esempio sotto.
"""

data[["total sulfur dioxide", "chlorides"]]

"""- Effettueremo la standardizzazione dopo aver applicato il metodo *Hold-Out* per suddividere i dati in training e validation set in preparazione all'addestramento dei modelli.

- Prima però proviamo ad addestrare un modello d'esempio senza standardizzare le feature.
"""

from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.linear_model import Lasso
from sklearn.linear_model import Perceptron
from sklearn.model_selection import train_test_split

y = data["quality"] #variabile da predire
X = data.drop(columns=['quality']) #variabili predittive

"""- Applichiamo il metodo *Hold-Out* per partizionare casualmente i dati, riservandone i 2/3 al training set e il restante al validation set. Per partizionare definiremo il parametro *random_state* (ovvero un seed per generare divisioni casuali del set) pari a 42, un numero convenzionale.

- Evitiamo per ora di utilizzare una *pipeline*, per analizzare più agevolmente il modello.
"""

X_train, X_val, y_train, y_val = \
    train_test_split(X, y, test_size=1/3, random_state=42)

"""- Addestriamo un modello **Perceptron** d'esempio."""

model = Perceptron(random_state=42)
model.fit(X_train, y_train)
model.score(X_val, y_val)

"""- Il nostro modello, che utilizza variabili non standardizzate, produce scarsi risultati con un tasso di accuratezza che arriva a 37,1%. Ciò accade perchè le variabili presenti presentano ordini di grandezza troppo distanti.

- Proviamo adesso a **standardizzare** le feature e vedere come ciò influisca sui nostri risultati:
"""

scaler = StandardScaler()
Xn_train = scaler.fit_transform(X_train)
Xn_val = scaler.transform(X_val)

"""- Addestriamo nuovamente un modello **Perceptron**."""

model = Perceptron(random_state=42)
model.fit(Xn_train, y_train)

"""- Una volta addestrato il modello, possiamo trovare i valori dei pesi $\mathbf{w}$ e del bias $b$ rispettivamente negli attributi `coef_[0]` e `intercept_[0]`"""

model.coef_[0]

"""- Possiamo anche capire meglio quale variabile abbia piú peso rispetto ad un'altra."""

pd.Series(model.coef_[0], index=X_train.columns)

model.intercept_[0]

"""- Valutiamo nuovamente la bontá del nostro modello:"""

model.score(Xn_val, y_val)

"""- Dopo aver normalizzato le variabili, il nostro modello produce ottimi risultati con un accuratezza del 70%. Calcoliamo ora il numero di istanze corrette e quelle errate."""

correct_class = model.predict(Xn_val) == y_val
correct_class.value_counts()

"""- Il numero di istanze che il nostro modello interpreta correttamente é di 1505 contro le 642 errate.

#Regolarizzazione dei dati ed indagine delle feature piú rilveanti

Per individuare le feature meno importanti utilizziamo la regolarizzazione **LASSO o L1** per poi andarle eventualmente ad eliminarle.
"""

model_l1 = Perceptron(penalty='l1', alpha=0.0001, random_state=42)
model_l1.fit(Xn_train, y_train)

pd.Series(model_l1.coef_[0], index=X_train.columns)

model_l1.score(Xn_val, y_val)

"""- Come si puó notare dalla cella qui sopra, applicando una regolarizzazione LASSO il modello perde precisione. Come ulteriore dimostrazione creiamo un modello Princeptron utilizzando peró solo le feature non annullate."""

X = data[["volatile acidity", "citric acid", "chlorides", "free sulfur dioxide", "density", "pH", "sulphates", "alcohol"]]

X8_train, X8_val, y8_train, y8_val = \
    train_test_split(X, y, test_size=1/3, random_state=42)

scaler = StandardScaler()
X8n_train = scaler.fit_transform(X8_train)
X8n_val = scaler.transform(X8_val)

model_8f = Perceptron(random_state=42)
model_8f.fit(X8n_train, y8_train)

model_8f.score(X8n_val, y8_val)

"""- E' evidente quindi che per i nostri scopi andranno considerate tutte le feature. Sará comunque importante apportare una regolarizzazione a tutte le nostre feature.

Applichiamo ora la regolarizzazione **RIDGE O L2** al nostro modello d'esempio.
"""

model_l2 = Perceptron(penalty='l2', alpha=0.0001, random_state=42)
model_l2.fit(Xn_train, y_train)

pd.Series(model_l2.coef_[0], index=X_train.columns)

model_l2.score(Xn_val, y_val)

"""- Come si vede, diminuendo il peso delle feature lo score, ovvero l'accuratezza, del modello aumenta.

- E' peró risaputo che la suddivisone di istanze tra classi sbilanciate porti a molti errori di classificazione e quindi alla generazione di un modello inadeguato: il nostro dataset soffre di questo problema. Bisogna quindi sottolineare l'inaffidabilitá degli score di accuratezza ottenuti fin'ora.

- La matrice di confusione sottostante conferma quanto appena descritto: di 1379 vini buoni presenti nel validation set solo 815 sono state etichettate correttamente.
"""

from sklearn.metrics import confusion_matrix
import seaborn as sns

plt.figure(figsize=(8, 6))
plt.title('Matrice di confusione sul modello con L1', size=16)
sns.heatmap(pd.DataFrame(confusion_matrix(y_val, model_l1.predict(Xn_val)), columns=model_l1.classes_, index = model_l1.classes_), annot=True, cmap='Blues', fmt='d');

"""#Generazione di modelli di learning

##Risoluzione problema classi sbilanciate

Come dimostrato nella fase precedente c'é un forte sbilanciamento delle classi.
E' necessario quindi attuare un meccanismo di **oversampling**, ovvero aumentare il numero delle istanze nella classe meno rappresentata. Come spiegato a lezione, si utilizzerà la tecnica **SMOTE** o **Synthetic Minority Oversampling Technique** che fa parte del modulo imblearn.

- Procediamo quindi a fare l'oversampling del training set e del validation set. Come da documentazione si userà il metodo *fit_resample(X,y)*.
"""

from imblearn.over_sampling import SMOTE
import warnings
warnings.filterwarnings('ignore')

sm = SMOTE(random_state=42)

X_bal, y_bal = sm.fit_resample(X, y)
X_bal = pd.DataFrame(X_bal, columns=X.columns)
y_bal = pd.Series(y_bal, dtype="category", name="quality")

X_train_bal, y_train_bal = sm.fit_resample(X_train, y_train)
X_train_bal = pd.DataFrame(X_train_bal, columns=X_train.columns)
y_train_bal = pd.Series(y_train_bal, dtype="category", name="quality")

"""Prima del resampling:"""

pd.Series(y_train).value_counts()

pd.Series(y_train).value_counts().plot.pie(autopct='%1.1f%%',
        shadow=True, startangle=90);

"""Dopo il resampling:"""

y_train_bal.value_counts()

y_train_bal.value_counts().plot.pie(autopct='%1.1f%%',
        shadow=True, startangle=90);

"""- Le classi sono ora bilanciate.

- Creiamo prima dei dizionari per salvarci alcuni dati sui modelli che andremo a creare:
"""

scores = {}
recall = {}
precision = {}
f1 = {}
models = {}
confusion_matrices = {}

"""## Perceptron

La fase di modellazione ha inizio provando un modello basato su Perceptron. Il _Perceptron_ è un algoritmo di apprendimento molto semplice, concettualmente simile alla discesa gradiente. Per trovare i migliori parametri possibili usiamo _GridSearch_: la grid search testerá tutte le combinazioni di iperparametri possibili. Per questo modello Perceptron cerchiamo di ottenere il massimo rendimento possibile dai seguenti iperparametri:
- standardizzazione o meno delle feature
- tipo di regolarizzazione del modello
- peso della regolarizzazione
- stima o meno dell'intercetta
"""

from sklearn.model_selection import KFold
kf = KFold(5, shuffle=True, random_state=42)

per_mod = Pipeline([
    ("scaler", StandardScaler()),
    ("per", Perceptron(n_jobs=-1, random_state=42)) #n_jobs=-1 per utilizzare tutti i processori
])

#In generale, tutti i parametri disponibili per l'ottimizzazione in GridSearchCV sono disponibili tramite
#per_mod.get_params().keys()

#ATTENZIONE! Nel dizionario per_grid (per cercare la migliore ottimizzazione
#dei vari parametri) la chiave di ogni possibile valore del parametro dovrá essere
#esattamente la stessa contenuta in per_mod.get_params().keys()

per_grid = {
    "scaler": [None, StandardScaler()],
    "per__penalty": ["l2", "l1", "elasticnet"],
    "per__alpha": [0.0001, 0.001, 0.01, 1, 10],
    "per__fit_intercept": [False, True]
}

from sklearn.model_selection import GridSearchCV

#se voglio dentro cv_results_ i training score --> return_train_score=True
per_gs = GridSearchCV(estimator=per_mod, param_grid=per_grid, cv=kf, n_jobs=-1)

per_gs.fit(X_train_bal, y_train_bal)

print('CV score: {:.5f}'.format(per_gs.score(X_val, y_val)))
print('\n')
print('Best parameter:',per_gs.best_params_)
print('\n')
print('Best estimator:',per_gs.best_estimator_)

"""Oltre all'accuratezza come percentuale di classificazioni corrette, é importante tenere conto anche di altri parametri che indicano quanto il nostro modello sia preciso.

Oltre alla matrice di confusione vista prima ora utilizzeremo come metro di giudizio per la precisione del modello anche:
- F1 Score, misura unica della performance di un modello è la **_F1-measure_**, ovvero la media armonica tra precision e recall.
- Precision Score, la **_precision_** indica la percentuale di esempi classificati come "Good" che sono realmente tali.
- Recall score, la **_recall_** indica la percentuale di esempi realmente di classe "Bad" che sono stati rilevati essere tali dal modello.
"""

from sklearn.metrics import precision_score, recall_score, f1_score

#col parametro average si indica se si vuole la media(macro) o un vettore (None) o altro
#col parametro pos_label si indica la classe di riferimento

print('F1 score: {:.5f}'.format(f1_score(y_val, per_gs.predict(X_val), pos_label='Good')))
print('\n')
print('Precision score: {:.5f}'.format(precision_score(y_val, per_gs.predict(X_val), pos_label='Good')))
print('\n')
print('Recall score: {:.5f}'.format(recall_score(y_val, per_gs.predict(X_val), pos_label='Good')))
print('\n')
plt.figure(figsize=(8, 6))
plt.title('Matrice di confusione sul modello con Perceptron', size=16)
sns.heatmap(pd.DataFrame(confusion_matrix(y_val, per_gs.predict(X_val)), columns=per_gs.classes_, index = per_gs.classes_),
            annot=True, cmap='Oranges', fmt='d');

"""Salviamo le informazioni nei vari dizionari, che più tardi ci serviranno."""

scores["per"] = per_gs.score(X_val, y_val)
recall["per"] = recall_score(y_val, per_gs.predict(X_val), pos_label='Good')
precision["per"] = precision_score(y_val, per_gs.predict(X_val), pos_label='Good')
f1["per"] = f1_score(y_val, per_gs.predict(X_val), pos_label='Good')
models["per"] = per_gs.best_estimator_
confusion_matrices["per"] = pd.DataFrame(confusion_matrix(y_val, per_gs.predict(X_val)), columns=per_gs.classes_, index = per_gs.classes_)

"""Dalla matrice di confusione e dagli score si puó notare come il modello appare ancora piuttosto sbilanciato e leggermente preciso.

Testiamo ora un modello basato sulla **Regressione Logistica**.

## Regressione Logistica

La fase di modellazione prosegue testando un modello basato sulla Regresione Logistica. La  _Regressione Logistica_ è un  un modello di classificazione binaria basato sulla regressione lineare. Per trovare i migliori parametri possibili usiamo nuovamente _GridSearch_.

Per questo modello cerchiamo i seguenti migliori iperparametri in due casi:\
*Il primo, con regolarizzazione L2 e L1 senza specificare l1_ratio:*
- standardizzazione o meno delle feature
- tipo di regolarizzazione del modello
- peso della regolarizzazione
- stima o meno dell'intercetta

*Il secondo, con elasticnet specificando i valori possibili di l1_ratio:*
- standardizzazione o meno delle feature
- tipo di regolarizzazione del modello
- peso della regolarizzazione
- peso l1_ratio
- stima o meno dell'intercetta

P.S: in Python il secondo iperparametro che la Regressione Elastic Net introduce (α) si identifica con _l1_ratio_
"""

from sklearn.linear_model import LogisticRegression
#random_state indica il seed per la casualità usato nell'addestramento
#solver indica quale implementazione utilizzare tra per l'addestramento, "saga"
#è l'unica che supporta tutte le opzioni per la regolarizzazione

log_mod = Pipeline([
    ("scaler", StandardScaler()),
    ("log", LogisticRegression(solver='saga', random_state=42))
])

#np.logspace come in matlab, estremi e poi numero di passi tra gli estremi
log_grid = [
    {
        "scaler": [None, StandardScaler()],
        "log__penalty": ["l2", "l1"],
        "log__C": np.logspace(-3, 3, num=10),
        "log__fit_intercept": [False, True]
    },
    {
        "scaler": [None, StandardScaler()],
        "log__penalty": ["elasticnet"],
        "log__C": np.logspace(-3, 3, num=10),
        "log__l1_ratio": np.logspace(0,1,num=10),
        "log__fit_intercept": [False, True]
    }
]

log_gs = GridSearchCV(estimator=log_mod,param_grid=log_grid, cv=kf, n_jobs=-1)

log_gs.fit(X_train_bal, y_train_bal)

print('CV score: {:.5f}'.format(log_gs.score(X_val, y_val)))
print('\n')
print('Best parameter:',log_gs.best_params_)
print('\n')
print('Best estimator:',log_gs.best_estimator_)

print('F1 score: {:.5f}'.format(f1_score(y_val, log_gs.predict(X_val), pos_label='Good')))
print('\n')
print('Precision score: {:.5f}'.format(precision_score(y_val, log_gs.predict(X_val), pos_label='Good')))
print('\n')
print('Recall score: {:.5f}'.format(recall_score(y_val, log_gs.predict(X_val), pos_label='Good')))
print('\n')
plt.figure(figsize=(8, 6))
plt.title('Matrice di confusione sul modello con Regressione Logistica', size=16)
sns.heatmap(pd.DataFrame(confusion_matrix(y_val, log_gs.predict(X_val)), columns=log_gs.classes_, index = log_gs.classes_),
            annot=True, cmap='Reds', fmt='d');

scores["log"] = log_gs.score(X_val, y_val);
recall["log"] = recall_score(y_val, log_gs.predict(X_val), pos_label='Good')
precision["log"] = precision_score(y_val, log_gs.predict(X_val), pos_label='Good')
f1["log"] = f1_score(y_val, log_gs.predict(X_val), pos_label='Good')
models["log"] = log_gs.best_estimator_
confusion_matrices["log"] = pd.DataFrame(confusion_matrix(y_val, log_gs.predict(X_val)), columns=log_gs.classes_, index = log_gs.classes_)

"""Come si vede sopra il modello appena generato é piú preciso sotto tutti i punti di vista rispetto al classico Perceptron. Rimane peró leggermente alta la classificazione di falsi positivi. La modellazione prosegue con un modello di tipo **SVM** o **Support Vector Machines**.

## Support Vector Machines

Continuiamo i nostri test utilizzando ora un modello basato sul Support Vector Machines o SVM.  Il  _SVM o Support Vector Machines_ è un  un modello, definito in origine per problemi di classificazione binaria basato, con l'obiettivo di inviduare la separazione lineare ottimale tra le istanze delle due classi. Per trovare i migliori parametri possibili usiamo nuovamente _GridSearch_.

Per questo modello cerchiamo i seguenti migliori iperparametri in due casi:\
*Il primo, che utilizza la funzione Kernel Gaussiana RBF:*
- tipologia funzione Kernel
- peso delle classi
- peso della regolarizzazione
- peso del coefficente Kernel gamma (γ)

*Il secondo, SVC lineare:*
- tipologia funzione Kernel
- peso delle classi
- peso della regolarizzazione

Nota: é stato tentato anche il caso con la Funzione Kernel Polinomiale ma con scarsi risultati ed eccessivi tempi di addestramento. A tal proposito non verrá riportato in questo documento.
"""

from sklearn.svm import SVC

svc_mod = Pipeline([
    ("scaler", StandardScaler()),
    ("svc", SVC(random_state=42))
])

#svc_mod.get_params().keys()
#gamma si poteva impostare in maniera automatica con gamma='scale' quando si crea l'SVC

svc_grid =[ {
    "svc__kernel": ["rbf"],
    "svc__class_weight": [None, "balanced"],
    "svc__C": [0.01, 0.1, 1, 10,100],
    'svc__gamma': [0.01, 0.1, 1, 10,100]
}, {
    "svc__kernel": ["linear"],
    "svc__class_weight": [None, "balanced"],
    "svc__C": [0.01, 0.1, 1, 10,100],
}
]

svc_gs = GridSearchCV(estimator=svc_mod,param_grid=svc_grid, cv=kf, n_jobs=-1)

svc_gs.fit(X_train_bal, y_train_bal)

print('CV score: {:.5f}'.format(svc_gs.score(X_val, y_val)))
print('\n')
print('Best parameter:',svc_gs.best_params_)
print('\n')
print('Best estimator:',svc_gs.best_estimator_)

print('F1 score: {:.5f}'.format(f1_score(y_val, svc_gs.predict(X_val), pos_label='Good')))
print('\n')
print('Precision score: {:.5f}'.format(precision_score(y_val, svc_gs.predict(X_val), pos_label='Good')))
print('\n')
print('Recall score: {:.5f}'.format(recall_score(y_val, svc_gs.predict(X_val), pos_label='Good')))
print('\n')
plt.figure(figsize=(8, 6))
plt.title('Matrice di confusione sul modello con SVC', size=16)
sns.heatmap(pd.DataFrame(confusion_matrix(y_val, svc_gs.predict(X_val)), columns=svc_gs.classes_, index = svc_gs.classes_),
            annot=True, cmap='Greens', fmt='d');

scores["svc"] = svc_gs.score(X_val, y_val);
recall["svc"] = recall_score(y_val, svc_gs.predict(X_val), pos_label='Good')
precision["svc"] = precision_score(y_val, svc_gs.predict(X_val), pos_label='Good')
f1["svc"] = f1_score(y_val, svc_gs.predict(X_val), pos_label='Good')
models["svc"] = svc_gs.best_estimator_
confusion_matrices["svc"] = pd.DataFrame(confusion_matrix(y_val, svc_gs.predict(X_val)), columns=svc_gs.classes_, index = svc_gs.classes_)

"""## Albero Decisionale

Il prossimo modello che vedremo é generato attraverso un Albero Decisionale o Decision Tree. I modelli di classificazione visti finora si basano su iperpiani descritti da equazioni (lineari o non) su tutte le variabili. L' **Albero Decisionale** o **Decision Tree** costituisce un approccio differente: la classificazione avviene in base ad una serie di decisioni "semplici", basate ciascuna su una sola variabile. L'obiettivo é quindi:
- creare una struttura ad albero dove ogni nodo intermedio contiene un predicato con una variabile di input
- il predicato divide le istanze di input in 2 sottoinsiemi a cui corrispondono 2 nodi figli con relativi sotto alberi e così via ricorsivamente
"""

from sklearn.tree import DecisionTreeClassifier

features_num = X.columns.size

tree_mod = Pipeline([
    ("scaler", StandardScaler()),
    ("tree", DecisionTreeClassifier(class_weight="balanced", random_state=42))
])

#print(tree_mod.get_params().keys())

tree_grid = {'scaler': [None, StandardScaler()],
             'tree__min_samples_split': range(2, 6),
             'tree__min_samples_leaf': range(1, 6),
             'tree__max_depth': range(2,6),
             'tree__max_features': range(2, features_num)}

tree_gs = GridSearchCV(estimator=tree_mod,param_grid=tree_grid, cv=kf, n_jobs=-1)

tree_gs.fit(X_train_bal, y_train_bal)

print('CV score: {:.5f}'.format(tree_gs.score(X_val, y_val)))
print('\n')
print('Best parameter:',tree_gs.best_params_)
print('\n')
print('Best estimator:',tree_gs.best_estimator_)

print('F1 score: {:.5f}'.format(f1_score(y_val, tree_gs.predict(X_val), pos_label='Good')))
print('\n')
print('Precision score: {:.5f}'.format(precision_score(y_val, tree_gs.predict(X_val), pos_label='Good')))
print('\n')
print('Recall score: {:.5f}'.format(recall_score(y_val, tree_gs.predict(X_val), pos_label='Good')))
print('\n')
plt.figure(figsize=(8, 6))
plt.title('Matrice di confusione sul modello con Albero Decisionale', size=16)
sns.heatmap(pd.DataFrame(confusion_matrix(y_val, tree_gs.predict(X_val)), columns=tree_gs.classes_, index = tree_gs.classes_),
            annot=True, cmap='Purples', fmt='d');

"""Per visualizzare maggiormente il modello possiamo stampare _testualmente_ e _graficamente_ il tutto. Per una maggior comprensione la profonditá é stata limitata a 2."""

from sklearn.tree import export_text, plot_tree
features_name = ['fixed acidity', 'volatile acidity', 'citric acid', 'residual sugar',
       'chlorides', 'free sulfur dioxide', 'total sulfur dioxide', 'density',
       'pH', 'sulphates', 'alcohol']
print(export_text(tree_gs.best_estimator_[1], feature_names=features_name, max_depth=3))

plt.figure(figsize=(18, 12))
plot_tree(tree_gs.best_estimator_[1], feature_names=features_name,  max_depth=2, filled=True);

"""La rappresentazione testuale mostra come il modello classifichi ciascuna richiesta:
- se l'alcohol non é maggiore di 10,79 allora considera l'aciditá volatile
  - se non superiore a 0,20 allora considera il diossido di solfuro libero altrimenti ricontrolla l'aciditá volatile.
- E cosí via fino a quando non classifica la richiesta in una maniera o in un'altra.
"""

scores["tree"] = tree_gs.score(X_val, y_val);
recall["tree"] = recall_score(y_val, tree_gs.predict(X_val), pos_label='Good')
precision["tree"] = precision_score(y_val, tree_gs.predict(X_val), pos_label='Good')
f1["tree"] = f1_score(y_val, tree_gs.predict(X_val), pos_label='Good')
models["tree"] = tree_gs.best_estimator_
confusion_matrices["tree"] = pd.DataFrame(confusion_matrix(y_val, tree_gs.predict(X_val)), columns=tree_gs.classes_, index = tree_gs.classes_)

"""## Generazione Modelli Opzionali

Per interesse personale sono stati generati alcuni modelli solamente accennati o brevemente studiati durante il corso di Data Intensive. Per tali ragioni non verrano presi in considerazione nella decisione per il miglior modello generato.

Alcune difficoltà incontrate nello sviluppo di questi modelli sono state superate informandoci anche online tra cui dalla seguente fonte https://machinelearningmastery.com/ .

### KNeighborsClassifier

*KNN* è un algoritmo di apprendimento automatico supervisionato che può essere utilizzato per risolvere problemi di classificazione e regressione. Il principio di KNN è che il valore o la classe di un dato è determinato dai dati attorno a questo valore.

L'algoritmo di previsione calcola la distanza dal punto sconosciuto x, a tutti i punti in input. I punti vengono quindi ordinati aumentando la distanza da x. La previsione viene effettuata prevedendo l'etichetta di maggioranza dai punti più vicini "K".

La scelta di una K influenzerà la classe a cui verrà assegnato un nuovo punto.

P.S: Il valore passato a knc_n_neighbors rappresenta il valore K.
"""

from sklearn.neighbors import KNeighborsClassifier

knc_mod = Pipeline([
    ('scaler', StandardScaler()),
    ('knc', KNeighborsClassifier(n_jobs=-1))
])

knc_grid = {"scaler": [None, StandardScaler()],
            'knc__n_neighbors': range(1, 10, 1),
            'knc__weights': ['uniform', 'distance']}

knc_gs = GridSearchCV(estimator=knc_mod,param_grid=knc_grid, cv=kf, n_jobs=-1)

knc_gs.fit(X_train_bal, y_train_bal)

print('Best CV score: {:.5f}'.format(knc_gs.best_score_))
print('\n')
print('Best parameter:',knc_gs.best_params_)
print('\n')
print('Best estimator:',knc_gs.best_estimator_)

print('F1 score: {:.5f}'.format(f1_score(y_val, knc_gs.predict(X_val), pos_label='Good')))
print('\n')
print('Precision score: {:.5f}'.format(precision_score(y_val, knc_gs.predict(X_val), pos_label='Good')))
print('\n')
print('Recall score: {:.5f}'.format(recall_score(y_val, knc_gs.predict(X_val), pos_label='Good')))
print('\n')
plt.figure(figsize=(8, 6))
plt.title('Matrice di confusione sul modello con KNeighborsClassifier', size=16)
sns.heatmap(pd.DataFrame(confusion_matrix(y_val, knc_gs.predict(X_val)), columns=knc_gs.classes_, index = knc_gs.classes_),
            annot=True, cmap='flare', fmt='d');

"""### XGBoost

L’algoritmo *XGBoost*, acronimo di eXtreme Gradient Boosting, è un’implementazione specifica del metodo Gradient Boosting che utilizza approssimazioni più accurate per trovare il miglior modello ad albero.

La formazione di un XGBoost è una procedura iterativa che calcola ad ogni passo la migliore suddivisione possibile per il k-esimo albero elencando tutte le possibili strutture ancora disponibili in quel punto del percorso. Per questo modello XGBoost cerchiamo di ottenere il massimo rendimento possibile dai seguenti iperparametri:
- xgb__eta o learning rate
- xgb__max_depth o numero massimo di nodi consentiti dalla radice alla foglia più lontana di un albero
- xgb__n_estimators o numero di alberi in un modello XGBoost
- xgb__alpha o il peso della regolarizzazione
"""

from xgboost import XGBClassifier

xgb_mod = Pipeline([
    ('scaler', StandardScaler()),
    ('xgb', XGBClassifier(n_jobs=8, random_state=42, objective='binary:logistic'))
])

xgb_grid = {
    'xgb__eta': [0.002, 0.1, 0.5],
    'xgb__max_depth': [6],
    'xgb__n_estimators': [150, 300],
    'xgb__alpha': [0.0001, 0.001]
}

xgb_gs = GridSearchCV(estimator=xgb_mod,param_grid=xgb_grid, cv=kf, n_jobs=-1)

xgb_gs.fit(X_train_bal, y_train_bal)

print('Best CV score: {:.5f}'.format(xgb_gs.best_score_))
print('\n')
print('Best parameter:',xgb_gs.best_params_)
print('\n')
print('Best estimator:',xgb_gs.best_estimator_)

print('F1 score: {:.5f}'.format(f1_score(y_val, xgb_gs.predict(X_val), pos_label='Good')))
print('\n')
print('Precision score: {:.5f}'.format(precision_score(y_val, xgb_gs.predict(X_val), pos_label='Good')))
print('\n')
print('Recall score: {:.5f}'.format(recall_score(y_val, xgb_gs.predict(X_val), pos_label='Good')))
print('\n')
plt.figure(figsize=(8, 6))
plt.title('Matrice di confusione sul modello con XGBoost', size=16)
sns.heatmap(pd.DataFrame(confusion_matrix(y_val, xgb_gs.predict(X_val)), columns=xgb_gs.classes_, index = xgb_gs.classes_),
            annot=True, cmap='YlOrBr', fmt='d');

"""### Random Forest

Ogni modello utilizzato dalla previsione del *Random Forest* è un albero decisionale.

Una Random Forest combina molti alberi decisionali in un unico modello. Individualmente, le previsioni fatte dai singoli alberi decisionali potrebbero non essere accurate, ma combinate insieme, le previsioni saranno in media più vicine al risultato. Per questo modello Random Forest cerchiamo di ottenere il massimo rendimento possibile dai seguenti iperparametri:
- rf__n_estimators o numero di alberi in un modello Random Forest ovvero nella foresta
- rf__max_depth o numero massimo di nodi consentiti dalla radice alla foglia più lontana di un albero
- rf__min_samples_leaf o numero minimo di campioni richiesti per trovarsi in un nodo foglia.
- rf__min_samples_split o numero minimo di campioni richiesti per dividere un nodo interno
"""

from sklearn.ensemble import RandomForestClassifier

rf_mod = Pipeline([
    ('scaler', StandardScaler()),
    ('rf', RandomForestClassifier(n_jobs=-1, random_state=42))
])

rf_grid = {
    'rf__n_estimators': [100, 200, 300],
    'rf__max_depth': [2, 4, 6, 8, 10],
    'rf__min_samples_leaf': [1, 2, 4],
    'rf__min_samples_split': [2, 5, 10],
}

rf_gs = GridSearchCV(estimator=rf_mod,param_grid=rf_grid, cv=kf, n_jobs=-1)

rf_gs.fit(X_train_bal, y_train_bal)

print('Best CV score: {:.5f}'.format(rf_gs.best_score_))
print('\n')
print('Best parameter:',rf_gs.best_params_)
print('\n')
print('Best estimator:',rf_gs.best_estimator_)

print('F1 score: {:.5f}'.format(f1_score(y_val, rf_gs.predict(X_val), pos_label='Good')))
print('\n')
print('Precision score: {:.5f}'.format(precision_score(y_val, rf_gs.predict(X_val), pos_label='Good')))
print('\n')
print('Recall score: {:.5f}'.format(recall_score(y_val, rf_gs.predict(X_val), pos_label='Good')))
print('\n')
plt.figure(figsize=(8, 6))
plt.title('Matrice di confusione sul modello con Random Forest', size=16)
sns.heatmap(pd.DataFrame(confusion_matrix(y_val, rf_gs.predict(X_val)), columns=rf_gs.classes_, index = rf_gs.classes_),
            annot=True, cmap='Purples', fmt='d');

"""### MLPClassifier

Lo scopo di un classificatore lineare è arrivare ad avere un iperpiano, determinato dalla combinazione lineare di $n$ variabili in input, in grado di dividere le istanze in classi. Se introducessimo delle variabili nascoste (hidden) in modo che l'iperpiano individuato dal modello sia combinazione lineare delle variabili nascoste e ciascuna nascosta sia combinazione lineare di quelle in input, allora otterremmo un *Multi Layer Perceptron* in quanto si tratta di piu modelli lineari disposti a strati. Per questo modello Multi Layer Perceptron cerchiamo di ottenere il massimo rendimento possibile dai seguenti iperparametri:
  - mlp__hidden_layer_sizes o numero di variabili nascoste da introdurre
  - mlp__batch_size o batch size (dimensione del lotto) cioè il numero di istanze che prenderà ad ogni iterazione dell'addestramento
"""

from sklearn.neural_network import MLPClassifier
MLP_mod = Pipeline([
    ("scaler", StandardScaler()),
    ("mlp", MLPClassifier(random_state=42, activation="relu"))
])

MLP_grid = {
    "mlp__hidden_layer_sizes": [16, 32, (16, 8)],
    "mlp__batch_size": [100, 200]
}

MLP_gs = GridSearchCV(estimator=MLP_mod,param_grid=MLP_grid, cv=kf, n_jobs=-1)

MLP_gs.fit(X_train_bal, y_train_bal)

print('Best CV score: {:.5f}'.format(MLP_gs.best_score_))
print('\n')
print('Best parameter:',MLP_gs.best_params_)
print('\n')
print('Best estimator:',MLP_gs.best_estimator_)

print('F1 score: {:.5f}'.format(f1_score(y_val, MLP_gs.predict(X_val), pos_label='Good')))
print('\n')
print('Precision score: {:.5f}'.format(precision_score(y_val, MLP_gs.predict(X_val), pos_label='Good')))
print('\n')
print('Recall score: {:.5f}'.format(recall_score(y_val, MLP_gs.predict(X_val), pos_label='Good')))
print('\n')
plt.figure(figsize=(8, 6))
plt.title('Matrice di confusione sul modello con MLP', size=16)
sns.heatmap(pd.DataFrame(confusion_matrix(y_val, MLP_gs.predict(X_val)), columns=MLP_gs.classes_, index = MLP_gs.classes_),
            annot=True, cmap='Purples', fmt='d');

"""# Valutazione modelli generati

## Valutazione modelli attraverso metriche standard

Iniziamo ora a valutare i modelli fin'ora generati sulla base di alcuni informazioni estrapolate da essi.

La prima metrica che utilizzeremo per confrontare la bontà dei modelli è lo **scarto quadratico medio dell'errore dei relativi iperparametri trovati con la GridSearch**: esso corrisponde alla media dei quadrati delle differenze tra ciascun valore reale e la corrispondente predizione.
"""

from sklearn.model_selection import cross_val_score

name_model = {
    'per': 'Perceptron',
    'log': 'LogisticRegression',
    'svc': 'Support Vector Machines',
    'tree': 'Decision Tree'
}

for model in models.values():
  print(name_model[list(model.named_steps.keys())[1]] +
        ": %0.6f di accuratezza con deviazione standard di %0.6f" %
        (cross_val_score(models[list(model.named_steps.keys())[1]],X_val, y_val,cv=kf).mean(),
         cross_val_score(models[list(model.named_steps.keys())[1]],X_val, y_val,cv=kf).std()))
  print('\n')

"""Come possiamo notare in tutti i metodi i valori degli score sono molto vicini alla media (deviazioni standard piccole).

Confrontiamo i valori  di altri modi per valutare l'accuratezza di un classificatore:
- **(R$^2$) Score**
- **Recall**
- **Precision**
- **F1 Score**
"""

pd.DataFrame.from_dict(scores, orient="index", columns=["(R^2) Score"])

pd.DataFrame.from_dict(recall, orient="index", columns=["Recall"])

pd.DataFrame.from_dict(precision, orient="index", columns=["Precision"])

pd.DataFrame.from_dict(f1, orient="index", columns=["F1 Score"])

"""Dai risultati emerge come tutti i modelli generati diano buoni risultati, riuscendo a classificare correttamente la maggior parte delle istanze. In particolare è chiaro come il modello basato sul SVM sia il migliore (seguito dal Decision Tree) mentre quello istanziato attraverso Perceptron sia il peggiore: ragionevole in quanto il Perceptron è l'algoritmo più semplice e vecchio tra questi testati.

A conferma di quanto appena detto osserviamo le matrici di confusione di questi:
"""

plt.figure(figsize=(4, 3))
plt.title('Matrice di confusione sul modello con Perceptron', size=10)
sns.heatmap(pd.DataFrame(confusion_matrix(y_val, per_gs.predict(X_val)), columns=per_gs.classes_, index = per_gs.classes_),
            annot=True, cmap='Oranges', fmt='d');

plt.figure(figsize=(4, 3))
plt.title('Matrice di confusione sul modello con Regressione Logistica', size=10)
sns.heatmap(pd.DataFrame(confusion_matrix(y_val, log_gs.predict(X_val)), columns=log_gs.classes_, index = log_gs.classes_),
            annot=True, cmap='Reds', fmt='d');

plt.figure(figsize=(4, 3))
plt.title('Matrice di confusione sul modello con SVC', size=10)
sns.heatmap(pd.DataFrame(confusion_matrix(y_val, svc_gs.predict(X_val)), columns=svc_gs.classes_, index = svc_gs.classes_),
            annot=True, cmap='Greens', fmt='d');

plt.figure(figsize=(4, 3))
plt.title('Matrice di confusione sul modello con Albero Decisionale', size=10)
sns.heatmap(pd.DataFrame(confusion_matrix(y_val, tree_gs.predict(X_val)), columns=tree_gs.classes_, index = tree_gs.classes_),
            annot=True, cmap='Purples', fmt='d');

"""## Valutazione modelli con Intervallo di confidenza

Produciamo ora un confronto con intervallo di confidenza per quantificare l'accuratezza dei modelli. Fissiamo quindi la percentuale di certezza al 95%.
"""

def confidence(acc, N, Z=1.96):
    den = (2*(N+Z**2))
    var = (Z*np.sqrt(Z**2+4*N*acc-4*N*acc**2)) / den
    a = (2*N*acc+Z**2) / den
    inf = a - var
    sup = a + var
    return (inf, sup)

def calculate_accuracy(conf_matrix):
    return np.diag(conf_matrix).sum() / conf_matrix.sum().sum()

name_model = {
    'per': 'Perceptron',
    'log': 'LogisticRegression',
    'svc': 'Support Vector Machines',
    'tree': 'Decision Tree'
}

for model in models.values():
  print(name_model[list(model.named_steps.keys())[1]] +
        ": Estremo inferiore %0.6f" %
         confidence(calculate_accuracy(confusion_matrices[list(model.named_steps.keys())[1]]), len(X_val))[0])

  print(name_model[list(model.named_steps.keys())[1]] +
        ": Estremo superiore %0.6f" %
         confidence(calculate_accuracy(confusion_matrices[list(model.named_steps.keys())[1]]), len(X_val))[1])

  print('\n')

"""## Valutazione per confronto con Modello casuale

Creiamo ora un modello casuale, utilizzando **DummyClassifier**. Questo è uno classificatore molto semplice che mette a disposizione sklearn e che può essere utile per confrontare le sue previsioni con quelle reali.
"""

from sklearn.dummy import DummyClassifier

#strategy="uniform" garantisce uniformità randomica nella predizione
rand_mod = DummyClassifier(random_state=42,strategy="uniform")
rand_mod.fit(X_train_bal, y_train_bal)

print('CV score: {:.5f}'.format(rand_mod.score(X_val, y_val)))
print('\n')
print('F1 score: {:.5f}'.format(f1_score(y_val, rand_mod.predict(X_val), pos_label='Good')))
print('\n')
print('Precision score: {:.5f}'.format(precision_score(y_val, rand_mod.predict(X_val), pos_label='Good')))
print('\n')
print('Recall score: {:.5f}'.format(recall_score(y_val, rand_mod.predict(X_val), pos_label='Good')))
print('\n')
plt.figure(figsize=(8, 6))
plt.title('Matrice di confusione sul modello Casuale', size=16)
sns.heatmap(pd.DataFrame(confusion_matrix(y_val, rand_mod.predict(X_val)), columns=rand_mod.classes_, index = rand_mod.classes_),
            annot=True, cmap='Pastel1', fmt='d');

scores["rand"] = rand_mod.score(X_val, y_val);
recall["rand"] = recall_score(y_val, rand_mod.predict(X_val), pos_label='Good')
precision["rand"] = precision_score(y_val, rand_mod.predict(X_val), pos_label='Good')
f1["rand"] = f1_score(y_val, rand_mod.predict(X_val), pos_label='Good')
#models["rand"] = rand_mod
confusion_matrices["rand"] = pd.DataFrame(confusion_matrix(y_val, rand_mod.predict(X_val)), columns=rand_mod.classes_, index = rand_mod.classes_)

from scipy import stats

def difference_between_two_models(error1, error2, confidence):
    z_half_alfa = stats.norm.ppf(confidence)
    variance = (((1 - error1) * error1) / len(y_val)) + (((1 - error2) * error2) / len(y_val))
    d_minus = abs(error1 - error2) - z_half_alfa * (pow(variance, 0.5))
    d_plus = abs(error1 - error2) + z_half_alfa * (pow(variance, 0.5))
    print("Valore minimo: {}\nValore massimo: {}\n".format(d_minus, d_plus))

for model in models.values():
  print(name_model[list(model.named_steps.keys())[1]] + " VS Modello Random")
  print("\n")
  difference_between_two_models(1 - f1[list(model.named_steps.keys())[1]], 1 - f1["rand"], 0.99)

"""Ora abbiamo la certezza al 99% che i diversi modelli addestrati siano migliori di un modello casuale: nessuno dei valori che possiamo visualizzare contiene lo 0.

Risulta ormai evidente come il modello migliore sia quello addestrato con **Support Vector Machines**.

# Confronto finale e Conclusioni

Individuato il modello SVM come quello migliore analizziamolo.

Per individuare l'iperpiano di separazione migliore tra le due classi, SVM si basa sulle istanze più vicine ad esso, i cosiddetti **Support Vector**. Vediamo quindi quali features sono più positivamente/negativamente correlate e in che misura con la variabile *quality*.

Vediamo prima quante istanze son stato necessarie per individuare l'iperpiano:
"""

len(models['svc'][1].support_)

"""Come notiamo le istanze utilizzate son molteplici (come si nota sotto, circa il 64%), questo indica come una divisione netta e precisa tra le due classi sia difficile da individuare. Ancora più complicato per modelli lineare come Perceptron: un altro motivo per la sua scarsa precisione."""

len(models['svc'][1].support_)  / len(X_train_bal) * 100

"""Per visualizzare meglio il tutto possiamo definire una funzione che ci permetta di capire quanti punti siano stati usati per definire l'iperpiano. Plottando i punti che abbiamo usato per il training e quelli usati dal support vector, vedremo come molti combacino"""

def plot_data_and_support(X, model):

    #in Xs metto tutte le righe e colonne di X che hanno indice preso da model.support_
    Xs = X.iloc[model.support_]
    #print(Xs)

    plt.scatter(X.iloc[:, 3], X.iloc[:, 6], s=4, c='black')
    plt.scatter(Xs.iloc[:, 3], Xs.iloc[:, 6], s=4, c='yellow')

plot_data_and_support(X_train_bal, models['svc'][1])

"""Infine scopriamo quale feature son state più o meno rilevanti per i nostri scopi di classificazione."""

support_dataframe = pd.DataFrame(models['svc'][1].support_vectors_, columns=X_train_bal.columns)
#support_dataframe.describe()

mean = pd.Series([0.061903,0.035754,0.043108,-0.031508,0.083313,-0.021621,-0.036785,0.026544,0.031376,0.069159,0.016729], index= X_train_bal.columns)
mean.plot.bar()
plt.grid(axis="y")

"""Come si era previsto nella fase esplorativa dei dati, tutte le features hanno un peso più o meno significativo sul risultato del modello, è quindi importante considerarle tutte per determinare se un vino è buono o cattivo.

In definitiva considerando i valori ottenuti dal modello, ci riteniamo molto soddisfatti di questa nostra prima esperienza nel mondo del Data Science. Le valutazioni affrontate nel finale del documento ci mostrano come il modello SVM possa essere un ottimo strumento di supporto per la classificazione di un vino buono o cattivo.

E' stato scelto questo dataset dal momento che la famiglia di un membro del team di questo progetto (Romagnoli Giacomo) è proprietaria di un'azienda vinicola. Ci immaginiamo ad esempio che un'azienda come quella possa utilizzare questo dispositivo per scremare in un primo momento i vini buoni da quelli cattivi per poi farli assaggiare ad un ipotetico somelier che ne giudicherà l'effettiva bontà.
"""

reason: import numpy,from scipy


text:
import numpy as np
import openslide
import sys
import os
from PIL import Image
import datetime
import time
import cv2
from shutil import copyfile as cp
import multiprocessing as mp
import random

slide_name = sys.argv[2] + '/' + sys.argv[1];
output_folder = sys.argv[3] + '/' + sys.argv[1];
patch_size_40X = 2100;      # 350*8 --> 8*8 smaller non-overlap patches (350x350)
level = 0

if not os.path.exists(sys.argv[3]): os.mkdir(sys.argv[3])

#start = time.time()
#time.sleep(random.randint(100, 1000)/100.0)  # wait for 1 --> 10s to avoid concurrency

fdone = '{}/extraction_done.txt'.format(output_folder);
if os.path.isfile(fdone):
    print('fdone {} exist, skipping'.format(fdone));
    exit(0);

print('extracting {}'.format(output_folder));

if not os.path.exists(output_folder): os.mkdir(output_folder);

try:
    oslide = openslide.OpenSlide(slide_name);
    if openslide.PROPERTY_NAME_MPP_X in oslide.properties:
        mag = 10.0 / float(oslide.properties[openslide.PROPERTY_NAME_MPP_X]);
    elif "XResolution" in oslide.properties:
        mag = 10.0 / float(oslide.properties["XResolution"]);
    elif 'tiff.XResolution' in oslide.properties:   # for Multiplex IHC WSIs, .tiff images
        Xres = float(oslide.properties["tiff.XResolution"])
        if Xres < 10:
            mag = 10.0 / Xres;
        else:
            mag = 10.0 / (10000/Xres)       # SEER PRAD
    else:
        print('[WARNING] mpp value not found. Assuming it is 40X with mpp=0.254!', slide_name);
        mag = 10.0 / float(0.254);

    pw = int(patch_size_40X * mag / 40);  # scale patch size from 40X to 'mag'

    width = oslide.dimensions[0];
    height = oslide.dimensions[1];
except:
    print('{}: exception caught'.format(slide_name));
    exit(1);

print('height/width: {}/{}'.format(height, width))

corrs = []
for x in range(1, width, pw):
    for y in range(1, height, pw):
        if x + pw > width:
            pw_x = width - x;
        else:
            pw_x = pw;

        if y + pw > height:
            pw_y = height - y;
        else:
            pw_y = pw;

        if int(patch_size_40X * pw_x / pw) > 50 and int(patch_size_40X * pw_y / pw) > 50:
            corrs.append((x, y, pw_x, pw_y))

def extract_patch(corr):
    x, y, pw_x, pw_y = corr

    fname = '{}/{}_{}_{}_{}.png'.format(output_folder, x, y, pw, patch_size_40X)

    patch = oslide.read_region((x, y), 0, (pw_x, pw_y));
    #shahira: skip where the alpha channel is zero
    patch_arr = np.array(patch);
    if(patch_arr[:,:,3].max() == 0):    # this is blank regions
        return
    patch = patch.resize((int(patch_size_40X * pw_x / pw), int(patch_size_40X * pw_y / pw)), Image.ANTIALIAS);
    patch.save(fname);

start = time.time()
print(slide_name, len(corrs))
pool = mp.Pool(processes=4)
pool.map(extract_patch, corrs)

print('Elapsed time: ', (time.time() - start)/60.0)


text:
r"""Algorithms for partitioning the non-dominated space into rectangles.

Code is adapted from botorch.

References

.. [Couckuyt2012]
    I. Couckuyt, D. Deschrijver and T. Dhaene, "Towards Efficient
    Multiobjective Optimization: Multiobjective statistical criterions,"
    2012 IEEE Congress on Evolutionary Computation, Brisbane, QLD, 2012,
    pp. 1-8.

"""

from typing import Optional
import numpy as np

from litebo.utils.multi_objective.pareto import is_non_dominated

class NondominatedPartitioning(object):
    r"""A class for partitioning the non-dominated space into hyper-cells.

    Note: this assumes minimization.

    Note: it is only feasible to use this algorithm to compute an exact
    decomposition of the non-dominated space for `m<5` objectives (alpha=0.0).

    The alpha parameter can be increased to obtain an approximate partitioning
    faster. The `alpha` is a fraction of the total hypervolume encapsuling the
    entire pareto set. When a hypercell's volume divided by the total hypervolume
    is less than `alpha`, we discard the hypercell. See Figure 2 in
    [Couckuyt2012]_ for a visual representation.

    This numpy implementation is adapted from botorch, which is adapted from
    https://github.com/GPflow/GPflowOpt/blob/master/gpflowopt/pareto.py.
    """

    def __init__(
        self,
        num_objs: int,
        Y: Optional[np.ndarray] = None,
        alpha: float = 0.0,
        eps: Optional[float] = None,
    ) -> None:
        """Initialize NondominatedPartitioning.

        Args:
            num_objs: The number of objective functions
            Y: A `n x m`-dim array
            alpha: a thresold fraction of total volume used in an approximate
                decomposition.
            eps: a small value for numerical stability
        """
        self.alpha = alpha
        self.num_objs = num_objs
        self._eps = eps
        if Y is not None:
            self.update(Y=Y)

    @property
    def eps(self) -> float:
        if self._eps is not None:
            return self._eps
        try:
            return 1e-6 if self._pareto_Y.dtype == np.float32 else 1e-8
        except AttributeError:
            return 1e-6

    @property
    def pareto_Y(self) -> np.ndarray:
        r"""This returns the non-dominated set assuming minimization.

        Returns:
            A `n_pareto x m`-dim array of outcomes.
        """
        if not hasattr(self, "_pareto_Y"):
            raise Exception("pareto_Y has not been initialized")
        return self._pareto_Y

    def _update_pareto_Y(self) -> bool:
        r"""Update the non-dominated front."""
        non_dominated_mask = is_non_dominated(self.Y)
        pf = self.Y[non_dominated_mask]
        # sort by first objective
        new_pareto_Y = pf[np.argsort(pf[:, 0])]
        if not hasattr(self, "_pareto_Y") or not np.equal(
            new_pareto_Y, self._pareto_Y
        ):
            self._pareto_Y = new_pareto_Y
            return True
        return False

    def update(self, Y: np.ndarray) -> None:
        r"""Update non-dominated front and decomposition.

        Args:
            Y: A `n x m`-dim array of outcomes.
        """
        self.Y = Y
        is_new_pareto = self._update_pareto_Y()
        # Update decomposition if the pareto front changed
        if is_new_pareto:
            if self.num_objs > 2:
                self.binary_partition_non_dominated_space()
            else:
                self.partition_non_dominated_space_2d()

    def binary_partition_non_dominated_space(self):
        r"""Partition the non-dominated space into disjoint hypercells.

        This method works for an arbitrary number of outcomes, but is
        less efficient than `partition_non_dominated_space_2d` for the
        2-outcome case.
        """
        # Extend pareto front with the ideal and anti-ideal point
        ideal_point = self._pareto_Y.min(axis=0, keepdims=True) - 1
        anti_ideal_point = self._pareto_Y.max(axis=0, keepdims=True) + 1

        aug_pareto_Y = np.concatenate([ideal_point, self._pareto_Y, anti_ideal_point], axis=0)
        # The binary parititoning algorithm uses indices the augmented pareto front.
        aug_pareto_Y_idcs = self._get_augmented_pareto_front_indices()

        # Initialize one cell over entire pareto front
        cell = np.zeros((2, self.num_objs), dtype=int)
        cell[1] = aug_pareto_Y_idcs.shape[0] - 1
        stack = [cell]
        total_volume = (anti_ideal_point - ideal_point).prod()

        # hypercells contains the indices of the (augmented) pareto front
        # that specify that bounds of the each hypercell.
        # It is a `2 x num_cells x num_objs`-dim array
        self.hypercells = np.empty((2, 0, self.num_objs))
        outcome_idxr = np.arange(self.num_objs)

        # Use binary partitioning
        while len(stack) > 0:
            cell = stack.pop()
            cell_bounds_pareto_idcs = aug_pareto_Y_idcs[cell, outcome_idxr]
            cell_bounds_pareto_values = aug_pareto_Y[
                cell_bounds_pareto_idcs, outcome_idxr
            ]
            # Check cell bounds
            # - if cell upper bound is better than pareto front on all outcomes:
            #   - accept the cell
            # - elif cell lower bound is better than pareto front on all outcomes:
            #   - this means the cell overlaps the pareto front. Divide the cell along
            #   - its longest edge.
            if (
                ((cell_bounds_pareto_values[1] - self.eps) < self._pareto_Y)
                .any(axis=1)
                .all()
            ):
                # Cell is entirely non-dominated
                self.hypercells = np.concatenate(
                    [self.hypercells, np.expand_dims(cell_bounds_pareto_idcs, 1)], axis=1
                )
            elif (
                ((cell_bounds_pareto_values[0] + self.eps) < self._pareto_Y)
                .any(axis=1)
                .all()
            ):
                # The cell overlaps the pareto front
                # compute the distance (in integer indices)
                idx_dist = cell[1] - cell[0]
                any_not_adjacent = (idx_dist > 1).any()
                cell_volume = (
                    (cell_bounds_pareto_values[1] - cell_bounds_pareto_values[0])
                    .prod(axis=-1)
                )

                # Only divide a cell when it is not composed of adjacent indices
                # and the fraction of total volume is above the approximation
                # threshold fraction
                if (
                    any_not_adjacent
                    and ((cell_volume / total_volume) > self.alpha).all()
                ):
                    # Divide the test cell over its largest dimension
                    # largest (by index length)
                    longest_dim = np.argmax(idx_dist)
                    length = idx_dist[longest_dim]

                    new_length1 = int(round(length / 2.0))
                    new_length2 = length - new_length1

                    # Store divided cells
                    # cell 1: subtract new_length1 from the upper bound of the cell
                    # cell 2: add new_length2 to the lower bound of the cell
                    for bound_idx, length_delta in (
                        (1, -new_length1),
                        (0, new_length2),
                    ):
                        new_cell = cell.copy()
                        new_cell[bound_idx, longest_dim] += length_delta
                        stack.append(new_cell)

    def partition_non_dominated_space_2d(self) -> None:
        r"""Partition the non-dominated space into disjoint hypercells.

        This direct method works for `m=2` outcomes.
        """
        if self.num_objs != 2:
            raise Exception(
                "partition_non_dominated_space_2d requires 2 outputs, "
                "but num_objs=%d" % self.num_objs
            )
        pf_ext_idx = self._get_augmented_pareto_front_indices()
        range_pf_plus1 = np.arange(
            self._pareto_Y.shape[0] + 1
        )
        lower = np.stack([range_pf_plus1, np.zeros_like(range_pf_plus1)], axis=-1)
        upper = np.stack(
            [range_pf_plus1 + 1, pf_ext_idx[-range_pf_plus1 - 1, -1]], axis=-1
        )
        self.hypercells = np.stack([lower, upper], axis=0)

    def _get_augmented_pareto_front_indices(self) -> np.ndarray:
        r"""Get indices of augmented pareto front."""
        pf_idx = np.argsort(self._pareto_Y, axis=0)
        return np.concatenate(
            [
                np.zeros(
                    (1, self.num_objs), dtype=int
                ),
                # Add 1 because index zero is used for the ideal point
                pf_idx + 1,
                np.full(
                    (1, self.num_objs),
                    self._pareto_Y.shape[0] + 1,
                    dtype=int
                ),
            ],
            axis=0,
        )

    def get_hypercell_bounds(self, ref_point: np.ndarray) -> np.ndarray:
        r"""Get the bounds of each hypercell in the decomposition.

        Args:
            ref_point: A `m`-dim array containing the reference point.

        Returns:
            A `2 x num_cells x num_objs`-dim array containing the
                lower and upper vertices bounding each hypercell.
        """
        aug_pareto_Y = np.concatenate(
            [
                # -inf is the lower bound of the non-dominated space
                np.full(
                    (1, self.num_objs),
                    float("-inf")
                ),
                self._pareto_Y,
                np.expand_dims(ref_point, 0),
            ],
            axis=0,
        )
        return self._get_hypercell_bounds(aug_pareto_Y=aug_pareto_Y)

    def _get_hypercell_bounds(self, aug_pareto_Y: np.ndarray) -> np.ndarray:
        r"""Get the bounds of each hypercell in the decomposition.

        Args:
            aug_pareto_Y: A `n_pareto + 2 x m`-dim array containing
            the augmented pareto front.

        Returns:
            A `2 x num_cells x num_objs`-dim array containing the
                lower and upper vertices bounding each hypercell.
        """
        num_cells = self.hypercells.shape[1]
        outcome_idxr = np.tile(np.arange(self.num_objs), num_cells)

        # this array is 2 x (num_cells *num_objs) x 2
        # the batch dim corresponds to lower/upper bound
        cell_bounds_idxr = np.stack(
            [self.hypercells.reshape(2, -1),
             np.tile(outcome_idxr[None, :], (2, 1))],
            axis=-1,
        ).astype(int)

        def chunk(arr, chunks, axis):
            num_to_chunk = arr.shape[axis]
            chunk_size = np.ceil(num_to_chunk / chunks).astype(int)
            real_chunks = np.ceil(num_to_chunk / chunk_size).astype(int)
            return tuple(np.take(arr, indices=range(i*chunk_size, (i+1)*chunk_size), axis=axis) for i in range(real_chunks))

        indexers = chunk(cell_bounds_idxr, self.num_objs, axis=-1)
        cell_bounds_values = aug_pareto_Y[indexers].reshape(2, -1, self.num_objs)
        return cell_bounds_values

    def compute_hypervolume(self, ref_point: np.ndarray) -> float:
        r"""Compute the hypervolume for the given reference point.

        Note: This assumes minimization.

        This method computes the hypervolume of the non-dominated space
        and computes the difference between the hypervolume between the
        ideal point and hypervolume of the non-dominated space.

        Note there are much more efficient alternatives for computing
        hypervolume when m > 2 (which do not require partitioning the
        non-dominated space). Given such a partitioning, this method
        is quite fast.

        Args:
            ref_point: A `m`-dim array containing the reference point.

        Returns:
            The dominated hypervolume.
        """
        if (self._pareto_Y <= ref_point).any():
            raise ValueError(
                "The reference point must be smaller than all pareto_Y values."
            )
        ideal_point = self._pareto_Y.min(axis=0, keepdims=True).values
        ref_point = np.expand_dims(ref_point, 0)
        aug_pareto_Y = np.concatenate([ideal_point, self._pareto_Y, ref_point], axis=0)
        cell_bounds_values = self._get_hypercell_bounds(aug_pareto_Y=aug_pareto_Y)
        total_volume = (ref_point - ideal_point).prod()
        non_dom_volume = (
            (cell_bounds_values[1] - cell_bounds_values[0]).prod(axis=-1).sum()
        )
        return total_volume - non_dom_volume

max_stars_repo_name: usualwitch/lite-bo

max_issues_repo_name: usualwitch/lite-bo

max_forks_repo_name: usualwitch/lite-bo


text:
"""
pdb.py
readers for pdb file formats, for molssi best practices workshop

"""

import numpy as np

def open_pdb(file_location):

    # This function reads in a pdb file and returns the atom names and coordinates.
    with open(file_location) as f:
        data = f.readlines()

    coordinates = []
    symbols = []
    for line in data:
        if 'ATOM' in line[0:6] or 'HETATM' in line[0:6]:
            symbols.append(line[76:79].strip())
            atom_coordinates = [float(x) for x in line[30:55].split()]
            coordinates.append(atom_coordinates)

    # Convert coordinate list to numpy array
    coordinates = np.array(coordinates)
    return symbols, coordinates

max_stars_repo_name: jsperhac/molecool

max_issues_repo_name: jsperhac/molecool

max_forks_repo_name: jsperhac/molecool


text:
from __future__ import division
import numpy as np
import scipy as sp

def construct_A(n):
    A = np.diag([2] * n)

    A += np.diag([-1] * (n - 1), k=1)
    A += np.diag([-1] * (n - 1), k=-1)

    condn_no = np.linalg.cond(A, p=2)
    return A, condn_no

if __name__ == '__main__':
    A, condn_no = construct_A(512)
    print("A = \n", A)
    print("Condition number of A  = ", condn_no)

max_stars_repo_name: justachetan/scientific-computing

max_issues_repo_name: justachetan/scientific-computing

max_forks_repo_name: justachetan/scientific-computing


text:
import numpy as np
import EpidemiologicalModels.CellManagement as CellManagement

class CellSpaceConfiguration:

    def __init__(self, nRows, nColumns, xnRows = -1, xnColumns = -1, unRows = 0, unColumns = 0):
        self.__nRows = nRows; self.__nColumns = nColumns  # Dimensiones de la región interna inicial
        self.__xnRows = xnRows; self.__xnColumns = xnColumns  # Ubicación de la región interna inicial
        self.__unRows = unRows; self.__unColumns = unColumns  # Dimensiones del área externa
        self.nRows = max(self.__nRows, self.__xnRows)
        self.nColumns = max(self.__nColumns, self.__xnColumns)
        self.__createInternalSpace()

    def __validate(self):
        if self.__nRows < 0:
            print("Debe definir un número de filas mayor a cero.")
            return False
        elif self.__nColumns < 0:
            print("Debe definir un número de columnas mayor a cero.")
            return False
        else:
            return True

    def basicParameters(self):
        """Parámetros con los que se definió el sistema inicialmente."""
        return (self.__nRows, self.__nColumns, self.__xnRows, self.__xnColumns, self.__unRows, self.__unColumns)

    # Forma del espacio de células

    def __createInternalSpace(self):
        if self.__validate():
            if self.__xnRows == -1 or self.__nColumns == -1:
                self.system = np.zeros((self.__nRows, self.__nColumns))
            else:
                self.system = -np.ones((self.__xnRows, self.__xnColumns))
                self.system = self.rectangularBoundary(self.__nRows, self.__nColumns, self.__unRows, self.__unColumns)

    def __boundaryDefinition(self, boundaryConditions):
        """Definición de la estructura del sistema dadas las condiciones de frontera
        boundaryConditions : Lista con las posiciones con individuos dentro del sistema"""
        for condition in range(len(boundaryConditions)):
            self.system[boundaryConditions[condition][0],boundaryConditions[condition][1]] = 0
        return self.system

    def rectangularBoundary(self, rectangleRows, rectangleColumns, rowPosition, columnPosition):
        """Ubica una matriz nula de tamaño rectangleRows*rectangleColumns en la posición a,b del sistema"""
        boundaryConditions = []
        for row in range(rectangleRows):
            for column in range(rectangleColumns):
                boundaryConditions.append((rowPosition + row,columnPosition + column))
        return self.__boundaryDefinition(boundaryConditions)

    # Condición inicial para aplicar el modelo epidemiológico

    def initialLocationOfInfected(self, initialPercentageInfected, position = "random", percentageOfInfectedMisplaced = 0):
        """ubicación inicial de infectados
        position : random
                   northwest  north   northeast
                   west       center  east
                   southwest  south   southeast"""
        if position == "random":
            return CellManagement.CellManagement(self).InitialCondition(initialPercentageInfected)
        else:
            # Se divide la zona rectángular en 9 bloques
            a = int(self.nRows/3); b = int(self.nColumns/3)
            system = CellManagement.CellManagement(self).InitialCondition(initialPercentageInfected*percentageOfInfectedMisplaced).system
            infectedBlock = CellManagement.CellManagement(CellSpaceConfiguration(a,b)).InitialCondition(initialPercentageInfected-percentageOfInfectedMisplaced,).system
            if position == "northwest":
                for i in range(a):
                    for j in range(b):
                        system[i][j] = infectedBlock[i][j]
            elif position == "north":
                for i in range(a):
                    for j in range(b,2*b):
                        system[i][j] = infectedBlock[i][j-b]
            elif position == "northeast":
                for i in range(a):
                    for j in range(2*b,3*b):
                        system[i][j] = infectedBlock[i][j-2*b]
            elif position == "west":
                for i in range(a,a*2):
                    for j in range(b):
                        system[i][j] = infectedBlock[i-a][j]
            elif position == "center":
                for i in range(a,a*2):
                    for j in range(b,2*b):
                        system[i][j] = infectedBlock[i-a][j-b]
            elif position == "east":
                for i in range(a,a*2):
                    for j in range(2*b,3*b):
                        system[i][j]=infectedBlock[i-a][j-2*b]
            elif position == "southwest":
                for i in range(2*a,3*a):
                    for j in range(b):
                        system[i][j] = infectedBlock[i-2*a][j]
            elif position == "south":
                for i in range(2*a,3*a):
                    for j in range(b,2*b):
                        system[i][j] = infectedBlock[i-2*a][j-b]
            elif position == "southeast":
                for i in range(2*a,3*a):
                    for j in range(2*b,3*b):
                        system[i][j] = infectedBlock[i-2*a][j-2*b]
            return system


text:
from collections import Counter
from math import ceil, log2

from sympy.combinatorics.graycode import GrayCode, gray_to_bin

def get_gray_length(n):
    """
    :param n: The length of the raw data
    :return: The length of the gray code required for indices of this data
    """
    return ceil(log2(2 * n + 1))

def get_encoded_length(data_len):
    """
    :param data_len: The length of the raw data
    :return: The length of the encoded data
    """
    return data_len + get_gray_length(data_len) + 1

def flip_bit(bit):
    """
    :param bit: '0' or '1'
    :return: '1' or '0', respectively
    """
    return str(1 - int(bit))

def flip_string(string):
    """
    Flips the given binary string - '0' to '1' and vice versa.
    :param string: String to flip.
    :return: Flipped string.
    """
    return ''.join([flip_bit(bit) for bit in string])

def get_partially_flipped_string(data, i):
    """
    Get the partially flipped data corresponding to the index i from the (i-1)th partially-flipped string.
    :param data: Data to transform.
    :param i: The index of the flip.
    :return: Partially flipped data for the i'th index.
    """
    if 0 < i <= len(data):
        data = data[:i - 1] + flip_bit(data[i - 1]) + data[i:]
    elif i > len(data):
        data = data[:i - len(data) - 1] + flip_bit(data[i - len(data) - 1]) + data[i - len(data):]
    return data

def add_balance_bit(encoded_string, counter):
    """
    Adds a balance bit (a bit that's used to compensate for the 2-jumps in the encoded word's weight throughout
    the encoding process)

    :param encoded_string: The encoded (data + gray index) string.
    :param counter: Counter for the amount of 0's and 1's in encoded_string.
    :return: encoded_string with a balance bit as a suffix, if it's possible to balance it, otherwise None.
    """
    if 0 <= (counter['1'] - counter['0']) <= 2:
        return encoded_string + '0'
    elif 1 <= (counter['0'] - counter['1']) <= 2:
        return encoded_string + '1'

def encode_gray_knuth(string):
    """
    Encodes the given binary string using GrayKnuth® algorithm.
    :param string: Data to encode (represented as a binary string).
    :type string: str
    :return: Encoded data.
    :rtype: str
    """
    validate_input(string)
    gray_len = get_gray_length(len(string))
    gray_code = GrayCode(gray_len)
    flipped_data = string
    for i, suffix in zip(range(2 * len(string) + 1), gray_code.generate_gray()):
        flipped_data = get_partially_flipped_string(flipped_data, i)
        encoded_string = flipped_data + suffix
        counter = Counter(encoded_string)
        balanced_encoded_string = add_balance_bit(encoded_string, counter)
        if balanced_encoded_string:
            return balanced_encoded_string

def find_decoded_length(encoded_length):
    """
    Finds the length of the decoded data according to the length of the encoded data.
    Implemented using a binary search, relying on the fact that the encoded length as function of the decoded length
    is increasing.
    :param encoded_length: Length of the encoded data. (as a number and not a binary string)
    :return: Length of the decoded data. (as a number and not a binary string)
    """
    start = 0
    end = encoded_length

    while start <= end:
        mid = (start + end) // 2
        mid_encoded_length = get_encoded_length(mid)
        if encoded_length == mid_encoded_length:
            return mid

        if encoded_length < mid_encoded_length:
            end = mid - 1
        else:
            start = mid + 1

def decode_gray_knuth(string, decoded_length=0):
    """
    Decodes the given binary string using GrayKnuth® algorithm.
    :param string: Data to decode. (as a GrayKnuth binary string)
    :type string: str
    :param decoded_length: Optional. The expected length of the decoded data. If 0 - it is determined by the length
    of the encoded data.
    :type decoded_length: int
    :return: Decoded data. (as a binary string)
    :rtype: str
    """
    validate_input(string)
    if decoded_length == 0:
        decoded_length = find_decoded_length(len(string))
    gray_index = string[decoded_length:-1]
    index = int(gray_to_bin(gray_index), 2)
    if index <= decoded_length:
        return flip_string(string[:index]) + string[index:decoded_length]
    else:
        index = index - decoded_length
        return string[:index] + flip_string(string[index:decoded_length])

def validate_input(input_string):
    for c in input_string:
        if c not in ['0', '1']:
            raise Exception("The input should be binary string, the char {} is illegal".format(str(c)))

def encode_gray_knuth_from_file(filename):
    """
    Encodes binary string from a file, using GrayKnuth® algorithm.
    :param filename: a path to file that contains the Data to encode (a number represented in binary format).
    :type filename: str
    :return: Encoded data.
    :rtype: str
    """
    with open(filename, 'r') as f:
        return encode_gray_knuth(f.read())

def encode_gray_knuth_to_file(string, filename):
    """
    Encodes binary string to a file, using GrayKnuth® algorithm.
    :param string: Data to encode (represented as a binary string).
    :type string: str
    :param filename: a path to create the file that will contain the encoded string.
    :type filename: str
    :return: None
    """
    with open(filename, 'w') as f:
        f.write(encode_gray_knuth(string))

def decode_gray_knuth_from_file(filename):
    """
    Decodes the given binary string from a file using GrayKnuth® algorithm.
    :param filename: a path to the file that contains the encoded string.
    :type filename: str
    :return: None
    """
    with open(filename, 'r') as f:
        return decode_gray_knuth(f.read())

def decode_gray_knuth_to_file(string, filename, decoded_length=0):
    """
    Decodes the given binary string to a file using GrayKnuth® algorithm.
    :param string: Data to decode. (as a GrayKnuth binary string)
    :type string: str
    :param decoded_length: Optional. The expected length of the decoded data. If 0 - it is determined by the length
    of the encoded data.
    :type decoded_length: int
    :param filename: a path to create the file that will contain the decoded string.
    :type filename: str
    :return: None
    """
    decoded_word = decode_gray_knuth(string, decoded_length)
    with open(filename, 'w') as f:
        f.write(decoded_word)

max_stars_repo_name: itay-nakash/GrayCodeKnuth

max_issues_repo_name: itay-nakash/GrayCodeKnuth

max_forks_repo_name: itay-nakash/GrayCodeKnuth

reason: from sympy


text:
# Copyright 2019-2020 ETH Zurich and the DaCe authors. All rights reserved.
import dace
import numpy as np

@dace.program
def internal(out: dace.float64[2]):
    with dace.tasklet:
        o1 >> out[0]
        o2 >> out[1]
        o1 = 5
        o2 = 3

@dace.program
def inline_noinput(A: dace.float64[2]):
    for i in dace.map[0:1]:
        internal(A)

def test():
    A = np.random.rand(2)

    inline_noinput(A)

    diff = np.linalg.norm(A - np.array([5., 3.]))
    print('Difference:', diff)
    assert diff < 1e-5

if __name__ == "__main__":
    test()

max_stars_repo_name: tobiasholenstein/dace

max_issues_repo_name: tobiasholenstein/dace

max_forks_repo_name: tobiasholenstein/dace


text:
"""

"""
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

import cv2 as cv
import numpy as np
import skimage
from PIL import Image
import tensorflow as tf
from nets.dce_net import enhance_np
data = skimage.data
#
# class tfboard_loss_callback(tf.keras.callbacks.Callback):
#     """
#
#     """
#     def __init__(self, log_dir='./logs/loss/', losser=None, session=None):
#         super(tfboard_loss_callback, self).__init__()
#         self.losser = losser
#         self.session = session
#         self.writer = tf.summary.FileWriter(log_dir)
#         self.sumary_loss1 = tf.summary.scalar('loss_1', self.losser[7])
#         self.sumary_loss2 = tf.summary.scalar('loss_2', self.losser[8])
#         self.sumary_loss3 = tf.summary.scalar('loss_3', self.losser[9])
#         self.sumary_loss4 = tf.summary.scalar('loss_4', self.losser[10])
#
#
#     def on_epoch_end(self, epoch, logs=None):
#         """
#
#         :param epoch:
#         :param logs:
#         :return:
#         """
#         loss1, loss2, loss3, loss4 = self.session.run([self.sumary_loss1,
#                                                        self.sumary_loss2,
#                                                        self.sumary_loss3,
#                                                        self.sumary_loss4])
#         self.writer.add_summary(loss1, epoch)
#         self.writer.add_summary(loss2, epoch)
#         self.writer.add_summary(loss3, epoch)
#         self.writer.add_summary(loss4, epoch)

def colormap_jet(img):
    min = img.min()
    if min < 0:
        img -= min
    img *= (255.0/img.max())
    temp = cv.applyColorMap(np.uint8(img), 2)
    out = temp # cv.cvtColor(temp, cv.COLOR_BGR2RGB)
    return out

class conv2d_callback(tf.keras.callbacks.Callback):
    def __init__(self, log_dir='./logs/tmp/', feed_inputs_display=None):
        super(conv2d_callback, self).__init__()
        self.seen = 0
        self.feed_inputs_display = feed_inputs_display
        self.writer = tf.summary.FileWriter(log_dir)

    def custom_set_feed_input_to_display(self, feed_inputs_display):
        self.feed_inputs_display = feed_inputs_display

    # copied from the above answers;
    def make_image(self, numpy_img):
        height, width, channel = numpy_img.shape
        image = Image.fromarray(numpy_img)
        import io
        output = io.BytesIO()
        image.save(output, format='PNG')
        image_string = output.getvalue()
        output.close()
        return tf.Summary.Image(height=height, width=width, colorspace=channel, encoded_image_string=image_string)

    # A callback has access to its associated model through the class property self.model.
    def on_epoch_end(self, epoch, logs=None):
        logs = logs or {}
        self.seen += 1
        if epoch % 1 == 0:  # every 200 iterations or batches, plot the costumed images using TensorBorad;
            summary_str = []
            for i in range(len(self.feed_inputs_display)):
                feature, _, _ = self.feed_inputs_display[i]

                summary_str.append(tf.Summary.Value(tag='Compare{}/lowlight_img'.format(i), image=self.make_image(
                    feature)))  # function colormap_jet(), defined above;

                lowlight_img = np.expand_dims(feature, axis=0)/255.0

                output_np = tf.keras.backend.get_session().run(
                            self.model.output, feed_dict={self.model.input: lowlight_img}
                        )

                enhance_img = enhance_np([lowlight_img, output_np])
                enhance_img = np.squeeze(enhance_img, axis=0)

                #
                min = enhance_img.min()
                if min < 0:
                    enhance_img -= min
                enhance_img *= (255.0 / enhance_img.max())
                enhance_img = enhance_img.astype(np.uint8)
                # enhance_img = cv.cvtColor(enhance_img, cv.COLOR_BGR2RGB)

                summary_str.append(tf.Summary.Value(tag='Compare{}/enhanced_img'.format(i), image=self.make_image(
                    enhance_img)))
                #

                layer_names = [layer.name for layer in self.model.layers if 'conv' in layer.name]
                for name in layer_names:
                    disp_pred, output_np = tf.keras.backend.get_session().run(
                            [self.model.get_layer(name).output, self.model.output], feed_dict={self.model.input: lowlight_img}
                        )
                    disp_pred = np.squeeze(disp_pred, axis=0)
                    # disp_pred = np.squeeze(self.model.predict_on_batch(feature), axis = 0)

                    disp_size = disp_pred.shape[-1]
                    for j in range(disp_size):
                        gray_img = disp_pred[:, :, j]
                        gray_img = np.expand_dims(gray_img, axis=-1)

                        summary_str.append(
                            tf.Summary.Value(tag='{}/img{}/out{}'.format(name, i, j), image=self.make_image(
                                colormap_jet(gray_img))))

            self.writer.add_summary(tf.Summary(value=summary_str), global_step=self.seen)

class save_ckpt_callback(tf.keras.callbacks.Callback):
    def __init__(self, saver, session, save_path):
        super(save_ckpt_callback, self).__init__()
        self.saver = saver
        self.session = session
        self.save_path = save_path

    def on_epoch_end(self, epoch, logs=None):
        self.saver.save(self.session, self.save_path, global_step=epoch)


text:
import math
import numpy as np
from suitcase.structure import Structure
from suitcase.fields import BaseField
from suitcase.fields import Magic
from suitcase.fields import CRCField
from suitcase.fields import SLInt8
from suitcase.fields import ULInt8

PACKET_START_CHAR = "!"
PROTOCOL_FLOAT_LENGTH = 7
PROTOCOL_UINT16_LENGTH = 4
CHECKSUM_LENGTH = 2
TERMINATOR_LENGTH = 2

MPU_INIT_FAILURE = "mpu_init failed!"
MPU_FIRMWARE_LOAD_ERROR = "Firmware Load ERROR"

UPDATE_RATE_HZ = 100

MSGID_YPR_UPDATE = 'y'
YPR_UPDATE_MESSAGE_LENGTH = 34
YPR_UPDATE_YAW_VALUE_INDEX = 2
YPR_UPDATE_PITCH_VALUE_INDEX = 9
YPR_UPDATE_ROLL_VALUE_INDEX = 16
YPR_UPDATE_COMPASS_VALUE_INDEX = 23
YPR_UPDATE_CHECKSUM_INDEX = 30
YPR_UPDATE_TERMINATOR_INDEX = 32

MSGID_QUATERNION_UPDATE = 'q'
QUATERNION_UPDATE_MESSAGE_LENGTH = 53
QUATERNION_UPDATE_QUAT1_VALUE_INDEX = 2
QUATERNION_UPDATE_QUAT2_VALUE_INDEX = 6
QUATERNION_UPDATE_QUAT3_VALUE_INDEX = 10
QUATERNION_UPDATE_QUAT4_VALUE_INDEX = 14
QUATERNION_UPDATE_ACCEL_X_VALUE_INDEX = 18
QUATERNION_UPDATE_ACCEL_Y_VALUE_INDEX = 22
QUATERNION_UPDATE_ACCEL_Z_VALUE_INDEX = 26
QUATERNION_UPDATE_MAG_X_VALUE_INDEX = 30
QUATERNION_UPDATE_MAG_Y_VALUE_INDEX = 34
QUATERNION_UPDATE_MAG_Z_VALUE_INDEX = 38
QUATERNION_UPDATE_TEMP_VALUE_INDEX = 42
QUATERNION_UPDATE_CHECKSUM_INDEX = 49
QUATERNION_UPDATE_TERMINATOR_INDEX = 51

MSGID_GYRO_UPDATE = 'g'
GYRO_UPDATE_MESSAGE_LENGTH = 46
GYRO_UPDATE_GYRO_X_VALUE_INDEX = 2
GYRO_UPDATE_GYRO_Y_VALUE_INDEX = 6
GYRO_UPDATE_GYRO_Z_VALUE_INDEX = 10
GYRO_UPDATE_ACCEL_X_VALUE_INDEX = 14
GYRO_UPDATE_ACCEL_Y_VALUE_INDEX = 18
GYRO_UPDATE_ACCEL_Z_VALUE_INDEX = 22
GYRO_UPDATE_MAG_X_VALUE_INDEX = 26
GYRO_UPDATE_MAG_Y_VALUE_INDEX = 30
GYRO_UPDATE_MAG_Z_VALUE_INDEX = 34
GYRO_UPDATE_TEMP_VALUE_INDEX = 38
GYRO_UPDATE_CHECKSUM_INDEX = 42
GYRO_UPDATE_TERMINATOR_INDEX = 44

MSGID_STREAM_CMD = 'S'
STREAM_CMD_MESSAGE_LENGTH = 9
STREAM_CMD_STREAM_TYPE_YPR = MSGID_YPR_UPDATE
STREAM_CMD_STREAM_TYPE_QUATERNION = MSGID_QUATERNION_UPDATE
STREAM_CMD_STREAM_TYPE_GYRO = MSGID_GYRO_UPDATE
STREAM_CMD_STREAM_TYPE_INDEX = 2
STREAM_CMD_UPDATE_RATE_HZ_INDEX = 3
STREAM_CMD_CHECKSUM_INDEX = 5
STREAM_CMD_TERMINATOR_INDEX = 7

MSG_ID_STREAM_RESPONSE = 's'
STREAM_RESPONSE_MESSAGE_LENGTH = 46
STREAM_RESPONSE_STREAM_TYPE_INDEX = 2
STREAM_RESPONSE_GYRO_FULL_SCALE_DPS_RANGE = 3
STREAM_RESPONSE_ACCEL_FULL_SCALE_G_RANGE = 7
STREAM_RESPONSE_UPDATE_RATE_HZ = 11
STREAM_RESPONSE_YAW_OFFSET_DEGREES = 15
STREAM_RESPONSE_QUAT1_OFFSET = 22
STREAM_RESPONSE_QUAT2_OFFSET = 26
STREAM_RESPONSE_QUAT3_OFFSET = 30
STREAM_RESPONSE_QUAT4_OFFSET = 34
STREAM_RESPONSE_FLAGS = 38
STREAM_RESPONSE_CHECKSUM_INDEX = 42
STREAM_RESPONSE_TERMINATOR_INDEX = 44

def valid_message(data):
    if type(data) != str:
        return False
    if len(data) < 2:
        return False
    if data[0] != PACKET_START_CHAR:
        return False

    if data[1] is MSGID_YPR_UPDATE:
        return len(data) == YPR_UPDATE_MESSAGE_LENGTH
    elif data[1] is MSGID_QUATERNION_UPDATE:
        return len(data) == QUATERNION_UPDATE_MESSAGE_LENGTH
    elif data[1] is MSGID_GYRO_UPDATE:
        return len(data) == GYRO_UPDATE_MESSAGE_LENGTH
    elif data[1] is MSGID_STREAM_CMD:
        return len(data) == STREAM_CMD_MESSAGE_LENGTH
    elif data[1] is MSG_ID_STREAM_RESPONSE:
        return len(data) == STREAM_RESPONSE_MESSAGE_LENGTH

    return False

def checksum(data):
    ret = 0
    for char in data:
        ret += ord(char)
    return "%02X" % (ret % 256)

class ProtocolChecksum(BaseField):
    def __init__(self, **kwargs):
        BaseField.__init__(self, **kwargs)
        self.bytes_required = 2

    def pack(self, stream):
        if type(self._value) is not str:
            raise RuntimeError, "Not a string"
        stream.write(self._value)

    def unpack(self, data):
        self._value = str(data)

class ProtocolUInt16(BaseField):
    def __init__(self, **kwargs):
        BaseField.__init__(self, **kwargs)
        self.bytes_required = PROTOCOL_UINT16_LENGTH

    def pack(self, stream):
        if self._value is not None:
            stream.write("%04X" % np.uint16(self._value))
        else:
            stream.write("0000")

    def unpack(self, data):
        # Expecting a hexadecimal string

        self._value = np.uint16(int(data, 16))

class ProtocolSInt16(BaseField):
    def __init__(self, **kwargs):
        BaseField.__init__(self, **kwargs)
        self.bytes_required = PROTOCOL_UINT16_LENGTH

    def pack(self, stream):
        if self._value is not None:
            stream.write("%04X" % np.uint16(self._value))
        else:
            stream.write("0000")

    def unpack(self, data):
        # Expecting a hexadecimal string
        self._value = np.int16(int(data, 16))

class ProtocolUInt8(BaseField):
    def __init__(self, **kwargs):
        BaseField.__init__(self, **kwargs)
        self.bytes_required = 1

    def pack(self, stream):
        if self._value is not None:
            stream.write("%02X" % np.uint8(self._value))
        else:
            stream.write("00")

    def unpack(self, data):
        self._value = np.int8(ord(data))

class ProtocolFloat(BaseField):
    def __init__(self, **kwargs):
        BaseField.__init__(self, **kwargs)
        self.bytes_required = PROTOCOL_FLOAT_LENGTH

    def pack(self, stream):
        if self._value is not None:
            f = float(self._value)
            val = ""
            if f < 0:
                val += "-"
            else:
                val += " "

            f = abs(f)
            val += "%03d." % f
            val += "%02d" % int(round((f - math.floor(f)) * 100))
            stream.write(val)
        else:
            stream.write(" 000.00")

    def unpack(self, data):
        #print data
        self._value = float(data)

class YPRUpdate(Structure):
    start = Magic("!")
    type = Magic("y")
    yaw = ProtocolFloat()
    pitch = ProtocolFloat()
    roll = ProtocolFloat()
    compass_heading = ProtocolFloat()
    checksum = CRCField(ProtocolChecksum(), checksum, 0, -4)
    termination = Magic("\r\n")

class QuaternionUpdate(Structure):
    start = Magic("!")
    type = Magic("q")
    q1 = ProtocolUInt16()
    q2 = ProtocolUInt16()
    q3 = ProtocolUInt16()
    q4 = ProtocolUInt16()
    accel_x = ProtocolSInt16()
    accel_y = ProtocolSInt16()
    accel_z = ProtocolSInt16()
    mag_x = ProtocolSInt16()
    mag_y = ProtocolSInt16()
    mag_z = ProtocolSInt16()
    temp_c = ProtocolFloat()
    checksum = CRCField(ProtocolChecksum(), checksum, 0, -4)
    termination = Magic("\r\n")

class StreamCommand(Structure):
    start = Magic("!")
    type = Magic(MSGID_STREAM_CMD)
    stream_type = ULInt8()
    update_rate_hz = ProtocolUInt8()
    checksum = CRCField(ProtocolChecksum(), checksum, 0, -4)
    termination = Magic("\r\n")

class StreamResponse(Structure):
    start = Magic("!")
    type = Magic("s")

    stream_type = ProtocolUInt8()

    gyro_fsr_dps = ProtocolUInt16()
    accel_fsr_g = ProtocolUInt16()
    update_rate_hz = ProtocolUInt16()
    yaw_offset_degrees = ProtocolFloat()
    q1_offset = ProtocolUInt16()
    q2_offset = ProtocolUInt16()
    q3_offset = ProtocolUInt16()
    q4_offset = ProtocolUInt16()
    flags = ProtocolUInt16()

    checksum = CRCField(ProtocolChecksum(), checksum, 0, -4)
    termination = Magic("\r\n")

def is_stream_response(data):
    return len(data) >= 2 and data[1] == MSG_ID_STREAM_RESPONSE

def is_ypr_update(data):
    return len(data) >= 2 and data[1] == MSGID_YPR_UPDATE

def is_quaternion_update(data):
    return len(data) >= 2 and data[1] == MSGID_QUATERNION_UPDATE

def is_gyro_update(data):
    return len(data) >= 2 and data[1] == MSGID_GYRO_UPDATE

def is_mpu_init_failure(data):
    if MPU_INIT_FAILURE in data:
        return True
    elif MPU_FIRMWARE_LOAD_ERROR in data:
        return True
    return False

def make_quaternion_cmd_packet():
    ret = StreamCommand()
    ret.update_rate_hz = UPDATE_RATE_HZ
    ret.stream_type = ord(MSGID_QUATERNION_UPDATE)
    return ret

def fix_quaternion(q):
    q2 = map(lambda x:  x/16384.0, q)
    for i in xrange(0, 3):
        if q2[i] >= 2:
            q2[i] -= 4
    return q2

if __name__ == "__main__":
    # y = YPRUpdate()
    # y.yaw = -222.11
    # y.pitch = 1
    # y.roll = 2
    # y.compas_heading = 3

    # print y.pack()
    # print float(" 052.00")
    # ypr = YPRUpdate.from_data("!y-000.09 006.18-005.11 147.26DF\r\n")
    # print y.yaw
    # print (make_quaternion_cmd_packet().pack())
    # q = QuaternionUpdate.from_data("!q22C8DF8A2AC8FE90D2DADAB2E5D6FEDD021F004A 023.554C\r\n")
    # print "22C8", int("22C8", 16), np.uint16(int("22C8", 16))
    # print q.temp_c

    s = StreamResponse.from_data("!sq07D000020064-081.38228BDF782AEAFE80000250\r\n")
    print s