---
title: "Lab9"
author: "Michał Ciach"
date: "November 25, 2017"
output: 
  html_document:
    toggle_code: True
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(ggplot2)
library(reshape2)
library(readr)
library(stringr)
```

## Hidden Markov Models 

## Text correction


The following function performs basic formatting of a text and encodes it in a factor of phonemes:

```{r}
phonemes <- str_split("a ą b c ć d e ę f g h i j k l ł m n ń o ó p r s ś t u ú v w x y z ź ż", pattern=" ")[[1]]
phonemes <- c(phonemes, c("rz", "ch", "ci", "zi", "dż", "dz", "dź", "sz", "cz", "si", " "))
factorize.text <- function(str, phnms=phonemes){
  # Converts a string to a factor of phonemes
  str <- str_replace_all(str, '[^\\s\\w]*', '')  # replace all non-word non-whitespace characters with an empty character
  str <- str_replace_all(str, '[0-9]*', '')  # replace all numbers with an empty character
  str <- str_replace_all(str, '[\\s]+', ' ')  # replace a sequence of at least one whitespace character with a single space
  str <- str_to_lower(str)
  L <- str_length(str)  # need to use str_length instead of length because from R's perspective str is a vector of length one
  fctr <- factor(character(L), levels=phnms)  # allocate memory
  str_i = 1  # position in the original string
  fctr_i = 1  # position in the factor
  print.thr <- round(L/10)  # used to print percentage of string processed
  while(str_i < L){
    if(str_i %% print.thr == 0){
      print(paste0(10*str_i/print.thr, "% processed"))
    }
    s <- str_sub(str, str_i, str_i+1)  # take a pair of letters
    if(s %in% phnms){ 
      # if the pair is in the phoneme list, add it to the factor
      fctr[fctr_i] <- s
      fctr_i <- fctr_i + 1  
      str_i <- str_i + 2
    } else {
      # otherwise, add the first letter only
      fctr[fctr_i] <- str_sub(s, 1, 1)
      fctr[fctr_i+1] <- str_sub(s, 2, 2)
      fctr_i <- fctr_i + 1
      str_i <- str_i + 1
    }
  }
  fctr <- fctr[1:fctr_i]  # the final factor will be shorter because we've merged some of the letters
  return(fctr)
}
```

```{r, include=F, eval=F}
# For tutorial preparation; would be better to introduce errors to text not factor
mutate <- function(fctr, p=0.01){
  # Simulates ortographical errors; p = probablity of error
  for(i in 1:length(fctr)){
    if(fctr[i] == "u"){
      fctr[i] <- sample(c("u", "ó"), 1, prob=c(1-p, p))
    }
    if(fctr[i] == "ó"){
      fctr[i] <- sample(c("ó", "u"), 1, prob=c(1-p, p))
    }
    if(fctr[i] == "ż"){
      fctr[i] <- sample(c("ż", "rz"), 1, prob=c(1-p, p))
    }
    if(fctr[i] == "rz"){
      fctr[i] <- sample(c("rz", "ż"), 1, prob=c(1-p, p))
    }
    if(fctr[i] == "ć"){
      fctr[i] <- sample(c("ć", "ci"), 1, prob=c(1-p, p))
    }
    if(fctr[i] == "ci"){
      fctr[i] <- sample(c("ci", "ć"), 1, prob=c(1-p, p))
    }
    if(fctr[i] == "ś"){
      fctr[i] <- sample(c("ś", "si"), 1, prob=c(1-p, p))
    }
    if(fctr[i] == "si"){
      fctr[i] <- sample(c("si", "ś"), 1, prob=c(1-p, p))
    }
    if(fctr[i] == "ź"){
      fctr[i] <- sample(c("ź", "zi"), 1, prob=c(1-p, p))
    }
    if(fctr[i] == "zi"){
      fctr[i] <- sample(c("zi", "ź"), 1, prob=c(1-p, p))
    }
    if(fctr[i] == "ch"){
      fctr[i] <- sample(c("ch", "h"), 1, prob=c(1-p, p))
    }
    if(fctr[i] == "h"){
      fctr[i] <- sample(c("h", "ch"), 1, prob=c(1-p, p))
    }
  }
  return(fctr)
}

#mutalice <- mutate(alice, p=0.1)
# # mutalice.str <- str_c(mutalice, collapse="")  # that's not a good way to do it because we'd need to linearize 
# # write_file(mutalice, "Datasets/Ortografia/Mutalicja.txt")
#mutalice2 <- as.character(mutalice)
#write(mutalice2, file="Datasets/Ortografia/Mutalicja.txt")
```


Load the original text:

```{r}
alice <- read_file("Datasets/Ortografia/Alicja.txt")
```

Factorize the original text (it may take a while - loops in R are very slow): 

```{r}
alice <- factorize.text(alice)
```

The text with errors is already factorized and written in a convenient format. Read it and split over newline character. This procedure will add a blank character at the end, because the file ends with a newline character. Remove the last element and convert the character vector to a factor. You need to specify the factor levels explicitly to make it consistent with the original text.

```{r}
mutalice <- read_file("Datasets/Ortografia/Mutalicja.txt")
mutalice <- str_split(mutalice, "\n")[[1]]
mutalice <- mutalice[-length(mutalice)]
mutalice <- factor(mutalice, levels=levels(alice))
```

Inspect the number of errors.

```{r}
print(paste("The text contains", round(100*mean(alice != mutalice), 2), "percent of ortographical errors"))
```

Compute the average occurences of each letter in the original text.

```{r}
letter.frequency <- tapply(alice, alice, length)
letter.frequency <- letter.frequency/sum(letter.frequency)
round(letter.frequency, 2)
```

We will use the texts for supervised training of the Hidden Markov Model. The transition matrix needs to be trained from the original text. Use both texts to train the observation matrix.

We will train the matrix using a ML estimate of transition probabilities. Despite an advanced mathematical derivation, it's the most natural way you could imagine: just count the numbers of transitions and normalize the rows to sum to 1. Note that if you use apply to normalize the rows, you need to transpose the resulting matrix - apply stacks the vectors column by column. A transposition is implemented in the function `t()`. The same needs to be done with the observation matrix. We'll also do a trick - we'll add one observation to each field of the observation matrix. Even though this adds errors, it actually makes the algorithm work better. Try to comment the trick, see what happens when we correct the text and think why. 

```{r}
M <- matrix(0, nrow=length(phonemes), ncol=length(phonemes), dimnames=list(phonemes, phonemes))
Y <- matrix(0, nrow=length(phonemes), ncol=length(phonemes), dimnames=list(phonemes, phonemes))
for(i in seq(1, length(alice), by=2)){
  M[alice[i], alice[i+1]] = M[alice[i], alice[i+1]] + 1
  Y[alice[i], mutalice[i]] = Y[alice[i], mutalice[i]] + 1
  Y[alice[i+1], mutalice[i+1]] = Y[alice[i+1], mutalice[i+1]] + 1
}
M <- t(apply(M, 1, function(x) x/sum(x)))
Y <- Y + 1  # the trick
Y <- t(apply(Y, 1, function(x) x/sum(x)))
```

We can now check the probabilities of some errors:

```{r}
round(Y['ć', ], 3)
```

We will use the trained matrices to correct the errors in the text. What we want to do is find the most probable vector of *states*, i.e. the correct letters, given the vector of *observations*. Since we have trained the matrices, we can use the Viterbi algorithm. We'll use the algorithm implemented in the library `HMM`. We will assume that the probability distribution of the first letter is given by the occurences in the original text.

```{r}
library(HMM)
alice.hmm <- initHMM(States=levels(mutalice), # set of states of Markov process
                     Symbols=levels(mutalice),  # set of possible observations, in our case it's the set of levels
                     startProbs=letter.frequency,
                     transProbs=M,
                     emissionProbs=Y)
```

We have initialized our Hidden Markov Model. Now, run the Viterbi algorithm. Correcting the whole text would take a lot of time; We'll just correct a fragment. 

```{r}
L = 4000
coralice <- viterbi(alice.hmm, mutalice[1:L])  # Corrected mutalice 
```

Let's look at a fragment the corrected text.

```{r}
print(c("Mutated text:", str_c(mutalice[1:100], collapse='')))
print(c("Corrected text:", str_c(coralice[1:100], collapse='')))
```

```{r}
print(c("Mutated text:", str_c(mutalice[3673:4000], collapse='')))
print(c("Corrected text:", str_c(coralice[3673:4000], collapse='')))
```

Our approach has corrected quite many errors. Note that it didn't correct the word "morze" to "może" - but the former is a fully correct word, and our method does not look at the context. Our approach has also introduced an error, by changing "ówdzie" to "uwdzie". That's reasonable, because "ó" is very rare at the beginning of a word. 

Notice also that our approach doesn't handle foreign words very well (which is also reasonable, because it was trained on a Polish text):

```{r}
str_c(mutalice[193:248], collapse='')
str_c(coralice[193:248], collapse='')
```

```{r}
str_c(mutalice[263:324], collapse='')
str_c(coralice[263:324], collapse='')
```

```{r}
print(paste("Before correction:", sum(alice[1:L] != mutalice[1:L]), "errors"))
print(paste("After correction:", sum(alice[1:L] != coralice[1:L]), "errors (including foreign words)"))
```

## Unsupervised learning

Our next task is separating stems from inflectional suffixes. For example, in the word "kraina", the stem is the "krain" and the inflectional suffix is "a". We want to make an algorithm which will automatically detect suffixes. But, we have a problem: we have no training set! That is why we need to use *unsupervised learning* - our HMM needs to figure the transition matrix and emmission probabilites all by itselt. 

The Baum-Welch algorithm, covered during the lecture, performs an unsupervised learning of a Hidden Markov Model. We'll use it to train our model. First, we need to state the initial transition matrix and emmission probabilities. We'll use some basic knowledge about stems and suffixes. First, if a letter belongs to a stem, then the next letter will rather still belong to a stem than to a suffix. Second, there are no whitespaces in the stem, but a suffix always contains one, so there is a pretty high probability of emmitting a whitespace at the stem. Also, the text will always start from a stem. 

In general, it's not a good idea to start from any matrices. In particular, avoid starting your algorithm on uniform probabilities, because usually that's a local minimum. 

```{r}
stem.M <- matrix(c(0.9, 0.1, 0.4, 0.6), byrow=T, ncol=2, dimnames=list(c("Stem", "Suffix"), c("Stem", "Suffix")))  # used as a starting point for BW algorithm
stem.M
```

```{r}
stem.Start <- c(1, 0)
stem.Y <- matrix(1, nrow=2, ncol=length(phonemes), dimnames=list(c("Stem", "Suffix"), phonemes)) 
stem.Y[1, " "] <- 0 # no whitespace at stem
stem.Y[2, " "] <- 3 # slightly higher probability of emmission of whitespace at the suffix
stem.Y[1, ] <- stem.Y[1,]/sum(stem.Y[1, ])
stem.Y[2, ] <- stem.Y[2, ]/sum(stem.Y[2, ])
```

We're all set with our starting point. Now, create the HMM. 

```{r}
stem.hmm <- initHMM(States=c("Stem", "Suffix"),
                    Symbols=phonemes,
                    startProbs = stem.Start,
                    transProbs=stem.M,
                    emissionProbs=stem.Y)
```

We have initialized our model with our basic knowledge on stems and suffixes. Now, to train it, use `baumWelsch` function from the `HMM` library. The first word in our text is "Lewis", which is not a Polish word, and may disturb our results. That's why we won't use it for training our algorithm. Be warned: BW algorithm may take a considerable amount of time! That's why we'll take a very short fragment of the text. Below, I took around 10 000 characters, because I can leave my laptop overnight and let it compute; during tutorial, you should take at most 1 000. It's rather important that your character vector starts with the beginning of a new word. The ending is less important. 

```{r}
fragmentalice <- alice[10025:15000]
stem.trained <- baumWelch(stem.hmm, fragmentalice)$hmm
```

The Baum Welch algorithm has now trained our model. Take a look at the inferred transition probabilities:

```{r}
stem.trained$transProbs
```

To get the most probable vector of states, we need to use the Viterbi algorithm. 

```{r}
stemalice <- viterbi(stem.trained, fragmentalice)
```

We can now visualize our results. We'll highlight the stems by capital letters.

```{r}
capitalice <- ifelse(stemalice=="Stem", str_to_upper(as.character(fragmentalice)), as.character(fragmentalice)) 
str_c(capitalice, collapse="")
```

Our approach has failed spectacularly. There are at least two reasons for this: First, we've used a very short fragment of the text. Second, there are words which do not have suffixes. Another possible reason is that there are some other features of words which have a higher influence on the letter distribution. The BW algorithm infers *some* states which influence the letter distribution - but they may not be the stem and suffix! That's why we can never fully trust our results, and always need to check if they agree with our predictions.

Let's inspect reason 2: words with no suffixes. We'll add another state to our model, called "Other". Initialize the matrix: 

```{r}
stem2.M <- matrix(c(0.9, 0.1, 0.0, 0.2, 0.6, 0.2, 0.4, 0.0, 0.6), byrow=T, ncol=3, dimnames=list(c("Stem", "Suffix", "Other"), c("Stem", "Suffix", "Other")))  # used as a starting point for BW algorithm
stem2.M
```

Initialize starting point and emmission probabilities:

```{r}
stem2.Start <- c(1, 0, 0)
stem2.Y <- matrix(1, nrow=3, ncol=length(phonemes), dimnames=list(c("Stem", "Suffix", "Other"), phonemes)) 
stem2.Y[1, " "] <- 0 # no whitespace at stem
stem2.Y[2, " "] <- 3 # slightly higher probability of emmission of whitespace at the suffix
stem2.Y[1, ] <- stem2.Y[1,]/sum(stem2.Y[1, ])
stem2.Y[2, ] <- stem2.Y[2, ]/sum(stem2.Y[2, ])
stem2.Y[3, ] <- stem2.Y[3, ]/sum(stem2.Y[3, ])
```

Initialize the model:

```{r}
stem2.hmm <- initHMM(States=c("Stem", "Suffix", "Other"),
                     Symbols=phonemes,
                     startProbs = stem2.Start,
                     transProbs=stem2.M,
                     emissionProbs=stem2.Y)
```

Train the model.

```{r}
stem2.trained <- baumWelch(stem2.hmm, fragmentalice)$hmm
```

Predict the states: 

```{r}
stemalice2 <- viterbi(stem2.trained, fragmentalice)
```

Now, highlight the "Suffix" state. But, be careful: Just because we've named it "suffix", doesn't mean that it actually corresponds to suffixes. The names of states don't matter for BW algorithm. Since the algorithm doesn't know what "Suffix" or "Other" mean, it might as well swap those classes. That's why we need to inspect all the classes.

```{r}
capitalice2 <- ifelse(stemalice2=="Suffix", str_to_upper(as.character(fragmentalice)), as.character(fragmentalice)) 
str_c(capitalice2, collapse="")
```

Now, it's much better! In many words the "Suffix" class actually corresponds to the suffixes. But, the words with no suffixes, like "że" or "się", are not that well predicted:

```{r}
capitalice2 <- ifelse(stemalice2=="Other", str_to_upper(as.character(fragmentalice)), as.character(fragmentalice)) 
str_c(capitalice2, collapse="")
```

The "Stem" class is also rather poorly predicted, but it agrees in some words. Using a larger training set (e.g. 200 000 letters) could potentially improve the results, but we would wait for hours for the algorithm to finish.

```{r}
capitalice2 <- ifelse(stemalice2=="Stem", str_to_upper(as.character(fragmentalice)), as.character(fragmentalice)) 
str_c(capitalice2, collapse="")
```

Let's predict the stems for the beginning of the text. Note that we don't need to train the model again - it's already trained on the other part of the text. Just use the Viterbi algorithm.

```{r}
stemalice3 <- viterbi(stem2.trained, alice[50030:51000])
```

Highlight the suffixes: 

```{r}
capitalice3 <- ifelse(stemalice3=="Suffix", str_to_upper(as.character(alice[50030:51000])), as.character(alice[50030:51000])) 
str_c(capitalice3, collapse="")
```

Stems:

```{r}
capitalice3 <- ifelse(stemalice3=="Stem", str_to_upper(as.character(alice[50030:51000])), as.character(alice[50030:51000])) 
str_c(capitalice3, collapse="")
```

Others:

```{r}
capitalice3 <- ifelse(stemalice3=="Other", str_to_upper(as.character(alice[50030:51000])), as.character(alice[50030:51000])) 
str_c(capitalice3, collapse="")
```

Another spectacular failure! The model works at first, but then classifies everything as stems. One of the reasons for this is what we've been dealing with before: *overfitting*. We've performed an extensive training on a small dataset, which has some particular features. In particular, there is no letter "h" in the training set. That's why the probability of emmiting this letter is always 0 in the trained model.

```{r}
unique(fragmentalice)
```

```{r}
stem2.trained$emissionProbs[, 7:15]
```

The letter "h" is present in the test set - and this is exactly the beginning of the spectacular failure. After this letter, each sequence of states has an equal probability: zero. The algorithm doesn't know what to do and just assumes that everything is a stem.

In this situation, one of the ways to improve our model is to make it worse. We need to force it to learn less from the training set. We will do this by specifying a maximum number of iterations of the BW algorithm. This will also allow us to use a larger training set.

```{r}
longeralice <- alice[10025:40000]
stem2.trained <- baumWelch(stem2.hmm, longeralice, maxIterations=10)$hmm
```

Let's check the new emission matrix. 

```{r}
stem2.trained$emissionProbs[, 7:15]
```

```{r}
stem2.trained$transProbs
```

The probability of emmission of the letter "h" is now non-zero. Let's check how our model performs now:

```{r}
stemalice3 <- viterbi(stem2.trained, alice[50030:51000])
```

Highlight the suffixes: 

```{r}
capitalice3 <- ifelse(stemalice3=="Suffix", str_to_upper(as.character(alice[50030:51000])), as.character(alice[50030:51000])) 
str_c(capitalice3, collapse="")
```

Stems:

```{r}
capitalice3 <- ifelse(stemalice3=="Stem", str_to_upper(as.character(alice[50030:51000])), as.character(alice[50030:51000])) 
str_c(capitalice3, collapse="")
```

Others:

```{r}
capitalice3 <- ifelse(stemalice3=="Other", str_to_upper(as.character(alice[50030:51000])), as.character(alice[50030:51000])) 
str_c(capitalice3, collapse="")
```

Now, it's much better (even though still not great). To get the best results, we would need to experiment some more with our parameters, initial conditions, and the model itself. 

The message of this tutorial: Hidden Markov Models are an art rather than a recipe. There is no simple theorem on how to make a good model and train it properly. We usually need to experiment with several solutions and know some tricks. 

