# Učitaj osnovne biblioteke...
import sklearn
import codecs
import mlutils
import matplotlib.pyplot as plt
import pgmpy as pgm
%pylab inline

"""
Explanation: Sveučilište u Zagrebu
Fakultet elektrotehnike i računarstva  
Strojno učenje 2018/2019
http://www.fer.unizg.hr/predmet/su

Laboratorijska vježba 5: Probabilistički grafički modeli, naivni Bayes, grupiranje i vrednovanje klasifikatora
Verzija: 1.4
Zadnji put ažurirano: 11. siječnja 2019.
(c) 2015-2019 Jan Šnajder, Domagoj Alagić  
Objavljeno: 11. siječnja 2019.
Rok za predaju: 21. siječnja 2019. u 07:00h

Upute
Peta laboratorijska vježba sastoji se od tri zadatka. U nastavku slijedite upute navedene u ćelijama s tekstom. Rješavanje vježbe svodi se na dopunjavanje ove bilježnice: umetanja ćelije ili više njih ispod teksta zadatka, pisanja odgovarajućeg kôda te evaluiranja ćelija. 
Osigurajte da u potpunosti razumijete kôd koji ste napisali. Kod predaje vježbe, morate biti u stanju na zahtjev asistenta (ili demonstratora) preinačiti i ponovno evaluirati Vaš kôd. Nadalje, morate razumjeti teorijske osnove onoga što radite, u okvirima onoga što smo obradili na predavanju. Ispod nekih zadataka možete naći i pitanja koja služe kao smjernice za bolje razumijevanje gradiva (nemojte pisati odgovore na pitanja u bilježnicu). Stoga se nemojte ograničiti samo na to da riješite zadatak, nego slobodno eksperimentirajte. To upravo i jest svrha ovih vježbi.
Vježbe trebate raditi samostalno. Možete se konzultirati s drugima o načelnom načinu rješavanja, ali u konačnici morate sami odraditi vježbu. U protivnome vježba nema smisla.
End of explanation
"""


from pgmpy.models import BayesianModel
from pgmpy.factors.discrete.CPD import TabularCPD
from pgmpy.inference import VariableElimination

model = BayesianModel([('C', 'S'), ('C', 'R'), ('S', 'W'), ('R', 'W')])

cpd_c = TabularCPD(variable='C', variable_card=2, values=[[0.5, 0.5]])
cpd_s = TabularCPD(variable='S', evidence=['C'], evidence_card=[2], 
                   variable_card=2,
                   values=[[0.9, 0.5], 
                           [0.1, 0.5]])
cpd_r = TabularCPD(variable='R', evidence=['C'], evidence_card=[2],
                   variable_card=2,
                   values=[[0.2, 0.8],
                           [0.8, 0.2]])
cpd_w = TabularCPD(variable='W', evidence=['S', 'R'], evidence_card=[2,2],
                  variable_card=2,
                  values=[[1, 0.1, 0.1, 0.01],
                          [0, 0.9, 0.9, 0.99]])

model.add_cpds(cpd_c, cpd_r, cpd_s, cpd_w)
model.check_model()

infer = VariableElimination(model)
print(infer.query(['W'])['W'])
print(infer.query(['S'], evidence={'W': 1})['S'])
print(infer.query(['R'], evidence={'W': 1})['R'])
print(infer.query(['C'], evidence={'S': 1, 'R': 1})['C'])
print(infer.query(['C'])['C'])

"""
Explanation: 1. Probabilistički grafički modeli -- Bayesove mreže
Ovaj zadatak bavit će se Bayesovim mrežama, jednim od poznatijih probabilističkih grafičkih modela (probabilistic graphical models; PGM). Za lakše eksperimentiranje koristit ćemo programski paket pgmpy. Molimo Vas da provjerite imate li ovaj paket te da ga instalirate ako ga nemate.
(a)
Prvo ćemo pogledati udžbenički primjer s prskalicom. U ovom primjeru razmatramo Bayesovu mrežu koja modelira zavisnosti između oblačnosti (slučajna varijabla $C$), kiše ($R$), prskalice ($S$) i mokre trave ($W$). U ovom primjeru također pretpostavljamo da već imamo parametre vjerojatnosnih distribucija svih čvorova. Ova mreža prikazana je na sljedećoj slici:

Koristeći paket pgmpy, konstruirajte Bayesovu mrežu iz gornjeg primjera. Zatim, koristeći egzaktno zaključivanje, postavite sljedeće posteriorne upite: $P(w=1)$, $P(s=1|w=1)$, $P(r=1|w=1)$, $P(c=1|s=1, r=1)$ i $P(c=1)$. Provedite zaključivanje na papiru i uvjerite se da ste ispravno konstruirali mrežu. Pomoći će vam službena dokumentacija te primjeri korištenja (npr. ovaj).
End of explanation
"""



print(infer.query(['S'], evidence={'W': 1, 'R': 1})['S'])
print(infer.query(['S'], evidence={'W': 1, 'R': 0})['S'])
print(infer.query(['R'], evidence={'W': 1, 'S': 1})['R'])
print(infer.query(['R'], evidence={'W': 1, 'S': 0})['R'])

"""
Explanation: Q: Koju zajedničku vjerojatnosnu razdiobu ova mreža modelira? Kako tu informaciju očitati iz mreže?
Q: U zadatku koristimo egzaktno zaključivanje. Kako ono radi?
Q: Koja je razlika između posteriornog upita i MAP-upita?
Q: Zašto je vjerojatnost $P(c=1)$ drugačija od $P(c=1|s=1,r=1)$ ako znamo da čvorovi $S$ i $R$ nisu roditelji čvora $C$?
(b)
Efekt objašnjavanja (engl. explaining away) zanimljiv je fenomen u kojem se događa da se dvije varijable "natječu" za objašnjavanje treće. Ovaj fenomen može se primijetiti na gornjoj mreži. U tom se slučaju varijable prskalice ($S$) i kiše ($R$) "natječu" za objašnjavanje mokre trave ($W$). Vaš zadatak je pokazati da se fenomen zaista događa.
End of explanation
"""


model.is_active_trail('C','W')

"""
Explanation: Q: Kako biste svojim riječima opisali ovaj fenomen, koristeći se ovim primjerom?
(c)
Koristeći BayesianModel.is_active_trail provjerite jesu li varijable oblačnosti ($C$) i mokre trave ($W$) uvjetno nezavisne. Što mora vrijediti kako bi te dvije varijable bile uvjetno nezavisne? Provjerite korištenjem iste funkcije.
End of explanation
"""


from sklearn.model_selection import train_test_split
spam_X, spam_y = mlutils.load_SMS_dataset('./spam.csv')

spam_X_train, spam_X_test, spam_y_train, spam_y_test = \
    train_test_split(spam_X, spam_y, train_size=0.7, test_size=0.3, random_state=69)

"""
Explanation: Q: Kako možemo na temelju grafa saznati koje dvije varijable su, uz neka opažanja, uvjetno nezavisne?
Q: Zašto bismo uopće htjeli znati koje su varijable u mreži uvjetno nezavisne?
2. Vrednovanje modela (klasifikatora)
Kako bismo se uvjerili koliko naš naučeni model zapravo dobro radi, nužno je provesti evaluaciju modela. Ovaj korak od presudne je važnosti u svim primjenama strojnog učenja, pa je stoga bitno znati provesti evaluaciju na ispravan način.
Vrednovat ćemo modele na stvarnom skupu podataka SMS Spam Collection [1], koji se sastoji od 5,574 SMS-poruka klasificiranih u dvije klase: spam (oznaka: spam) i ne-spam (oznaka: ham). Ako već niste, preuzmite skup podataka s poveznice ili sa stranice kolegija i stavite ga u radni direktorij (otpakirajte arhivu i preimenujte datoteku u spam.csv po potrebi). Sljedeći komad kôda učitava skup podataka i dijeli ga na podskupove za učenje i testiranje.
[1] Almeida, T.A., GÃmez Hidalgo, J.M., Yamakami, A. Contributions to the Study of SMS Spam Filtering: New Collection and Results. Proceedings of the 2011 ACM Symposium on Document Engineering (DOCENG'11), Mountain View, CA, USA, 2011.
End of explanation
"""


from sklearn.pipeline import Pipeline
from sklearn.feature_extraction.text import TfidfVectorizer
from sklearn.decomposition import TruncatedSVD
from sklearn.preprocessing import Normalizer
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import accuracy_score

"""
Explanation: (a)
Prije nego što krenemo u vrednovanje modela za klasifikaciju spama, upoznat ćete se s jednostavnijom apstrakcijom cjelokupnog procesa učenja modela u biblioteci scikit-learn. Ovo je korisno zato što se učenje modela često sastoji od mnoštva koraka prije sâmog pozivanja magične funkcije fit: ekstrakcije podataka, ekstrakcije značajki, standardizacije, skaliranja, nadopunjavanjem nedostajućih vrijednosti i slično. 
U "standardnom pristupu", ovo se svodi na pozamašan broj linija kôda u kojoj konstantno proslijeđujemo podatke iz jednog koraka u sljedeći, tvoreći pritom cjevovod izvođenja. Osim nepreglednosti, ovakav pristup je često i sklon pogreškama, s obzirom na to da je dosta jednostavno proslijediti pogrešan skup podataka i ne dobiti pogrešku pri izvođenju kôda. Stoga je u biblioteci scikit-learn uveden razred pipeline.Pipeline. Kroz ovaj razred, svi potrebni koraci učenja mogu se apstrahirati iza jednog cjevovoda, koji je opet zapravo model s fit i predict funkcijama.
U ovom zadatku ćete napraviti samo jednostavni cjevovod modela za klasifikaciju teksta, koji se sastoji od pretvorbe teksta u vektorsku reprezentaciju vreće riječi s TF-IDF-težinama, redukcije dimenzionalnosti pomoću krnje dekompozicije singularnih vrijednosti, normalizacije, te konačno logističke regresije.
NB: Nije sasvim nužno znati kako rade ovi razredi pomoću kojih dolazimo do konačnih značajki, ali preporučamo da ih proučite ako vas zanima (posebice ako vas zanima obrada prirodnog jezika).
End of explanation
"""


# TF-IDF
vectorizer = TfidfVectorizer(stop_words="english", ngram_range=(1, 2), max_features=500)
spam_X_feat_train = vectorizer.fit_transform(spam_X_train)
# Smanjenje dimenzionalnosti
reducer = TruncatedSVD(n_components=300, random_state=69)
spam_X_feat_train = reducer.fit_transform(spam_X_feat_train)
# Normaliziranje
normalizer = Normalizer()
spam_X_feat_train = normalizer.fit_transform(spam_X_feat_train)
# NB
clf = LogisticRegression(solver='lbfgs')
clf.fit(spam_X_feat_train, spam_y_train)

# I sada ponovno sve ovo za testne podatke.
spam_X_feat_test = vectorizer.transform(spam_X_test)
spam_X_feat_test = reducer.transform(spam_X_feat_test)
spam_X_feat_test = normalizer.transform(spam_X_feat_test)

print(accuracy_score(spam_y_test, clf.predict(spam_X_feat_test)))

x_test123 = ["You were selected for a green card, apply here for only 50 USD!!!",
         "Hey, what are you doing later? Want to grab a cup of coffee?"]
x_test = vectorizer.transform(x_test123)
x_test = reducer.transform(x_test)
x_test = normalizer.transform(x_test)
print(clf.predict(x_test))

"""
Explanation: Prvo, prilažemo kôd koji to radi "standardnim pristupom":
End of explanation
"""



vectorizer = TfidfVectorizer(stop_words="english", ngram_range=(1, 2), max_features=500)
reducer = TruncatedSVD(n_components=300, random_state=69)
normalizer = Normalizer()
clf = LogisticRegression(solver='lbfgs')

pipeline = Pipeline([('vectorizer', vectorizer), ('reducer', reducer), ('normalizer', normalizer), ('clf', clf)])

pipeline.fit(spam_X_train, spam_y_train)
print(accuracy_score(spam_y_test, pipeline.predict(spam_X_test)))
print(pipeline.predict(x_test123))

"""
Explanation: Vaš zadatak izvesti je dani kôd korištenjem cjevovoda. Proučite razred pipeline.Pipeline.
NB Ne treba vam više od svega nekoliko naredbi.
End of explanation
"""


from sklearn.metrics import classification_report, accuracy_score

print(classification_report(y_pred=pipeline.predict(spam_X_test), y_true=spam_y_test))

"""
Explanation: (b)
U prošlom smo podzadatku ispisali točnost našeg modela. Ako želimo vidjeti koliko je naš model dobar po ostalim metrikama, možemo iskoristiti bilo koju funkciju iz paketa metrics. Poslužite se funkcijom metrics.classification_report, koja ispisuje vrijednosti najčešćih metrika. (Obavezno koristite naredbu print kako ne biste izgubili format izlaza funkcije.) Ispišite ponovno točnost za usporedbu.
End of explanation
"""


from sklearn.dummy import DummyClassifier
rando = DummyClassifier(strategy='uniform')

pipeline = Pipeline([('vectorizer', vectorizer), ('reducer', reducer), ('normalizer', normalizer), ('clf', rando)])

pipeline.fit(spam_X_train, spam_y_train)
print(classification_report(spam_y_test, pipeline.predict(spam_X_test)))

mfc = DummyClassifier(strategy='most_frequent')

pipeline = Pipeline([('vectorizer', vectorizer), ('reducer', reducer), ('normalizer', normalizer), ('clf', mfc)])

pipeline.fit(spam_X_train, spam_y_train)
print(classification_report(spam_y_test, pipeline.predict(spam_X_test)))


"""
Explanation: Potreba za drugim metrikama osim točnosti može se vidjeti pri korištenju nekih osnovnih modela (engl. baselines). Možda najjednostavniji model takvog tipa je model koji svrstava sve primjere u većinsku klasu (engl. most frequent class; MFC) ili označuje testne primjere nasumično (engl. random). Proučite razred dummy.DummyClassifier i pomoću njega stvorite spomenute osnovne klasifikatore. Opet ćete trebati iskoristiti cjevovod kako biste došli do vektorskog oblika ulaznih primjera, makar ovi osnovni klasifikatori koriste samo oznake pri predikciji.
End of explanation
"""


from sklearn.model_selection import cross_val_score, KFold

kf = KFold(n_splits=5)


vectorizer = TfidfVectorizer(stop_words="english", ngram_range=(1, 2), max_features=500)
reducer = TruncatedSVD(n_components=300, random_state=69)
normalizer = Normalizer()
clf = LogisticRegression(solver='lbfgs')

pipeline = Pipeline([('vectorizer', vectorizer), ('reducer', reducer), ('normalizer', normalizer), ('clf', clf)])

for train_index, test_index in kf.split(spam_X):
    X_train, X_test = spam_X[train_index], spam_X[test_index]
    y_train, y_test = spam_y[train_index], spam_y[test_index]
    
    pipeline.fit(X_train, y_train)

    print()
    print(cross_val_score(estimator=pipeline, X=X_test, y=y_test, cv=5))
    print(accuracy_score(spam_y_test, pipeline.predict(spam_X_test)))

    

"""
Explanation: Q: Na temelju ovog primjera objasnite zašto točnost nije uvijek prikladna metrika.
Q: Zašto koristimo F1-mjeru?
(c)
Međutim, provjera za kakvom smo posegli u prošlom podzadatku nije robusna. Stoga se u strojnom učenju obično koristi k-struka unakrsna provjera. Proučite razred model_selection.KFold i funkciju model_selection.cross_val_score te izračunajte procjenu pogreške na cijelom skupu podataka koristeći peterostruku unakrsnu provjeru. 
NB: Vaš model je sada cjevovod koji sadrži čitavo pretprocesiranje. Također, u nastavku ćemo se ograničiti na točnost, ali ovi postupci vrijede za sve metrike.
End of explanation
"""


from sklearn.model_selection import GridSearchCV

param_grid = {
    'vectorizer__max_features': [500, 1000],
    'reducer__n_components': [100, 200, 300]
}


vectorizer = TfidfVectorizer(stop_words="english", ngram_range=(1, 2), max_features=500)
reducer = TruncatedSVD()
normalizer = Normalizer()
clf = LogisticRegression(solver='lbfgs')

pipeline = Pipeline([('vectorizer', vectorizer), ('reducer', reducer), ('normalizer', normalizer), ('clf', clf)])


f = GridSearchCV(pipeline, param_grid, cv = 3).fit(spam_X_train)
print(f.best_estimator_.fit(spam_X_train, spam_y_train).predict(spam_X_test))

"""
Explanation: Q: Zašto "obična" unakrsna provjera nije dovoljno robusna?
Q: Što je to stratificirana k-struka unakrsna provjera? Zašto ju često koristimo?  
(d)
Gornja procjena pogreške je u redu ako imamo već imamo model (bez ili s fiksiranim hiperparametrima). Međutim, mi želimo koristiti model koji ima optimalne vrijednosti hiperparametara te ih je stoga potrebno optimirati korištenjem pretraživanja po rešetci (engl. grid search). Očekivano, biblioteka scikit-learn već ima ovu funkcionalnost u razredu model_selection.GridSearchCV. Jedina razlika vaše implementacije iz prošlih vježbi (npr. kod SVM-a) i ove jest ta da ova koristi k-struku unakrsnu provjeru.
Prije optimizacije vrijednosti hiperparametara, očigledno moramo definirati i samu rešetku vrijednosti hiperparametara. Proučite kako se definira ista kroz rječnik u primjeru. 
Proučite spomenuti razred te pomoću njega pronađite i ispišite najbolje vrijednosti hiperparametara cjevovoda iz podzadatka (a): max_features $\in {500, 1000}$ i n_components $\in { 100, 200, 300 }$ korištenjem pretraživanja po rešetci na skupu za učenje ($k=3$, kako bi išlo malo brže).
End of explanation
"""


from sklearn.model_selection import GridSearchCV, KFold

def nested_kfold_cv(clf, param_grid, X, y, k1=10, k2=3):
    err = []
    kfold = KFold(n_splits=k1, shuffle=True, random_state=42)
    
    for index_train, index_test in kfold.split(X):
        
        X_train, y_train, X_test, y_test = X[index_train], y[index_train], X[index_test], y[index_test]
        f = GridSearchCV(clf, param_grid, cv=k2).fit(X_train, y_train)
        h = f.best_estimator_.fit(X_train, y_train).predict(X_test)
        
        err.append(zero_one_loss(y_test, h))
        
    return err

"""
Explanation: Q: Koja se metrika optimira pri ovoj optimizaciji?
Q: Kako biste odredili broj preklopa $k$?  
(e)
Ako želimo procijeniti pogrešku, ali pritom i napraviti odabir modela, tada se okrećemo ugniježđenoj k-strukoj unakrsnoj provjeri (engl. nested k-fold cross validation). U ovom zadatku ćete ju sami implementirati.
Implementirajte funkciju nested_kfold_cv(clf, param_grid, X, y, k1, k2) koja provodi ugniježđenu k-struku unakrsnu provjeru. Argument clf predstavlja vaš klasifikator, param_grid rječnik vrijednosti hiperparametara (isto kao i u podzadatku (d)), X i y označeni skup podataka, a k1 i k2 broj preklopa u vanjskoj, odnosno unutarnjoj petlji. Poslužite se razredima model_selection.GridSearchCV i  model_selection.KFold.
Funkcija vraća listu pogrešaka kroz preklope vanjske petlje.
End of explanation
"""


np.random.seed(1337)
C1_scores_5folds = np.random.normal(78, 4, 5)
C2_scores_5folds = np.random.normal(81, 2, 5)

C1_scores_10folds = np.random.normal(78, 4, 10)
C2_scores_10folds = np.random.normal(81, 2, 10)

C1_scores_50folds = np.random.normal(78, 4, 50)
C2_scores_50folds = np.random.normal(81, 2, 50)

"""
Explanation: Q: Kako biste odabrali koji su hiperparametri generalno najbolji, a ne samo u svakoj pojedinačnoj unutarnjoj petlji?
Q: Čemu u konačnici odgovara procjena generalizacijske pogreške?
(f)
Scenarij koji nas najviše zanima jest usporedba dvaju klasifikatora, odnosno, je li jedan od njih zaista bolji od drugog. Jedini način kako to možemo zaista potvrditi jest statističkom testom, u našem slučaju uparenim t-testom. Njime ćemo se baviti u ovom zadatku.
Radi bržeg izvođenja, umjetno ćemo generirati podatke koji odgovaraju pogreškama kroz vanjske preklope dvaju klasifikatora (ono što bi vratila funkcija nested_kfold_cv):
End of explanation
"""


from scipy.stats import ttest_rel


print(ttest_rel(C1_scores_5folds, C2_scores_5folds).pvalue)
print(ttest_rel(C1_scores_10folds, C2_scores_10folds).pvalue)
print(ttest_rel(C1_scores_50folds, C2_scores_50folds).pvalue)


"""
Explanation: Iskoristite ugrađenu funkciju scipy.stats.ttest_rel za provedbu uparenog t-testa i provjerite koji od ova modela je bolji kada se koristi 5, 10 i 50 preklopa.
End of explanation
"""


from sklearn.datasets import make_blobs

Xp, yp = make_blobs(n_samples=300, n_features=2, centers=[[0, 0], [3, 2.5], [0, 4]], 
                    cluster_std=[0.45, 0.3, 0.45], random_state=96)
plt.scatter(Xp[:,0], Xp[:,1], c=yp, cmap=plt.get_cmap("cool"), s=20)

"""
Explanation: Q: Koju hipotezu $H_0$ i alternativnu hipotezu $H_1$ testiramo ovim testom?
Q: Koja pretpostavka na vjerojatnosnu razdiobu primjera je napravljena u gornjem testu? Je li ona opravdana?
Q: Koji je model u konačnici bolji i je li ta prednost značajna uz $\alpha = 0.05$?  
3. Grupiranje
U ovom zadatku ćete se upoznati s algoritmom k-sredina (engl. k-means), njegovim glavnim nedostatcima te pretpostavkama. Također ćete isprobati i drugi algoritam grupiranja: model Gaussovih mješavina (engl. Gaussian mixture model).
(a)
Jedan od nedostataka algoritma k-sredina jest taj što unaprijed zahtjeva broj grupa ($K$) u koje će grupirati podatke. Ta informacija nam često nije dostupna (kao što nam nisu dostupne ni oznake primjera) te je stoga potrebno nekako izabrati najbolju vrijednost hiperparametra $K$. Jedan od naivnijih pristupa jest metoda lakta/koljena (engl. elbow method) koju ćete isprobati u ovom zadatku.
U svojim rješenjima koristite ugrađenu implementaciju algoritma k-sredina, dostupnoj u razredu cluster.KMeans. 
NB: Kriterijska funkcija algoritma k-sredina još se i naziva inercijom (engl. inertia). Za naučeni model, vrijednost kriterijske funkcije $J$ dostupna je kroz razredni atribut inertia_.
End of explanation
"""


Ks = range(1,16)

from sklearn.cluster import KMeans

Js = []

for K in Ks:
    J = KMeans(n_clusters=K).fit(Xp).inertia_
    Js.append(J)
    
plot(Ks, Js)

"""
Explanation: Iskoristite skup podataka Xp dan gore. Isprobajte vrijednosti hiperparametra $K$ iz $[0,1,\ldots,15]$. Ne trebate dirati nikakve hiperparametre modela osim $K$. Iscrtajte krivulju od $J$ u ovisnosti o broju grupa $K$. Metodom lakta/koljena odredite vrijednost hiperparametra $K$.
End of explanation
"""


from mlutils import plot_silhouette

for K in [2,3,5]:
    plot_silhouette(K, Xp)

"""
Explanation: Q: Koju biste vrijednost hiperparametra $K$ izabrali na temelju ovog grafa? Zašto? Je li taj odabir optimalan? Kako to znate?  
Q: Je li ova metoda robusna?
Q: Možemo li izabrati onaj $K$ koji minimizira pogrešku $J$? Objasnite.
(b)
Odabir vrijednosti hiperparametra $K$ može se obaviti na mnoštvo načina. Pored metode lakta/koljena, moguće je isto ostvariti i analizom siluete (engl. silhouette analysis). Za to smo pripremili funkciju mlutils.plot_silhouette koja za dani broj grupa i podatke iscrtava prosječnu vrijednost koeficijenta siluete i vrijednost koeficijenta svakog primjera (kroz grupe). 
Vaš je zadatak isprobati različite vrijednosti hiperparametra $K$, $K \in {2, 3, 5}$ i na temelju dobivenih grafova odlučiti se za optimalan $K$.
End of explanation
"""


from sklearn.datasets import make_blobs

X1, y1 = make_blobs(n_samples=1000, n_features=2, centers=[[0, 0], [1.3, 1.3]], cluster_std=[0.15, 0.5], random_state=96)
plt.scatter(X1[:,0], X1[:,1], c=y1, cmap=plt.get_cmap("cool"), s=20)

"""
Explanation: Q: Kako biste se gledajući ove slike odlučili za $K$?
Q: Koji su problemi ovog pristupa?
(c)
U ovom i sljedećim podzadatcima fokusirat ćemo se na temeljne pretpostavke algoritma k-sredina te što se događa ako te pretpostavke nisu zadovoljene. Dodatno, isprobat ćemo i grupiranje modelom Gaussovih mješavina (engl. Gaussian Mixture Models; GMM) koji ne nema neke od tih pretpostavki.
Prvo, krenite od podataka X1, koji su generirani korištenjem funkcije datasets.make_blobs, koja stvara grupe podataka pomoću izotropskih Gaussovih distribucija.
End of explanation
"""


k_means = KMeans(n_clusters=2).fit(X1)

plt.scatter(X1[:,0], X1[:,1], c=k_means.predict(X1), cmap=plt.get_cmap("cool"), s=20);

"""
Explanation: Naučite model k-sredina (idealno pretpostavljajući $K=2$) na gornjim podatcima i prikažite dobiveno grupiranje (proučite funkciju scatter, posebice argument c).
End of explanation
"""


from sklearn.datasets import make_circles

X2, y2 = make_circles(n_samples=1000, noise=0.15, factor=0.05, random_state=96)
plt.scatter(X2[:,0], X2[:,1], c=y2, cmap=plt.get_cmap("cool"), s=20)

"""
Explanation: Q: Što se dogodilo? Koja je pretpostavka algoritma k-sredina ovdje narušena?
Q: Što biste morali osigurati kako bi algoritam pronašao ispravne grupe?
(d)
Isprobajte algoritam k-sredina na podatcima generiranim korištenjem funkcije datasets.make_circles, koja stvara dvije grupe podataka tako da je jedna unutar druge.
End of explanation
"""


k_means = KMeans(n_clusters=2).fit(X2)

plt.scatter(X2[:,0], X2[:,1], c=k_means.predict(X2), cmap=plt.get_cmap("cool"), s=20);

"""
Explanation: Ponovno, naučite model k-sredina (idealno pretpostavljajući $K=2$) na gornjim podatcima i prikažite dobiveno grupiranje (proučite funkciju scatter, posebice argument c).
End of explanation
"""


X31, y31 = make_blobs(n_samples=1000, n_features=2, centers=[[0, 0]], cluster_std=[0.2], random_state=69)
X32, y32 = make_blobs(n_samples=50, n_features=2, centers=[[0.7, 0.5]], cluster_std=[0.15], random_state=69)
X33, y33 = make_blobs(n_samples=600, n_features=2, centers=[[0.8, -0.4]], cluster_std=[0.2], random_state=69)
plt.scatter(X31[:,0], X31[:,1], c="#00FFFF", s=20)
plt.scatter(X32[:,0], X32[:,1], c="#F400F4", s=20)
plt.scatter(X33[:,0], X33[:,1], c="#8975FF", s=20)

# Just join all the groups in a single X.
X3 = np.vstack([X31, X32, X33])
y3 = np.hstack([y31, y32, y33])

"""
Explanation: Q: Što se dogodilo? Koja je pretpostavka algoritma k-sredina ovdje narušena?
Q: Što biste morali osigurati kako bi algoritam pronašao ispravne grupe?
(e)
Završno, isprobat ćemo algoritam na sljedećem umjetno stvorenom skupu podataka:
End of explanation
"""


k_means = KMeans(n_clusters=3).fit(X3)

plt.scatter(X3[:,0], X3[:,1], c=k_means.predict(X3), cmap=plt.get_cmap("cool"), s=20);

"""
Explanation: Ponovno, naučite model k-sredina (ovaj put idealno pretpostavljajući $K=3$) na gornjim podatcima i prikažite dobiveno grupiranje (proučite funkciju scatter, posebice argument c).
End of explanation
"""


from sklearn.mixture import GaussianMixture


gauss_X1 = GaussianMixture(n_components = 2).fit(X1)
gauss_X2 = GaussianMixture(n_components = 2).fit(X2)
gauss_X3 = GaussianMixture(n_components = 3).fit(X3)

h_X1 = gauss_X1.predict(X1)
h_X2 = gauss_X2.predict(X2)
h_X3 = gauss_X3.predict(X3)

figure(1)
plt.scatter(X1[:,0], X1[:,1], c=h_X1, cmap=plt.get_cmap("cool"), s=20);

figure(2)
plt.scatter(X2[:,0], X2[:,1], c=h_X2, cmap=plt.get_cmap("cool"), s=20);

figure(3)
plt.scatter(X3[:,0], X3[:,1], c=h_X3, cmap=plt.get_cmap("cool"), s=20);


"""
Explanation: Q: Što se dogodilo? Koja je pretpostavka algoritma k-sredina ovdje narušena?
Q: Što biste morali osigurati kako bi algoritam pronašao ispravne grupe?
(f)
Sada kada ste se upoznali s ograničenjima algoritma k-sredina, isprobat ćete grupiranje modelom mješavine Gaussa (Gaussian Mixture Models; GMM), koji je generalizacija algoritma k-sredina (odnosno, algoritam k-sredina specijalizacija je GMM-a). Implementacija ovog modela dostupna je u mixture.GaussianMixture. Isprobajte ovaj model (s istim pretpostavkama o broju grupa) na podacima iz podzadataka (c)-(e). Ne morate mijenjati nikakve hiperparametre ni postavke osim broja komponenti.
End of explanation
"""


import itertools as it
from scipy.special import comb

def rand_index_score(y_gold, y_predict):
    N = len(y_gold)
    grupa1 = ([ y for i, y in enumerate(y_gold) if y_predict[i] == 0])
    grupa2 = ([ y for i, y in enumerate(y_gold) if y_predict[i] == 1])
    grupa3 = ([ y for i, y in enumerate(y_gold) if y_predict[i] == 2])
    
    n = [[len([y for y in g if y == i])
            for i in [0,1,2]]
            for g in [grupa1, grupa2, grupa3]]
    
    a = sum([(comb(nnn, 2)) for nn in n for nnn in nn])
    
    b = n[0][0] * (n[1][1] + n[1][2] + n[2][1] + n[2][2]) + \
        n[0][1] * (n[1][0] + n[1][2] + n[2][0] + n[2][2]) + \
        n[0][2] * (n[1][0] + n[1][1] + n[2][0] + n[2][1]) + \
        n[1][0] * (n[2][1] + n[2][2]) + \
        n[1][1] * (n[2][0] + n[2][2]) + \
        n[1][2] * (n[2][0] + n[2][1])
        
    
    return (a+b) / comb(N,2)

"""
Explanation: (g)
Kako vrednovati točnost modela grupiranja ako imamo stvarne oznake svih primjera (a u našem slučaju imamo, jer smo mi ti koji smo generirali podatke)? Često korištena mjera jest Randov indeks koji je zapravo pandan točnosti u zadatcima klasifikacije. Implementirajte funkciju rand_index_score(y_gold, y_predict) koja ga računa. Funkcija prima dva argumenta: listu stvarnih grupa kojima primjeri pripadaju (y_gold) i listu predviđenih grupa (y_predict). Dobro će vam doći funkcija itertools.combinations.
End of explanation
"""


y_pred = KMeans(n_clusters=3).fit(Xp).predict(Xp)

rand_index_score(yp, y_pred)

"""
Explanation: Q: Zašto je Randov indeks pandan točnosti u klasifikacijskim problemima?
Q: Koji su glavni problemi ove metrike?
Q: Kako vrednovati kvalitetu grupiranja ako nenamo stvarne oznake primjera? Je li to uopće moguće?
End of explanation
"""


