movielens %>%
  dplyr::count(movieId) %>%
  ggplot(aes(n)) +
  geom_histogram(bins = 30, color = "black", fill = "cornflowerblue") +
  scale_x_log10() +
  ggtitle("Movies")

movielens %>% 
  dplyr::count(userId) %>% 
  ggplot(aes(n)) +
  geom_histogram(bins = 30, color = "black", fill = "cornflowerblue") +
  scale_x_log10() +
  ggtitle("Users")

movielens <- movielens %>%
  mutate(year = as.numeric(format(as.Date(as.POSIXct(timestamp, 
                                                     origin = "1970-01-01"), 
                                          format ="%Y/%m/%d"), "%Y")))

pattern <- "\\s[(][0-9]{4}[)]"

# Add Release Year and Years Since Release to the data
movielens <- movielens %>% 
  mutate(Release_Year = as.numeric(substr(str_extract(string = title, 
                                                     pattern = pattern),3,6)),
         years_since_release = year - Release_Year)

# Double check the above achieved the desired effect
movielens %>% head(n=10) %>% knitr::kable()

movielens %>% group_by(Release_Year) %>%
  summarise(mean_Release_Year_Rating = mean(rating)) %>%
  ggplot(aes(Release_Year,mean_Release_Year_Rating)) +
  geom_point() + geom_smooth(method = "lm")

movielens %>% group_by(years_since_release) %>%
  summarise(mean_rating_by_years_since_release = mean(rating)) %>%
  ggplot(aes(years_since_release, mean_rating_by_years_since_release)) + 
  geom_point() + geom_smooth(method = "lm")

naive_rmse <- RMSE(test_set$rating, mu_hat)
naive_rmse %>% knitr::kable()

RMSE_results <- tibble(Method = "The Average", RMSE = naive_rmse) 
RMSE_results %>% knitr::kable()

train_set %>%
  group_by(movieId) %>%
  summarize(avg_rating = mean(rating)) %>% 
  filter(n()>=100) %>%
  ggplot(aes(avg_rating)) +
  geom_histogram(bins = 30, color = "black", fill = "cornflowerblue")

# fit <- lm(rating ~ as.factor(userId), data = train_set)

# In order to keep the code as clean as possible, _hat notation has been dropped

mu <- mean(train_set$rating)

movie_avgs <- train_set %>% 
  group_by(movieId) %>% 
  summarize(b_i = mean(rating - mu))

predicted_ratings <- mu + test_set %>%
  left_join(movie_avgs, by='movieId') %>% 
  .$b_i

model_2_rmse <- RMSE(predicted_ratings, test_set$rating)
RMSE_results <- bind_rows(RMSE_results, 
                          tibble(Method = "Movie Effect Model", 
                                 RMSE = model_2_rmse ))
RMSE_results %>% knitr::kable()

train_set %>%
  group_by(userId) %>%
  summarize(avg_rating = mean(rating)) %>% 
  filter(n()>=100) %>%
  ggplot(aes(avg_rating)) +
  geom_histogram(bins = 30, color = "black", fill = "cornflowerblue")

user_avgs <- train_set %>%
  left_join(movie_avgs, by='movieId') %>%
  group_by(userId) %>%
  summarize(b_u = mean(rating - mu - b_i))

predicted_ratings <- test_set %>%
  left_join(movie_avgs, by='movieId') %>%
  left_join(user_avgs, by='userId') %>%
  mutate(pred = mu + b_i + b_u) %>% 
  .$pred

model_3_rmse <- RMSE(predicted_ratings, test_set$rating)

RMSE_results <- bind_rows(RMSE_results,
                          tibble(Method="Movie + User Effects Model",
                                 RMSE = model_3_rmse ))
RMSE_results %>% knitr::kable()

time_avgs <- train_set %>%
  left_join(movie_avgs, by='movieId') %>%
  left_join(user_avgs, by='userId') %>% 
  group_by(years_since_release) %>%
  summarize(b_t = mean(rating - mu - b_i - b_u))

genre_avgs <- train_set %>% 
  left_join(movie_avgs, by='movieId') %>%
  left_join(user_avgs, by='userId') %>%
  left_join(time_avgs, by='years_since_release') %>% 
  group_by(genres) %>% 
  summarize(b_g = mean(rating - mu - b_i - b_u - b_t))

predicted_ratings <- test_set %>%
  left_join(movie_avgs, by='movieId') %>%
  left_join(user_avgs, by='userId') %>%
  left_join(time_avgs, by='years_since_release') %>% 
  left_join(genre_avgs, by='genres') %>% 
  mutate(pred = mu + b_i + b_u + b_t + b_g) %>% 
  .$pred

model_4_rmse <- RMSE(predicted_ratings, test_set$rating)

RMSE_results <- bind_rows(RMSE_results,
                          tibble(Method=
                                   "Movie + User + Time + Genre Effects Model",
                                 RMSE = model_4_rmse ))
RMSE_results %>% knitr::kable()

train_set %>% 
  left_join(movie_avgs, by='movieId') %>%
  mutate(residual = rating - (mu + b_i)) %>%
  arrange(desc(abs(residual))) %>%
  select(title, residual) %>% slice(1:10) %>% knitr::kable()

movie_titles <- train_set %>%
  select(movieId, title) %>%
  distinct()

ten_best <- movie_avgs %>% 
  left_join(movie_titles, by="movieId") %>%
  arrange(desc(b_i)) %>%
  select(title, b_i) %>%
  slice(1:10,) %>% knitr::kable(name="Best")

ten_worst <- movie_avgs %>% 
  left_join(movie_titles, by="movieId") %>%
  arrange(b_i) %>%
  select(title, b_i) %>%
  slice(1:10) %>% knitr::kable(name = "Worst")

ten_best

ten_worst

train_set %>% 
  dplyr::count(movieId) %>%
  left_join(movie_avgs) %>%
  left_join(movie_titles, by="movieId") %>%
  arrange(desc(b_i)) %>%
  select(title, b_i, n) %>%
  slice(1:10) %>% knitr::kable()

train_set %>% 
  dplyr::count(movieId) %>%
  left_join(movie_avgs) %>%
  left_join(movie_titles, by="movieId") %>%
  arrange(b_i) %>%
  select(title, b_i, n) %>%
  slice(1:10) %>% knitr::kable()

lambda <- 4

mu <- mean(train_set$rating)

movie_reg_avgs <- train_set %>%
  group_by(movieId) %>%
  summarize(b_i = sum(rating - mu)/(n()+lambda), n_i = n())

tibble(original = movie_avgs$b_i,
       regularlized = movie_reg_avgs$b_i,
       n = movie_reg_avgs$n_i) %>%
  ggplot(aes(original, regularlized, size=sqrt(n))) +
  geom_point(shape=1, alpha=0.5)

train_set %>%
  dplyr::count(movieId) %>%
  left_join(movie_reg_avgs) %>%
  left_join(movie_titles, by="movieId") %>%
  arrange(desc(b_i)) %>%
  select(title, b_i, n) %>%
  slice(1:10) %>% knitr::kable()

train_set %>%
  dplyr::count(movieId) %>%
  left_join(movie_reg_avgs) %>%
  left_join(movie_titles, by="movieId") %>%
  arrange(b_i) %>%
  select(title, b_i, n) %>%
  slice(1:10) %>% knitr::kable()

predicted_ratings <- test_set %>%
  left_join(movie_reg_avgs, by='movieId') %>%
  mutate(pred = mu + b_i) %>%
  .$pred

model_3_rmse <- RMSE(predicted_ratings, test_set$rating)

RMSE_results <- bind_rows(RMSE_results,
                          tibble(Method="Regularized Movie Effect Model",
                                 RMSE = model_3_rmse ))
RMSE_results %>% knitr::kable()

lambdas <- seq(0, 10, 0.25)

rmses <- sapply(lambdas, function(l){
  mu <- mean(train_set$rating)
  
  b_i <- train_set %>%
    group_by(movieId) %>%
    summarize(b_i = sum(rating - mu)/(n()+l))

  b_u <- train_set %>%
    left_join(b_i, by="movieId") %>%
    group_by(userId) %>%
    summarize(b_u = sum(rating - b_i - mu)/(n()+l))
  
  b_t <- train_set %>%
    left_join(b_i, by='movieId') %>%
    left_join(b_u, by='userId') %>%
    group_by(years_since_release) %>%
    summarize(b_t = sum(rating - b_i - b_u - mu)/(n()+l))
  
  b_g <- train_set %>% 
    left_join(b_i, by='movieId') %>%
    left_join(b_u, by='userId') %>%
    left_join(b_t, by='years_since_release') %>% 
    group_by(genres) %>% 
    summarize(b_g = sum(rating - b_i - b_u - b_t - mu)/(n()+l))

  predicted_ratings <- test_set %>%
    left_join(b_i, by = "movieId") %>%
    left_join(b_u, by = "userId") %>%
    left_join(b_t, by = 'years_since_release') %>% 
    left_join(b_g, by = 'genres') %>% 
    mutate(pred = mu + b_i + b_u + b_t + b_g) %>%
    .$pred

return(RMSE(predicted_ratings, test_set$rating))
})

qplot(lambdas, rmses)

lambda <- lambdas[which.min(rmses)]
lambda %>% knitr::kable()

RMSE_results <- bind_rows(
  RMSE_results,
  tibble(Method="Regularized Movie + User + Time + Genre Model",
         RMSE = min(rmses)))
RMSE_results %>% knitr::kable()

# Sample the train set for movies that have received 600 or more 
# ratings and users who users who have rated 600 or more movies

train_small <- train_set %>%
  group_by(movieId) %>%
  filter(n() >= 600) %>% 
  ungroup() %>%
  group_by(userId) %>%
  filter(n() >= 600) %>% ungroup()

# Convert to matrix
y <- train_small %>%
  select(userId, movieId, rating) %>%
  spread(movieId, rating) %>%
  as.matrix()

# keep rownames
rownames(y)<- y[,1]

y <- y[,-1]

# select distinct titles
movie_titles <- train_set %>% 
  select(movieId, title) %>%
  distinct()

# add colnames
colnames(y) <- with(movie_titles, title[match(colnames(y), movieId)])

# calc vectorization
y <- sweep(y, 2, colMeans(y, na.rm=TRUE))
y <- sweep(y, 1, rowMeans(y, na.rm=TRUE))

# display correlation between times of movies
m_1 <- "Godfather, The (1972)"
m_2 <- "Godfather: Part II, The (1974)"
p1 <- qplot(y[ ,m_1], y[,m_2], xlab = m_1, ylab = m_2)

m_1 <- "Godfather, The (1972)"
m_3 <- "Goodfellas (1990)"
p2 <- qplot(y[ ,m_1], y[,m_3], xlab = m_1, ylab = m_3)

m_4 <- "You've Got Mail (1998)"
m_5 <- "Sleepless in Seattle (1993)"
p3 <- qplot(y[ ,m_4], y[,m_5], xlab = m_4, ylab = m_5)
gridExtra::grid.arrange(p1, p2 ,p3, ncol = 3)

cor(y[, c(m_1, m_2, m_3, m_4, m_5)], use="pairwise.complete") %>%
  knitr::kable()

y[is.na(y)] <- 0
y <- sweep(y, 1, rowMeans(y))
pca <- prcomp(y)

dim(pca$rotation) %>% knitr::kable()

dim(pca$x) %>% knitr::kable()

plot(pca$sdev)

var_explained <- cumsum(pca$sdev^2/sum(pca$sdev^2))
plot(var_explained)

pcs <- data.frame(pca$rotation, name = colnames(y))
pcs %>%  
  ggplot(aes(PC1, PC2)) + 
  geom_point() +
  geom_text_repel(aes(PC1, PC2, label=name),
                  data = filter(pcs,
                                PC1 < -.01 | 
                                  PC1 > .01 | 
                                  PC2 < -0.075 | 
                                  PC2 > .01))

pcs %>% select(name, PC1) %>% arrange(PC1) %>% slice(1:10) %>% tibble() %>%
  knitr::kable()

pcs %>% select(name, PC1) %>% arrange(desc(PC1)) %>% slice(1:10) %>% 
  tibble() %>% knitr::kable()

pcs %>% select(name, PC2) %>% arrange(PC2) %>% slice(1:10) %>% tibble() %>% 
  knitr::kable()

pcs %>% select(name, PC2) %>% arrange(desc(PC2)) %>% slice(1:10) %>% 
  tibble() %>% knitr::kable()

# A quick read of the help files states that I must have matrices in the 
# required realRatingMatrix type. So first up I will have to coerce the 
# movielens data to a matrix
y <- movielens %>%
  select(userId, movieId, rating) %>%
  spread(movieId, rating) %>%
  as.matrix()

# setting names properly
rownames(y)<- y[,1]

y <- y[,-1]

# Finally coerce to the recommenderlab specific "realRatingMatrix"
movielenseM <- as(y, "realRatingMatrix")

# split the data by a 80/20 split which mimics the second split of the big_train 
# set to train and test sets

e <- evaluationScheme(movielenseM, method="split", 
                      train=0.8, k=4, given=80, goodRating = 3)

# create a libmf recommender using training data
r <- Recommender(getData(e, "train"), "LIBMF")
                                                            
# create predictions for the test data using known ratings (the given param 
# from the evaluation scheme)
p <- predict(r, getData(e, "known"), type="ratings")

# calculate the average RMSE when User specific
acc <- calcPredictionAccuracy(p, getData(e, "unknown"), byUser=TRUE)
acc <- mean(acc[,1])

# Append to the RMSE results tibble
RMSE_results <- bind_rows(RMSE_results,
                    tibble(Method = "Matrix Factorization - Recommenderlab",
                           RMSE = acc))
# Show the RMSE improvement
RMSE_results %>%
  knitr::kable()

set.seed(1985, sample.kind = "Rounding") # This is a randomized algorithm

# This takes four minutes to run and should not be a concern for anyone on
# a modern laptop. 

# Convert the train and test sets into recosystem input format
train_data <- with(train_set, data_memory(user_index = userId,
                                          item_index = movieId,
                                          rating = rating))
test_data <- with(test_set, data_memory(user_index = userId,
                                        item_index = movieId,
                                        rating = rating))


# Create the model object
r <- recosystem::Reco()

# Select the best tuning parameters using cross-validation
opts <- r$tune(train_data, opts = list(dim = c(10, 20, 30),
                                       costp_l2 = c(0.01, 0.1),         
                                       costq_l2 = c(0.01, 0.1),         
                                       costp_l1 = 0,                    
                                       costq_l1 = 0,                    
                                       lrate    = c(0.01, 0.1),         
                                       nthread  = 8, 
                                       niter    = 10,                   
                                       verbose  = FALSE)) 
# Train the algorithm
r$train(train_data, opts = c(opts$min, 
                             nthread = 8, 
                             niter = 100, 
                             verbose = FALSE))

# Calculate the predicted values
reco_pred <- r$predict(test_data, out_memory())

RMSE_results <- bind_rows(RMSE_results,
                    tibble(Method = 
                             "Parallel Matrix Factorization - Recosystem",
                           RMSE = RMSE(test_set$rating, reco_pred)))
# Show the RMSE improvement
RMSE_results %>%
  knitr::kable()

# I use the lambda that reduced the RMSE the most in the train / test model to
# avoid over training of the final validation model. It is stored as lambda.

mu <- mean(big_train$rating)
  
b_i <- big_train %>%
  group_by(movieId) %>%
  summarize(b_i = sum(rating - mu)/(n()+lambda))

b_u <- big_train %>%
  left_join(b_i, by="movieId") %>%
  group_by(userId) %>%
  summarize(b_u = sum(rating - b_i - mu)/(n()+lambda))
  
b_t <- big_train %>%
  left_join(b_i, by='movieId') %>%
  left_join(b_u, by='userId') %>%
  group_by(years_since_release) %>%
  summarize(b_t = sum(rating - b_i - b_u - mu)/(n()+lambda))
  
b_g <- big_train %>% 
  left_join(b_i, by='movieId') %>%
  left_join(b_u, by='userId') %>%
  left_join(b_t, by='years_since_release') %>% 
  group_by(genres) %>% 
  summarize(b_g = sum(rating - b_i - b_u - b_t - mu)/(n()+lambda))

predicted_ratings <- validation %>%
  left_join(b_i, by = "movieId") %>%
  left_join(b_u, by = "userId") %>%
  left_join(b_t, by = 'years_since_release') %>% 
  left_join(b_g, by = 'genres') %>% 
  mutate(pred = mu + b_i + b_u + b_t + b_g) %>%
  .$pred

RMSE_Reg <- RMSE(predicted_ratings, validation$rating)

Final_Validation_RMSE_results_table <- tibble(Method = 
"Regularized Movie + User + Time + Genre Model", RMSE = RMSE_Reg) 

Final_Validation_RMSE_results_table %>% knitr::kable()

# I have already coerced the movielens data into the required format and it is
# stored as movielenseM which is a realRatingMatrix so there is no need to 
# repeat those steps. However, I do need to repeat the split. The recommenderlab 
# way is evaluationScheme and the code below shows a 90/10 split which 
# is the same proportion as the original big_train / validation split, which is
# also 90/10.
e <- evaluationScheme(movielenseM, method="split", 
                      train=0.9, k=1, given=90)

# create a libmf recommender using training data
r <- Recommender(getData(e, "train"), "LIBMF")

# create predictions for the test data using known ratings (the given param 
# from the evaluation scheme)
p <- predict(r, getData(e, "known"), type="ratings")

# calculate the average RMSE when User specific
acc <- calcPredictionAccuracy(p, getData(e, "unknown"), byUser=TRUE)
acc <- mean(acc[,1])

# Append to the RMSE results tibble
Final_Validation_RMSE_results_table <- bind_rows(
  Final_Validation_RMSE_results_table,
  tibble(Method = "Final Validation - Recommenderlab - Matrix Factorization",
         RMSE = acc))

# Show the Results
Final_Validation_RMSE_results_table %>% knitr::kable()

# Final Validation

set.seed(1986, sample.kind="Rounding")

train_data <- with(big_train, data_memory(user_index = userId, #big_train
                                    item_index = movieId,
                                    rating = rating))
test_data <- with(validation, data_memory(user_index = userId, #VALIDATION
                                          item_index = movieId,
                                          rating = rating))

r <- recosystem::Reco()

# Select the best tuning parameters using cross-validation
opts <- r$tune(train_data, opts = list(dim = c(10, 20, 30),
                                       costp_l2 = c(0.01, 0.1),         
                                       costq_l2 = c(0.01, 0.1),         
                                       costp_l1 = 0,                    
                                       costq_l1 = 0,                    
                                       lrate    = c(0.01, 0.1),         
                                       nthread  = 8, 
                                       niter    = 10,                   
                                       verbose  = FALSE)) 
# Train the algorithm
r$train(train_data, opts = c(opts$min, 
                             nthread = 8, 
                             niter = 100, 
                             verbose = FALSE))

# Calculate the predicted values
reco_final_pred <- r$predict(test_data, out_memory())

Final_Validation_RMSE_results_table <- bind_rows(
  Final_Validation_RMSE_results_table,
  tibble(Method = 
           "Final Validation - Recosystem - True Parallel Matrix Factorization",
         RMSE = RMSE(validation$rating, reco_final_pred)))
Final_Validation_RMSE_results_table %>% knitr::kable()