library(keras)
library(tibble)
# Data
boston_housing <- dataset_boston_housing()

c(train_data, train_labels) %<-% boston_housing$train
c(test_data, test_labels) %<-% boston_housing$test
paste0("Training entries: ", length(train_data), ", labels: ", length(train_labels))
train_data[1, ] # Display sample features, notice the different scales

column_names <- c('CRIM', 'ZN', 'INDUS', 'CHAS', 'NOX', 'RM', 'AGE',
                  'DIS', 'RAD', 'TAX', 'PTRATIO', 'B', 'LSTAT')
train_df <- as_tibble(train_data)
colnames(train_df) <- column_names
train_df
train_labels[1:10] # Display first 10 entries


############################################
# Feature engineering

# Test data is *not* used when calculating the mean and std.
# Normalize training data
train_data <- scale(train_data)

# Use means and standard deviations from training set to normalize test set
col_means_train <- attr(train_data, "scaled:center")
col_stddevs_train <- attr(train_data, "scaled:scale")
test_data <- scale(test_data, center = col_means_train, scale = col_stddevs_train)

train_data[1, ] # First training sample, normalized

######################################
# OLS
olstrain = data.frame(cbind(train_labels, train_data))
olstest  = data.frame(cbind(test_labels, test_data))
colnames(olstrain) = c('y', 'CRIM', 'ZN', 'INDUS', 'CHAS', 'NOX', 'RM', 'AGE',
                       'DIS', 'RAD', 'TAX', 'PTRATIO', 'B', 'LSTAT')
colnames(olstest) = c('y', 'CRIM', 'ZN', 'INDUS', 'CHAS', 'NOX', 'RM', 'AGE',
                      'DIS', 'RAD', 'TAX', 'PTRATIO', 'B', 'LSTAT')
ols = lm(y~.,data=olstrain)
summary(ols)
olspred = predict(ols, newdata=olstest)

library(Metrics)
mae(test_labels, olspred)

##################################
# Keras

#lrate=0.02
lrate=0.1

# OLS-like results are sensitive to choice of LR and EPOCHS (both must match)

epochs = 100
build_model <- function() {
  model <- keras_model_sequential() %>%

    # OLSisch
    layer_dense(units = 1, input_shape = dim(train_data)[2])

    # Simple NN
    #layer_dense(units = 64, activation = "relu", input_shape = dim(train_data)[2]) %>%
    #layer_dense(units = 1)

    model %>% compile(
      loss = "mse",
      optimizer = optimizer_rmsprop(lr=lrate),
      metrics = list("mean_absolute_error")
    )

  model
}

model <- build_model()
#model %>% summary()

# Fit the model and store training stats
history <- model %>% fit(
  train_data,
  train_labels,
  epochs = epochs,
  validation_split = 0.2,
  verbose = 1
)

c(loss, kmae) %<-% (model %>% evaluate(test_data, test_labels, verbose = 0))

print("======================================================")
print(paste0("Learning rate: ", lrate))
print(paste0("Keras: Mean absolute error on test set: $", sprintf("%.2f", kmae*1000)))
print(paste0("Keras mean val_mae:                     $", sprintf("%.2f", mean(history$metrics$val_mean_absolute_error[30:100])*1000)))
print(paste0("O L S: Mean absolute error on test set: $", sprintf("%.2f", mae(test_labels, olspred)*1000)))
print(paste0("Diff.                                 : $", sprintf("%.2f", kmae*1000 - mae(test_labels, olspred)*1000)))