{"cells":[{"metadata":{},"cell_type":"markdown","source":"**Hello World.**\n\nI published a similar kernel a year ago when we were asked to predict formation and bandgap energy in the [Nomad2018 Predicting Transparent Conductors](https://www.kaggle.com/c/nomad2018-predict-transparent-conductors/data) challenge. As in the past competition, I hope that the subsequent code might be a useful resource to gain some interesting insights.\n\n**What is this kernel about?**\n\nThe [PCA](https://en.wikipedia.org/wiki/Principal_component_analysis) algorithm is leveraged to project the 3D points that are contained in the *xyz*-files onto the 2D plane. While in the 3D space some patterns might remain unrevealed, the projection along their greatest variance might provide a more accessible view to the point configurations. Moreover, I included connecting lines between points that represent the *scalar coupling constant*."},{"metadata":{},"cell_type":"markdown","source":"Loading the libraries:"},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# General libraries\nimport pandas as pd\nimport numpy as np\n\n# Plotting and Visualization Library\n%matplotlib inline\nfrom mpl_toolkits.mplot3d import Axes3D\nimport matplotlib.pyplot as plt\n\n# Math Libraries\nfrom sklearn.decomposition import PCA\nfrom scipy.spatial import ConvexHull","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Loading the train set and saving the unique atom configurations to an array:"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# Load the main train and test data\ntrain = pd.read_csv(\"../input/train.csv\")\nmolecule_name = train.molecule_name.unique()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Function to load the xyz-files:"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def get_xyz_data(filepath):\n    \"\"\"A function to load the 3D points contained in the xyz-Files to a pandas DataFrame.\n\n    Args:\n        filename (string): Path to the xyz-file\n\n    Returns:\n        pandas DataFrame: A table in which the first three columns correspond to the point's\n                          xyz-coordinates and the fourth column the atom abbreviation\n\n    Note:\n        This function is adapted from Tony Y: https://www.kaggle.com/tonyyy\n\n    \"\"\"\n    A = pd.DataFrame(columns=list('ABCD'))\n    with open(filepath) as f:\n        k = 0\n        for line in f.readlines():\n            x = line.split()\n            if k > 1:\n                newrowA = pd.DataFrame([[x[1],x[2],x[3],x[0]]], columns=list('ABCD'))\n                A = A.append(newrowA)\n            k = k + 1\n    return A","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Function to project and to draw the points:"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def plot_pca(index, showConvexHull=False):\n    \"\"\"A function that projects the 3D points onto its principal components using PCA algorithm.\n       Further, the points (in 3D and 2D) are drawn while the target scalar_coupling_constant is\n       represented as the thickness of the linking line between the points.\n       \n    Args:\n        index (int): An index that picks from the molecule_name array a molecule name\n        showConvexHull(bool): An optional value. If set to true the 2D plots will include the convex hull\n                              around the point configuration.\n    \"\"\"\n    \n    fn = \"../input/structures/{}.xyz\".format(molecule_name[index])\n    train_xyz = get_xyz_data(fn)\n    temp = train[train.molecule_name == molecule_name[index]]\n    color_dict = { 'C':'black', 'H':'blue', 'O':'red', 'N':'green' }\n    \n    minimalValue = np.min(temp.scalar_coupling_constant)\n    maximalValue = np.max(temp.scalar_coupling_constant)\n    # TODO: One can probably find a much better scaling for the thickness of the lines\n    thickness = 3*temp.scalar_coupling_constant.values / maximalValue\n    \n    matrix = train_xyz\n    colour = matrix[\"D\"]\n    matrix = matrix[[\"A\",\"B\",\"C\"]].values\n    matrix = matrix.astype(float)\n    \n    pca = PCA(n_components=3)\n    X_r = pca.fit(matrix).transform(matrix)\n    \n    df_ = pd.DataFrame(np.round(X_r,2))\n        \n    x = np.array(matrix[:,0])\n    y = np.array(matrix[:,1])\n    z = np.array(matrix[:,2])\n\n    fig = plt.figure(figsize=(20,20))\n    \n    ax = fig.add_subplot(221, projection='3d')\n    ax.scatter(x, y, z, c=[color_dict[i] for i in colour], marker='o', s=70)\n    ax.set_xlabel('X Coordinate')\n    ax.set_ylabel('Y Coordinate')\n    ax.set_zlabel('Z Coordinate')\n    ax.set_title(molecule_name[index])\n    \n    ax = fig.add_subplot(222)\n    plt.scatter(X_r[:, 0], X_r[:, 1], color=[color_dict[i] for i in colour], alpha=.8, lw=1, s=70)\n    if(showConvexHull):\n        hull = ConvexHull(X_r[:,[0,1]])\n        volume_1 = hull.volume\n        plt.plot(X_r[hull.vertices,0], X_r[hull.vertices,1], 'r--', lw=1)\n    for k in range(0,len(temp)):\n        a = temp.iloc[k,2]\n        b = temp.iloc[k,3]\n        plt.plot(X_r[[a, b],0], X_r[[a, b],1], 'b', lw=thickness[k])\n    plt.title(molecule_name[index])\n    ax.set_xlabel('First Principal Component')\n    ax.set_ylabel('Second Principal Component')\n    ax.set_title(molecule_name[index])\n    \n    ax = fig.add_subplot(223)\n    plt.scatter(X_r[:, 0], X_r[:, 2], color=[color_dict[i] for i in colour], alpha=.8, lw=1, s=70)\n    if(showConvexHull):\n        hull = ConvexHull(X_r[:,[0,2]])\n        volume_2 = hull.volume\n        plt.plot(X_r[hull.vertices,0], X_r[hull.vertices,2], 'r--', lw=1)\n    for k in range(0,len(temp)):\n        a = temp.iloc[k,2] #'atom_index_0']\n        b = temp.iloc[k,3] #'atom_index_1']\n        plt.plot(X_r[[a, b],0], X_r[[a, b],2], 'b', lw=thickness[k])\n    plt.title(molecule_name[index])\n    ax.set_xlabel('First Principal Component')\n    ax.set_ylabel('Third Principal Component')\n    ax.set_title(molecule_name[index])\n    \n    ax = fig.add_subplot(224)\n    plt.scatter(X_r[:, 1], X_r[:, 2], color=[color_dict[i] for i in colour], alpha=.8, lw=1, s=70)\n    if(showConvexHull):\n        hull = ConvexHull(X_r[:,[1,2]])\n        volume_3 = hull.volume\n        plt.plot(X_r[hull.vertices,1], X_r[hull.vertices,2], 'r--', lw=1)\n    for k in range(0,len(temp)):\n        a = temp.iloc[k,2]\n        b = temp.iloc[k,3]\n        plt.plot(X_r[[a, b],1], X_r[[a, b],2], 'b', lw=thickness[k])\n    plt.title(molecule_name[index])  \n    ax.set_xlabel('Second Principal Component')\n    ax.set_ylabel('Third Principal Component')\n    ax.set_title(molecule_name[index])\n    \n    plt.show()\n    \n    print(\"On the first PC are approx. \" + str(len(df_[0].unique())) + \" distinct coordinates with atoms\")\n    print(\"On the second PC are approx. \" + str(len(df_[1].unique())) + \" distinct coordinates with atoms\")\n    print(\"On the third PC are approx. \" + str(len(df_[2].unique())) + \" distinct coordinates with atoms\")\n    \n    if(showConvexHull):\n        print(\"\")\n        print(\"Area covered by the first and second principal component: \" + str(volume_1))\n        print(\"Area covered by the first and third principal component: \" + str(volume_2))\n        print(\"Area covered by the second and third principal component: \" + str(volume_3))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now Let's execute the code. The color configuration for the atoms is set as follows:\n- C : black\n- H : blue\n- O : red\n- N : green"},{"metadata":{"trusted":true},"cell_type":"code","source":"index = [345,654,1000,1337,2789,10000]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_pca(index=index[0])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_pca(index=index[1])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_pca(index=index[2])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_pca(index=index[3])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_pca(index=index[4])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_pca(index=index[5])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"I hope that the visualizations contribute to discovering interesting patterns and that it will help you to proceed in the competition. Good luck to everyone!\n\nBest, Max"}],"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}