creating-bond-centered-2d-images.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "_kg_hide-input": true
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import pickle\n",
    "import time\n",
    "\n",
    "structures = pd.read_csv('../Data/structures.csv')\n",
    "train = pd.read_csv('../Data/test.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "_kg_hide-input": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "selecting 24195 bonds from 1JHN\n"
     ]
    }
   ],
   "source": [
    "symetry = '1JHN'\n",
    "file_prefix = '../Data/Images/'+symetry+'/'\n",
    "train = train[train['type'] == symetry]\n",
    "import os\n",
    "if not os.path.exists(file_prefix):\n",
    "    os.makedirs(file_prefix)\n",
    "print(f'selecting {train.shape[0]} bonds from {symetry}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "_kg_hide-input": true
   },
   "outputs": [],
   "source": [
    "def map_atom_info(df, atom_idx):\n",
    "    df = pd.merge(df, structures, how = 'left',\n",
    "                  left_on  = ['molecule_name', f'atom_index_{atom_idx}'],\n",
    "                  right_on = ['molecule_name',  'atom_index'])\n",
    "    \n",
    "    df = df.drop('atom_index', axis=1)\n",
    "    df = df.rename(columns={'atom': f'atom_{atom_idx}',\n",
    "                            'x': f'x_{atom_idx}',\n",
    "                            'y': f'y_{atom_idx}',\n",
    "                            'z': f'z_{atom_idx}'})\n",
    "    return df\n",
    "\n",
    "train = map_atom_info(train, 0)\n",
    "train = map_atom_info(train, 1)\n",
    "\n",
    "train_p_0 = train[['x_0', 'y_0', 'z_0']].values\n",
    "train_p_1 = train[['x_1', 'y_1', 'z_1']].values\n",
    "d_bond = train_p_0 - train_p_1\n",
    "nbonds = len(d_bond)\n",
    "\n",
    "train_p_0 = train[['x_0', 'y_0', 'z_0']].values\n",
    "train_p_1 = train[['x_1', 'y_1', 'z_1']].values\n",
    "center = (train_p_0 + train_p_1) / 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "we define a function that performs the integration along the circles:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def density_integration(df, nb_angles):\n",
    "    density_integrale = []\n",
    "    for icouple in range(df.shape[0]): \n",
    "\n",
    "        step_length = length[icouple]\n",
    "        p = df[icouple, :, :].copy()   \n",
    "        density = np.zeros(nb_angles)\n",
    "\n",
    "        for iatom in range(len(atoms)):\n",
    "            val = p - atoms[iatom]\n",
    "            dist = np.linalg.norm(val, axis=1)\n",
    "            density += 1 / dist**2\n",
    "\n",
    "        integrale = density.sum() * step_length\n",
    "        density_integrale.append(integrale)\n",
    "    return density_integrale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We define images of size 35 x 35:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "nb_r = 35\n",
    "nb_R = 35\n",
    "nb_angles = 45\n",
    "\n",
    "values = [(r, R) for r in np.linspace(5, 0.01, nb_r) for R in np.linspace(-5, 5, nb_R)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5.314911842346191"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "start_time = time.time()\n",
    "\n",
    "image_file = 'images.pkl'\n",
    "bond_file = 'bond.pkl'\n",
    "\n",
    "#bond_list = []\n",
    "#image_set = []\n",
    "\n",
    "ids = train['id'].values\n",
    "for ibond in range(5): \n",
    "    #__________________\n",
    "    d = d_bond[ibond]\n",
    "    orig = center[ibond]\n",
    "\n",
    "    molec = train.iloc[ibond]['molecule_name']\n",
    "    #bond_list.append(molec+'_'+str(ids[ibond]))\n",
    "\n",
    "    norm = np.linalg.norm(d)\n",
    "    d = d / norm\n",
    "    #______________________\n",
    "    if d[2] != 0:\n",
    "        u = np.ones(3)\n",
    "        u[2] = -(d[0]*u[0] + d[1]*u[1]) / d[2]\n",
    "        norm = np.linalg.norm(u)\n",
    "        u = u / norm\n",
    "    elif d[0] != 0:\n",
    "        u = np.ones(3)\n",
    "        u[0] = - d[1]*u[1] / d[0]\n",
    "        norm = np.linalg.norm(u)\n",
    "        u = u / norm\n",
    "    else:\n",
    "        u = np.ones(3)\n",
    "        u[1] = 0\n",
    "        norm = np.linalg.norm(u)\n",
    "        u = u / norm\n",
    "    #____________________________________\n",
    "    v = np.cross(d, u)\n",
    "    norm = np.linalg.norm(v)\n",
    "    v = v / norm\n",
    "    #_______________________________\n",
    "    new_data = []\n",
    "    length = []\n",
    "    for r, R in values:\n",
    "        center_circle = R * d  + orig \n",
    "        #______________________________________________________________\n",
    "        val = []\n",
    "        for theta in np.linspace(0, 2 * np.pi, nb_angles):\n",
    "            xp = r * np.cos(theta)\n",
    "            yp = r * np.sin(theta)\n",
    "            val.append([xp, yp])\n",
    "    #         print(theta)\n",
    "        #________________________________\n",
    "        new_pt = []\n",
    "        for pt in val:\n",
    "            new_pt.append(center_circle + (pt[0] * u + pt[1] * v ))\n",
    "        new_data.append(np.array(new_pt))\n",
    "\n",
    "        length.append(r * (2 * np.pi) / (nb_angles - 1))\n",
    "    #__________________________________\n",
    "    new_data = np.array(new_data)\n",
    "    #__________________________________\n",
    "    image = []\n",
    "    for atom in ['C', 'N', 'O', 'H', 'F']:\n",
    "        atoms = structures[(structures['molecule_name'] == molec) &\n",
    "                           (structures['atom'] == atom)][['x', 'y', 'z']].values\n",
    "\n",
    "        temp = np.sum(1/np.sum((new_data[:,:,np.newaxis,:] - atoms)**2,3),(1,2)) * length\n",
    "        image.append(temp.reshape(nb_r, nb_R))\n",
    "\n",
    "    #image_set.append(image)\n",
    "\n",
    "    with open(file_prefix+molec+'_'+str(ids[ibond])+'.pkl', 'wb') as fp:\n",
    "        pickle.dump(image, fp)\n",
    "\n",
    "# with open(bond_file, 'wb') as fp:\n",
    "#     pickle.dump(bond_list, fp)\n",
    "\n",
    "time.time() - start_time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(5, 3)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "atoms.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1225, 45, 3)"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "new_data.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "9.95 ms ± 77.1 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "temp = np.sum(1/np.sum((new_data[:,:,np.newaxis,:] - atoms)**2,3),(1,2)) * length"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([3.88533578, 4.12935821, 4.38592007, ..., 0.01181474, 0.01056351,\n",
       "       0.00950516])"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "temp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "92 ms ± 4.23 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "data = density_integration(new_data, nb_angles)\n",
    "temp2 = np.array(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}