# Univerzitet u Novom Sadu, Fakultet tehničkih nauka, Novi Sad, Srbija
# Studijski program OAS Informacioni inženjering
# Predmet Metode i tehnike nauke o podacima

# Pomoćni sadržaj

# Korišćeni podaci
#
# Skup podataka "Annual Social and Economic Supplements"
# (segment za 2023. godinu)
# - Internet sajt Popisnog biroa SAD
# - podaci iz Godišnjeg društvenog i ekonomskog dodatka za 2023. godinu 
#   iz Tekućeg istraživanja stanovništva (SAD)
#   - pružalac podataka: Popisni biro SAD
#   - vremenska odrednica za segment podataka: 2023.
#   - datum najnovijeg povezanog ažuriranja: 4. 9. 2024.
#   - format podataka: CSV i drugi
#   - podskup: podaci o domaćinstvima (datoteka "hhpub23.csv" među datotekama
#     unutar arhivske datoteke "asecpub23csv.zip")
#
# Internet strana:
# https://www.census.gov/data/datasets/time-series/demo/cps/cps-asec.html


# %% Biblioteke

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from keras.utils import to_categorical
from keras.models import Sequential
from keras.layers import Dense
from keras.optimizers import SGD, Adam
from sklearn.metrics import confusion_matrix, ConfusionMatrixDisplay


# %% Generator vrednosti

gv = np.random.default_rng(83)


# %% Učitavanje i odabir podataka

domaćinstva = pd.read_csv("hhpub23.csv")

print("Izvorne kolone:")
print(domaćinstva.columns.values)

domaćinstva = domaćinstva.loc[:, ["GEDIV", "GTCBSASZ", "H_LIVQRT", "HH5TO18", 
                                  "HEARNVAL", "HINC_FR", "HINC_SE", "HINC_WS", 
                                  "HANN_YN", "HCSP_YN", "HDIV_YN", "HFIN_YN", 
                                  "HINC_UC", "HFOODSP", "NOW_HPUB", 
                                  "NOW_HPRIV"]]

domaćinstva.rename(columns={
    "GEDIV": "područje", 
    "GTCBSASZ": "veličina_metropolitanskog_područja", 
    "H_LIVQRT": "vrsta_stambenog_prostora", 
    "HH5TO18": "broj_osoba_od_5_do_18_godina", 
    "HEARNVAL": "ukupna_zarada", 
    "HINC_FR": "samozaposlenje_poljoprivreda", 
    "HINC_SE": "samozaposlenje_poslovanje", 
    "HINC_WS": "plata", 
    "HANN_YN": "anuitet", 
    "HCSP_YN": "alimentacija", 
    "HDIV_YN": "akcije", 
    "HFIN_YN": "finansijska_pomoć", 
    "HINC_UC": "naknada_za_nezaposlene", 
    "HFOODSP": "bonovi_za_hranu", 
    "NOW_HPUB": "javno_osiguranje", 
    "NOW_HPRIV": "privatno_osiguranje"
    }, inplace=True)

domaćinstva = domaćinstva.loc[domaćinstva["ukupna_zarada"] > 0, :]

print("Odabrani podaci:")
print(domaćinstva)

print("Odabrani kolone:")
print(domaćinstva.dtypes)


# %% Analiza podataka

ukupna_zarada_min = domaćinstva["ukupna_zarada"].min()
ukupna_zarada_medijana = domaćinstva["ukupna_zarada"].median()
ukupna_zarada_prosek = domaćinstva["ukupna_zarada"].mean()
ukupna_zarada_maks = domaćinstva["ukupna_zarada"].max()
print("Medijalna ukupna zarada: {:.2f}".format(ukupna_zarada_medijana))
print("Prosečna ukupna zarada: {:.2f}".format(ukupna_zarada_prosek))

granice = [33, 67]
percentili = np.percentile(domaćinstva["ukupna_zarada"], granice)
print("{}. percentil za ukupnu zaradu: {:.2f}".format(granice[0], 
                                                      percentili[0]))
print("{}. percentil za ukupnu zaradu: {:.2f}".format(granice[1], 
                                                      percentili[1]))

sekvenca = np.arange(ukupna_zarada_min, ukupna_zarada_maks + 10000, 10000)

plt.figure(figsize=(15, 6))
plt.hist(domaćinstva["ukupna_zarada"], bins=sekvenca, color="orange")

plt.axvline(percentili[0], color="blue", linestyle="--", 
            linewidth=1, label="{}. percentil".format(granice[0]))
plt.axvline(percentili[1], color="red", linestyle="--", 
            linewidth=1, label="{}. percentil".format(granice[1]))

plt.title("Ukupna zarada domaćinstva")
plt.xlabel("Ukupna zarada")
plt.ylabel("Učestalost")
plt.legend()


# %% Pripremanje podataka za postupke obučavanja i validacije neuronskih mreža

domaćinstva["ukupna_zarada_klasa"] = domaćinstva["ukupna_zarada"].apply(
    lambda x: 0 if x <= percentili[0] else 1 if x <= percentili[1] else 2)

prediktori = ["područje", "veličina_metropolitanskog_područja", 
              "vrsta_stambenog_prostora", "broj_osoba_od_5_do_18_godina", 
              "samozaposlenje_poljoprivreda", "samozaposlenje_poslovanje",
              "plata", "anuitet", "alimentacija", "akcije", 
              "finansijska_pomoć", "naknada_za_nezaposlene", 
              "bonovi_za_hranu", "javno_osiguranje", "privatno_osiguranje"]

cilj = "ukupna_zarada_klasa"
cilj_brojnost = domaćinstva[cilj].nunique()

podaci_obučavanje = domaćinstva.sample(frac=0.8, random_state=gv.bit_generator)
podaci_obučavanje_ulaz = podaci_obučavanje.get(prediktori)
podaci_obučavanje_izlaz = to_categorical(podaci_obučavanje.get(cilj), 
                                         cilj_brojnost)

podaci_validacija = domaćinstva.drop(index=podaci_obučavanje.index)
podaci_validacija_ulaz = podaci_validacija.get(prediktori)
podaci_validacija_izlaz = to_categorical(podaci_validacija.get(cilj),
                                         cilj_brojnost)

stvarne_klase = podaci_validacija.get(cilj)


# %% Obučavanje i validacija neuronskih mreža

broj_epoha = 20

print("Model A")
nm_a = Sequential()
nm_a.add(Dense(2, activation="tanh"))
nm_a.add(Dense(cilj_brojnost, activation="softmax"))
nm_a.compile(loss="categorical_crossentropy", 
             optimizer=SGD(), 
             metrics=["categorical_accuracy"])
rezultati_a = nm_a.fit(podaci_obučavanje_ulaz, podaci_obučavanje_izlaz, 
                       batch_size=256, epochs=broj_epoha, verbose=1, 
                       validation_data=(podaci_validacija_ulaz, 
                                        podaci_validacija_izlaz))
procene_a = nm_a.predict(podaci_validacija_ulaz)
procene_klase_a = np.argmax(procene_a, axis=1) 

print("Model B")
nm_b = Sequential()
nm_b.add(Dense(24, activation="sigmoid"))
nm_b.add(Dense(24, activation="sigmoid"))
nm_b.add(Dense(cilj_brojnost, activation="softmax"))
nm_b.compile(loss="categorical_crossentropy", 
             optimizer=Adam(), 
             metrics=["categorical_accuracy"])
rezultati_b = nm_b.fit(podaci_obučavanje_ulaz, podaci_obučavanje_izlaz, 
                       batch_size=64, epochs=broj_epoha, verbose=1, 
                       validation_data=(podaci_validacija_ulaz, 
                                        podaci_validacija_izlaz))
procene_b = nm_b.predict(podaci_validacija_ulaz)
procene_klase_b = np.argmax(procene_b, axis=1) 


# %% Ispitivanje performansi neuronskih mreža

matrica_konfuzije_a = confusion_matrix(stvarne_klase, procene_klase_a)
print("Matrica konfuzije za model A")
print(matrica_konfuzije_a)

matrica_konfuzije_b = confusion_matrix(stvarne_klase, procene_klase_b)
print("Matrica konfuzije za model B")
print(matrica_konfuzije_b)

prikaz_a = ConfusionMatrixDisplay(matrica_konfuzije_a)
prikaz_b = ConfusionMatrixDisplay(matrica_konfuzije_b)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 6))
fig.suptitle("Matrice konfuzije")
prikaz_a.plot(cmap=plt.cm.Purples, ax=ax1)
ax1.set(xlabel="Procenjena klasa", ylabel="Stvarna klasa")
ax1.set_title("Model A")
prikaz_b.plot(cmap=plt.cm.Greens, ax=ax2)
ax2.set(xlabel="Procenjena klasa", ylabel="Stvarna klasa")
ax2.set_title("Model B")

plt.figure(figsize=(8, 8))
plt.plot(range(1, broj_epoha + 1), 
         rezultati_a.history["categorical_accuracy"], 
         ":", color="purple", label="model A: obučavanje")
plt.plot(range(1, broj_epoha + 1), 
         rezultati_a.history["val_categorical_accuracy"], 
         "--", color="purple", label="model A: validacija")
plt.plot(range(1, broj_epoha + 1), 
         rezultati_b.history["categorical_accuracy"], 
         ":", color="green", label="model B: obučavanje")
plt.plot(range(1, broj_epoha + 1), 
         rezultati_b.history["val_categorical_accuracy"], 
         "--", color="green", label="model B: validacija")
plt.xlim(0)
plt.title("Performanse neuronskih mreža")
plt.xticks([i for i in range(1, broj_epoha + 1) if i % 5 == 0])
plt.xlabel("Epoha")
plt.ylabel("Kategorijska tačnost")
plt.legend()