{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Exploring Additional Data in Predicting Molecular Properties\n\nIn this competition we are asked to **predict magnetic interactions between a pair of atoms**. As we all know, literally everything is made of atoms. So if we could get deeper insight into the structure and dynamics of molecules, that could advance many fields of science including environmental science, pharmaceutical science, and materials science. \n\nI like this kind of science competition, as I can defintely gain a reward. Winning a kaggle competition is wonderful, but getting to know a new field of science is surely another reward regardless of your final result may be.\n\nHere I explore **the additional data, which are consisted of 5 csv files (dipole_moments.csv, magnetic_shielding_tensors.csv, mulliken_charges.csv, potential_energy.csv, scalar_coupling_contributions.csv)**. They are provided only for training data, so these infomration may not be super useful for our predictions. But those data are given. Why don't we explore what they are? At least we can learn new things and potentially get better understandings about training & test data! \n\nI am not a chemist or anything, so your feedbacks are very welcome:D\n\n## Files"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load in \n\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\n\nimport os\nprint(os.listdir(\"../input\"))\n\n# Any results you write to the current directory are saved as output.","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Libraries"},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"# libraries\nimport seaborn as sns\nimport matplotlib.pyplot as plt\nfrom mpl_toolkits.mplot3d import Axes3D\n\nfrom scipy import stats\nimport gc\nimport warnings\nwarnings.filterwarnings(\"ignore\")\n\nprint(\"Libraries were loaded.\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Loading data"},{"metadata":{"trusted":true},"cell_type":"code","source":"# load data\ntrain_df = pd.read_csv('../input/train.csv')\npotential_energy_df = pd.read_csv('../input/potential_energy.csv')\nmulliken_charges_df = pd.read_csv('../input/mulliken_charges.csv')\nscalar_coupling_contributions_df = pd.read_csv('../input/scalar_coupling_contributions.csv')\nmagnetic_shielding_tensors_df = pd.read_csv('../input/magnetic_shielding_tensors.csv')\ndipole_moments_df = pd.read_csv('../input/dipole_moments.csv')\nstructure_df = pd.read_csv('../input/structures.csv')\ntest_df = pd.read_csv('../input/test.csv')\n\n\nprint(\"All the data were loaded.\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Basic information about each file\nWhat does each file look like?"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# What are inside those files?\ndfs = [train_df, potential_energy_df, mulliken_charges_df, \n       scalar_coupling_contributions_df, magnetic_shielding_tensors_df, \n       dipole_moments_df, structure_df, test_df]\nnames = [\"train_df\", \"potential_energy_df\", \"mulliken_charges_df\", \n       \"scalar_coupling_contributions_df\", \"magnetic_shielding_tensors_df\", \n       \"dipole_moments_df\", \"structure_df\", \"test_df\"]\n\n# display info about a DataFrame\ndef dispDF(df, name):\n    print(\"========== \" + name + \" ==========\")\n    print(\"SHAPE ----------------------\")\n    print(df.shape)\n    print('')\n    print(\"HEAD ----------------------\")\n    print(df.head(5))\n    print('')\n    print(\"DATA TYPE ----------------\")\n    print(df.dtypes)\n    print('')\n    print(\"UNIQUES -------------------\")\n    print(df.nunique())\n    print('')\n    print(\"======================================\")\n\npd.set_option('display.expand_frame_repr', False)\nfor df, name in zip(dfs, names):\n    dispDF(df, name)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now we look at each file one by one:) I try to visualize the content of each file."},{"metadata":{"trusted":true,"_kg_hide-input":false,"_kg_hide-output":false},"cell_type":"code","source":"# colors\ncolors = sns.color_palette(\"cubehelix\", 8)\nsns.set()\n\nsubsample = 100","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Dipole moments\n\n### What are the dipole moments?\nDipole moments quantify **a separation of charge between atoms**. \n\nWhen a pair of different atoms forms an ionic bond or a covalent bond, an atom with the stronger **electronegativity (F > O > Cl > N > Br > I > S > C > H > metals)** pulls shared electrons away from the other. For example, in CO2, O has stronger electronegativity than C. That means that O pulls electrons from C, generating the dipole moments between a pair of the two atoms. \n\nHowever, CO2 as a molecule itself does not have dipole moments because in its structure (O = C = O), Os pull electrons away from C equally but in the opposite direction, which leads to overall zero separation of charge between atoms.\n\nH2O is another story. Unlike CO2, H2O has the asymmetric structure. Therefore, although O pulls electrons from Hs, a separation of charge still remains.\n\n![](https://i.ibb.co/4YJbS0Z/dipole-Moments.png)\n\nReference:\n- CHEMISTRY LibreTexts: Dipole Moments\nhttps://chem.libretexts.org/Bookshelves/Physical_and_Theoretical_Chemistry_Textbook_Maps/Supplemental_Modules_(Physical_and_Theoretical_Chemistry)/Physical_Properties_of_Matter/Atomic_and_Molecular_Properties/Dipole_Moments\n\n- CHEMISTRY LibreTexts: Electronegativity\nhttps://chem.libretexts.org/Bookshelves/Physical_and_Theoretical_Chemistry_Textbook_Maps/Supplemental_Modules_(Physical_and_Theoretical_Chemistry)/Physical_Properties_of_Matter/Atomic_and_Molecular_Properties/Electronegativity\n\nLet's see how contents of the file look like.\n\n### File description\n> dipole_moments.csv - contains the molecular electric dipole moments. These 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"},{"metadata":{"trusted":true},"cell_type":"code","source":"# display info about the data frame\ndispDF(dipole_moments_df, \"dipole moments\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's look at the charge distribution. Note that I sampled data to reduce the computational burden."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"fig = plt.figure()\nax = Axes3D(fig)\n# ax = fig.add_subplot(111, projection='3d')\n\nscatter_colors = sns.color_palette(\"husl\", 85003)\n\n# 3D scatter\nax.scatter(dipole_moments_df['X'][::subsample], dipole_moments_df['Y'][::subsample],\n           dipole_moments_df['Z'][::subsample], s=30, alpha=0.5, c=scatter_colors[::subsample])\n\nax.set_xlabel('X')\nax.set_ylabel('Y')\nax.set_zlabel('Z')\nax.set_title('Dipole Moment')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Each molecule has a different color. There seems to be no cluster or anything, but we can find some outliers.\n\nHere is a histogram of the distances of molecules in this space."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"dipole_moments_df['dm_distance'] = np.asarray([x**2 + y**2 + z**2 for x, y, z in zip(dipole_moments_df['X'],dipole_moments_df['Y'], dipole_moments_df['Z'])])\n\nfig, ax = plt.subplots(1, 2, figsize=(12, 6))\nax = ax.flatten()\n\n# original distribution\nsns.distplot(dipole_moments_df['dm_distance'], color=colors[0], kde=False, norm_hist=False, ax=ax[0])\nax[0].set_xlabel('distance')\nax[0].set_title('dipole moment')\n\n# in log\nsns.distplot(np.log(dipole_moments_df['dm_distance'] + 0.00001), color=colors[0], kde=False, norm_hist=False, ax=ax[1])\nax[1].set_xlabel('log distance')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As you see, the original histogram (left) is highly right-scewed. Hence it becomes approximately normal after the log-transform (right).\n\nMaybe it is worthwhile making a list of \"outlier\" molecules in terms of the distance in the dipole-moment-space. Here I set the threshold to be 100. This list may be useful for our predictions."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"dipole_moments_df['dm_outliers'] = np.zeros(dipole_moments_df.shape[0]).astype(int)\ndipole_moments_df.loc[dipole_moments_df['dm_distance'] > 100, 'dm_outliers'] = int(1)\nprint(\"outliers (dipole moments): \" + str(np.sum(dipole_moments_df['dm_outliers'] == 1)) + \" molecules\")\ndipole_moments_df.head(7)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Magnetic shielding tensors\n### What are the magnetic shielding tensors?\nMagnetic shielding tensors describe the magnetic field surrounding the nucleus.\n\nThe **NMR (Nuclear Magnetic Resonance)** puts a strong magnetic field upon the molecule. Since the electrons orbit the nucleus (which by the way consists of neutrons and protons), they slightly influence with so-called shielding, or how the external magnetic alters the energy levels of charged particles. \n\nAs the dynamics of theses electrons are influenced by other atoms due to electronegativity, the shielding is inevitably influenced by the local structure or composition of atoms in the molecule. In other words, we can estimate them by NMR. **The magnetic shielding tensors are resulting 3-by-3 tensor (matrix) to represent the measured magnetic field surrounding the nucleus**.\n\n![](https://i.ibb.co/k00CK8g/PCA-tensors.png)\n\nHere the red solid lines indicate the principal components of the shift tensor (σxx, σyy, σzz). Connecting these principal components results in the black ellipses, which represents the electron field surrounding the nucleus.\n\nReference:\n- CHEMISTRY LibreTexts: Chemical Shift (Shielding)\nhttps://chem.libretexts.org/Bookshelves/Physical_and_Theoretical_Chemistry_Textbook_Maps/Supplemental_Modules_(Physical_and_Theoretical_Chemistry)/Spectroscopy/Magnetic_Resonance_Spectroscopies/Nuclear_Magnetic_Resonance/NMR%3A_Theory/NMR_Interactions/Chemical_Shift_(Shielding)\n\n### Data description\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."},{"metadata":{"trusted":true,"_kg_hide-input":false},"cell_type":"code","source":"# display info about the data frame\ndispDF(magnetic_shielding_tensors_df, \"magnetic shielding tensors\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This file stores the 3-by-3 matrix (tensor). I only use the principal components (XX, YY, and ZZ) to visualize their distributions."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"fig = plt.figure()\nax = fig.add_subplot(111, projection='3d')\n\nscatter_colors = sns.color_palette(\"husl\", 29)\n\n# 3D scatter\nfor i in range(29):\n    xx = magnetic_shielding_tensors_df.loc[magnetic_shielding_tensors_df['atom_index']==i, 'XX']\n    yy = magnetic_shielding_tensors_df.loc[magnetic_shielding_tensors_df['atom_index']==i, 'YY']\n    zz = magnetic_shielding_tensors_df.loc[magnetic_shielding_tensors_df['atom_index']==i, 'ZZ']\n    ax.scatter(xx[::subsample*10], yy[::subsample*10], zz[::subsample*10], s=30, alpha=0.5, c=scatter_colors[i])\n\nax.set_xlabel('X')\nax.set_ylabel('Y')\nax.set_zlabel('Z')\nax.set_title('Magnetic shielding tensors')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"It is hard to see but there are 29 colors representing the 29 unique atom index. Are there ... clusters? I mean, small (left side) and big (right side) ones?"},{"metadata":{},"cell_type":"markdown","source":"# Potential energy\n### What is the potential energy?\nPotential energy is relatively easy to understand in the physics. It is **the energy held by an object relative to others**. For example, in the figure below, you can easily imagine what would happen if the string was released. The elastic potential energy of the bow is then transformed into kinetic energy.\n\n![](https://i.ibb.co/Y39t4xx/2019-06-24-12-35-28.png)\n\nThis example is essentially similar to the case with molecules. Atoms form a certain bond (e.g. ionic) to be a molecule such that the overall energy becomes lower (more stable). In other words, **the bonds between atoms have potential energy, which is converted into kinetic energy once the bonds are broken**.\n\nReference:\n- CHEMISTRY LibreTexts: Potential energy https://chem.libretexts.org/Bookshelves/Physical_and_Theoretical_Chemistry_Textbook_Maps/Supplemental_Modules_(Physical_and_Theoretical_Chemistry)/Thermodynamics/Energies_and_Potentials/Potential_Energy\n- [Wikipedia: Potential energy](https://en.wikipedia.org/wiki/Potential_energy)\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."},{"metadata":{"trusted":true},"cell_type":"code","source":"# display info about the data frame\ndispDF(potential_energy_df, \"potential energy\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# potential energy\nfig, ax = plt.subplots(1, 1, figsize=(8, 6))\n\n# histogram\nsns.distplot(potential_energy_df['potential_energy'], \n             kde=False, color=colors[2], ax=ax)\n\n# median\npe_median = potential_energy_df['potential_energy'].median()\nax.axvline(pe_median, color='r', linestyle='--', lw=4)\nax.text(pe_median + 30, 15000, 'median = ' + str(pe_median), fontsize=12, color='r') ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can see that molecules are bimodally distributed with respect to the potential energy. Let's use the median split to have two clusters."},{"metadata":{"trusted":true},"cell_type":"code","source":"# median split\nhighPE_molecules = potential_energy_df.loc[potential_energy_df['potential_energy'] >= pe_median]\nlowPE_molecules = potential_energy_df.loc[potential_energy_df['potential_energy'] < pe_median]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(str(highPE_molecules.shape[0]) + \" high potential energy molecules:\")\nhighPE_molecules.head(7)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(str(lowPE_molecules.shape[0]) + \" low potential energy molecules:\")\nlowPE_molecules.head(7)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# low (0) and high (0) potential energy\npotential_energy_df['energy_class'] = np.zeros(potential_energy_df.shape[0]).astype(int)\npotential_energy_df.loc[potential_energy_df['potential_energy'] >= pe_median, 'energy_class'] = int(1)\npotential_energy_df.head(7)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Mulliken charges\n### What are the Mulliken charges?\nThe Mulliken charges characterize **the electric charge distribution in a molecule**. \n\nA molecular orbital consists of multiple atomic orbitals, so we can estimate the electric charge distribution by a linear combination of these oribitals. The Mulliken charge is assigned to each atom as the sum of the linear terms over all the orbitals that belong to the atom.\n\n![](https://i.ibb.co/PwvdC2w/mulliken.png)\n\nHere net Mulliken charge distribution on H, C and N atoms in histidines are shown as an example (Alia et al., 2004).\n\nReference\n- [Mulliken population analysis](https://en.wikipedia.org/wiki/Mulliken_population_analysis)\n- [Alia et al (2004) Heteronuclear 2D ( 1 H- 13 C) MAS NMR Resolves the Electronic Structure of Coordinated Histidines in Light-Harvesting Complex II: Assessment of Charge Transfer and Electronic Delocalization Effect](https://www.researchgate.net/publication/8893081_Heteronuclear_2D_1_H-_13_C_MAS_NMR_Resolves_the_Electronic_Structure_of_Coordinated_Histidines_in_Light-Harvesting_Complex_II_Assessment_of_Charge_Transfer_and_Electronic_Delocalization_Effect)\n\n### Data description\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."},{"metadata":{"trusted":true},"cell_type":"code","source":"# display info about the data frame\ndispDF(mulliken_charges_df, \"mulliken charges\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# distribution of mulliken_charge\nmulliken_charges_df[\"mulliken_charge\"].hist()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Mulliken charges seems to be slightly above 0 in many cases."},{"metadata":{},"cell_type":"markdown","source":"# Structure\nThis structure.csv is served for us to visualize molecules in 3D (e.g. [How To: Easy Visualization of Molecules.](https://www.kaggle.com/borisdee/how-to-easy-visualization-of-molecules)). I wonder if this file can be useful beyound the visualization.\n\n### Data description\n> structure.zip - folder containing molecular structure (xyz) files, where the first line is the number of atoms in the molecule, followed by a blank line, and then a line for every atom, where the first column contains the atomic element (H for hydrogen, C for carbon etc.) and the remaining columns contain the X, Y and Z cartesian coordinates (a standard format for chemists and molecular visualization programs)."},{"metadata":{"trusted":true},"cell_type":"code","source":"# display info about the data frame\ndispDF(structure_df, \"structure\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# # Let's visualize one molecule anyway\n# !pip install ase","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# import ase\n# from ase import Atoms\n# import ase.visualize\n\n# positions = structure_df.loc[structure_df['molecule_name'] == 'dsgdb9nsd_000001', ['x', 'y', 'z']]\n# symbols = structure_df.loc[structure_df['molecule_name'] == 'dsgdb9nsd_000001', 'atom']\n# ase.visualize.view(Atoms(positions=positions, symbols=symbols), viewer=\"x3d\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"~~This molecule (dsgdb9nsd_000001) is apparently CH4 (Methane). Looks pretty.~~\n\nSince structure.csv contains atom name in a column, it may be a good idea to put another column about electronegativity. Do you not remember electronegativity of each atom? Me neither. This is exactly what internet is for.\n\n![electronegativity](https://i.ibb.co/wWnW8Vs/electronegativity.png)\n\ntaken from [What would cause an atom to have a low electronegativity value? ](https://socratic.org/questions/what-would-cause-an-atom-to-have-a-low-electronegativity-value)\n\nAccording to this table, electronegativities are 2.20 for H, 2.55 for C, 3.04 for N, 3.44 for O, and 3.98 for F."},{"metadata":{"trusted":true},"cell_type":"code","source":"# add electronegativity to df\nstructure_df['electronegativity'] = structure_df['atom']\nstructure_df['electronegativity'] = structure_df['electronegativity'].map({'H': 2.20, 'C': 2.55, 'N': 3.04, 'O': 3.44, 'F': 3.98})\nstructure_df.head(12)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Scaler coupling contributions\n### What is the scaler coupling?\nAs electrons orbit the nucleus (neutrons + protons), the nucleus becomes a magnetic dipole and spins to the random orientations. The scaler coupling (J-coupling in an isotropic liquid, an indirect dipole-dipole coupling) arises from **an indirect interaction between two nuclear spins**. \n\nAs the following data description says, the scaler coupling tells a lot about the molecule such as **the Fermi contact interaction (interaction between an electron and an atomic nucleus)**, **the spin-dipolar interaction**, **the paramagnetic (there are unpaired electrons in an orbital) or the diamagnetic ( electrons are paired or their total spin is 0) spin-orbit interaction (interaction of a particle's spin with its motion inside a potential)**.\n\nReference\n- [Wikipedia: Fermi contact interaction](https://en.wikipedia.org/wiki/Fermi_contact_interaction)\n- [Wikipedia: Spin-orbit interaction](https://en.wikipedia.org/wiki/Spin%E2%80%93orbit_interaction)\n- [Diamagnetism and Paramagnetism](https://courses.lumenlearning.com/introchem/chapter/diamagnetism-and-paramagnetism/)\n\n### Data description\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. "},{"metadata":{"trusted":true},"cell_type":"code","source":"# display info about the data frame\ndispDF(scalar_coupling_contributions_df, \"scalar coupling contributions\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"fig, ax = plt.subplots(2,1, figsize=(12, 8))\nax = ax.flatten()\n\n# proportions of fc, sd, pso, dso in scaler-coupling constant\nfc = scalar_coupling_contributions_df['fc']\nsd = scalar_coupling_contributions_df['sd']\npso = scalar_coupling_contributions_df['pso']\ndso = scalar_coupling_contributions_df['dso']\n\ncolumns = ['fc', 'sd', 'pso', 'dso']\ncontribution_colors = sns.color_palette(\"Set2\", 4)\nnrows = scalar_coupling_contributions_df.shape[0]\n\nfor i, c in enumerate(columns):\n    contributions = 100 * scalar_coupling_contributions_df[c].abs() \\\n                / (np.abs(fc) + np.abs(sd) + np.abs(pso) + np.abs(dso))\n    ax[0].plot(np.arange(0, nrows, subsample*100), contributions[::subsample*100], c=contribution_colors[i], label=c)\n    \nax[0].set_xlabel('molecule')\nax[0].set_ylabel('% to scaler coupling constant')\nax[0].legend()\n\n# unique counts of molecular type\ncounts = np.zeros(scalar_coupling_contributions_df['type'].nunique())\nfor i, u in enumerate(scalar_coupling_contributions_df['type'].unique()):\n    counts[i] = np.sum(scalar_coupling_contributions_df['type'].values == u)\n    \nsns.barplot(x=scalar_coupling_contributions_df['type'].unique(), y=counts, ax=ax[1])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"From the top panel we can know that our target, scaler coupling constant, is in most cases determined by the Fermi Contact, which tells us about an interaction between an electron and an atomic nucleus. \n\nFrom the bottom panel we can see that the type of coupling is not uniformly distributed. We have **pretty much data about 3JHC, 2JHC, and 1JHC** but **limited amount for 1JHN, 2JHN, and 3JHN**. \n\nLet's see if these types of molecules are related to the values of interactions (fc, sd, pso, & dso). The following snipets were suggested by [Dmitrii Borkin](https://www.kaggle.com/hedgehoginfog). Thanks a lot!"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# Values of each interaction vs molecular types  \nfig, ax = plt.subplots(4, 1, figsize=(20, 14)) \nax = ax.flatten()\nfor i, col in enumerate(columns): \n    means = scalar_coupling_contributions_df[[\"type\", col]].groupby(by='type')[col].mean()\n    SDs = scalar_coupling_contributions_df[[\"type\", col]].groupby(by='type')[col].std() \n    means.plot(kind=\"bar\", yerr=SDs, ax=ax[i])\n    ax[i].set_ylabel(col)\n    \nplt.tight_layout()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Interesting, **values of each interaction type varies a lot across molecular types**. For example, the molecular type 2JHH gives rise to predominantly positive values of sd (the spin-dipolar interaction) and pso (the paramagnetic spin-orbit interaction) but negative values of fc (the Fermi contact interaction) and dso (the diamagnetic spin-orbit interaction)."},{"metadata":{},"cell_type":"markdown","source":"# Combine all the files into one\nIt seems possible to combine all the additional files as one. Such a big file may be useful for our feature engineering.\n\nNote that \n\n**\"dipole_moments\" and \"potential_energy\"** files are about **molecule**.\n\n**\"magnetic_shielding_tensor\", \"mulliken_charge\", and \"structure\"** files are about **an atom in a molecule**. Don't forget that only \"structure\" contains data for both train and test. \n\n**\"scalar_coupling_contributions\"** file is about **\"a pair of atoms in a molecule\"**.\n\nSo we got to be careful not to loose infomration when combining these files."},{"metadata":{"trusted":true},"cell_type":"code","source":"# combine \"dipole_moments_df\" and \"potential_energy_df\" (The both have 85003 rows, information per molecule)\nDM_PE_df = pd.merge(dipole_moments_df, potential_energy_df, on='molecule_name')\n\ndel dipole_moments_df, potential_energy_df, train_df, test_df\ngc.collect()\n\nprint(\"There are {} rows and {} columns in DM_PE_df.\".format(DM_PE_df.shape[0], DM_PE_df.shape[1]))\nDM_PE_df.head(12)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# combine \"magnetic_shielding_tensors_df\" and \"mulliken_charges_df\" (The both have 1533537 rows, information per atom in a molecule)\nMST_MC_df = pd.merge(magnetic_shielding_tensors_df, mulliken_charges_df, on=['molecule_name', 'atom_index'])\n\ndel magnetic_shielding_tensors_df, mulliken_charges_df\ngc.collect()\n\nprint(\"There are {} rows and {} columns in DM_PE_df.\".format(MST_MC_df.shape[0], MST_MC_df.shape[1]))\nMST_MC_df.head(12)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# combine these two\nMST_MC_DM_PE_df = pd.merge(MST_MC_df, DM_PE_df, on='molecule_name', how='left')\n\ndel MST_MC_df, DM_PE_df\ngc.collect()\n\nprint(\"There are {} rows and {} columns in DM_PE_df.\".format(MST_MC_DM_PE_df.shape[0], MST_MC_DM_PE_df.shape[1]))\nMST_MC_DM_PE_df.head(12)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### reduce memory burden\nThanks to https://www.kaggle.com/speedwagon/permutation-importance"},{"metadata":{"trusted":true},"cell_type":"code","source":"# lighter structure\ndef reduce_mem_usage(df, verbose=True):\n    numerics = ['int16', 'int32', 'int64', 'float16', 'float32', 'float64']\n    start_mem = df.memory_usage().sum() / 1024**2\n    for col in df.columns:\n        col_type = df[col].dtypes\n        \n        if col_type in numerics:\n            c_min = df[col].min()\n            c_max = df[col].max()\n            if str(col_type)[:3] == 'int':\n                if c_min > np.iinfo(np.int8).min and c_max < np.iinfo(np.int8).max:\n                    df[col] = df[col].astype(np.int8)\n                elif c_min > np.iinfo(np.int16).min and c_max < np.iinfo(np.int16).max:\n                    df[col] = df[col].astype(np.int16)\n                elif c_min > np.iinfo(np.int32).min and c_max < np.iinfo(np.int32).max:\n                    df[col] = df[col].astype(np.int32)\n                elif c_min > np.iinfo(np.int64).min and c_max < np.iinfo(np.int64).max:\n                    df[col] = df[col].astype(np.int64)\n            else:\n                c_prec = df[col].apply(lambda x: np.finfo(x).precision).max()\n                if c_min > np.finfo(np.float16).min and c_max < np.finfo(np.float16).max and c_prec == np.finfo(np.float16).precision:\n                    df[col] = df[col].astype(np.float16)\n                elif c_min > np.finfo(np.float32).min and c_max < np.finfo(np.float32).max and c_prec == np.finfo(np.float32).precision:\n                    df[col] = df[col].astype(np.float32)\n                else:\n                    df[col] = df[col].astype(np.float64)\n    end_mem = df.memory_usage().sum() / 1024**2\n    if verbose: print('Mem. usage decreased to {:5.2f} Mb ({:.1f}% reduction)'.format(end_mem, 100 * (start_mem - end_mem) / start_mem))\n    return df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"MST_MC_DM_PE_df = reduce_mem_usage(MST_MC_DM_PE_df)\nscalar_coupling_contributions_df = reduce_mem_usage(scalar_coupling_contributions_df)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# combine it with \"scaler_coupling_contributions_df\" (information per a pair of atoms in a molecule)\ncombined_df0 = pd.merge(scalar_coupling_contributions_df, MST_MC_DM_PE_df, \n                           left_on=['molecule_name','atom_index_0'], right_on=['molecule_name','atom_index'], how='left')\nprint(\"There are {} rows and {} columns in combined_df0.\".format(combined_df0.shape[0], combined_df0.shape[1]))\ncombined_df0.head(12)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Here the combined file is based on \"atom_index_0\", but of course you can do the same using \"atom_index_1\" like the following. You can even concatenate these two big files together to have a complete information about a pair of molecules. But I could not manage due to the memory limit of this kernel. "},{"metadata":{"trusted":true},"cell_type":"code","source":"# # combine it with \"scaler_coupling_contributions_df\" (information per a pair of atoms in a molecule)\n# combined_df1 = pd.merge(scalar_coupling_contributions_df, MST_MC_DM_PE_df, \n#                            left_on=['molecule_name','atom_index_1'], right_on=['molecule_name','atom_index'], how='left')\n\n# del scalar_coupling_contributions_df, MST_MC_DM_PE_df\n# gc.collect()\n\n# print(\"There are {} rows and {} columns in combined_df0.\".format(combined_df1.shape[0], combined_df1.shape[1]))\n# combined_df1.head(12)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# # combine these two\n# combined_df = pd.merge(combined_df0, combined_df1, on=['molecule_name'])\n\n# del combined_df0, combined_df1\n# gc.collect()\n\n# print(\"There are {} rows and {} columns in combined_df.\".format(combined_df.shape[0], combined_df.shape[1]))\n# combined_df.head(12)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"TO BE UPDATED..."}],"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}