# Učitaj osnovne biblioteke...
import numpy as np
import sklearn
import matplotlib.pyplot as plt
import scipy as sp
%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 1: Regresija
Verzija: 1.1
Zadnji put ažurirano: 12. listopada 2018.
(c) 2015-2018 Jan Šnajder, Domagoj Alagić, Mladen Karan 
Objavljeno: 12. listopada 2018.
Rok za predaju: 22. listopada 2018. u 07:00h

Upute
Prva laboratorijska vježba sastoji se od osam zadataka. 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
"""


X = np.array([[0],[1],[2],[4]])
y = np.array([4,1,2,5])

"""
Explanation: Zadatci
1. Jednostavna regresija
Zadan je skup primjera $\mathcal{D}={(x^{(i)},y^{(i)})}_{i=1}^4 = {(0,4),(1,1),(2,2),(4,5)}$. Primjere predstavite matrixom $\mathbf{X}$ dimenzija $N\times n$ (u ovom slučaju $4\times 1$) i vektorom oznaka $\textbf{y}$, dimenzija $N\times 1$ (u ovom slučaju $4\times 1$), na sljedeći način:
End of explanation
"""


from sklearn.preprocessing import PolynomialFeatures
Phi = PolynomialFeatures(1, False, True).fit_transform(X)
print(Phi)

"""
Explanation: (a)
Proučite funkciju PolynomialFeatures iz biblioteke sklearn i upotrijebite je za generiranje matrice dizajna $\mathbf{\Phi}$ koja ne koristi preslikavanje u prostor više dimenzije (samo će svakom primjeru biti dodane dummy jedinice; $m=n+1$).
End of explanation
"""


from numpy import linalg

w = np.dot(np.dot(np.linalg.inv(np.dot(np.transpose(Phi), Phi)), np.transpose(Phi)), y)
print(w)

w2 = np.dot(np.linalg.pinv(Phi), y)
print(w2)

"""
Explanation: (b)
Upoznajte se s modulom linalg. Izračunajte težine $\mathbf{w}$ modela linearne regresije kao $\mathbf{w}=(\mathbf{\Phi}^\intercal\mathbf{\Phi})^{-1}\mathbf{\Phi}^\intercal\mathbf{y}$. Zatim se uvjerite da isti rezultat možete dobiti izračunom pseudoinverza $\mathbf{\Phi}^+$ matrice dizajna, tj. $\mathbf{w}=\mathbf{\Phi}^+\mathbf{y}$, korištenjem funkcije pinv.
End of explanation
"""


from sklearn.metrics import mean_squared_error

h = np.dot(Phi, w)
print (h)

error = mean_squared_error(y, h)
print (error)

plt.plot(X, y, '+', X, h, linewidth = 1)
plt.axis([-3, 6, -1, 7])

"""
Explanation: Radi jasnoće, u nastavku je vektor $\mathbf{x}$ s dodanom dummy jedinicom $x_0=1$ označen kao $\tilde{\mathbf{x}}$.
(c)
Prikažite primjere iz $\mathcal{D}$ i funkciju $h(\tilde{\mathbf{x}})=\mathbf{w}^\intercal\tilde{\mathbf{x}}$. Izračunajte pogrešku učenja prema izrazu $E(h|\mathcal{D})=\frac{1}{2}\sum_{i=1}^N(\tilde{\mathbf{y}}^{(i)} - h(\tilde{\mathbf{x}}))^2$. Možete koristiti funkciju srednje kvadratne pogreške mean_squared_error iz modula sklearn.metrics.
Q: Gore definirana funkcija pogreške $E(h|\mathcal{D})$ i funkcija srednje kvadratne pogreške nisu posve identične. U čemu je razlika? Koja je "realnija"?
End of explanation
"""


try:
    np.dot(np.linalg.inv(Phi), y)
except LinAlgError as err:
    print(err)

"""
Explanation: (d)
Uvjerite se da za primjere iz $\mathcal{D}$ težine $\mathbf{w}$ ne možemo naći rješavanjem sustava $\mathbf{w}=\mathbf{\Phi}^{-1}\mathbf{y}$, već da nam doista treba pseudoinverz.
Q: Zašto je to slučaj? Bi li se problem mogao riješiti preslikavanjem primjera u višu dimenziju? Ako da, bi li to uvijek funkcioniralo, neovisno o skupu primjera $\mathcal{D}$? Pokažite na primjeru.
End of explanation
"""


from sklearn.linear_model import LinearRegression

lr = LinearRegression().fit(Phi, y)

w2 = [lr.intercept_, lr.coef_[1]]
h2 = lr.predict(Phi)
error2 = mean_squared_error(y, h)

print ('staro: ')
print (w)
print (h)
print (error)

print('novo: ')
print (w2)
print (h2)
print (error2)

"""
Explanation: (e)
Proučite klasu LinearRegression iz modula sklearn.linear_model. Uvjerite se da su težine koje izračunava ta funkcija (dostupne pomoću atributa coef_ i intercept_) jednake onima koje ste izračunali gore. Izračunajte predikcije modela (metoda predict) i uvjerite se da je pogreška učenja identična onoj koju ste ranije izračunali.
End of explanation
"""


from numpy.random import normal
def make_labels(X, f, noise=0) :
    return map(lambda x : f(x) + (normal(0,noise) if noise>0 else 0), X)

def make_instances(x1, x2, N) :
    return sp.array([np.array([x]) for x in np.linspace(x1,x2,N)])

N = 50
sigma = 200
fun = lambda x :5 + x - 2*x**2 - 5*x**3
x = make_instances(-5, 5, N)
y = list(make_labels(x, fun, sigma))

y6a = y
x6a = x

"""
Explanation: 2. Polinomijalna regresija i utjecaj šuma
(a)
Razmotrimo sada regresiju na većem broju primjera. Koristite funkciju make_labels(X, f, noise=0) koja uzima matricu neoznačenih primjera $\mathbf{X}{N\times n}$ te generira vektor njihovih oznaka $\mathbf{y}{N\times 1}$. Oznake se generiraju kao $y^{(i)} = f(x^{(i)})+\mathcal{N}(0,\sigma^2)$, gdje je $f:\mathbb{R}^n\to\mathbb{R}$ stvarna funkcija koja je generirala podatke (koja nam je u stvarnosti nepoznata), a $\sigma$ je standardna devijacija Gaussovog šuma, definirana parametrom noise. Za generiranje šuma koristi se funkcija numpy.random.normal. 
Generirajte skup za učenje od $N=50$ primjera uniformno distribuiranih u intervalu $[-5,5]$ pomoću funkcije $f(x) = 5 + x -2 x^2 -5 x^3$ uz šum  $\sigma=200$:
End of explanation
"""


plt.figure(figsize=(10, 5))
plt.plot(x, fun(x), 'r', linewidth = 1)
plt.scatter(x, y)


"""
Explanation: Prikažite taj skup funkcijom scatter.
End of explanation
"""


from sklearn.preprocessing import PolynomialFeatures

Phi = PolynomialFeatures(3).fit_transform(x.reshape(-1, 1))
w = np.dot(np.linalg.pinv(Phi), y)
h = np.dot(Phi, w)
error = mean_squared_error(y, h)
print(error)

plt.figure(figsize=(10,5))
plt.scatter(x, y)
plt.plot(x, h, 'r', linewidth=1)

"""
Explanation: (b)
Trenirajte model polinomijalne regresije stupnja $d=3$. Na istom grafikonu prikažite naučeni model $h(\mathbf{x})=\mathbf{w}^\intercal\tilde{\mathbf{x}}$ i primjere za učenje. Izračunajte pogrešku učenja modela.
End of explanation
"""


Phi_d = []; 
w_d = []; 
h_d = [];
err_d = [];

d = [1, 3, 5, 10, 20]

for i in d:
    Phi_d.append(PolynomialFeatures(i).fit_transform(x.reshape(-1,1)))
    
for i in range(0, len(d)):
    w_d.insert(i, np.dot(np.linalg.pinv(Phi_d[i]), y))
    h_d.insert(i, np.dot(Phi_d[i], w_d[i]))

for i in range(0, len(d)):
    err_d.insert(i, mean_squared_error(y, h_d[i]))
    print (str(d[i]) + ': ' + str(err_d[i]))


fig = plt.figure(figsize=(15, 20))
fig.subplots_adjust(wspace=0.2) 

for i in range(0, len(d)):
    ax = fig.add_subplot(5, 2, i+1)
    ax.scatter(x, y);
    ax.plot(x, h_d[i], 'r', linewidth = 1)
                   

"""
Explanation: 3. Odabir modela
(a)
Na skupu podataka iz zadatka 2 trenirajte pet modela linearne regresije $\mathcal{H}_d$ različite složenosti, gdje je $d$ stupanj polinoma, $d\in{1,3,5,10,20}$. Prikažite na istome grafikonu skup za učenje i funkcije $h_d(\mathbf{x})$ za svih pet modela (preporučujemo koristiti plot unutar for petlje). Izračunajte pogrešku učenja svakog od modela.
Q: Koji model ima najmanju pogrešku učenja i zašto?
End of explanation
"""


from sklearn.model_selection import train_test_split

X_train, X_test, y_train, y_test = train_test_split(x, y, test_size = 0.5)
err_train = [];
err_test = [];
d = range(0, 20)

for i in d:
    Phi_train = PolynomialFeatures(i).fit_transform(X_train.reshape(-1, 1))
    Phi_test = PolynomialFeatures(i).fit_transform(X_test.reshape(-1, 1))
    w_train = np.dot(np.linalg.pinv(Phi_train), y_train)
    h_train = np.dot(Phi_train, w_train)
    h_test = np.dot(Phi_test, w_train)
    
    err_train.insert(i, np.log(mean_squared_error(y_train, h_train)))
    err_test.insert(i, np.log(mean_squared_error(y_test, h_test)))

plt.figure(figsize=(10,5))
plt.plot(d, err_train, d, err_test)
plt.grid()

"""
Explanation: (b)
Razdvojite skup primjera iz zadatka 2 pomoću funkcije cross_validation.train_test_split na skup za učenja i skup za ispitivanje u omjeru 1:1. Prikažite na jednom grafikonu pogrešku učenja i ispitnu pogrešku za modele polinomijalne regresije $\mathcal{H}_d$, sa stupnjem polinoma $d$ u rasponu $d\in [1,2,\ldots,20]$. Radi preciznosti, funkcije $h(\mathbf{x})$ iscrtajte na cijelom skupu primjera (ali pogrešku generalizacije računajte, naravno, samo na ispitnome skupu). Budući da kvadratna pogreška brzo raste za veće stupnjeve polinoma, umjesto da iscrtate izravno iznose pogrešaka, iscrtajte njihove logaritme.
NB: Podjela na skupa za učenje i skup za ispitivanje mora za svih pet modela biti identična.
Q: Je li rezultat u skladu s očekivanjima? Koji biste model odabrali i zašto?
Q: Pokrenite iscrtavanje više puta. U čemu je problem? Bi li problem bio jednako izražen kad bismo imali više primjera? Zašto?
End of explanation
"""


N2 = [100, 200, 1000];
sigma = [100, 200, 500];

X_train4c_temp = [];
X_test4c_temp = [];
y_train4c_temp = [];
y_test4c_temp = [];

x_tmp = np.linspace(-5, 5, 1000);
X_train, X_test = train_test_split(x_tmp, test_size = 0.5)

for i in range(0, 3):
    
    y_tmp_train = list(make_labels(X_train, fun, sigma[i]))
    y_tmp_test = list(make_labels(X_test, fun, sigma[i]))
    for j in range(0,3):   
        X_train4c_temp.append(X_train[0:int(N2[j]/2)])
        X_test4c_temp.append(X_test[0:int(N2[j]/2)])
        y_train4c_temp.append(y_tmp_train[0:int(N2[j]/2)])
        y_test4c_temp.append(y_tmp_test[0:int(N2[j]/2)])

    
err_tr = [];
err_tst = [];

for i in range(0, 9):
    X_train4c = X_train4c_temp[i]
    X_test4c = X_test4c_temp[i]
    y_train4c = y_train4c_temp[i]
    y_test4c = y_test4c_temp[i]
    
    err_train4c = [];
    err_test4c = [];
    d4c = range(0, 20)

    for j in d4c:
        Phi_train4c = PolynomialFeatures(j).fit_transform(X_train4c.reshape(-1, 1))
        Phi_test4c = PolynomialFeatures(j).fit_transform(X_test4c.reshape(-1, 1))
        w_train4c = np.dot(np.linalg.pinv(Phi_train4c), y_train4c)
        h_train4c = np.dot(Phi_train4c, w_train4c)
        h_test4c = np.dot(Phi_test4c, w_train4c)
        err_train4c.insert(j, np.log(mean_squared_error(y_train4c, h_train4c)))
        err_test4c.insert(j, np.log(mean_squared_error(y_test4c, h_test4c)))
        
    err_tr.append(err_train4c);
    err_tst.append(err_test4c);
    
    
fig = plt.figure(figsize=(15, 10))
fig.subplots_adjust(wspace=0.2, hspace = 0.35) 

Nn = [100, 200, 1000, 100, 200, 1000, 100, 200, 1000]
sgm = [100, 100, 100, 200, 200, 200, 500, 500, 500]

for i in range(0, 9):    
    ax = fig.add_subplot(3, 3, i+1)
    plt.plot(d, err_tr[i], d, err_tst[i]); grid;
    ax.grid();

"""
Explanation: (c)
Točnost modela ovisi o (1) njegovoj složenosti (stupanj $d$ polinoma), (2) broju primjera $N$, i (3) količini šuma. Kako biste to analizirali, nacrtajte grafikone pogrešaka kao u 3b, ali za sve kombinacija broja primjera $N\in{100,200,1000}$ i količine šuma $\sigma\in{100,200,500}$ (ukupno 9 grafikona). Upotrijebite funkciju subplots kako biste pregledno posložili grafikone u tablicu $3\times 3$. Podatci se generiraju na isti način kao u zadatku 2.
NB: Pobrinite se da svi grafikoni budu generirani nad usporedivim skupovima podataka, na sljedeći način. Generirajte najprije svih 1000 primjera, podijelite ih na skupove za učenje i skupove za ispitivanje (dva skupa od po 500 primjera). Zatim i od skupa za učenje i od skupa za ispitivanje načinite tri različite verzije, svaka s drugačijom količinom šuma (ukupno 2x3=6 verzija podataka). Kako bi simulirali veličinu skupa podataka, od tih dobivenih 6 skupova podataka uzorkujte trećinu, dvije trećine i sve podatke. Time ste dobili 18 skupova podataka -- skup za učenje i za testiranje za svaki od devet grafova.
Q: Jesu li rezultati očekivani? Obrazložite.
End of explanation
"""


lam = [0, 1, 10]
y = np.array([4,1,2,5])

Phi4a = PolynomialFeatures(3).fit_transform(X)
w_L2 = [];

def w_reg(lam): 
    t1 = np.dot(Phi4a.T, Phi4a) + np.dot(lam, np.eye(4))
    t2 = np.dot(np.linalg.inv(t1), Phi4a.T)
    return np.dot(t2, y)

for i in range(0, 3):
    w_L2.insert(i, w_reg(lam[i]))
    print (w_reg(lam[i]))


"""
Explanation: 4. Regularizirana regresija
(a)
U gornjim eksperimentima nismo koristili regularizaciju. Vratimo se najprije na primjer iz zadatka 1. Na primjerima iz tog zadatka izračunajte težine $\mathbf{w}$ za polinomijalni regresijski model stupnja $d=3$ uz L2-regularizaciju (tzv. ridge regression), prema izrazu $\mathbf{w}=(\mathbf{\Phi}^\intercal\mathbf{\Phi}+\lambda\mathbf{I})^{-1}\mathbf{\Phi}^\intercal\mathbf{y}$. Napravite izračun težina za regularizacijske faktore $\lambda=0$, $\lambda=1$ i $\lambda=10$ te usporedite dobivene težine.
Q: Kojih je dimenzija matrica koju treba invertirati?
Q: Po čemu se razlikuju dobivene težine i je li ta razlika očekivana? Obrazložite.
End of explanation
"""


from sklearn.linear_model import Ridge

for i in lam:
    w = [];
    w_L22 = Ridge(alpha = i).fit(Phi4a, y)
    
    w.append(w_L22.intercept_)
    for i in range(0, len(w_L22.coef_[1:])):
        w.append(w_L22.coef_[i])
    
    print (w)

"""
Explanation: (b)
Proučite klasu Ridge iz modula sklearn.linear_model, koja implementira L2-regularizirani regresijski model. Parametar $\alpha$ odgovara parametru $\lambda$. Primijenite model na istim primjerima kao u prethodnom zadatku i ispišite težine $\mathbf{w}$ (atributi coef_ i intercept_).
Q: Jesu li težine identične onima iz zadatka 4a? Ako nisu, objasnite zašto je to tako i kako biste to popravili.
End of explanation
"""


x5a = linspace(-5, 5, 50);
f = (5 + x5a - 2*(x5a**2) - 5*(x5a**3));
y5a = f + normal(0, 200, 50);

lamd = [0, 100]
dd = [2, 10]
h5a = []

for i in lamd:
    for j in dd:
        Phi5a = PolynomialFeatures(j).fit_transform(x5a.reshape(-1,1))
        w_5a = np.dot(np.dot(np.linalg.inv(np.dot(Phi5a.T, Phi5a) + np.dot(i, np.eye(j+1))), Phi5a.T), y5a);
        h_5a = np.dot(Phi5a, w_5a)
        h5a.append(h_5a)
       

lamdd = [0, 0, 100, 100]
ddd = [2, 10, 2, 10]

fig = plt.figure(figsize=(15, 10))
fig.subplots_adjust(wspace=0.2, hspace = 0.2) 

for i in range(0, len(lamdd)):    
    ax = fig.add_subplot(2, 2, i+1)
    plt.plot(x5a, h5a[i], 'r', linewidth = 2)
    plt.scatter(x5a, y5a);

"""
Explanation: 5. Regularizirana polinomijalna regresija
(a)
Vratimo se na slučaj $N=50$ slučajno generiranih primjera iz zadatka 2. Trenirajte modele polinomijalne regresije $\mathcal{H}_{\lambda,d}$ za $\lambda\in{0,100}$ i $d\in{2,10}$ (ukupno četiri modela). Skicirajte pripadne funkcije $h(\mathbf{x})$ i primjere (na jednom grafikonu; preporučujemo koristiti plot unutar for petlje).
Q: Jesu li rezultati očekivani? Obrazložite.
End of explanation
"""


X5a_train, X5a_test, y5a_train, y5a_test = train_test_split(x5a, y5a, test_size = 0.5)
err5a_train = [];
err5a_test = [];
d = 20;
lambda5a = range(0, 50)

for i in lambda5a:
    Phi5a_train = PolynomialFeatures(d).fit_transform(X5a_train.reshape(-1, 1))
    Phi5a_test = PolynomialFeatures(d).fit_transform(X5a_test.reshape(-1, 1))
    w5a_train = np.dot(np.dot(np.linalg.inv(np.dot(Phi5a_train.T, Phi5a_train) + np.dot(i, np.eye(d+1))), Phi5a_train.T), y5a_train);
    h5a_train = np.dot(Phi5a_train, w5a_train)
    h5a_test = np.dot(Phi5a_test, w5a_train)
    
    err5a_train.insert(i, np.log(mean_squared_error(y5a_train, h5a_train)))
    err5a_test.insert(i, np.log(mean_squared_error(y5a_test, h5a_test)))

plt.figure(figsize=(8,4))
plt.plot(lambda5a, err5a_train, lambda5a, err5a_test);
plt.grid(), plt.xlabel('$\lambda$'), plt.ylabel('err');
plt.legend(['ucenje', 'ispitna'], loc='best');

"""
Explanation: (b)
Kao u zadataku 3b, razdvojite primjere na skup za učenje i skup za ispitivanje u omjeru 1:1. Prikažite krivulje logaritama pogreške učenja i ispitne pogreške u ovisnosti za model $\mathcal{H}_{d=20,\lambda}$, podešavajući faktor regularizacije $\lambda$ u rasponu $\lambda\in{0,1,\dots,50}$.
Q: Kojoj strani na grafikonu odgovara područje prenaučenosti, a kojoj podnaučenosti? Zašto?
Q: Koju biste vrijednosti za $\lambda$ izabrali na temelju ovih grafikona i zašto?
End of explanation
"""


def nonzeroes(coef, tol=1e-6): 
    return len(coef) - len(coef[sp.isclose(0, coef, atol=tol)])

"""
Explanation: 6. L1-regularizacija i L2-regularizacija
Svrha regularizacije jest potiskivanje težina modela $\mathbf{w}$ prema nuli, kako bi model bio što jednostavniji. Složenost modela može se okarakterizirati normom pripadnog vektora težina $\mathbf{w}$, i to tipično L2-normom ili L1-normom. Za jednom trenirani model možemo izračunati i broj ne-nul značajki, ili L0-normu, pomoću sljedeće funkcije:
End of explanation
"""


from sklearn.linear_model import Ridge
from sklearn.linear_model import Lasso

lambda6a = range(1,100)
d6a = 20
X6a_train, X6a_test, y6a_train, y6a_test = train_test_split(x6a, y6a, test_size = 0.5)
Phi6a_train = PolynomialFeatures(d6a).fit_transform(X6a_train.reshape(-1,1))

L0 = [];
L1 = [];
L2 = [];

L1_norm = lambda w: sum(abs(w));
L2_norm = lambda w: math.sqrt(np.dot(w.T, w));

for i in lambda6a:
    w6a = np.dot(np.dot(np.linalg.inv(np.dot(Phi6a_train.T, Phi6a_train) + np.dot(i, np.eye(d6a+1))), Phi6a_train.T), y6a_train);
    
    L0.append(nonzeroes(w6a))
    L1.append(L1_norm(w6a))
    L2.append(L2_norm(w6a))
    
plot(lambda6a, L0, lambda6a, L1, lambda6a, L2, linewidth = 1)
grid()

"""
Explanation: (a)
Za ovaj zadatak upotrijebite skup za učenje i skup za testiranje iz zadatka 3b. Trenirajte modele L2-regularizirane polinomijalne regresije stupnja $d=20$, mijenjajući hiperparametar $\lambda$ u rasponu ${1,2,\dots,100}$. Za svaki od treniranih modela izračunajte L{0,1,2}-norme vektora težina $\mathbf{w}$ te ih prikažite kao funkciju od $\lambda$.
Q: Objasnite oblik obiju krivulja. Hoće li krivulja za $\|\mathbf{w}\|_2$ doseći nulu? Zašto? Je li to problem? Zašto?
Q: Za $\lambda=100$, koliki je postotak težina modela jednak nuli, odnosno koliko je model rijedak?
End of explanation
"""



L0 = [];
L1 = [];
L2 = [];

for i in lambda6a:
    lass = Lasso(alpha = i, tol = 0.115).fit(Phi6a_train, y6a_train)
    w6b = lass.coef_
    
    L0.append(nonzeroes(w6b))
    L1.append(L1_norm(w6b))
    L2.append(L2_norm(w6b))
    
plot(lambda6a, L0, lambda6a, L1, lambda6a, L2, linewidth = 1)
legend(['L0', 'L1', 'L2'], loc = 'best')
grid()

"""
Explanation: (b)
Glavna prednost L1-regularizirane regresije (ili LASSO regression) nad L2-regulariziranom regresijom jest u tome što L1-regularizirana regresija rezultira rijetkim modelima (engl. sparse models), odnosno modelima kod kojih su mnoge težine pritegnute na nulu. Pokažite da je to doista tako, ponovivši gornji eksperiment s L1-regulariziranom regresijom, implementiranom u klasi  Lasso u modulu sklearn.linear_model.
End of explanation
"""


n_data_points = 500
np.random.seed(69)

# Generiraj podatke o bodovima na prijamnom ispitu koristeći normalnu razdiobu i ograniči ih na interval [1, 3000].
exam_score = np.random.normal(loc=1500.0, scale = 500.0, size = n_data_points) 
exam_score = np.round(exam_score)
exam_score[exam_score > 3000] = 3000
exam_score[exam_score < 0] = 0

# Generiraj podatke o ocjenama iz srednje škole koristeći normalnu razdiobu i ograniči ih na interval [1, 5].
grade_in_highschool = np.random.normal(loc=3, scale = 2.0, size = n_data_points)
grade_in_highschool[grade_in_highschool > 5] = 5
grade_in_highschool[grade_in_highschool < 1] = 1

# Matrica dizajna.
grades_X = np.array([exam_score,grade_in_highschool]).T

# Završno, generiraj izlazne vrijednosti.
rand_noise = np.random.normal(loc=0.0, scale = 0.5, size = n_data_points)
exam_influence = 0.9
grades_y = ((exam_score / 3000.0) * (exam_influence) + (grade_in_highschool / 5.0) \
            * (1.0 - exam_influence)) * 5.0 + rand_noise
grades_y[grades_y < 1] = 1
grades_y[grades_y > 5] = 5

"""
Explanation: 7. Značajke različitih skala
Često se u praksi možemo susreti sa podatcima u kojima sve značajke nisu jednakih magnituda. Primjer jednog takvog skupa je regresijski skup podataka grades u kojem se predviđa prosjek ocjena studenta na studiju (1--5) na temelju dvije značajke: bodova na prijamnom ispitu (1--3000) i prosjeka ocjena u srednjoj školi. Prosjek ocjena na studiju izračunat je kao težinska suma ove dvije značajke uz dodani šum.
Koristite sljedeći kôd kako biste generirali ovaj skup podataka.
End of explanation
"""


plt.figure()
plot(exam_score, grades_y, 'r+')
grid()

plt.figure()
plot(grade_in_highschool, grades_y, 'g+')
grid()

"""
Explanation: a)
Iscrtajte ovisnost ciljne vrijednosti (y-os) o prvoj i o drugoj značajki (x-os). Iscrtajte dva odvojena grafa.
End of explanation
"""



    w = Ridge(alpha = 0.01).fit(grades_X, grades_y).coef_
    print(w)

"""
Explanation: b)
Naučite model L2-regularizirane regresije ($\lambda = 0.01$), na podacima grades_X i grades_y:
End of explanation
"""


from sklearn.preprocessing import StandardScaler

#grades_y.reshape(-1, 1)

scaler = StandardScaler()
scaler.fit(grades_X)
grades_X_fixed = scaler.transform(grades_X)

scaler2 = StandardScaler()
scaler2.fit(grades_y.reshape(-1, 1))
grades_y_fixed = scaler2.transform(grades_y.reshape(-1, 1))

"""
Explanation: Sada ponovite gornji eksperiment, ali prvo skalirajte podatke grades_X i grades_y i spremite ih u varijable grades_X_fixed i grades_y_fixed. Za tu svrhu, koristite StandardScaler.
End of explanation
"""


grades_X_fixed_colinear = [[g[0],g[1],g[1]] for g in grades_X_fixed]



"""
Explanation: Q: Gledajući grafikone iz podzadatka (a), koja značajka bi trebala imati veću magnitudu, odnosno važnost pri predikciji prosjeka na studiju? Odgovaraju li težine Vašoj intuiciji? Objasnite.  
8. Multikolinearnost i kondicija matrice
a)
Izradite skup podataka grades_X_fixed_colinear tako što ćete u skupu grades_X_fixed iz
zadatka 7b duplicirati zadnji stupac (ocjenu iz srednje škole). Time smo efektivno uveli savršenu multikolinearnost.
End of explanation
"""



w = Ridge(alpha = 0.01).fit(grades_X_fixed_colinear, grades_y_fixed).coef_
print(w)

"""
Explanation: Ponovno, naučite na ovom skupu L2-regularizirani model regresije ($\lambda = 0.01$).
End of explanation
"""


w_001s = []
w_1000s = []

for i in range(10):
    X_001, X_1000, y_001, y_1000 = train_test_split(grades_X_fixed_colinear, grades_y_fixed, test_size = 0.5)
    w_001 = Ridge(alpha = 0.01).fit(X_001, y_001).coef_
    w_1000 = Ridge(alpha = 0.01).fit(X_1000, y_1000).coef_
    
    w_001s.append(w_001[0])
    w_1000s.append(w_1000[0])
    #print(w_001)
    #print(w_1000)
    #print()
    
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
scaler.fit(w_001s)
print(scaler.mean_)

scaler = StandardScaler()
scaler.fit(w_1000s)
print(scaler.mean_)

"""
Explanation: Q: Usporedite iznose težina s onima koje ste dobili u zadatku 7b. Što se dogodilo?
b)
Slučajno uzorkujte 50% elemenata iz skupa grades_X_fixed_colinear i naučite dva modela L2-regularizirane regresije, jedan s $\lambda=0.01$, a jedan s $\lambda=1000$. Ponovite ovaj pokus 10 puta (svaki put s drugim podskupom od 50% elemenata).  Za svaki model, ispišite dobiveni vektor težina u svih 10 ponavljanja te ispišite standardnu devijaciju vrijednosti svake od težina (ukupno šest standardnih devijacija, svaka dobivena nad 10 vrijednosti).
End of explanation
"""


lam = 0.01

phi = grades_X_fixed_colinear
s = np.add(np.dot(np.transpose(phi), phi), lam * np.identity(len(a)))
print(np.linalg.cond(s))

lam = 10

phi = grades_X_fixed_colinear
s = np.add(np.dot(np.transpose(phi), phi), lam * np.identity(len(a)))
print(np.linalg.cond(s))

"""
Explanation: Q: Kako regularizacija utječe na stabilnost težina?
Q: Jesu li koeficijenti jednakih magnituda kao u prethodnom pokusu? Objasnite zašto.
c)
Koristeći numpy.linalg.cond izračunajte kondicijski broj matrice $\mathbf{\Phi}^\intercal\mathbf{\Phi}+\lambda\mathbf{I}$, gdje je $\mathbf{\Phi}$ matrica dizajna (grades_fixed_X_colinear). Ponovite i za $\lambda=0.01$ i za $\lambda=10$.
End of explanation
"""


