---
title: "Classification report"
output: html_document
---

```{r, echo=F, message=F}
library(dplyr)
library(reshape2)
library(knitr)
library(ggplot2)
library(pheatmap)
library(RColorBrewer)
library(png)

PATH = "/home/hpc_tanel/deep-yeast/make-data/data/HOwt/"
VAL_PRED_PATH = "../random-forest/rf_preds_val.csv"
TEST_PRED_PATH = "../random-forest/rf_preds_test.csv"
MODEL_NAME = "Random forest"
```

### Model: `r I(MODEL_NAME)`

```{r, echo=F}
doc = read.table(paste0(PATH, "HOwt_doc.txt"), sep=";", header=T, as.is=T)

val = read.table(paste0(PATH, "HOwt_val.txt"), as.is=T)
names(val) = c("src", "y")
val = left_join(val, doc, by=c("y"="label_idx"))
val$pred = read.table(VAL_PRED_PATH, as.is=T)[, 1]
val$pred_label = doc$label[val$pred + 1]

test = read.table(paste0(PATH, "HOwt_test.txt"), as.is=T)
names(test) = c("src", "y")
test = left_join(test, doc, by=c("y"="label_idx"))
test$pred = read.table(TEST_PRED_PATH, as.is=T)[, 1]
test$pred_label = doc$label[test$pred + 1]

```


### Accuracy
```{r, echo=F}
get_accuracy = function(data, n=2){
  round(100*mean(data$pred == data$y), n)
}

```

`r MODEL_NAME` accuracy on validation data was `r get_accuracy(val, 1)/100`, and on test data `r get_accuracy(test, 1)/100`.

### Accuracy per category

```{r, echo=F, fig.width=10, fig.height=7, message=F}
recall_tbl = test %>%
  group_by(label) %>%
  summarise(num_imgs=n(),
            num_correct=sum(y==pred),
            recall=mean(y==pred))

ggplot(recall_tbl, aes(x=factor(label), y=recall)) + 
  geom_histogram(stat="identity", fill="#53cfff") +
  geom_text(aes(label=round(recall, 2)), vjust=1.2, color="white") + 
  theme_light(base_size=20) +
  xlab("") +
  ylab("recall") + 
  ggtitle("Recall per group - test data\n") +
  theme(plot.title=element_text(size=18),
        axis.text.x = element_text(angle = 90, hjust = 1)) +
  ylim(0,1)

precision_tbl = test %>%
  group_by(pred_label) %>%
  summarise(precision=mean(y==pred))

ggplot(precision_tbl, aes(x=factor(pred_label), y=precision)) + 
  geom_histogram(stat="identity", fill="#53cfff") +
  geom_text(aes(label=round(precision, 2)), vjust=1.2, color="white") + 
  theme_light(base_size=20) +
  xlab("") +
  ylab("precision") + 
  ggtitle("Precision per group - test data\n") +
  theme(plot.title=element_text(size=18),
        axis.text.x = element_text(angle = 90, hjust = 1)) + 
  ylim(0,1)

metrics_tbl = left_join(precision_tbl, recall_tbl, by=c("pred_label"="label"))
metrics_tbl.m = melt(metrics_tbl[, c("pred_label", "precision", "recall")])

ggplot(metrics_tbl.m, aes(x=factor(pred_label), y=value, fill=variable)) + 
  geom_bar(stat="identity", position="dodge") +
  theme_minimal(base_size=20) +
  xlab("") +
  ylab("metric") + 
  ggtitle("Metrics for test data\n") +
  theme(plot.title=element_text(size=18),
        axis.text.x = element_text(angle = 90, hjust = 1)) +
  scale_fill_brewer(palette="Set1") +
  ylim(0,1)

get_per_category_accuracy = function(x){
  preds = test$pred_label
  preds = ifelse(preds==x, 1, 0)
  trues = test$label
  trues = ifelse(trues==x, 1, 0)
  mean(preds==trues)
}

accuracy_per_group = sapply(metrics_tbl$pred_label, get_per_category_accuracy)
metrics_tbl$accuracy = accuracy_per_group

ggplot(metrics_tbl, aes(x=factor(pred_label), y=accuracy)) + 
  geom_histogram(stat="identity", fill="#53cfff") +
  geom_text(aes(label=round(accuracy, 2)), vjust=1.2, color="white") + 
  theme_light(base_size=20) +
  xlab("") +
  ylab("accuracy") + 
  ggtitle("Accuracy per group - test data\n") +
  theme(plot.title=element_text(size=18),
        axis.text.x = element_text(angle = 90, hjust = 1)) + 
  ylim(0,1)

metrics_tbl_print = metrics_tbl
metrics_tbl_print$precision = round(metrics_tbl_print$precision, 3)
metrics_tbl_print$recall = round(metrics_tbl_print$recall, 3)
metrics_tbl_print$accuracy = round(metrics_tbl_print$accuracy, 3)
colnames(metrics_tbl_print) = c("Class", "Precision", "Number of cell images", "Correct calls", "Recall", "Accuracy")
kable(metrics_tbl_print)

#print(paste0("Mean precision over classes: ", round(mean(metrics_tbl$precision), 2)))
#print(paste0("Median precision over classes: ", round(median(metrics_tbl$precision), 2)))
#print(paste0("Mean recall over classes: ", round(mean(metrics_tbl$recall), 2)))
#print(paste0("Median recall over classes: ", round(median(metrics_tbl$recall), 2)))
#print(paste0("Mean accuracy over classes: ", round(mean(metrics_tbl$accuracy), 2)))
#print(paste0("Median accuracy over classes: ", round(median(metrics_tbl$accuracy), 2)))

write.table(metrics_tbl, file="RF_metrics.csv", quote=FALSE, sep=";", row.names=F)
```

Across classes, the mean and median precision were `r round(mean(metrics_tbl$precision), 2)` and `r round(median(metrics_tbl$precision), 2)` respectively, recall `r round(mean(metrics_tbl$recall), 2)` and `r round(median(metrics_tbl$recall), 2)`, and accuracy `r round(mean(metrics_tbl$accuracy), 2)` and `r round(median(metrics_tbl$accuracy), 2)`.

### Confusion matrix

Each element $M_{ij}$ in the confusion matrix $M$ shows the number of observations from class $i$ which were predicted as class $j$. The confusion matrix in the figure below has been normalized by class support size (each row has been divided by the number of elements in that class).

```{r, echo=F, fig.width=8}
#write.table(test, file="RF_test_matrix.csv", quote=FALSE, sep=";", row.names=F)
cm = table(test$y, test$pred)
colnames(cm) = doc$label
rownames(cm) = doc$label

tmp_mat = prop.table(as.matrix(cm), 1)
pheatmap(tmp_mat, cluster_rows = F, cluster_cols = F,
         color = colorRampPalette((brewer.pal(n = 7, name = "Blues")))(100))

kable(cm)
```

To see misclassification patterns more clearly, the diagonal of the confusion matrix has been excluded in the following figure.

```{r, echo=F, fig.width=8}
diag(tmp_mat) = NA
pheatmap(tmp_mat, cluster_rows = F, cluster_cols = F,
         color = colorRampPalette((brewer.pal(n = 7, name = "Blues")))(100))

```


### Protein classification accuracy

```{r, echo=F}
extract_code = function(x){
  # input: "plate01/001002000-0-5-1.jpg"
  # output: "plate01/001002000"
  return(strsplit(x, "-")[[1]][1])
}
test$code = sapply(test$src, extract_code)

most_frequent = function(x){
  names(sort(table(x), decreasing = T))[1]
}

df = test %>%
  group_by(code) %>%
  summarise(label=most_frequent(label),
            class=most_frequent(y),
            pred=most_frequent(pred))

```

Across proteins, the mean and median accuracy were `r sprintf("%.3f", mean(df$class == df$pred))` and `r sprintf("%.3f", median(df$class == df$pred))`, respectively.

```{r, echo=FALSE}
df2 = df %>%
  group_by(label) %>%
  summarise(num_imgs=n(),
            num_correct=sum(class==pred),
            accuracy=num_correct/num_imgs)
colnames(df2) = c("Class", "Number of proteins", "Correct calls", "Accuracy")
df2$Accuracy = round(df2$Accuracy, 3)
kable(df2)
```


### Protein classification accuracy (at least 10 cell images)

```{r, echo=F}
df = test %>%
  group_by(code) %>%
  summarise(label=most_frequent(label),
            class=most_frequent(y),
            pred=most_frequent(pred),
            num=n()) %>%
  filter(num >= 10) 

```

Across proteins with at least 10 cell images, the mean and median accuracy were `r sprintf("%.3f", mean(df$class == df$pred))` and `r sprintf("%.3f", median(df$class == df$pred))`, respectively.

```{r, echo=F}
df2 = df %>%
  group_by(label) %>%
  summarise(num_imgs=n(),
            num_correct=sum(class==pred),
            accuracy=num_correct/num_imgs)
colnames(df2) = c("Class", "Number of proteins", "Correct calls", "Accuracy")
df2$Accuracy = round(df2$Accuracy, 3)
kable(df2)
```


```{r, eval=FALSE, echo=F}
df = test %>%
  group_by(code) %>%
  summarise(label=most_frequent(label),
            class=most_frequent(y),
            pred=most_frequent(pred),
            num=n()) %>%
  group_by(num) %>%
  summarise(num_imgs=n(),
            num_correct=sum(class==pred),
            accuracy=num_correct/num_imgs)
write.table(df, file="RF_num_imgs_metrics.csv", quote=FALSE, sep=";", row.names=F)
```

### Misclassifications - false negatives

Original and scaled cell images, respectively.

```{r, echo=F, fig.width=10, fig.height=2}
set.seed(123)
DATA_PATH = "/storage/hpc_tanel/images-png-merged-patches/HOwt/"
DATA_PATH_SCALED = "/storage/hpc_tanel/images-png-merged-patches-scaled/HOwt/"

for(i in 1:12){
  class = doc$label[i]
  sub = test[(test$label == class) & (test$y != test$pred), ]
  sample_src = sample(sub$src, min(20, nrow(sub)))
  
  fnr = 1 - metrics_tbl$recall[metrics_tbl$pred_label == class]
  fnr = round(fnr, 3)
  print(paste0(class, ", false negative rate: ", fnr))
  par(mfrow=c(2, 10), mar=c(0.1, 0.1, 0.1, 0.1)/10)
  for(j in 1:length(sample_src)){
    path = paste0(DATA_PATH, sample_src[j])
    img = readPNG(path, native=TRUE)
    plot.new()
    rasterImage(img,-0.03,-0.03,1.03,1.03,interpolate=FALSE)
  }
  par(mfrow=c(2, 10), mar=c(0.1, 0.1, 0.1, 0.1)/10)
  for(j in 1:length(sample_src)){
    path = paste0(DATA_PATH_SCALED, sample_src[j])
    img = readPNG(path, native=TRUE)
    plot.new()
    rasterImage(img,-0.03,-0.03,1.03,1.03,interpolate=FALSE)
  }
}
```

### Misclassifications - false positives

```{r, echo=F, fig.width=10, fig.height=2}
DATA_PATH = "/storage/hpc_tanel/images-png-merged-patches/HOwt/"
DATA_PATH_SCALED = "/storage/hpc_tanel/images-png-merged-patches-scaled/HOwt/"

for(i in 1:12){
  class = doc$label_idx[i]
  sub = test[(test$pred == class) & (test$y != test$pred), ]
  sample_src = sample(sub$src, min(20, nrow(sub)))
  
  print(doc$label[i])
  par(mfrow=c(2, 10), mar=c(0.1, 0.1, 0.1, 0.1)/10)
  for(j in 1:length(sample_src)){
    path = paste0(DATA_PATH, sample_src[j])
    img = readPNG(path, native=TRUE)
    plot.new()
    rasterImage(img,-0.03,-0.03,1.03,1.03,interpolate=FALSE)
  }
  par(mfrow=c(2, 10), mar=c(0.1, 0.1, 0.1, 0.1)/10)
  for(j in 1:length(sample_src)){
    path = paste0(DATA_PATH_SCALED, sample_src[j])
    img = readPNG(path, native=TRUE)
    plot.new()
    rasterImage(img,-0.03,-0.03,1.03,1.03,interpolate=FALSE)
  }
}
```

### Common mistakes (top 5)

```{r, echo=F, fig.width=10, fig.height=2}
DATA_PATH = "/storage/hpc_tanel/images-png-merged-patches/HOwt/"
DATA_PATH_SCALED = "/storage/hpc_tanel/images-png-merged-patches-scaled/HOwt/"

cm = table(test$label, test$pred_label)
tmp_mat = prop.table(as.matrix(cm), 1)
diag(tmp_mat) = 0

max_vals = sort(tmp_mat, decreasing = T)[1:5]

for(i in 1:length(max_vals)){
  val = max_vals[i]
  idx = which(tmp_mat==val, arr.ind=T)
  idx_row = idx[1]
  idx_col = idx[2]
  
  true_class = rownames(tmp_mat)[idx_row]
  pred_class = colnames(tmp_mat)[idx_col]
  num_imgs = metrics_tbl$num_imgs[metrics_tbl$pred_label == true_class]
  count = cm[idx_row, idx_col]
  print(paste0("True class: ", true_class, 
               " (total images ", num_imgs,
               "), predicted as: ", pred_class,
               ", count: ", count,
               " (", sprintf("%.1f", count/num_imgs*100), "% error)"))
  
  sub = test[(test$label == true_class) & (test$pred_label == pred_class), ]
  sample_src = sample(sub$src, min(20, nrow(sub)))
  
  par(mfrow=c(2, 10), mar=c(0.1, 0.1, 0.1, 0.1)/10)
  for(j in 1:length(sample_src)){
    path = paste0(DATA_PATH, sample_src[j])
    img = readPNG(path, native=TRUE)
    plot.new()
    rasterImage(img,-0.03,-0.03,1.03,1.03,interpolate=FALSE)
  }
  par(mfrow=c(2, 10), mar=c(0.1, 0.1, 0.1, 0.1)/10)
  for(j in 1:length(sample_src)){
    path = paste0(DATA_PATH_SCALED, sample_src[j])
    img = readPNG(path, native=TRUE)
    plot.new()
    rasterImage(img,-0.03,-0.03,1.03,1.03,interpolate=FALSE)
  }
}
```
