# Studijny program: "Pravdepodobnost a matematicka statistika"
# Vyucujuci: Radoslav Harman, KAMS FMFI UK Bratislava

# Klasicke metricke skalovanie (multidimensional scaling; MDS)

### Vzdialenosti 39 Slovenskych miest
### Data su zo stranky http://www.slovakiasite.com/sk/vzdialenosti.php

mesta <- read.table("http://www.iam.fmph.uniba.sk/ospm/Harman/data/mesta_SK.txt", header = TRUE) 
print(mesta)
n <- 39

# Vypocet MDS 
help(cmdscale) # Matematicky zaujimavy parameter: add!

mesta.mds <- cmdscale(mesta, k = n - 1, eig = TRUE)
eig <- mesta.mds$eig
plot(eig); grid()

# Kontrola suhlasu s teoriou
A <- -0.5*as.matrix(mesta*mesta)
P <- diag(n) - (rep(1, n) %*% t(rep(1, n)))/n
B <- P %*% A %*% P
eigen(B)$values

# Miera fitu v 2D
sum(abs(eig[1:2]))/sum(abs(eig))

# Zobrazenie v 2D
x <- mesta.mds$points[, 1]
y <- mesta.mds$points[, 2]
plot(x, y, xlim = range(x)*1.2, ylim = range(y)*1.2, cex = 8, asp = 1)
text(x,y + 0.4, labels = colnames(mesta), col = "red", cex = 0.6)


### Pribuznost vybranych europskych krajin podla priebehu epidemie koronavirusu
### Data su prevzate zo stranky https://ourworldindata.org/coronavirus

# Nacitanie a zakladny prehlad o datach
cvd <- read.csv("c:/covid.csv")
dim(cvd)

# Pozrime si nahodny riadok (kvoli citatelnosti transponovany)
t(cvd[sample(1:nrow(cvd), 1), ])

# Preco si treba davat velky pozor:
summary(cvd$new_deaths)
cvd[!is.na(cvd$new_deaths) & (cvd$new_deaths < 0), c(3, 4, 9)]
cvd[!is.na(cvd$new_deaths) & (cvd$new_deaths > 5000), c(3, 4, 9)]

# Vyberme len niektore europske krajiny
vyber.krajin <- c("Austria", "Belarus", "Belgium", "Bulgaria", "Croatia", "Czechia",
           "Denmark", "Estonia", "Finland", "France", "Germany", "Greece",
           "Hungary", "Ireland", "Italy", "Latvia", "Lithuania", "Netherlands",       
           "Norway", "Poland", "Portugal", "Romania", "Russia", "Serbia", "Slovakia",
           "Slovenia", "Spain", "Sweden", "Switzerland", "Ukraine", "United Kingdom") 
n.krajin <- length(vyber.krajin)
cvd <- cvd[cvd$location %in% vyber.krajin, ]; dim(cvd)
cvd$location <- droplevels(cvd$location)
as.matrix(table(cvd$location))

as.matrix(table(cvd$date))
# Vyberme len datumy, v ktorych ma kazda zo zvolenych krajin zaznam
cvd <- cvd[as.character(cvd$date) > "2020-03-07", ]; dim(cvd)
cvd$date <- droplevels(cvd$date)
as.matrix(table(cvd$date))
vyber.dni <- levels(cvd$date)
n.dni <- length(vyber.dni)

cvd.mat <- matrix(0, nrow = length(vyber.krajin), ncol = length(vyber.dni))
cvd.clean <- as.data.frame(cvd.mat)
names(cvd.clean) <- vyber.dni
row.names(cvd.clean) <- vyber.krajin

# Vytvorime maticu, v ktorej su objekty vybrane krajiny a premenne new_deaths_smoothed pre dni
# Toto je naprogramovane velmi neefektivne. Du: najst efektivnejsi sposob :)
for (i.krajina in 1:n.krajin) {
  for (i.den in 1:n.dni) {
    riadok <- cvd[(cvd$location == vyber.krajin[i.krajina]) &
                  (cvd$date == vyber.dni[i.den]), ]
    if (is.na(riadok$new_deaths_smoothed_per_million)) riadok$new_deaths_smoothed_per_million <- 0
    cvd.clean[i.krajina, i.den] <- riadok$new_deaths_smoothed_per_million 
  }
}

paleta <- rainbow(n.krajin)
plot(c(0, 1), c(1, n.krajin), axes = FALSE, type = "n", xlab = "", ylab = "")
for (i.krajina in 1:n.krajin) 
  text(0.5, i.krajina, paste(i.krajina, vyber.krajin[i.krajina]),
       col = paleta[i.krajina])

plot(c(1, ncol(cvd.clean)), c(min(cvd.clean), max(cvd.clean)), type = "n"); grid()
for (i.krajina in seq_along(vyber.krajin)) {
  lines(t(cvd.clean[i.krajina,]), col = paleta[i.krajina])
}
# Ta krajina so zapornym peakom je Spanielsko

# Zobrazenie krajin podla "pribuznosti vyvoja v case"
# zalozime na matici nepodobnosti vyvoja poctov umrti
# vypocitanej cez robustny korelacny koeficient "Kendallovo tau"

# Vypocet matice korelacii
C <- matrix(1, ncol = n.krajin, nrow = n.krajin)
row.names(C) <- colnames(C) <- vyber.krajin
for (i1 in 1:(n.krajin - 1)) {
  for (i2 in (i1 + 1):n.krajin) {
    C[i2, i1] <- C[i1, i2] <- cor(as.numeric(cvd.clean[i1, ]),
                                  as.numeric(cvd.clean[i2, ]), method = "kendall")
  }
}

# Pozrime sa na korelacnu maticu C
library(corrplot)
corrplot(C)

# Transformacia na maticu nepodobnosti
D <- -log(C) # Da sa aj jednoducho D <- 1 - C

# Vypocitame vizualizaciu pomocou mds
cvd.mds <- cmdscale(D, k = nrow(D) - 1, eig = TRUE)
eig <- cvd.mds$eig; plot(eig)
sum(abs(eig[1:2]))/sum(abs(eig))

x <- cvd.mds$points[, 1]
y <- cvd.mds$points[, 2]

# Zobrazenie v 2D
plot(x, y, xlim = range(x)*1.2, ylim = range(y)*1.2, type = "n", asp = 1)
text(x, y, labels = vyber.krajin, cex = 1, col = paleta)

# Detailny pohlad na stred obrazku
plot(x, y, xlim = range(x)*0.4, ylim = range(y)*0.4, type = "n", asp = 1)
text(x, y, labels = vyber.krajin, cex = 1, col = paleta)