{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Summary\n\n[Raphaël CC's kernel](https://www.kaggle.com/daijin12/coulomb-interaction-high-perf-no-loop-almost) can finish the computation of 5 nearest atoms' Coulomb interaction in just under 8 hours for all molecules and for `H`, `C`, `O`, `N`, `F` atoms. \n\nHowever, if you wanna compute 8 or 9 nearest atoms, this kernel will time out since there are more than 10 `H` and `C` atoms in some molecules. The biggest reason is that `pandas` uses only 1 CPU core when manipulating dataframes, which is not very economical considering there are 4 cores at our disposal. So in this a kernel parallel version of the above computation is implemented, where even more complicated computation like Yukawa potentials for all molecules in under 4 hours.\n\nTo use the output file, simple concat or merge it with the original structures, then merge it with the train and test dataframes, and you are good to go.\n\n### Reference:\n\n* [coulomb_interaction - speed up!](https://www.kaggle.com/rio114/coulomb-interaction/notebook)\n* [coulomb_interaction - Parallelized](https://www.kaggle.com/brandenkmurray/coulomb-interaction-parallelized/notebook)\n* [Coulomb interaction - High perf, no loop (almost)](https://www.kaggle.com/daijin12/coulomb-interaction-high-perf-no-loop-almost)"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np\nfrom tqdm import tqdm_notebook as tqdm\nimport pandas as pd \nimport multiprocessing as mp\nimport warnings\nwarnings.filterwarnings(\"ignore\")","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"structures = pd.read_csv('../input/structures.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"nuclear_charge = {'H':1.0, 'C':6, 'N':7, 'O':8, 'F':9}\nstructures['nuclear_charge'] = [nuclear_charge[x] for x in structures['atom'].values]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Yukawa potential\n\nI am not a chemist, but from [the Wikipedia entry](https://en.wikipedia.org/wiki/Yukawa_potential):\n$$\n{\\displaystyle V_{\\text{Yukawa}}(r)=-\\frac{q_{1}q_{2}}{4\\pi \\epsilon_{0}}{\\frac {e^{-\\alpha mr}}{r}},}\n$$\n$r$ is the radial distance to the particle (atom), and $\\alpha$ is another scaling constant, so that ${\\displaystyle 1/\\alpha m}$ is the range. So here I make a simplification that let the range be the 0.5 times the maximum distance among all other atoms in a molecule to the atom of interest."},{"metadata":{"trusted":true},"cell_type":"code","source":"def compute_all_yukawa(x):   \n    #Apply compute_all_dist2 to each atom \n    return x.apply(compute_yukawa_matrix,axis=1,x2=x)\n\ndef compute_yukawa_matrix(x,x2):\n    # atoms in the molecule which are not the processed one\n    notatom = x2[(x2.atom_index != x[\"atom_index\"])].reset_index(drop=True) \n    # processed atom\n    atom = x[[\"x\",\"y\",\"z\"]]\n    charge = x[['nuclear_charge']]\n    \n    # compute distance from to processed atom to each other\n    notatom['dist'] = ((notatom[[\"x\",\"y\",\"z\"]].values - atom.values)**2).sum(axis=1)\n    notatom['dist'] = np.sqrt(notatom['dist'].astype(np.float32))\n    notatom['dist'] = charge.values*notatom[['nuclear_charge']].values.reshape(-1)\\\n                    *np.exp(-2*notatom['dist']/notatom['dist'].max())/notatom['dist']\n\n    # sort atom per the smallest distance (highest 1/r**2) per group of C/H/N... \n    s = notatom.groupby(\"atom\")[\"dist\"].transform(lambda x : x.sort_values(ascending=False))\n    \n    # keep only the five nearest atoms per group of C/H/N...\n    index0, index1=[],[]\n    for i in notatom.atom.unique():\n        for j in range(notatom[notatom.atom == i].shape[0]):\n            if j < 5:\n                index1.append(\"dist_\" + i + \"_\" + str(j))\n            index0.append(j)\n    s.index = index0\n    s = s[s.index < 5]\n    s.index = index1\n    \n    return s","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Benchmark using first 100 molecules"},{"metadata":{"trusted":true},"cell_type":"code","source":"small_idx = structures.molecule_name.isin(structures.molecule_name.unique()[:100])\n_smallstruct = structures[small_idx]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Using the current [fastest way](https://www.kaggle.com/daijin12/coulomb-interaction-high-perf-no-loop-almost) to compute takes about 20 seconds."},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nsmallstruct1 = _smallstruct.groupby(\"molecule_name\").apply(compute_all_yukawa)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"smallstruct1.head(10)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Multiprocessing\n\nThe following approach makes use of the `groupby` to get an iterator, so that we can use `multiprocessing` which saves more than 50% of the computation time.\n\nReference: [Parallel operations over a Pandas DF](https://www.kaggle.com/gvyshnya/parallel-operations-over-a-pandas-df)"},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nchunk_iter = _smallstruct.groupby(['molecule_name'])\npool = mp.Pool(4) # use 4 processes\n\nfunclist = []\nfor df in tqdm(chunk_iter):\n    # process each data frame\n    f = pool.apply_async(compute_all_yukawa,[df[1]])\n    funclist.append(f)\n\nresult = []\nfor f in tqdm(funclist):\n    result.append(f.get(timeout=120)) # timeout in 120 seconds = 2 mins\n\n# combine chunks with transformed data into a single structure file\nsmallstruct2 = pd.concat(result)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"smallstruct2.head(10)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Just to make sure we are getting the same thing by two methods."},{"metadata":{"trusted":true},"cell_type":"code","source":"np.allclose(smallstruct2.fillna(0), smallstruct1.fillna(0))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Compute Yukawa interaction for all molecules\n\nWithout the parallelization, it takes about 11 hours to run."},{"metadata":{"trusted":true},"cell_type":"code","source":"chunk_iter = structures.groupby(['molecule_name'])\npool = mp.Pool(4) # use 4 CPU cores\n\nfunclist = []\nfor df in tqdm(chunk_iter):\n    # process each data frame\n    f = pool.apply_async(compute_all_yukawa,[df[1]])\n    funclist.append(f)\n\nresult = []\nfor f in tqdm(funclist):\n    result.append(f.get()) \n\n# combine chunks with transformed data into a single training set\nstructures_yukawa = pd.concat(result)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"structures_yukawa.to_csv('structures_yukawa.csv',index=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## How to use this file in CHAMPS competition\n\nSimply concatenate it with the existing `structures` dataframe:\n```python\nstructures = pd.concat([structures, structures_yukawa], axis=1)\n```\nThen merge it with train or test dataframe using a modified version of Andrew's routine\n\n```python\ndef map_atom_info(df_1, df_2, atom_idx):\n    df = pd.merge(df_1, df_2, 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    return df\n\nfor atom_idx in [0,1]:\n    train = map_atom_info(train, structures, atom_idx)\n    train = train.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```"}],"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}