# Load Packages: library(data.table) library(Matrix) library(xgboost) # Load Data: #list.files("../input/champs-scalar-coupling") TrainData <- fread("../input/champs-scalar-coupling/train.csv") TestData <- fread("../input/champs-scalar-coupling/test.csv") InfoStructure <- fread("../input/champs-scalar-coupling/structures.csv") trainID <- TrainData$id testID <- TestData$id atom_properties <- data.table(atom = c("H", "C", "N", "O", "F"), electroneg = c(2.2, 2.55, 3.04, 3.44, 3.98)) BondingTypes <- levels(factor(c(TrainData$type, TestData$type) )) InfoStructure <- merge(InfoStructure, atom_properties) setindex(InfoStructure, molecule_name) # Summarize information on structure InfoStructure_1 <- InfoStructure names(InfoStructure_1)[c(1,3:7)] <- paste0(names(InfoStructure_1)[c(1,3:7)],"_1") InfoStructure_cross <- merge(InfoStructure, InfoStructure_1, by = c("molecule_name"), allow.cartesian = TRUE) names(InfoStructure_cross)[c(2:7)] <- paste0(names(InfoStructure_cross)[c(2:7)],"_0") rm(InfoStructure_1) InfoStructure_cross <- InfoStructure_cross[atom_index_0 != atom_index_1] InfoStructure_cross[, Dist := sqrt((x_1-x_0)^2+(y_1-y_0)^2+(z_1-z_0)^2)] InfoStructure_cross[, DistInv := 1/Dist] InfoStructure_cross[, DistInv3 := 1/Dist^3] InfoStructure_cross[, DifElNeg := electroneg_1-electroneg_0] InfoStructure_cross[, ElNegByDist := DifElNeg/Dist] MinDistAtomPairs <- InfoStructure_cross[, list(MinPairDist = min(Dist)), by=c("atom_0", "atom_1")] MinDistAtomPairs[atom_0 == "H" & atom_1 == "H", MinPairDist := 0] MinDistAtomPairs[MinPairDist > 2, MinPairDist := 0] InfoStructure_cross <- merge(InfoStructure_cross, MinDistAtomPairs, by = c("atom_0", "atom_1")) paste("dimensions of full blown structure information:", paste(dim(InfoStructure_cross), collapse = " x ")) # Summarydata by Atom: AtomInfo <- InfoStructure_cross[, list(MinDist = min(Dist), MaxDist = max(Dist), MeanDist = mean(Dist), MedDist = median(Dist), MinDistInv = min(DistInv), MaxDistInv = max(DistInv), MeanDistInv = mean(DistInv), MedDistInv = median(DistInv), MinDistInv3 = min(DistInv3), MaxDistInv3 = max(DistInv3), MeanDistInv3 = mean(DistInv3), MedDistInv3 = median(DistInv3)), by=c("molecule_name", "atom_index_0")] AtomInfo2 <- InfoStructure_cross[Dist <= MinPairDist*1.4, list(Neighbors = length(Dist), NeighborH = sum(atom_1=="H"), NeighborC = sum(atom_1=="C"), NeighborN = sum(atom_1=="N"), NeighborO = sum(atom_1=="O"), NeighborF = sum(atom_1=="F"), MeanDifElNeg = mean(DifElNeg), MinDifElNeg = min(DifElNeg), MaxDifElNeg = max(DifElNeg), MeanElNegDist = mean(ElNegByDist), MinElNegDist = min(ElNegByDist), MaxElNegDist = max(ElNegByDist), MeanDistN = mean(Dist), MinDistN = min(Dist), MaxDistN = max(Dist), MeanDistInvN = mean(DistInv), MinDistInvN = min(DistInv), MaxDistInvN = max(DistInv), MeanDistInv3N = mean(DistInv3), MinDistInv3N = min(DistInv3), MaxDistInv3N = max(DistInv3), MinRelDist = min(Dist/MinPairDist), MaxRelDist = max(Dist/MinPairDist), MeanRelDist = mean(Dist/MinPairDist), MedRelDist = median(Dist/MinPairDist)), by=c("molecule_name", "atom_index_0")] AtomInfo <- merge(AtomInfo, AtomInfo2, by=c("molecule_name", "atom_index_0")) paste("number of Atom-wise parameters:", ncol(AtomInfo)) rm(AtomInfo2, MinDistAtomPairs) # Summarydata by molecule: ToDo: number of C, H, N... evtl. auch in AtomInfo2 MolInfo <- InfoStructure[, list(NumberAtoms = length(unique(atom_index)), MinElNeg = min(electroneg), MeanElNeg = mean(electroneg), MaxElNeg = max(electroneg), SdElNeg = sd(electroneg)), by=c("molecule_name")] MolInfo <- merge(MolInfo, InfoStructure_cross[, list(MinDistMol = min(Dist), MaxDistMol = max(Dist), MeanDistMol = mean(Dist), MedDistMol = median(Dist)), by=c("molecule_name")], by=c("molecule_name")) MolInfo <- merge(MolInfo, AtomInfo[, list(MinMinDifElNeg = min(MinDifElNeg), MaxMinDifElNeg = max(MinDifElNeg), MinMaxDifElNeg = min(MaxDifElNeg), MaxMaxDifElNeg = max(MaxDifElNeg), MinMinElNegDist = min(MinElNegDist), MaxMinElNegDist = max(MinElNegDist), MinMaxElNegDist = min(MaxElNegDist), MaxMaxElNegDist = max(MaxElNegDist)), by=c("molecule_name")], by=c("molecule_name")) InputData <- rbind(TrainData[, 1:5], TestData) setindex(InputData, molecule_name) paste("number of molecule-wise parameters:", ncol(MolInfo)) # Aggregate data by bonding, Atom and molecule to Input dataset: InputData <- merge(InputData, InfoStructure_cross[,c("molecule_name", "atom_index_0", "atom_index_1","Dist","DistInv","DistInv3","DifElNeg","ElNegByDist")], by=c("molecule_name", "atom_index_0", "atom_index_1") ) InputData <- merge(InputData, MolInfo, by = "molecule_name" ) FirstNewCol <- ncol(InputData)+1 InputData <- merge(InputData, AtomInfo, by.x = c("molecule_name", "atom_index_1"), by.y = c("molecule_name", "atom_index_0") ) names(InputData)[FirstNewCol:ncol(InputData)] <- paste0(names(InputData)[FirstNewCol:ncol(InputData)], "_1") FirstNewCol <- ncol(InputData)+1 InputData <- merge(InputData, AtomInfo, by = c("molecule_name", "atom_index_0") ) names(InputData)[FirstNewCol:ncol(InputData)] <- paste0(names(InputData)[FirstNewCol:ncol(InputData)], "_0") setorder(InputData, id) rm(InfoStructure_cross, AtomInfo, MolInfo); gc() # Plot the 4 most highly correlated features to the scalar coupling constant for each bonding type par(mfrow = c(2, 2)) for (i in 1:length(BondingTypes)) { tmpCorData <- merge(InputData[id %in% trainID & type == BondingTypes[i], c(4, 6:101)], TrainData[id %in% trainID & type == BondingTypes[i], list(id, scalar_coupling_constant)], by = "id") corMatrix <- cor(tmpCorData) for (j in 1:4) { corCol <- order(corMatrix[98, 1:97], decreasing = TRUE)[j] plot(unlist(tmpCorData[, 98])~as.data.frame(tmpCorData)[, corCol], xlab = names(tmpCorData)[corCol], ylab = names(tmpCorData)[98], main = paste0("#", j, " correlation")) } mtext(paste("bonding type:", BondingTypes[i]), outer = TRUE, line=-1.5) }