{"cells":[{"metadata":{},"cell_type":"markdown","source":"**Contents:**\n* [1. Problem Statement](#first-bullet)\n* [2. Data Load](#second-bullet)\n* [3. Exploratory Data Analysis](#third-bullet)\n* [4. Data Preparation](#fourth-bullet)\n* [5. Clustering](#fifth-bullet)"},{"metadata":{},"cell_type":"markdown","source":"----\n<a class=\"anchor\" id=\"first-bullet\"></a>\n## 1. Problem Statement\n\n\nThe main target of this competition is to develop an algorithm that can predict the magnetic interaction between two atoms in a molecule (i.e., the **scalar coupling constant**). For more information about this constant, this [link](http://nmrwiki.org/wiki/index.php?title=Scalar_coupling) provides a great explanation.\n\nFor this competition, we shouldn't predict all the atom pairs in each molecule rather, only the pairs that are explicitly listed in the train and test files. For example, some molecules contain Fluorine (F), but we must not predict the scalar coupling constant for any pair that includes F.\n\nThe training and test splits are by molecule, so that no molecule in the training data is found in the test data.\n\n"},{"metadata":{},"cell_type":"markdown","source":"### 1.1. Available data\n \n* train.csv - the training set, consisting of 4658147 observations and 6 variables. The first column (molecule_name) is the name of the molecule where the coupling constant originates, the second (atom_index_0) and third column (atom_index_1) are the atom indices of the atom-pair creating the coupling. The fourth column (scalar_coupling_constant) is the scalar coupling constant that we want to be able to predict.\n\n* test.csv - the test set; same info as train, without the target variable. It consists of 2505542 observations and 5 variables.\n\n* sample_submission.csv - a sample submission file in the correct format.\n\n* structures.csv - this file contains the same information as the individual xyz structure files, but in a single file. It has 2358657 observations and 6 variables.\n\n\nThere are also some additional files that can be used to get more information about the Training set.\nThose are the following:\n* dipole_moments.csv - has 85003 observations and 4 variables. Contains the molecular electric dipole moments, which are three dimensional vectors that indicate the charge distribution in the molecule. The first column (molecule_name) are the names of the molecule, the second to fourth column are the X, Y and Z components respectively of the dipole moment. \n\n* magnetic_shielding_tensors.csv - contains the magnetic shielding tensors for all atoms in the molecules. The first column (molecule_name) contains the molecule name, the second column (atom_index) contains the index of the atom in the molecule, the third to eleventh columns contain the XX, YX, ZX, XY, YY, ZY, XZ, YZ and ZZ elements of the tensor/matrix respectively.\n\n* mulliken_charges.csv - contains the mulliken charges for all atoms in the molecules. The first column (molecule_name) contains the name of the molecule, the second column (atom_index) contains the index of the atom in the molecule, the third column (mulliken_charge) contains the mulliken charge of the atom.\n\n* potential_energy.csv - contains the potential energy of the molecules. The first column (molecule_name) contains the name of the molecule, the second column (potential_energy) contains the potential energy of the molecule.\n\n* scalar_coupling_contributions.csv - The scalar coupling constants in train.csv (or corresponding files) are a sum of four terms. scalar_coupling_contributions.csv contain all these terms. The first column (molecule_name) are the name of the molecule, the second (atom_index_0) and third column (atom_index_1) are the atom indices of the atom-pair, the fourth column indicates the type of coupling, the fifth column (fc) is the Fermi Contact contribution, the sixth column (sd) is the Spin-dipolar contribution, the seventh column (pso) is the Paramagnetic spin-orbit contribution and the eighth column (dso) is the Diamagnetic spin-orbit contribution.\n"},{"metadata":{},"cell_type":"markdown","source":"----\n<a class=\"anchor\" id=\"second-bullet\"></a>\n## 2. Data load"},{"metadata":{},"cell_type":"markdown","source":"### 2.1 Python libraries"},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":" ## Import Python libraries\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the \"../input/\" directory.\n# For example, running this (by clicking run or pressing Shift+Enter) will list the files in the input directory\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport matplotlib.pyplot as plt\n%matplotlib inline\nimport seaborn as sns\nsns.set(color_codes=True)\n\nimport scipy as sp #collection of functions for scientific computing and advance mathematics\nfrom scipy import stats\nfrom scipy.stats import norm, skew\nfrom scipy.special import boxcox1p\nfrom scipy.stats import boxcox_normmax\n\n\nfrom sklearn.cluster import KMeans\n\n\nfrom sklearn.model_selection import GroupKFold\nfrom sklearn.preprocessing import LabelEncoder\nfrom sklearn.ensemble import RandomForestRegressor\n\nimport os\nprint(os.listdir(\"../input\"))\n\n\n# warnings mute\nimport warnings\nwarnings.simplefilter('ignore')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### 2.2 Data Import"},{"metadata":{"trusted":true},"cell_type":"code","source":"train = pd.read_csv(\"../input/train.csv\")\ntest = pd.read_csv(\"../input/test.csv\")\nsub = pd.read_csv(\"../input/sample_submission.csv\")\nstructures = pd.read_csv(\"../input/structures.csv\")\ndipole = pd.read_csv(\"../input/dipole_moments.csv\")\npotential = pd.read_csv(\"../input/potential_energy.csv\")\nmulliken = pd.read_csv(\"../input/mulliken_charges.csv\")\ncontributions = pd.read_csv(\"../input/scalar_coupling_contributions.csv\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"----\n<a class=\"anchor\" id=\"third-bullet\"></a>\n## 3. Exploratory Data Analysis"},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"train.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"structures.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"dipole.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"potential.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"mulliken.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"contributions.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"----\n<a class=\"anchor\" id=\"fourth-bullet\"></a>\n## 4. Data preparation"},{"metadata":{},"cell_type":"markdown","source":"### 4.1 Table joins\n\nCreate a dataframe joining the information from different datasets."},{"metadata":{"trusted":true},"cell_type":"code","source":"df1 = pd.merge(train,  dipole)\ndf = pd.merge(df1, potential)\n\n# reorder columns: locate the target variable at the end\ndf = df[[\"id\", \"molecule_name\", \"atom_index_0\", \"atom_index_1\", \n         \"type\", \"X\", \"Y\", \"Z\", \"potential_energy\", \"scalar_coupling_constant\"]]\n\ndf.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# check if some null values were introduced\ndf.isnull().any()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Note:\nN = 3000\nprint(\"The original dataset has {} records, which is quite a lot.\\nIn some cases, we will select just {} rows to speed up the computation.\". format(df.shape[0], N))\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### 4.2 Groupping"},{"metadata":{"trusted":true},"cell_type":"code","source":"df.groupby(\"type\").median()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"structures.groupby(\"atom\").median()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#different molecules\nstructures[\"molecule_name\"].nunique()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df.groupby('type').count()['molecule_name'].sort_values().plot(kind='barh',\n                                                                color='steelblue',\n                                                               figsize=(15, 8),\n                                                               title='Count of Coupling Type in Train Set')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"structures.groupby('atom').count()['molecule_name'].sort_values().plot(kind='barh',\n                                                                color='steelblue',\n                                                               figsize=(15, 5),\n                                                               title='Count of Atom Type in Structures Set')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_plot = df.drop('id', 1)\ndf_plot.hist(figsize=(12,12))\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.pairplot(df_plot)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### 4.3 Normality check and correlations\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Testing for normal distribution hypothesis in numerical features\ntest_normality = lambda x: stats.shapiro(x.fillna(0))[1] < 0.01\nnumerical_features = [f for f in df.columns if df.dtypes[f] != 'object']\nnormal = pd.DataFrame(df[numerical_features])\nnormal = normal.apply(test_normality)\nprint(normal)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Calculate correlations\ncorr = df.corr(method='spearman')\nmask = np.zeros_like(corr)\nmask[np.triu_indices_from(mask)] = True\n# Heatmap\nplt.figure(figsize=(12, 7))\nsns.heatmap(corr,\n            vmax=.5,\n            mask=mask,\n            #annot=True, \n            fmt='.2f',\n            linewidths=.2, cmap=\"YlGnBu\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### 4.4 Study of the target variable"},{"metadata":{"trusted":true},"cell_type":"code","source":"target = df[\"scalar_coupling_constant\"]\n\n# let's get some stats on the 'scalar_coupling_constant' variable\nprint(\"Statistics for the scalar_coupling_constant training dataset:\\n\")\nprint(\"Minimum value: {:,.2f}\".format(np.min(target)))\nprint(\"Maximum value: {:,.2f}\".format(np.max(target)))\nprint(\"Mean value: {:,.2f}\".format(np.mean(target)))\nprint(\"Median value {:,.2f}\".format(np.median(target)))\nprint(\"Standard deviation: {:,.2f}\".format(np.std(target)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#  To get a visual of the outliers, let's create a box plot.\nsns.boxplot(y = target)\nplt.ylabel('target')\nplt.title(' ');\n\n# count number of outliers in the variable\nQ1 = target.quantile(0.25)\nQ3 = target.quantile(0.75)\nIQR = Q3 - Q1\nprint(\"IQR value: {}\\n# of outliers: {}\".format(\n    IQR,\n    ((target < (Q1 - 1.5 * IQR)) | (target > (Q3 + 1.5 * IQR))).sum()))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# let's plot a histogram with the fitted parameters used by the function\nsns.distplot(target , fit=norm);\n(mu, sigma) = norm.fit(target)\nplt.legend(['Normal dist. ($\\mu=$ {:.2f} and $\\sigma=$ {:.2f} )'.format(mu, sigma)],\n            loc='best')\nplt.title('Price (Log)');\nprint(\"Skewness: %f\" % target.skew())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_sampled = df.sample(n = N, random_state = 0)\n\ntype_1JHC = df_sampled[df_sampled[\"type\"] == \"1JHC\"]\ntype_1JHN = df_sampled[df_sampled[\"type\"] == \"1JHN\"]\ntype_2JHC = df_sampled[df_sampled[\"type\"] == \"2JHC\"]\ntype_2JHH = df_sampled[df_sampled[\"type\"] == \"2JHH\"]\ntype_2JHN = df_sampled[df_sampled[\"type\"] == \"2JHN\"]\ntype_3JHC = df_sampled[df_sampled[\"type\"] == \"3JHC\"]\ntype_3JHH = df_sampled[df_sampled[\"type\"] == \"3JHH\"]\ntype_3JHN = df_sampled[df_sampled[\"type\"] == \"3JHN\"]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.set(rc={'figure.figsize':(17.7,8.27)})\n\nax = sns.kdeplot(type_1JHC[\"scalar_coupling_constant\"], label = \"type 1JHC\")\nsns.kdeplot(type_1JHN[\"scalar_coupling_constant\"], label = \"type 1JHN\")\n\nsns.kdeplot(type_2JHC[\"scalar_coupling_constant\"], label = \"type 2JHC\")\nsns.kdeplot(type_2JHH[\"scalar_coupling_constant\"], label = \"type 2JHH\")\nsns.kdeplot(type_2JHN[\"scalar_coupling_constant\"], label = \"type 2JHN\")\n\nsns.kdeplot(type_3JHC[\"scalar_coupling_constant\"], label = \"type 3JHC\")\nsns.kdeplot(type_3JHH[\"scalar_coupling_constant\"], label = \"type 3JHH\")\nsns.kdeplot(type_3JHN[\"scalar_coupling_constant\"], label = \"type 3JHN\")\n\nax.set(xlabel='Scalar coupling constant', ylabel='Density', title='Distribution of the scalar coupling constant by type')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We see that the type 1JHC molecules have in general higher values of the scalar coupling constant. The type 1JHN is also quite different from the rest. The molecules of type 2JHH are the ones with the lower coupling constant values.\n\nWe can make the same visualization using a ridge plot with the seaborn library (check this [link](https://seaborn.pydata.org/examples/kde_ridgeplot.html))."},{"metadata":{"trusted":true},"cell_type":"code","source":"df_sampled = df.sample(n = N, random_state = 0)\n\nsns.set(style=\"white\", rc={\"axes.facecolor\": (0, 0, 0, 0)})\n\n# Initialize the FacetGrid object\npal = sns.cubehelix_palette(10, rot=-.25, light=.7)\ng = sns.FacetGrid(df_sampled, row =\"type\", hue=\"type\", aspect=15, height=.5, palette=pal)\n# Draw the densities in a few steps\ng.map(sns.kdeplot, \"scalar_coupling_constant\", clip_on=False, shade=True, alpha=1, lw=1.5, bw=.2)\ng.map(sns.kdeplot, \"scalar_coupling_constant\", clip_on=False, color=\"w\", lw=1.6, bw=.2)\ng.map(plt.axhline, y=0, lw=1, clip_on=False)\n\n# Define and use a simple function to label the plot in axes coordinates\ndef label(x, color, label):\n    ax = plt.gca()\n    ax.text(0, .2, label, fontweight=\"bold\", color=color,\n            ha=\"left\", va=\"center\", transform=ax.transAxes)\n\n\ng.map(label, \"scalar_coupling_constant\")\n\n# Set the subplots to overlap\ng.fig.subplots_adjust(hspace=-.25)\n\n# Remove axes details that don't play well with overlap\ng.set_titles(\"\")\ng.set(yticks=[])\ng.despine(bottom=True, left=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_sampled = df.sample(n = N, random_state = 0)\nsns.relplot(x=\"scalar_coupling_constant\", y=\"potential_energy\", hue=\"type\", data=df_sampled);","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_sampled = df.sample(n = N, random_state = 0)\nsns.catplot(x=\"type\", y=\"scalar_coupling_constant\", kind=\"swarm\", data=df_sampled);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"----\n<a class=\"anchor\" id=\"fifth-bullet\"></a>\n## 5. Clustering\n\nLet's see if we can find some groups of molecules with the same characteristics."},{"metadata":{},"cell_type":"markdown","source":"### 5.1 Number of clusters\n\nTo select the most adecuate number of clusters to choose, we will make use of the elbow rule."},{"metadata":{"trusted":true},"cell_type":"code","source":"#clusters with the numerical variables considered\n\n\n\ndf_numeric = df[[\"atom_index_0\", \"atom_index_1\", \"X\", \"Y\", \"Z\", \"potential_energy\", \"scalar_coupling_constant\"]]\n\nsse = {}\nfor k in range(1, 10):\n    kmeans = KMeans(n_clusters=k, max_iter=1000, random_state=0).fit(df_numeric)\n    df_numeric[\"clusters\"] = kmeans.labels_\n    sse[k] = kmeans.inertia_ # Inertia: Sum of distances of samples to their closest cluster center\nplt.figure()\nplt.plot(list(sse.keys()), list(sse.values()))\nplt.title(\"Elbow Rule\")\nplt.xlabel(\"Number of clusters\")\nplt.ylabel(\"SSE\")\nplt.show()\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We will select k = 3 clusters."},{"metadata":{},"cell_type":"markdown","source":"### 5.2 Run the algorithm"},{"metadata":{"trusted":true},"cell_type":"code","source":"df_numeric = df[[\"atom_index_0\", \"atom_index_1\", \"X\", \"Y\", \"Z\", \"potential_energy\", \"scalar_coupling_constant\"]]\nk = 3\n\n# --- subselect data to speed up computations --- if needed just delete the following lines -------\nN = 2000\ndf_numeric = df_numeric.sample(n = N, random_state = 0)\n## ------------------------------------------------------------------------------------------------\n\n\n# build the clustering model\nkmeans = KMeans(n_clusters=k)\nkmeans.fit(df_numeric)\n\n# create the clusters\nclusters = kmeans.predict(df_numeric)\n\ndf_numeric[\"clusters\"] = clusters\ndf_numeric.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### 5.3 Visualization of the clusters\n\nWe will perform t-SNE on the data to be able to visualize the clusters in two dimensions."},{"metadata":{"trusted":true},"cell_type":"code","source":"# apply t-SNE to the data to visualize it in 2-Dimensions\n\n# standardize the data -------\nfrom sklearn.preprocessing import StandardScaler\nscaler = StandardScaler()\n\nscaler.fit(df_numeric) # fit the scaler to the data\nscaled_data = scaler.transform(df_numeric) # scale the data\n\n\n## t-SNE -------\nfrom sklearn.manifold import TSNE\n# call t-SNE object\ntsne = TSNE(random_state=0)\n# fit it to the data: scaled data previously obtained: for t-SNE data should be also scaled\nTSNE_data = tsne.fit_transform(scaled_data)\n\n## plot ------\nplt.figure(figsize=(7,7))\n# loop through all the digits and add the points to the graph\nfor i in range(0,k):\n    plt.scatter(TSNE_data[clusters == i,0], TSNE_data[clusters == i,1], alpha = 1, label = \"Cluster {}\".format(i))\nplt.xlabel(\"First component\")\nplt.ylabel(\"Second Component\")\nplt.title(\"T-SNE\")\nplt.legend(loc='best')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Model\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"df_model = df.drop('molecule_name', axis = 1)\n# create validation set\nfrom sklearn.model_selection import train_test_split\n\nX_train = df_model.drop('scalar_coupling_constant', 1)\ny_train = df_model[\"scalar_coupling_constant\"]\nX_train, X_val, y_train, y_val = train_test_split(X_train, y_train, test_size = 0.3, random_state=0)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Label Encoding\nlbl = LabelEncoder()\nlbl.fit(list(train['type'].values) + list(train['type'].values))\ntrain['type'] = lbl.transform(list(train['type'].values))\ntest['type'] = lbl.transform(list(test['type'].values))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"molecules = train.pop('molecule_name')\ntest = test.drop('molecule_name', axis=1)\n\ny = train.pop('scalar_coupling_constant')\n\nyoof = np.zeros(len(train))\nyhat = np.zeros(len(test))\n\nn_splits = 3\ngkf = GroupKFold(n_splits=n_splits) # we're going to split folds by molecules\n\nfold = 0\nfor in_index, oof_index in gkf.split(train, y, groups=molecules):\n    fold += 1\n    print(f'fold {fold} of {n_splits}')\n    X_in, X_oof = train.values[in_index], train.values[oof_index]\n    y_in, y_oof = y.values[in_index], y.values[oof_index]\n    reg = RandomForestRegressor(n_estimators=250,\n                                max_depth=9,\n                                min_samples_leaf=3,\n                                n_jobs=-1)\n    reg.fit(X_in, y_in)\n    yoof[oof_index] = reg.predict(X_oof)\n    yhat += reg.predict(test)\n\nyhat /= n_splits","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sample_submission = pd.read_csv('../input/sample_submission.csv', index_col='id')\n\nbenchmark = sample_submission.copy()\nbenchmark['scalar_coupling_constant'] = yhat\nbenchmark.to_csv('atomic_distance_benchmark.csv')","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.4","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}