{"cells":[{"metadata":{},"cell_type":"markdown","source":"# This code generates and saves features from [Giba's R kernel (v8, -1.17 LB)](https://www.kaggle.com/titericz/giba-r-data-table-simple-features-1-17-lb)"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"rm(list=ls())\nlibrary(data.table)\nlibrary(tictoc)\nset.seed(0)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train <- fread(\"../input/train.csv\")\ntest  <- fread(\"../input/test.csv\" )\nsub <- fread('../input/sample_submission.csv')\n\n#Use ID as groupKfold = 5 folds\ntrain[, ID := as.integer(.GRP)%%4 , by=\"molecule_name\" ]\ntest[ , ID := as.integer(.GRP)%%4 , by=\"molecule_name\" ]\ndim(train)\ndim(test)\nhead(train,30)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Merge atom positions"},{"metadata":{"trusted":true},"cell_type":"code","source":"structures_0 <- fread(\"../input/structures.csv\")\nstructures_1 <- copy(structures_0)\n\nsetnames(structures_0, names(structures_0)[-1], paste0(names(structures_0)[-1], \"_0\"))\nsetnames(structures_1, names(structures_1)[-1], paste0(names(structures_1)[-1], \"_1\"))\nsetnames(structures_0, names(structures_0)[-c(1:2)], paste0(\"structure_\", names(structures_0)[-c(1:2)]))\nsetnames(structures_1, names(structures_1)[-c(1:2)], paste0(\"structure_\", names(structures_1)[-c(1:2)]))\n\ntrain <- merge(train, structures_0, all.x = T,\n               by = c(\"molecule_name\", \"atom_index_0\"))\ntrain <- merge(train, structures_1, all.x = T,\n               by = c(\"molecule_name\", \"atom_index_1\"))\ntest  <- merge(test , structures_0, all.x = T,\n               by = c(\"molecule_name\", \"atom_index_0\"))\ntest  <- merge(test , structures_1, all.x = T,\n               by = c(\"molecule_name\", \"atom_index_1\"))\n\ntrain_test <- rbind(train, test, use.names = T, fill = T)\n\nrm(train, test); gc()\ndim(train_test)\nhead(train_test)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Encode atoms"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_test[,typei:=.GRP, by=\"type\"]\n\ntrain_test[, N1 := length(unique(atom_index_1)) ,\n           by=c(\"molecule_name\",\"atom_index_0\")  ]\ntrain_test[, N2 := length(unique(atom_index_0)) ,\n           by=c(\"molecule_name\",\"atom_index_1\")  ]\n\ndt0 <- train_test[, paste(type,collapse=\" \"),\n                  by=c(\"molecule_name\",\"atom_index_0\") ]\ndt1 <- train_test[, paste(type,collapse=\" \"),\n                  by=c(\"molecule_name\",\"atom_index_1\") ]\nsetnames(dt0,colnames(dt0), c(\"molecule_name\",\"atom_index\",\"V1\") )\nsetnames(dt1,colnames(dt1), c(\"molecule_name\",\"atom_index\",\"V1\") )\n\ndt <- rbind( dt0, dt1  )\ndt[, link := .GRP , by=\"V1\" ]\ndt[, N := .N, by=\"link\"]\ndt <- dt[, list(group=link[1]) , keyby=\"V1\" ]\ndt0[, link := dt[J(dt0$V1)]$group ]\ndt1[, link := dt[J(dt1$V1)]$group ]\nsetkeyv(dt0, c(\"molecule_name\",\"atom_index\"))\nsetkeyv(dt1, c(\"molecule_name\",\"atom_index\"))\n\ntrain_test[, link0 := dt0[J(train_test$molecule_name, train_test$atom_index_0)]$link ]\ntrain_test[, link1 := dt1[J(train_test$molecule_name, train_test$atom_index_1)]$link ]\ntrain_test[, linkN := .N, by=c(\"link0\",\"link1\") ]\nhead(train_test)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Distances"},{"metadata":{"trusted":true},"cell_type":"code","source":"# difference and distance features\ntrain_test[, \":=\"(dist_xyz = sqrt((structure_x_0 - structure_x_1) ^ 2 +\n                                    (structure_y_0 - structure_y_1) ^ 2 +\n                                    (structure_z_0 - structure_z_1) ^ 2)\n)]\n\ntrain_test[, inv_dist0 := 1/sum(1/dist_xyz^3) ,\n           by=c(\"molecule_name\",\"atom_index_0\") ]\ntrain_test[, inv_dist1 := 1/sum(1/dist_xyz^3) ,\n           by=c(\"molecule_name\",\"atom_index_1\") ]\ntrain_test[, inv_distP := (inv_dist0*inv_dist1)/(inv_dist0+inv_dist1) ]\nhead(train_test)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Some distance/radius based features"},{"metadata":{"trusted":true},"cell_type":"code","source":"atomic_radius <- data.table(\n  atom = c('H','C','N','O','F'),\n  R = c(0.38,0.77,0.75,0.73,0.71),\n  E = c(2.2,2.55,3.04,3.44,3.98)\n  ,key=\"atom\"\n)\ntrain_test[, R0 := atomic_radius[J(train_test$structure_atom_0)]$R  ]\ntrain_test[, R1 := atomic_radius[J(train_test$structure_atom_1)]$R  ]\ntrain_test[, E0 := atomic_radius[J(train_test$structure_atom_0)]$E  ]\ntrain_test[, E1 := atomic_radius[J(train_test$structure_atom_1)]$E  ]\n\ntrain_test[, inv_dist0R := 1/sum(1/((dist_xyz-R0-R1)^2)) ,\n           by=c(\"molecule_name\",\"atom_index_0\") ]\ntrain_test[, inv_dist1R := 1/sum(1/((dist_xyz-R0-R1)^2)) ,\n           by=c(\"molecule_name\",\"atom_index_1\") ]\ntrain_test[, inv_distPR := (inv_dist0R*inv_dist1R)/(inv_dist0R+inv_dist1R) ]\n\ntrain_test[, inv_dist0E := 1/sum(1/((dist_xyz*(0.5*E0+0.5*E1))^2)) ,\n           by=c(\"molecule_name\",\"atom_index_0\") ]\ntrain_test[, inv_dist1E := 1/sum(1/((dist_xyz*(0.5*E0+0.5*E1))^2)) ,\n           by=c(\"molecule_name\",\"atom_index_1\") ]\ntrain_test[, inv_distPE := (inv_dist0E*inv_dist1E)/(inv_dist0E+inv_dist1E) ]\n\ntrain_test[, linkM0:= mean(inv_dist0), by=c(\"link0\") ]\ntrain_test[, linkM1:= mean(inv_dist0), by=c(\"link1\") ]\ntrain_test[, linkM0:= linkM0 - inv_dist0 ]\ntrain_test[, linkM1:= linkM1 - inv_dist1 ]\n\nhead(train_test)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Simple Aggregations"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_test <- merge(train_test, train_test[, \n                                  .(\n                                   min_molecule_atom_0_dist_xyz  = min(dist_xyz),\n                                   mean_molecule_atom_0_dist_xyz = mean(dist_xyz),\n                                   max_molecule_atom_0_dist_xyz  = max(dist_xyz),\n                                   sd_molecule_atom_0_dist_xyz   = sd(dist_xyz)\n), .(molecule_name, atom_index_0)], by = c(\"molecule_name\", \"atom_index_0\"))\n\ntrain_test <- merge(train_test, train_test[, \n                                  .(\n                                   min_molecule_atom_1_dist_xyz  = min(dist_xyz),\n                                   mean_molecule_atom_1_dist_xyz = mean(dist_xyz),\n                                   max_molecule_atom_1_dist_xyz  = max(dist_xyz),\n                                   sd_molecule_atom_1_dist_xyz   = sd(dist_xyz)\n), .(molecule_name, atom_index_1)], by = c(\"molecule_name\", \"atom_index_1\"))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Coulomb, Van der Walls & Yukawa interactions"},{"metadata":{"trusted":true},"cell_type":"code","source":"dstruct = fread('../input/structures.csv')\ndstruct[,i:=.I]\nlist_atoms = dstruct[,.(molecule_name,atom_index,i)]\nmd = as.matrix(dstruct[,.(x,y,z)])\natom_type =unlist(dstruct[,atom])\nrm(dstruct)\nsystem.time(\n    {\n     aij<-merge(list_atoms,list_atoms,by='molecule_name',allow.cartesian=TRUE)\n     aij[,atom_index.y:=NULL]\n     colnames(aij)<-c('molecule_name','atom_index','i','j')\n     aij<-aij[i!=j]\n     aij[,atj:=atom_type[j]]\n    }\n)\nrs = md[aij[,i],]-md[aij[,j],]\nsystem.time(aij[,d:=sqrt(rowSums(rs*rs))])\nsystem.time(daij<-aij[,.(coulomb=sum(1/d),vander=sum(1/d/d/d),yukawa=sum(exp(-d)/d)),by=.(molecule_name,atom_index,atj)])  #fixed Yukawa by /d\nsystem.time(daij2<-dcast(daij,molecule_name+atom_index~atj,value.var=c('coulomb','yukawa','vander')))\ndaij2[is.na(daij2)]=0\n\nrm(list_atoms,atom_type,aij,rs,daij)\n\ndataCY=merge(train_test[,.(molecule_name,atom_index_0,atom_index_1)],daij2,by.x=c('molecule_name','atom_index_0'),by.y=c('molecule_name','atom_index'),sort=F)\ndataCY=merge(dataCY,daij2,by.x=c('molecule_name','atom_index_1'),by.y=c('molecule_name','atom_index'),sort=F)\n\nrm(daij2)\n\ntrain_test=cbind(train_test,dataCY[,-3,with=F])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Counts"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_test[train_test$structure_atom_1=='C',\n           countC0:= .N,\n           by=c(\"molecule_name\",\"atom_index_0\")  ]\ntrain_test[train_test$structure_atom_1=='H',\n           countH0:= .N,\n           by=c(\"molecule_name\",\"atom_index_0\")  ]\ntrain_test[train_test$structure_atom_1=='N',\n           countN0:= .N,\n           by=c(\"molecule_name\",\"atom_index_0\")  ]\ntrain_test[train_test$structure_atom_1=='C',\n           countC1 := .N,\n           by=c(\"molecule_name\",\"atom_index_1\")  ]\ntrain_test[train_test$structure_atom_1=='H',\n           countH1 := .N,\n           by=c(\"molecule_name\",\"atom_index_1\")  ]\ntrain_test[train_test$structure_atom_1=='N',\n           countN1 := .N,\n           by=c(\"molecule_name\",\"atom_index_1\")  ]\n\ntrain_test[train_test$structure_atom_1=='C',\n           distC0 := sum(dist_xyz)/countC0  ,\n           by=c(\"molecule_name\",\"atom_index_0\")  ]\ntrain_test[train_test$structure_atom_1=='H',\n           distH0 := sum(dist_xyz)/countH0  ,\n           by=c(\"molecule_name\",\"atom_index_0\")  ]\ntrain_test[train_test$structure_atom_1=='N',\n           distN0 := sum(dist_xyz)/countN0  ,\n           by=c(\"molecule_name\",\"atom_index_0\")  ]\ntrain_test[train_test$structure_atom_1=='C',\n           distC1 := sum(dist_xyz)/countC1  ,\n           by=c(\"molecule_name\",\"atom_index_0\")  ]\ntrain_test[train_test$structure_atom_1=='H',\n           distH1 := sum(dist_xyz)/countH1  ,\n           by=c(\"molecule_name\",\"atom_index_0\")  ]\ntrain_test[train_test$structure_atom_1=='N',\n           distN1 := sum(dist_xyz)/countN1  ,\n           by=c(\"molecule_name\",\"atom_index_0\")  ]\n\ntrain_test[, countC0:= NULL  ]\ntrain_test[, countH0:= NULL  ]\ntrain_test[, countN0:= NULL  ]\ntrain_test[, countC1:= NULL  ]\ntrain_test[, countH1:= NULL  ]\ntrain_test[, countN1:= NULL  ]\nhead(train_test)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Distance features by atom type"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_test[, pos := 1:.N, by=c(\"molecule_name\",\"atom_index_0\") ]\n\nsort_dist <- function( par, typ, k=\"H\" ){\n  res <- matrix( 999, nrow=length(par), ncol=13 )\n  par[ typ!=k ] <- 999\n  for( i in 1:nrow(res) ){\n    res[i, 1:length(par) ] <- par\n    res[i, i ] <- 999\n    res[i,] <- res[i, order(res[i,],decreasing=F) ]\n  }\n  res[res==999] <- NA\n  return( list( res[,1],res[,2],res[,3],res[,4],res[,5] ) )\n}\n\ndt <- train_test[, sort_dist(par=dist_xyz,typ=structure_atom_1,k=\"H\") ,\n                 keyby=c(\"molecule_name\",\"atom_index_0\") ]\ndt[, pos := 1:.N, by=c(\"molecule_name\",\"atom_index_0\") ]\nsetkeyv( dt, c(\"molecule_name\",\"atom_index_0\",\"pos\")  )\n\ntrain_test[, adH1 := dt[J(train_test$molecule_name, train_test$atom_index_0, train_test$pos)]$V1 ]\ntrain_test[, adH2 := dt[J(train_test$molecule_name, train_test$atom_index_0, train_test$pos)]$V2 ]\ntrain_test[, adH3 := dt[J(train_test$molecule_name, train_test$atom_index_0, train_test$pos)]$V3 ]\ntrain_test[, adH4 := dt[J(train_test$molecule_name, train_test$atom_index_0, train_test$pos)]$V4 ]\nhead(train_test)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"dt <- train_test[, sort_dist(par=dist_xyz,typ=structure_atom_1,k=\"C\") ,\n                 keyby=c(\"molecule_name\",\"atom_index_0\") ]\ndt[, pos := 1:.N, by=c(\"molecule_name\",\"atom_index_0\") ]\nsetkeyv( dt, c(\"molecule_name\",\"atom_index_0\",\"pos\")  )\n\ntrain_test[, adC1 := dt[J(train_test$molecule_name, train_test$atom_index_0, train_test$pos)]$V1 ]\ntrain_test[, adC2 := dt[J(train_test$molecule_name, train_test$atom_index_0, train_test$pos)]$V2 ]\ntrain_test[, adC3 := dt[J(train_test$molecule_name, train_test$atom_index_0, train_test$pos)]$V3 ]\ntrain_test[, adC4 := dt[J(train_test$molecule_name, train_test$atom_index_0, train_test$pos)]$V4 ]\nhead(train_test)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"dt <- train_test[, sort_dist(par=dist_xyz,typ=structure_atom_1,k=\"N\") ,\n                 keyby=c(\"molecule_name\",\"atom_index_0\") ]\ndt[, pos := 1:.N, by=c(\"molecule_name\",\"atom_index_0\") ]\nsetkeyv( dt, c(\"molecule_name\",\"atom_index_0\",\"pos\")  )\n\ntrain_test[, adN1 := dt[J(train_test$molecule_name, train_test$atom_index_0, train_test$pos)]$V1 ]\ntrain_test[, adN2 := dt[J(train_test$molecule_name, train_test$atom_index_0, train_test$pos)]$V2 ]\ntrain_test[, adN3 := dt[J(train_test$molecule_name, train_test$atom_index_0, train_test$pos)]$V3 ]\ntrain_test[, adN4 := dt[J(train_test$molecule_name, train_test$atom_index_0, train_test$pos)]$V4 ]\nhead(train_test)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Atom type counts"},{"metadata":{"trusted":true},"cell_type":"code","source":"structure <- fread('../input/structures.csv')\ndt <- structure[, .N, keyby= c(\"molecule_name\",\"atom\") ]\ndt1 <- dt[ atom=='C' ]\ndt2 <- dt[ atom=='H' ]\ndt3 <- dt[ atom=='N' ]\ndt4 <- dt[ atom=='F' ]\ndt5 <- dt[ atom=='O' ]\n\ntrain_test[, NC := dt1[J(train_test$molecule_name)]$N ]\ntrain_test[, NH := dt2[J(train_test$molecule_name)]$N ]\ntrain_test[, NN := dt3[J(train_test$molecule_name)]$N ]\ntrain_test[, NF := dt4[J(train_test$molecule_name)]$N ]\ntrain_test[, NO := dt5[J(train_test$molecule_name)]$N ]\nhead(train_test)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Drop useless features"},{"metadata":{"trusted":true},"cell_type":"code","source":"features <- colnames(train_test)[c(1:ncol(train_test))]\nfeatures <- features[ !features %in% c(\"id\",\"molecule_name\",\n                                       \"atom_index_0\",\"atom_index_1\",\n                                       \"type\",\"scalar_coupling_constant\",\n                                       \"ID\",\n                                       \"structure_atom_0\",\"structure_atom_1\",\n                                       \"structure_x_0\",\"structure_y_0\",\"structure_z_0\",\n                                       \"structure_x_1\",\"structure_y_1\",\"structure_z_1\",\n                                       \"typei\",\"pos\",\n                                       \"R0\",\"R1\",\"E0\",\"E1\"\n                                       ) ]\nlength(features);features\nhead(train_test[,features,with=F])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Write train_test as .csv"},{"metadata":{"trusted":true},"cell_type":"code","source":"write.csv(train_test[,features,with=F], file=gzfile(\"train_test_giba_117.csv.gz\"))","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"R","language":"R","name":"ir"},"language_info":{"mimetype":"text/x-r-source","name":"R","pygments_lexer":"r","version":"3.4.2","file_extension":".r","codemirror_mode":"r"}},"nbformat":4,"nbformat_minor":1}