{"cells":[{"metadata":{},"cell_type":"markdown","source":"Quantum Machine 9 (QM9) dataset is now uploaded to Kaggle, [find here](https://www.kaggle.com/zaharch/quantum-machine-9-aka-qm9). \n\nNote that QM9 contains extra information both for the train and the test datasets of the competition. **This kernel extracts the features from the dataset and saves them into file data.covs.pickle for convenience of use**. I think you can even add this kernel to your pipeline of kernels to experiment with the features, or just download the output file (it takes a few hours to create it). The kernel also includes a simple LightGBM model on the extracted features with feature importance graph. Spoiler: mulliken partial charges for the two atoms are on top.\n\nNote: I am not sure how to extract the information from the list of frequencies correctly, currently I have just taken min/max and mean values of the list.\n\nDisclaimer: **the dataset is not allowed to use for your final submissions in this competition**. But we can still learn from it."},{"metadata":{},"cell_type":"markdown","source":"**Does QM9 contain the information from extra files given in the competition?**\n\n1. dipole_moments.csv contains X,Y,Z values per molecule and I found that sqrt(X^2+Y^2+Z^2)=mu where mu is given in QM9\n2. mulliken_charges.csv matches the mulliken charges from QM9\n3. scalar_coupling_contributions.csv I can't find this info in QM9\n4. magnetic_shielding_tensors.csv I can't find this info in QM9"},{"metadata":{},"cell_type":"markdown","source":"\nQM9 contains the structure information, and additionally the following information both for the train and the test:\n\n1. Mulliken partial charge for each atom\n2. Frequencies for degrees of freedom\n3. SMILES from GDB9 and for relaxed geometry\n4. InChI for GDB9 and for relaxed geometry\n\nand also the following 17 properties per molecule:\n\n`\nI. Property  Unit         Description\n 1  tag       -            gdb9; string constant to ease extraction via grep\n 2  index     -            Consecutive, 1-based integer identifier of molecule\n 3  A         GHz          Rotational constant A\n 4  B         GHz          Rotational constant B\n 5  C         GHz          Rotational constant C\n 6  mu        Debye        Dipole moment\n 7  alpha     Bohr^3       Isotropic polarizability\n 8  homo      Hartree      Energy of Highest occupied molecular orbital (HOMO)\n 9  lumo      Hartree      Energy of Lowest occupied molecular orbital (LUMO)\n10  gap       Hartree      Gap, difference between LUMO and HOMO\n11  r2        Bohr^2       Electronic spatial extent\n12  zpve      Hartree      Zero point vibrational energy\n13  U0        Hartree      Internal energy at 0 K\n14  U         Hartree      Internal energy at 298.15 K\n15  H         Hartree      Enthalpy at 298.15 K\n16  G         Hartree      Free energy at 298.15 K\n17  Cv        cal/(mol K)  Heat capacity at 298.15 K\n`"},{"metadata":{},"cell_type":"markdown","source":"Example of QM9 data format (for dsgdb9nsd_000001.xyz):\n\n`\n5\ngdb 1\t157.7118\t157.70997\t157.70699\t0.\t13.21\t-0.3877\t0.1171\t0.5048\t35.3641\t0.044749\t-40.47893\t-40.476062\t-40.475117\t-40.498597\t6.469\t\nC\t-0.0126981359\t 1.0858041578\t 0.0080009958\t-0.535689\nH\t 0.002150416\t-0.0060313176\t 0.0019761204\t 0.133921\nH\t 1.0117308433\t 1.4637511618\t 0.0002765748\t 0.133922\nH\t-0.540815069\t 1.4475266138\t-0.8766437152\t 0.133923\nH\t-0.5238136345\t 1.4379326443\t 0.9063972942\t 0.133923\n1341.307\t1341.3284\t1341.365\t1562.6731\t1562.7453\t3038.3205\t3151.6034\t3151.6788\t3151.7078\nC\tC\t\nInChI=1S/CH4/h1H4\tInChI=1S/CH4/h1H4\n`"},{"metadata":{"trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport scipy.io as sio\nimport os\nfrom pathlib import Path\nimport csv\nimport pickle\nfrom joblib import Parallel, delayed\nfrom tqdm import tqdm_notebook as tqdm\nimport pdb\nimport lightgbm as lgb\nimport xgboost as xgb\nimport random\nfrom sklearn.model_selection import GroupKFold, StratifiedKFold, KFold\nimport seaborn as sn\n\nPATH_QM9 = Path('../input/quantum-machine-9-aka-qm9/dsgdb9nsd.xyz')\nPATH_BASE = Path('../input/champs-scalar-coupling')\nPATH_WORKING = Path('../working')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train = pd.read_csv(PATH_BASE/'train.csv')\ntest = pd.read_csv(PATH_BASE/'test.csv')\n\nboth = pd.concat([train, test], axis=0, sort=False)\nboth = both.set_index('molecule_name',drop=False)\n\nboth.sort_index(inplace=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def processQM9_file(filename):\n    path = PATH_QM9/filename\n    molecule_name = filename[:-4]\n    \n    row_count = sum(1 for row in csv.reader(open(path)))\n    na = row_count-5\n    freqs = pd.read_csv(path,sep=' |\\t',engine='python',skiprows=row_count-3,nrows=1,header=None)\n    sz = freqs.shape[1]\n    is_linear = np.nan\n    if 3*na - 5 == sz:\n        is_linear = False\n    elif 3*na - 6 == sz:\n        is_linear = True\n    \n    stats = pd.read_csv(path,sep=' |\\t',engine='python',skiprows=1,nrows=1,header=None)\n    stats = stats.loc[:,2:]\n    stats.columns = ['rc_A','rc_B','rc_C','mu','alpha','homo','lumo','gap','r2','zpve','U0','U','H','G','Cv']\n    \n    stats['freqs_min'] = freqs.values[0].min()\n    stats['freqs_max'] = freqs.values[0].max()\n    stats['freqs_mean'] = freqs.values[0].mean()\n    stats['linear'] = is_linear\n    \n    mm = pd.read_csv(path,sep='\\t',engine='python', skiprows=2, skipfooter=3, names=range(5))[4]\n    if mm.dtype == 'O':\n        mm = mm.str.replace('*^','e',regex=False).astype(float)\n    stats['mulliken_min'] = mm.min()\n    stats['mulliken_max'] = mm.max()\n    stats['mulliken_mean'] = mm.mean()\n    \n    stats['molecule_name'] = molecule_name\n    \n    data = pd.merge(both.loc[[molecule_name],:].reset_index(drop=True), stats, how='left', on='molecule_name')\n    data['mulliken_atom_0'] = mm[data['atom_index_0'].values].values\n    data['mulliken_atom_1'] = mm[data['atom_index_1'].values].values\n    \n    return data\n\ndef processQM9_list(files):\n    df = pd.DataFrame()\n    for i,filename in enumerate(files):\n        stats = processQM9_file(filename)\n        df = pd.concat([df, stats], axis = 0)\n    return df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"all_files = os.listdir(PATH_BASE/'structures')\n\n%time result = Parallel(n_jobs=4, temp_folder=PATH_WORKING) \\\n    (delayed(processQM9_list)(all_files[100*idx:min(100*(idx+1), len(all_files))]) for idx in tqdm(range(int(np.ceil(len(all_files)/100)))))\n\ndata = pd.concat(result)\ndata = data.reset_index(drop=True)\ndata.to_pickle(PATH_WORKING/'data.covs.pickle')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"data = pd.read_pickle(PATH_WORKING/'data.covs.pickle')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Model"},{"metadata":{"trusted":true},"cell_type":"code","source":"def setSeeds(seed = 1234):\n    random.seed(seed)\n    np.random.seed(seed)\n    os.environ['PYTHONHASHSEED'] = str(seed)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"cat_names = ['type','linear']\nfor i in cat_names:\n    data[i] = pd.factorize(data[i])[0]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"data_train = data.loc[~data['scalar_coupling_constant'].isnull()].reset_index(drop=True)\ndata_test = data.loc[data['scalar_coupling_constant'].isnull()].reset_index(drop=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"params_lgb = {'application': 'regression_l1',\n              'metric': 'l1',\n              'num_leaves': 90,\n              'max_depth': 7,\n              'learning_rate': 1,\n              'bagging_freq' : 1,\n              'bagging_fraction': 0.9,\n              'feature_fraction': 0.9,\n              'min_split_gain': 0.02,\n              'min_child_samples': 50,\n              'min_child_weight': 0.01,\n              'lambda_l2': 0.05,\n              'lambda_l1': 0.01,\n              'verbosity': -1,\n              'data_random_seed': 17}\n\nFOLDS_VALID = 5\nmeta_cols = ['id', 'molecule_name', 'atom_index_0', 'atom_index_1', 'scalar_coupling_constant']\ntarget = 'scalar_coupling_constant'\n\nfi = pd.DataFrame()\nmodels = []\n\nsetSeeds(1234)\n\ndata_filt = data_train.drop(meta_cols, axis=1, errors='ignore')\npredictions = np.zeros(len(data_filt))\n\nkf = KFold(n_splits=FOLDS_VALID, random_state=1234, shuffle=True)\nfold_splits = kf.split(data_filt)\n\nfor i, (dev_index, val_index) in enumerate(fold_splits):\n\n    print('Fold', i)\n\n    Xt, Xv = data_filt.loc[dev_index, :], data_filt.loc[val_index, :]\n    yt, yv = data_train.loc[dev_index, target], data_train.loc[val_index, target]\n\n    d_train = lgb.Dataset(Xt, yt)\n    d_valid = lgb.Dataset(Xv, yv)\n\n    watchlist = [d_train, d_valid]\n    model = lgb.train(params_lgb,\n                      train_set=d_train,\n                      num_boost_round=500,\n                      valid_sets=watchlist,\n                      verbose_eval=100,\n                      early_stopping_rounds=100)\n\n    predictions[val_index] = model.predict(Xv, num_iteration = model.best_iteration)\n\n    fold_importance = pd.DataFrame()\n    fold_importance[\"feature\"] = Xt.columns\n    fold_importance[\"importance\"] = model.feature_importance()\n    fold_importance[\"fold\"] = i\n    fi = pd.concat([fi, fold_importance], axis=0)\n\n    models.append(model)\n\nprint('MAE', (predictions - data_train[target]).abs().mean())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"data_test_filt = data_test.drop(meta_cols, axis=1, errors='ignore')\n\npreds = np.zeros((FOLDS_VALID, len(data_test)))\nfor i in range(FOLDS_VALID):\n    print('Fold', i)\n    preds[i,:] = models[i].predict(data_test_filt, num_iteration=models[i].best_iteration)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub = pd.read_csv(PATH_BASE/'sample_submission.csv')\nsub['id'] = data_test['id']\nsub['scalar_coupling_constant'] = preds.mean(0)\nsub.to_csv(PATH_WORKING/'submission.csv', index=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fi['score'] = fi[[\"feature\", \"importance\"]].groupby('feature').transform('mean')\n\nplt.figure(figsize=(16, 12))\nsn.barplot(x=\"importance\", y=\"feature\", data=fi.sort_values(by=\"score\", ascending=False))\nplt.title('LGB Features (avg over folds)')","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}