Supporting code for HEA

---
title: "Supplementary_Code_HEA_Analysis"
author: "Ross Kunz"
date: "1/21/2020"
output: word_document
---
  
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

Copyright 2021, Battelle Energy Alliance, LLC

This document is included for a reproducible demonstration of all of the machine learning methods applied in the paper.

# Packages required:

```{r}
library(R6) # required for the class structure
library(ggplot2) # required for plotting
library(ggpubr) # required for plotting
library(mgcv) # required for the generalized linear model
library(neuralnet) # required for multi-layer perceptron
library(randomForest) # required for random forest
library(e1071) # required for support vector machine
```


# Preprocessing the matminer data based on the Wen Alloys.

```{r}
matminer_pre_process = function(file_path, splitSet = '10training90test'){
  all_terms = c("Yang.delta",                                            
                "Yang.omega",                                                  
                "APE.mean",                                                      
                "Radii.local.mismatch",                                    
                "Radii.gamma",                                                  
                "Configuration.entropy",                                         
                "Atomic.weight.mean",                                            
                "Lambda.entropy",                                               
                "Electronegtivity.delta",                                       
                "Electronegativity.local.mismatch",                             
                "VEC.mean",                                                   
                "Mixing.enthalpy.structure.fcc" ,                                
                "Mean.cohesive.energy",                                        
                "Interant.electrons",                                            
                "Interant.s.electrons",                                        
                "Shear.modulus.mean",                                            
                "Shear.modulus.delta",                                          
                "Shear.modulus.local.mismatch",                                 
                "Shear.modulus.strength.model",                                  
                "APE.mean_APE.mean",                                             
                "Configuration.entropy_Configuration.entropy",                  
                "Lambda.entropy_Lambda.entropy",                                 
                "VEC.mean_VEC.mean",                                        
                "Mixing.enthalpy.structure.fcc_Configuration.entropy" ,          
                "Mixing.enthalpy.structure.fcc_Electronegativity.local.mismatch",
                "Mixing.enthalpy.structure.fcc_VEC.mean",                       
                "Mixing.enthalpy.structure.fcc_Mixing.enthalpy.structure.fcc",  
                "Mean.cohesive.energy_Mixing.enthalpy.structure.fcc",            
                "Mean.cohesive.energy_Mean.cohesive.energy",                    
                "Interant.electrons_Radii.local.mismatch",                      
                "Interant.electrons_Configuration.entropy",                    
                "Interant.electrons_Electronegativity.local.mismatch",          
                "Interant.electrons_Mean.cohesive.energy",                       
                "Interant.electrons_Interant.electrons"  ,                       
                "Shear.modulus.mean_Electronegtivity.delta",                     
                "Shear.modulus.mean_Mean.cohesive.energy" ,                      
                "Shear.modulus.mean_Interant.s.electrons" ,                      
                "Shear.modulus.local.mismatch_Interant.electrons",               
                "Shear.modulus.strength.model_Shear.modulus.mean") 
  
  fiveElements = read.csv(file_path)
  total_samples = dim(fiveElements)[1]
  fiveElements = fiveElements[-which(names(fiveElements) %in% c("composition", 'Weight.Fraction', 'Atomic.Fraction'))]
  
  response = fiveElements$Elastic.modulus
  originalSet = fiveElements[,-which(names(fiveElements) %in% c("Elastic.modulus", 'Interant.f.electrons', 'Interant.p.electrons','Interant.d.electrons', 'Temperature','Weight.Fraction', 'Atomic.Fraction',
                                                                'Mixing.enthalpy', 'Mixing.enthalpy.structure.bcc', 'Mixing.enthalpy.structure.hcp', 'Total.weight'))]
  response = response[-which(originalSet$Radii.local.mismatch == 0)]
  originalSet = originalSet[-which(originalSet$Radii.local.mismatch == 0),]
  
  originalSet = originalSet[-which(response == min(response)),]
  response = response[-which(response == min(response))]
  
  originalSet$Yang.delta[which(is.infinite(originalSet$Yang.delta) == T)] = 0
  elementNames = which(names(originalSet) %in% c("Fe", "Cr", "Ni", "Cu", "Co"))
  
  elementMatrix = originalSet[, elementNames]
  justElements = as.matrix(elementMatrix)
  predElement = scale(elementMatrix)
  elementMatrix$elastic.modulus = response
  
  originalSet = originalSet[,-elementNames]
  if(splitSet == '70training30test'){
    set.seed(1234)
    tempDim = dim(originalSet)[1]
    trainSamples = sample(1:tempDim,size = round(.7 * tempDim), replace = F)
  }else if(splitSet == 'CrNi'){
    trainSamples = unique(c(which(elementMatrix$Cr == 0), which(elementMatrix$Ni == 0)))
  }else if(splitSet == '40training60test'){
    set.seed(1234)
    tempDim = dim(originalSet)[1]
    trainSamples = sample(1:tempDim,size = round(.4 * tempDim), replace = F)
  }else if(splitSet == '20training80test'){
    set.seed(1234)
    tempDim = dim(originalSet)[1]
    trainSamples = sample(1:tempDim,size = round(.2 * tempDim), replace = F)
  }else if(splitSet == '10training90test'){
    set.seed(1234)
    tempDim = dim(originalSet)[1]
    trainSamples = sample(1:tempDim,size = round(.1 * tempDim), replace = F)
  }
  
  trainingData = originalSet[trainSamples, ]
  trainingMean = apply(trainingData, 2, mean)
  trainingSD = apply(trainingData, 2, sd)
  trainingData = sweep(trainingData,2,trainingMean,'-')
  trainingData = as.matrix(sweep(trainingData,2,trainingSD,'/'))
  
  trainingResponse = response[trainSamples]
  trainingResponse_mean = mean(trainingResponse)
  trainingResponse_sd = sd(trainingResponse)
  trainingResponse = (trainingResponse - trainingResponse_mean) / trainingResponse_sd 
  
  trainingAsDF = as.data.frame(trainingData)
  trainingAsDF$response = (trainingResponse - trainingResponse_mean) / trainingResponse_sd
  
  testData = originalSet[-trainSamples,]
  testResponse = (response[-trainSamples] - trainingResponse_mean) / trainingResponse_sd
  testData = sweep(testData,2,trainingMean,'-')
  testData = as.matrix(sweep(testData,2,trainingSD,'/'))
  
  originalSet = sweep(originalSet,2,trainingMean,'-')
  originalSet = as.matrix(sweep(originalSet,2,trainingSD,'/'))
  
  originalSetWithElements = as.data.frame(cbind(originalSet, justElements))
  originalSetWithElements$response = (response - trainingResponse_mean) / trainingResponse_sd
  response = (response - mean(response)) / sd(response)
  
  elements_as_list = rep('a', dim(justElements)[1])
  for(i in 1:dim(elementMatrix)[1]){
    elements_as_list[i] = paste(colnames(justElements)[which(justElements[i,] != 0)], collapse = '')
  }
  num_elements = nchar(elements_as_list) / 2
  
  result = list('trainingData' = trainingData, 'testData' = testData,
                'trainingResponse' = trainingResponse, 'testResponse' = testResponse,
                'response' = response, 'justElements' = justElements,
                'AllData' = originalSetWithElements, 'elements_as_list' = elements_as_list,
                'num_elements' = num_elements, 'trainingSamples' = trainSamples, 
                'total_samples' = total_samples, 'training_mean' = trainingResponse_mean,
                'training_sd' = trainingResponse_sd)
}
```

# Class for applying regression to the HEA data

```{r}
HEA_Regression = R6::R6Class('HEA_Regression', list(
  # Training information
  train_predictors = NULL,
  train_response = NULL,
  train_model = NULL,
  train_fitted_values = NULL,
  train_residuals = NULL,
  train_rmse = NULL,
  train_r2 = NULL,
  # Testing information
  test_predictors = NULL,
  test_response = NULL,
  test_fitted_values = NULL,
  test_residuals = NULL,
  test_rmse = NULL,
  test_r2 = NULL,
  initialize = function(train_predictors, train_response, test_predictors = NULL, test_response = NULL){
    self$train_predictors = train_predictors
    self$train_response = train_response
    self$test_predictors = test_predictors
    self$test_response = test_response
  },
  fit = function(model){
    self$train_model = switch(
      model,
      'gam' = {
        require(mgcv)
        self$to_data_frame()
        gam(response~Miracle.Radius.mean + VEC.mean + 
              s(Mean.cohesive.energy, k =12) + 
              Shear.modulus.mean, 
            data = self$train_predictors, method = 'REML', select = TRUE)
      },
      'mlp' = {
        require(neuralnet)
        self$to_data_frame()
        neuralnet(response ~ .,data = self$train_predictors, hidden = 4, act.fct = 'tanh')
      },
      'ols' = {
        self$to_data_frame()
        lm(response ~ ., data = self$train_predictors)
      },
      'rf' = {
        require(randomForest)
        randomForest(self$train_predictors, self$train_response)
      },
      'svm' = {
        require(e1071)
        svm(self$train_predictors, self$train_response)
      }
    )
    self$train_fitted_values = as.numeric(predict(self$train_model, self$train_predictors))
    self$train_residuals = self$train_response - self$train_fitted_values
    self$train_rmse = sqrt(mean(self$train_residuals^2))
    self$train_r2 = 1 - sum(self$train_residuals^2) / sum((self$train_response - mean(self$train_response))^2)
    if(!is.null(self$test_predictors)){
      self$test_fitted_values = as.numeric(predict(self$train_model, self$test_predictors))
      self$test_residuals = self$test_response - self$test_fitted_values
      self$test_rmse = sqrt(mean(self$test_residuals^2))
      self$test_r2 = 1 - sum(self$test_residuals^2) / sum((self$test_response - mean(self$test_response))^2)
    }
    invisible(self)
  },
  plot_points = function(x = 'train_fitted_values', 
                         y = 'train_residuals', 
                         x_label = 'Fitted', 
                         y_label = 'Residuals', 
                         point_size = 1,
                         point_color = '#000000',
                         select_index = NULL,
                         save_plot_name = NULL, 
                         with_lm_line = T){
    
    plotDF = data.frame(X = self[[x]], Y = self[[y]])
    if(!is.null(select_index)){
      plotDF = plotDF[select_index, ]
    }
    g = ggplot(plotDF, aes(x = X, y = Y))+ 
      geom_point(size = point_size, colour = point_color)+ 
      ggplot2::theme(text = ggplot2::element_text(size = 24),
                     legend.position = "none",
                     panel.background = ggplot2::element_blank(), 
                     axis.line = ggplot2::element_line(colour = "black"),
                     panel.grid.minor = ggplot2::element_line(colour = "black", linetype = "dashed", size = .05),
                     panel.grid.major = ggplot2::element_line(colour = "black", linetype = "dashed", size = .1),
                     panel.border = ggplot2::element_rect(fill = "NA"),
                     legend.key.width = ggplot2::unit(1.5,"cm"),
                     legend.key =ggplot2:: element_rect(fill = "NA")
      ) +
      xlab(x_label) +
      ylab(y_label)
    
    if(with_lm_line){
      g = g + geom_smooth(method='lm', formula= y~x)
    }
    
    if(!is.null(save_plot_name)){
      ggsave(filename = save_plot_name, g, width = 10, height = 10)
    }
    return(g)
  },
  standardize_response = function(train_mean, train_sd){
    self$train_response = self$train_response * train_sd + train_mean
    self$train_fitted_values = self$train_fitted_values * train_sd + train_mean
    self$test_response = self$test_response * train_sd + train_mean
    self$test_fitted_values = self$test_fitted_values * train_sd + train_mean
  },
  to_data_frame = function(){
    self$train_predictors = as.data.frame(self$train_predictors)
    if(!is.null(self$test_predictors)){
      self$test_predictors = as.data.frame(self$test_predictors)
    }
    self$train_predictors$response = self$train_response
  }
))
```


# Analysis of the 111 plane

```{r}
clean_data = matminer_pre_process('data/matminer_elastic_modulus_111.csv')
trainingData = clean_data$trainingData
testData = clean_data$testData
trainingResponse = clean_data$trainingResponse
testResponse = clean_data$testResponse
response = clean_data$response
justElements = clean_data$justElements
elementMatrix = as.data.frame(justElements)
allData = clean_data$AllData
table(clean_data$num_elements[clean_data$trainingSamples])
```

## Fitting a linear model when only using the elements as predictors

```{r}
OLS = HEA_Regression$new(train_predictors = justElements, train_response = response)
OLS = OLS$fit(model = 'ols')
```

## Plotting the linear model results

```{r}
cbp1 <- c("#E69F00", "#56B4E9", "#009E73",
          "#D55E00", "#CC79A7")

temp_color = rep(0, dim(justElements)[1])
for(i in 1:length(temp_color)){
  temp_color[i] = length(which(justElements[i,] !=0))
}

g1 = OLS$plot_points(
  x = 'train_fitted_values',
  y = 'train_residuals',
  x_label = 'Fitted',
  y_label = 'Residuals', 
  select_index = which(temp_color == 2),
  point_size = 4,
  point_color = cbp1[1]
)

g2 = OLS$plot_points(
  x = 'train_fitted_values',
  y = 'train_residuals',
  x_label = 'Fitted',
  y_label = 'Residuals', 
  select_index = which(temp_color == 3),
  point_size = 4,
  point_color = cbp1[2]
)

g3 = OLS$plot_points(
  x = 'train_fitted_values',
  y = 'train_residuals',
  x_label = 'Fitted',
  y_label = 'Residuals', 
  select_index = which(temp_color == 4),
  point_size = 4,
  point_color = cbp1[3]
)

g4 = OLS$plot_points(
  x = 'train_fitted_values',
  y = 'train_residuals',
  x_label = 'Fitted',
  y_label = 'Residuals', 
  select_index = which(temp_color == 5),
  point_size = 4,
  point_color = cbp1[4]
)


g = ggarrange(
  g1,g2,g3,g4, 
  nrow=2, ncol=2, 
  labels = c("(a)", "(b)", "(c)", "(d)"), 
  font.label = list(size = 24), hjust =0
)
g
```

## Initializing a data frame for the results of all models

```{r}

result_df = data.frame('model' = c('rf', 'svm', 'mlp', 'gam'),
                       'train_rmse' = 0,
                       'test_rmse' = 0,
                       'test_R2' = 0)
all_residuals = 0
all_fitted = 0
all_methods = 'A'
```

## Plotting the results of random forest

```{r}
model_name = 'rf'
temp_obj = HEA_Regression$new(train_predictors = trainingData,
                        train_response = trainingResponse, 
                        test_predictors = testData,
                        test_response = testResponse)
temp_obj$fit(model_name)
model_position = which(result_df$model == model_name)
result_df$train_rmse[model_position] = temp_obj$train_rmse
result_df$test_rmse[model_position] = temp_obj$test_rmse
result_df$test_R2[model_position] = temp_obj$test_r2

temp_obj$plot_points(x='train_fitted_values', y = 'train_residuals',
               x_label = 'Fitted', y_label = 'Residuals')

temp_obj$standardize_response(clean_data$training_mean, clean_data$training_sd)

temp_obj$plot_points(x='test_fitted_values', y = 'test_response',
               x_label = 'Predicted', y_label = 'Observed')

all_residuals = c(all_residuals, temp_obj$train_residuals)
all_fitted = c(all_fitted, temp_obj$train_fitted_values)
all_methods = c(all_methods, rep(model_name, length(temp_obj$train_fitted_values)))
```

## Plotting the results of support vector machine
```{r}
model_name = 'svm'

temp_obj = HEA_Regression$new(train_predictors = trainingData,
                        train_response = trainingResponse, 
                        test_predictors = testData,
                        test_response = testResponse)
temp_obj$fit(model_name)
model_position = which(result_df$model == model_name)
result_df$train_rmse[model_position] = temp_obj$train_rmse
result_df$test_rmse[model_position] = temp_obj$test_rmse
result_df$test_R2[model_position] = temp_obj$test_r2

temp_obj$plot_points(x='train_fitted_values', y = 'train_residuals',
               x_label = 'Fitted', y_label = 'Residuals')

temp_obj$standardize_response(clean_data$training_mean, clean_data$training_sd)

temp_obj$plot_points(x='test_fitted_values', y = 'test_response',
               x_label = 'Predicted', y_label = 'Observed')

all_residuals = c(all_residuals, temp_obj$train_residuals)
all_fitted = c(all_fitted, temp_obj$train_fitted_values)
all_methods = c(all_methods, rep(model_name, length(temp_obj$train_fitted_values)))
```


## Plotting the results of multi-layer perceptron
```{r}
model_name = 'mlp'

temp_obj = HEA_Regression$new(train_predictors = trainingData,
                        train_response = trainingResponse, 
                        test_predictors = testData,
                        test_response = testResponse)
temp_obj$fit(model_name)
model_position = which(result_df$model == model_name)
result_df$train_rmse[model_position] = temp_obj$train_rmse
result_df$test_rmse[model_position] = temp_obj$test_rmse
result_df$test_R2[model_position] = temp_obj$test_r2

temp_obj$plot_points(x='train_fitted_values', y = 'train_residuals',
               x_label = 'Fitted', y_label = 'Residuals')

temp_obj$standardize_response(clean_data$training_mean, clean_data$training_sd)

temp_obj$plot_points(x='test_fitted_values', y = 'test_response',
               x_label = 'Predicted', y_label = 'Observed')

all_residuals = c(all_residuals, temp_obj$train_residuals)
all_fitted = c(all_fitted, temp_obj$train_fitted_values)
all_methods = c(all_methods, rep(model_name, length(temp_obj$train_fitted_values)))
```

## Plotting the results of generalized additive model
```{r}
model_name = 'gam'

temp_obj = HEA_Regression$new(train_predictors = trainingData,
                        train_response = trainingResponse, 
                        test_predictors = testData,
                        test_response = testResponse)
temp_obj$fit(model_name)
model_position = which(result_df$model == model_name)
result_df$train_rmse[model_position] = temp_obj$train_rmse
result_df$test_rmse[model_position] = temp_obj$test_rmse
result_df$test_R2[model_position] = temp_obj$test_r2

temp_obj$plot_points(x='train_fitted_values', y = 'train_residuals',
               x_label = 'Fitted', y_label = 'Residuals')

temp_obj$standardize_response(clean_data$training_mean, clean_data$training_sd)

temp_obj$plot_points(x='test_fitted_values', y = 'test_response',
               x_label = 'Predicted', y_label = 'Observed')

all_residuals = c(all_residuals, temp_obj$train_residuals)
all_fitted = c(all_fitted, temp_obj$train_fitted_values)
all_methods = c(all_methods, rep(model_name, length(temp_obj$train_fitted_values)))
```

# Comparison of all results on the 111 plane

```{r}
knitr::kable(result_df)
```


```{r}
plotDF = data.frame('Residuals' = all_residuals[-1],
                    'Fitted' = all_fitted[-1],
                    'Method' = all_methods[-1])

temp_obj = HEA_Regression$new(train_predictors = trainingData,
                        train_response = trainingResponse, 
                        test_predictors = allData,
                        test_response = allData$response)
temp_obj$fit('gam')
temp_obj$standardize_response(clean_data$training_mean, clean_data$training_sd)
predicted_val = temp_obj$test_fitted_values
all_resid = temp_obj$test_residuals
```

## Errors per composition

```{r, echo=F}
elements = c('Fe', 'Cr', 'Ni', 'Cu', 'Co')

names(elementMatrix) = elements
elements_as_list = rep('a', dim(elementMatrix)[1])

for(i in 1:dim(elementMatrix)[1]){
  elements_as_list[i] = paste(names(elementMatrix)[which(elementMatrix[i,] != 0)], collapse = '')
}

composition_results = data.frame(composition = sort(unique(elements_as_list)),
                                 rmse = 0,
                                 num_samp = 0)
for(i in 1:dim(composition_results)[1]){
  temp_loc = which(elements_as_list == composition_results$composition[i])
  composition_results$rmse[i] = round(sqrt(mean(all_resid[temp_loc]^2)), 4)
  composition_results$num_samp[i] = length(temp_loc)
}

names(elementMatrix) = elements
elements_as_list = rep('a', dim(elementMatrix)[1])

for(i in 1:dim(elementMatrix)[1]){
  elements_as_list[i] = paste(names(elementMatrix)[which(elementMatrix[i,] != 0)], collapse = '')
}

number_of_elements = (nchar(elements_as_list) / 2)

composition_results = composition_results[order(nchar(composition_results$composition)),]
composition_results

```

## Plotting fits per composition

```{r}

g1 = temp_obj$plot_points(
  x = 'test_fitted_values',
  y = 'test_response',
  x_label = 'Predicted',
  y_label = 'Observed', 
  select_index = which(number_of_elements == 2),
  point_size = 4,
  point_color = cbp1[1]
)

g2 = temp_obj$plot_points(
  x = 'test_fitted_values',
  y = 'test_response',
  x_label = 'Predicted',
  y_label = 'Observed', 
  select_index = which(number_of_elements == 3),
  point_size = 4,
  point_color = cbp1[2]
)


g3 = temp_obj$plot_points(
  x = 'test_fitted_values',
  y = 'test_response',
  x_label = 'Predicted',
  y_label = 'Observed', 
  select_index = which(number_of_elements == 4),
  point_size = 4,
  point_color = cbp1[3]
)

g4 = temp_obj$plot_points(
  x = 'test_fitted_values',
  y = 'test_response',
  x_label = 'Predicted',
  y_label = 'Observed', 
  select_index = which(number_of_elements == 5),
  point_size = 4,
  point_color = cbp1[4]
)

g = ggarrange(g1,g2,g3,g4, nrow=2, ncol=2, labels = c("(a)", "(b)", "(c)", "(d)"), font.label = list(size = 20), hjust =0)
g

```


## Verifying the GAM residuals are minimized by each element

```{r}

df_error = as.data.frame(justElements)
df_error$Residuals = temp_obj$test_residuals
df_error$num_components = number_of_elements

fit = lm(Residuals ~ Fe + Cr + Ni + Cu + Co - 1, data = df_error)
summary(fit)
fit$coefficients * clean_data$training_sd
```

## Errors by element

```{r}
element_names = colnames(justElements)

errors = rep(0, length(element_names))
error_by_element = list()
for(i in 1:length(element_names)){
  errors[i] = sqrt(mean(all_resid[grep(element_names[i], elements_as_list)]^2))
  error_by_element[[i]] = all_resid[grep(element_names[i], elements_as_list)]
}


error_by_composition = list()
for(i in 2:5){
  error_by_composition[[i -1]] = all_resid[which(nchar(elements_as_list)/2 == i)]
}

names(errors) = element_names
names(error_by_element) = element_names
errors 
```

## Testing differences in errors between Fe and all other elements

```{r}
temp_error = c(error_by_element$Cu^2, error_by_element$Cr^2, error_by_element$Ni^2, error_by_element$Co^2)
t.test(error_by_element$Fe^2, temp_error)
```

## Testing differences in errors between Fe and Ni

```{r}
t.test(error_by_element$Fe^2, error_by_element$Ni^2)
```

## Mean error per element

```{r}
lapply(error_by_element, mean)
```
## Testing differences in error from the 2 element system and all other compositions

```{r}
temp_error = c(error_by_composition[[2]], error_by_composition[[3]], error_by_composition[[4]])
t.test(error_by_composition[[1]], temp_error)
```


## Table of model results for all planes

```{r}
result_df = data.frame('model' = c('rf', 'svm', 'mlp', 'gam'),
                       'train_rmse_100' = 0,
                       'train_rmse_110' = 0,
                       'train_rmse_111' = 0,
                       'test_rmse_100' = 0,
                       'test_rmse_110' = 0,
                       'test_rmse_111' = 0,
                       'test_R2_100' = 0,
                       'test_R2_110' = 0,
                       'test_R2_111' = 0)

clean_data = matminer_pre_process('data/matminer_elastic_modulus_100.csv')
trainingData = clean_data$trainingData
testData = clean_data$testData
trainingResponse = clean_data$trainingResponse
testResponse = clean_data$testResponse
response = clean_data$response
justElements = clean_data$justElements
elementMatrix = as.data.frame(justElements)
allData = clean_data$AllData


model_name = 'rf'
temp_obj = HEA_Regression$new(train_predictors = trainingData,
                        train_response = trainingResponse, 
                        test_predictors = testData,
                        test_response = testResponse)
temp_obj$fit(model_name)
model_position = which(result_df$model == model_name)
result_df$train_rmse_100[model_position] = temp_obj$train_rmse
result_df$test_rmse_100[model_position] = temp_obj$test_rmse
result_df$test_R2_100[model_position] = temp_obj$test_r2

model_name = 'svm'
temp_obj = HEA_Regression$new(train_predictors = trainingData,
                        train_response = trainingResponse, 
                        test_predictors = testData,
                        test_response = testResponse)
temp_obj$fit(model_name)
model_position = which(result_df$model == model_name)
result_df$train_rmse_100[model_position] = temp_obj$train_rmse
result_df$test_rmse_100[model_position] = temp_obj$test_rmse
result_df$test_R2_100[model_position] = temp_obj$test_r2

model_name = 'mlp'
temp_obj = HEA_Regression$new(train_predictors = trainingData,
                        train_response = trainingResponse, 
                        test_predictors = testData,
                        test_response = testResponse)
temp_obj$fit(model_name)
model_position = which(result_df$model == model_name)
result_df$train_rmse_100[model_position] = temp_obj$train_rmse
result_df$test_rmse_100[model_position] = temp_obj$test_rmse
result_df$test_R2_100[model_position] = temp_obj$test_r2


model_name = 'gam'
temp_obj = HEA_Regression$new(train_predictors = trainingData,
                        train_response = trainingResponse, 
                        test_predictors = testData,
                        test_response = testResponse)
temp_obj$fit(model_name)
model_position = which(result_df$model == model_name)
result_df$train_rmse_100[model_position] = temp_obj$train_rmse
result_df$test_rmse_100[model_position] = temp_obj$test_rmse
result_df$test_R2_100[model_position] = temp_obj$test_r2




clean_data = matminer_pre_process('data/matminer_elastic_modulus_110.csv')
trainingData = clean_data$trainingData
testData = clean_data$testData
trainingResponse = clean_data$trainingResponse
testResponse = clean_data$testResponse
response = clean_data$response
justElements = clean_data$justElements
elementMatrix = as.data.frame(justElements)
allData = clean_data$AllData

model_name = 'rf'
temp_obj = HEA_Regression$new(train_predictors = trainingData,
                        train_response = trainingResponse, 
                        test_predictors = testData,
                        test_response = testResponse)
temp_obj$fit(model_name)
model_position = which(result_df$model == model_name)
result_df$train_rmse_110[model_position] = temp_obj$train_rmse
result_df$test_rmse_110[model_position] = temp_obj$test_rmse
result_df$test_R2_110[model_position] = temp_obj$test_r2

model_name = 'svm'
temp_obj = HEA_Regression$new(train_predictors = trainingData,
                        train_response = trainingResponse, 
                        test_predictors = testData,
                        test_response = testResponse)
temp_obj$fit(model_name)
model_position = which(result_df$model == model_name)
result_df$train_rmse_110[model_position] = temp_obj$train_rmse
result_df$test_rmse_110[model_position] = temp_obj$test_rmse
result_df$test_R2_110[model_position] = temp_obj$test_r2

model_name = 'mlp'
temp_obj = HEA_Regression$new(train_predictors = trainingData,
                        train_response = trainingResponse, 
                        test_predictors = testData,
                        test_response = testResponse)
temp_obj$fit(model_name)
model_position = which(result_df$model == model_name)
result_df$train_rmse_110[model_position] = temp_obj$train_rmse
result_df$test_rmse_110[model_position] = temp_obj$test_rmse
result_df$test_R2_110[model_position] = temp_obj$test_r2


model_name = 'gam'
temp_obj = HEA_Regression$new(train_predictors = trainingData,
                        train_response = trainingResponse, 
                        test_predictors = testData,
                        test_response = testResponse)
temp_obj$fit(model_name)
model_position = which(result_df$model == model_name)
result_df$train_rmse_110[model_position] = temp_obj$train_rmse
result_df$test_rmse_110[model_position] = temp_obj$test_rmse
result_df$test_R2_110[model_position] = temp_obj$test_r2

clean_data = matminer_pre_process('data/matminer_elastic_modulus_111.csv')
trainingData = clean_data$trainingData
testData = clean_data$testData
trainingResponse = clean_data$trainingResponse
testResponse = clean_data$testResponse
response = clean_data$response
justElements = clean_data$justElements
elementMatrix = as.data.frame(justElements)
allData = clean_data$AllData

model_name = 'rf'
temp_obj = HEA_Regression$new(train_predictors = trainingData,
                        train_response = trainingResponse, 
                        test_predictors = testData,
                        test_response = testResponse)
temp_obj$fit(model_name)
model_position = which(result_df$model == model_name)
result_df$train_rmse_111[model_position] = temp_obj$train_rmse
result_df$test_rmse_111[model_position] = temp_obj$test_rmse
result_df$test_R2_111[model_position] = temp_obj$test_r2

model_name = 'svm'
temp_obj = HEA_Regression$new(train_predictors = trainingData,
                        train_response = trainingResponse, 
                        test_predictors = testData,
                        test_response = testResponse)
temp_obj$fit(model_name)
model_position = which(result_df$model == model_name)
result_df$train_rmse_111[model_position] = temp_obj$train_rmse
result_df$test_rmse_111[model_position] = temp_obj$test_rmse
result_df$test_R2_111[model_position] = temp_obj$test_r2

model_name = 'mlp'
temp_obj = HEA_Regression$new(train_predictors = trainingData,
                        train_response = trainingResponse, 
                        test_predictors = testData,
                        test_response = testResponse)
temp_obj$fit(model_name)
model_position = which(result_df$model == model_name)
result_df$train_rmse_111[model_position] = temp_obj$train_rmse
result_df$test_rmse_111[model_position] = temp_obj$test_rmse
result_df$test_R2_111[model_position] = temp_obj$test_r2


model_name = 'gam'
temp_obj = HEA_Regression$new(train_predictors = trainingData,
                        train_response = trainingResponse, 
                        test_predictors = testData,
                        test_response = testResponse)
temp_obj$fit(model_name)
model_position = which(result_df$model == model_name)
result_df$train_rmse_111[model_position] = temp_obj$train_rmse
result_df$test_rmse_111[model_position] = temp_obj$test_rmse
result_df$test_R2_111[model_position] = temp_obj$test_r2

result_df[,-1] = round((result_df[,-1]), 2)
knitr::kable(result_df)
```

## Visual description of the differences in Youngs modulus per plane

```{r}
t1 = read.csv('data/matminer_elastic_modulus_100.csv')
t2 = read.csv('data/matminer_elastic_modulus_110.csv')
t3 = read.csv('data/matminer_elastic_modulus_111.csv')

fit1 = lm(t3$Elastic.modulus ~ t1$Elastic.modulus - 1)
summary(fit1)


fit2 = lm(t3$Elastic.modulus ~ t2$Elastic.modulus - 1)
summary(fit2)


temp_data = data.frame(elastic_modulus = c(sort(t1$Elastic.modulus)*fit1$coefficients,
                                           sort(t2$Elastic.modulus)*fit2$coefficients,
                                           sort(t3$Elastic.modulus)),
                       sample = rep(1:length(t1$Elastic.modulus),3),
                       Plane = rep(c('100','110', '111'), each = length(t1$Elastic.modulus)))
ggplot2::ggplot(temp_data, ggplot2::aes(sample, elastic_modulus, color = Plane)) +
  ggplot2::geom_point() + 
  ggplot2::xlab('Sample Index') +
  ggplot2::ylab("Young's Modulus") +
  ggplot2::theme(text = ggplot2::element_text(size = 24),
                 #legend.position = c(.2,.8),
                 panel.background = ggplot2::element_blank(), 
                 axis.line = ggplot2::element_line(colour = "black"),
                 panel.grid.minor = ggplot2::element_line(colour = "black", linetype = "dashed", size = .05),
                 panel.grid.major = ggplot2::element_line(colour = "black", linetype = "dashed", size = .1),
                 panel.border = ggplot2::element_rect(fill = "NA"),
                 legend.key.width = ggplot2::unit(1.5,"cm"),
                 legend.key = ggplot2::element_rect(fill = "NA")
  )


```