## if keras package has not been installed, install it.                         
if(!require(keras)) {install.packages("keras"); library(keras)}
#[1]-----
tmp <- read.csv("DatafingerprintSub.csv", header=TRUE)
## choose target                                                                
target1 <- scale(tmp[,"pKa"])
## remove the target variables             
data1_org <- as.matrix(tmp[, -1])
## remove variables with no variation                                           
is_nv <- apply(data1_org, 2, function(xxx){
  return(var(xxx) == 0)
})
data1_v <- data1_org[, !is_nv]
## scale each variable                                                          
data1 <- apply(data1_v, 2, scale)
nc <- ncol(data1); nr <- nrow(data1)
#[2]-----
model4 <- keras_model_sequential()
model4 %>%
  layer_dense(units=81, input_shape=nc) %>%
  layer_activation_leaky_relu(0.1) %>%
  layer_dense(units=27) %>%
  layer_activation_leaky_relu(0.1) %>%
  layer_dense(units=9) %>%
  layer_activation_leaky_relu(0.1) %>%
  layer_dense(units=3) %>%
  layer_activation_leaky_relu(0.1) %>%
  layer_dense(units=1)
pred4 <- predict(model4, data1)
model4 %>% compile(loss="mean_squared_error", optimizer=optimizer_adagrad(learning_rate=0.01))
#[3]-----
n_sample <- nrow(data1)
idx_tst <- sample(1:n_sample, n_sample*0.1)
idx_trn <- (1:n_sample)[-idx_tst]
data1_tst <- data1[idx_tst,]
data1_trn <- data1[idx_trn,]
target1_tst <- target1[idx_tst]
target1_trn <- target1[idx_trn]
#[4] training                                                                                                                                        
callbacks = list(callback_early_stopping(monitor = "val_loss",
                                         patience = 5, restore_best_weights = TRUE))
hist4 <- fit(model4, data1_trn, target1_trn, epochs=100,
             validation_data=list(data1_tst, target1_tst), verbose=1,
             callbacks = callbacks)
pred4_trn <- predict(model4, data1_trn)
pred4_tst <- predict(model4, data1_tst)
#[5]-----
matplot(cbind(hist4$metrics$loss, hist4$metrics$val_loss),
        type="l", xlab="Epoch", ylab="Loss", main="4 hidden layers")
legend("topright", legend=c("Train", "Test"), lty=1:2, col=1:2)
#[6]-----
plot(target1_trn, pred4_trn, xlab="Target", ylab="Predicted", main="4 hidden layers")
points(target1_tst, pred4_tst, col=2, pch=3)
abline(0, 1, lty=3, col=8)
legend("topleft", legend=c("Train", "Test"), pch=c(1, 3), col=1:2)