{"cells":[{"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_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"import scipy as sp #collection of functions for scientific computing and advance mathematics\nfrom scipy import stats\nfrom scipy.stats import norm, skew\nfrom scipy.special import boxcox1p\nfrom scipy.stats import boxcox_normmax\n\n\n# warnings mute\nimport warnings\nwarnings.simplefilter('ignore')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Data ingestion**"},{"metadata":{"trusted":true},"cell_type":"code","source":"def reduce_mem_usage(df):\n    \"\"\" iterate through all the columns of a dataframe and modify the data type\n        to reduce memory usage.        \n    \"\"\"\n    start_mem = df.memory_usage().sum() / 1024**2\n    print('Memory usage of dataframe is {:.2f} MB'.format(start_mem))\n    \n    for col in df.columns:\n        col_type = df[col].dtype\n        \n        if col_type != object:\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                if c_min > np.finfo(np.float16).min and c_max < np.finfo(np.float16).max:\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:\n                    df[col] = df[col].astype(np.float32)\n                else:\n                    df[col] = df[col].astype(np.float64)\n        else:\n            df[col] = df[col].astype('category')\n\n    end_mem = df.memory_usage().sum() / 1024**2\n    print('Memory usage after optimization is: {:.2f} MB'.format(end_mem))\n    print('Decreased by {:.1f}%'.format(100 * (start_mem - end_mem) / start_mem))\n    \n    return df\n\n\ndef import_data(file):\n    \"\"\"create a dataframe and optimize its memory usage\"\"\"\n    df = pd.read_csv(file, parse_dates=True, keep_date_col=True)\n    df = reduce_mem_usage(df)\n    return df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train = import_data(\"../input/train.csv\")\ntest =import_data(\"../input/test.csv\")\nsub = import_data(\"../input/sample_submission.csv\")\nstructures = import_data(\"../input/structures.csv\")\n#dipole = import_data(\"../input/dipole_moments.csv\")\n#potential = import_data(\"../input/potential_energy.csv\")\n#mulliken = import_data(\"../input/mulliken_charges.csv\")\ncontributions = import_data(\"../input/scalar_coupling_contributions.csv\")\n#This dataframe to be used for merging with all chemichal features ( Electronegativity , radius , isotope ..)\nstructures_notreal = structures","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\n# Importing gc module \nimport gc \n  \n# Returns the number of \n# objects it has collected \n# and deallocated \ncollected = gc.collect() \n  \n# Prints Garbage collector  \n# as 0 object \nprint(\"Garbage collector: collected\", \n          \"%d objects.\" % collected) ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Feature engineering**\n\nSince \"O\" & \"F\" are not present in train data set ,\nI've decided to fill properties related to those atoms by a random value for all properties . "},{"metadata":{},"cell_type":"markdown","source":"![](https://physicsworld.com/wp-content/uploads/2019/01/PWFeb19Crease-modern.jpg)"},{"metadata":{},"cell_type":"markdown","source":""},{"metadata":{},"cell_type":"markdown","source":"***Chemical properties***\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"#electronegativity pauling\nelectronegativity = {'H':2.2, 'C':2.55, 'N':3.04 , 'O':3.44, 'F':3.98 }\n#Type de charge \ncharge ={'H':0, 'C':1, 'N':1 , 'O':3.44, 'F':3.98 }\n#etat\netat ={'H':0, 'C':1, 'N':0, 'O':3.44, 'F':3.98 }\n#Masse kg/m^3\nmasse = {'H':76, 'C':3513, 'N':1026, 'O':3.44, 'F':3.98 }\n#volume cm^3/mole\nvolume = {'H':13.26, 'C':3.42, 'N':13.65, 'O':3.44, 'F':3.98 }\n#Rayon atomique (mesuré)\nrayon_am = {'H':25, 'C':70, 'N':65, 'O':3.44, 'F':3.98 }\n#Rayon atomique (calculé)\nrayon_ac = {'H':53, 'C':67, 'N':56 ,'O':0.73, 'F':0.71}\n#Rayon covalent\nrayon_c = {'H':38, 'C':77, 'N':75, 'O':3.44, 'F':3.98 }\n#Rayon ionique\nrayon_i = {'H':-3, 'C':4, 'N':-3, 'O':3.44, 'F':3.98 }\n#Rayon de Van der Waals\nrayon_vdw = {'H':120, 'C':170, 'N':155, 'O':3.44, 'F':3.98 }\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"***Physical properties***"},{"metadata":{"trusted":true},"cell_type":"code","source":"#Point de fusion\nfusion = {'H':-259.1, 'C':3546.9, 'N':-209.9, 'O':3.44, 'F':3.98 }\n#seuil d'ébulution minimal (celcius)\nebulution_min = {'H':-252.9, 'C':4826.9, 'N':-195.8, 'O':3.44, 'F':3.98 }\n#Enthalpie de fusion ΔHf (kj/mol)\nenthalpie_fusion = {'H':0.12, 'C':105, 'N':0.72, 'O':3.44, 'F':3.98 }\n#Enthalpie de vaporisation ΔHv\nenthalpie_vaporisation = {'H':0.46, 'C':710.9, 'N':5.58, 'O':3.44, 'F':3.98 }\n#Capacité thermique\ncapacite_thermique = {'H':14.3, 'C':0.71, 'N':1.04, 'O':3.44, 'F':3.98 }\n#Conductivité thermique\nconductivite_thermique = {'H':0.18, 'C':990, 'N':0.03, 'O':3.44, 'F':3.98 }\n\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Nuclear properties"},{"metadata":{"trusted":true},"cell_type":"code","source":"\n#Nb isotopes\nisotopes = {'H':3, 'C':12, 'N':12, 'O':3.44, 'F':3.98 }\n#Isotopes emeteurs\nisotopes_emeteurs = {'H':0, 'C':3, 'N':3, 'O':3.44, 'F':3.98 }\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Create dictionnary to store all properties related to atoms in structures dataframe"},{"metadata":{"trusted":true},"cell_type":"code","source":"dico_chemical_elements = {'electronegativity':electronegativity ,\n                         'charge':charge,\n                          'etat':etat,\n                          'masse':masse,\n                          'volume':volume,\n                          'rayon_am':rayon_am,\n                          'rayon_ac':rayon_ac,\n                          'rayon_c':rayon_c,\n                          'rayon_i':rayon_i,\n                          'rayon_vdw':rayon_vdw,\n                          'fusion':fusion,\n                          'ebulution_min':ebulution_min,\n                          'enthalpie_fusion':enthalpie_fusion,\n                          'enthalpie_vaporisation':enthalpie_vaporisation,\n                          'capacite_thermique':capacite_thermique,\n                          'conductivite_thermique':conductivite_thermique,\n                          'isotopes':isotopes,\n                          'isotopes_emeteurs':isotopes_emeteurs\n                         }","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Merge structures with elements in dico of properties"},{"metadata":{"trusted":true},"cell_type":"code","source":"def dico_todf(list_dicos, df):\n    for k,v in list_dicos.items():\n        df[k] = df['atom'].apply(lambda x : v[x])\n    return df \n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Mapping between structures & train/test"},{"metadata":{"trusted":true},"cell_type":"code","source":"\ndef 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","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":" Compute diheral angle"},{"metadata":{"trusted":true},"cell_type":"code","source":"def dihedral_angle(data): \n        \n    vals = np.array(data[:, 3:6], dtype=np.float64)\n    mol_names = np.array(data[:, 0], dtype=np.str)\n \n    result = np.zeros((data.shape[0], 2), dtype=object)\n    # use every 4 rows to compute the dihedral angle\n    for idx in range(0, vals.shape[0] - 4, 4):\n\n        a0 = vals[idx]\n        a1 = vals[idx + 1]\n        a2 = vals[idx + 2]\n        a3 = vals[idx + 3]\n        \n        b0 = a0 - a1\n        b1 = a2 - a1\n        b2 = a3 - a2\n        \n        # normalize b1 so that it does not influence magnitude of vector\n        # rejections that come next\n        b1 /= np.linalg.norm(b1)\n    \n        # vector rejections\n        # v = projection of b0 onto plane perpendicular to b1\n        #   = b0 minus component that aligns with b1\n        # w = projection of b2 onto plane perpendicular to b1\n        #   = b2 minus component that aligns with b1\n\n        v = b0 - np.dot(b0, b1) * b1\n        w = b2 - np.dot(b2, b1) * b1\n\n        # angle between v and w in a plane is the torsion angle\n        # v and w may not be normalized but that's fine since tan is y/x\n        x = np.dot(v, w)\n        y = np.dot(np.cross(b1, v), w)\n       \n        # We want all 4 first rows for every molecule to have the same value\n        # (in order to have the same length as the dataframe)\n        result[idx:idx + 4] = [mol_names[idx], np.degrees(np.arctan2(y, x))]\n        \n    return result","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from datetime import datetime\nstartTime = datetime.now()\ndihedral = dihedral_angle(structures[structures.groupby('molecule_name')['atom_index'].transform('count').ge(4)].groupby('molecule_name').head(4).values)\nprint('Time elapsed (hh:mm:ss.ms) {}'.format(datetime.now() - startTime))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"themap = {k:v for k, v in dihedral if k}\n# Add diehral and cos diehral angle to features \nstructures['dihedral'] = structures['molecule_name'].map(themap)\nstructures['cosdihedral'] = structures['dihedral'].map(np.cos)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Add electronegativity & radius to structures"},{"metadata":{"trusted":true},"cell_type":"code","source":"atoms = structures['atom'].values\natoms_en = [electronegativity[x] for x in (atoms)]\natoms_rad = [rayon_ac[x] for x in (atoms)]\n\nstructures['EN'] = atoms_en\nstructures['rad'] = atoms_rad\n\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Add bonds to features"},{"metadata":{"trusted":true},"cell_type":"code","source":"#Add bonds to features\ni_atom = structures['atom_index'].values\np = structures[['x', 'y', 'z']].values\np_compare = p\nm = structures['molecule_name'].values\nm_compare = m\nr = structures['rad'].values\nr_compare = r\n\nsource_row = np.arange(len(structures))\nmax_atoms = 28\n\nbonds = np.zeros((len(structures)+1, max_atoms+1), dtype=np.int8)\nbond_dists = np.zeros((len(structures)+1, max_atoms+1), dtype=np.float32)\n\n\nfor i in (range(max_atoms-1)):\n    p_compare = np.roll(p_compare, -1, axis=0)\n    m_compare = np.roll(m_compare, -1, axis=0)\n    r_compare = np.roll(r_compare, -1, axis=0)\n    \n    mask = np.where(m == m_compare, 1, 0) #Are we still comparing atoms in the same molecule?\n    dists = np.linalg.norm(p - p_compare, axis=1) * mask\n    r_bond = r + r_compare\n    \n    bond = np.where(np.logical_and(dists > 0.0001, dists < r_bond), 1, 0)\n    \n    source_row = source_row\n    target_row = source_row + i + 1 #Note: Will be out of bounds of bonds array for some values of i\n    target_row = np.where(np.logical_or(target_row > len(structures), mask==0), len(structures), target_row) #If invalid target, write to dummy row\n    \n    source_atom = i_atom\n    target_atom = i_atom + i + 1 #Note: Will be out of bounds of bonds array for some values of i\n    target_atom = np.where(np.logical_or(target_atom > max_atoms, mask==0), max_atoms, target_atom) #If invalid target, write to dummy col\n    \n    bonds[(source_row, target_atom)] = bond\n    bonds[(target_row, source_atom)] = bond\n    bond_dists[(source_row, target_atom)] = dists\n    bond_dists[(target_row, source_atom)] = dists\n\nbonds = np.delete(bonds, axis=0, obj=-1) #Delete dummy row\nbonds = np.delete(bonds, axis=1, obj=-1) #Delete dummy col\nbond_dists = np.delete(bond_dists, axis=0, obj=-1) #Delete dummy row\nbond_dists = np.delete(bond_dists, axis=1, obj=-1) #Delete dummy col\n\n\nbonds_numeric = [[i for i,x in enumerate(row) if x] for row in (bonds)]\nbond_lengths = [[dist for i,dist in enumerate(row) if i in bonds_numeric[j]] for j,row in enumerate((bond_dists))]\nbond_lengths_mean = [ np.mean(x) for x in bond_lengths]\nn_bonds = [len(x) for x in bonds_numeric]\n\nbond_data = {'n_bonds':n_bonds, 'bond_lengths_mean': bond_lengths_mean }\nbond_df = pd.DataFrame(bond_data)\nstructures = structures.join(bond_df)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Compute distinct bonds between atoms for each mollecule "},{"metadata":{"trusted":true},"cell_type":"code","source":"df_struct_aux= structures.groupby(['molecule_name'])['atom'].agg([('bonds_distc',lambda x : len(set(x.tolist())))]).reset_index()\nstructures=pd.merge(structures,df_struct_aux)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"[](http://)Fetch atoms properties from dictionnary "},{"metadata":{"trusted":true},"cell_type":"code","source":"structures=dico_todf(dico_chemical_elements,structures)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"structures","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Estimate the position of the mollecule barycentre ( Without taking into account the mass of the atoms ) "},{"metadata":{"trusted":true},"cell_type":"code","source":"structures_aux_xbary=structures.groupby(['molecule_name'])[['x','y','z']].agg({'avg':np.average}).reset_index()\nstructures_aux_xbary_=pd.DataFrame(structures_aux_xbary.get_values())\nstructures_aux_xbary_.columns=['molecule_name','x_bar','y_bar','z_bar']\nstructures_bary = pd.merge(structures,structures_aux_xbary_)\ndel structures_bary['atom']\nstructures=structures_bary","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Merge structures and train/test ( Create feature_x0 & feature_x1)"},{"metadata":{"trusted":true},"cell_type":"code","source":"train = map_atom_info(train, 0)\ntrain = map_atom_info(train, 1)\ntest = map_atom_info(test, 0)\ntest = map_atom_info(test, 1)\ndel train['x_bar_y']\ndel train['y_bar_y']\ndel train['z_bar_y']\ndel test['x_bar_y']\ndel test['y_bar_y']\ndel test['z_bar_y']\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.columns","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%env JOBLIB_TEMP_FOLDER=/tmp","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.columns","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_p_0 = train[['x_0', 'y_0', 'z_0']].values\ntrain_p_1 = train[['x_1', 'y_1', 'z_1']].values\ntest_p_0 = test[['x_0', 'y_0', 'z_0']].values\ntest_p_1 = test[['x_1', 'y_1', 'z_1']].values\n\ntrain['dist'] = np.linalg.norm(train_p_0 - train_p_1, axis=1)\ntest['dist'] = np.linalg.norm(test_p_0 - test_p_1, axis=1)\ntrain['dist_x'] = (train['x_0'] - train['x_1']) ** 2\ntest['dist_x'] = (test['x_0'] - test['x_1']) ** 2\ntrain['dist_y'] = (train['y_0'] - train['y_1']) ** 2\ntest['dist_y'] = (test['y_0'] - test['y_1']) ** 2\ntrain['dist_z'] = (train['z_0'] - train['z_1']) ** 2\ntest['dist_z'] = (test['z_0'] - test['z_1']) ** 2","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Format type of jointure "},{"metadata":{"trusted":true},"cell_type":"code","source":"train['type_0'] = train['type'].apply(lambda x: x[0])\ntest['type_0'] = test['type'].apply(lambda x: x[0])\ntrain['type_1'] = train['type'].apply(lambda x: x[1:])\ntest['type_1'] = test['type'].apply(lambda x: x[1:])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Compute mean of distances"},{"metadata":{"trusted":true},"cell_type":"code","source":"train['dist_to_type_mean'] = train['dist'] / train.groupby('type')['dist'].transform('mean')\ntest['dist_to_type_mean'] = test['dist'] / test.groupby('type')['dist'].transform('mean')\n\ntrain['dist_to_type_0_mean'] = train['dist'] / train.groupby('type_0')['dist'].transform('mean')\ntest['dist_to_type_0_mean'] = test['dist'] / test.groupby('type_0')['dist'].transform('mean')\n\ntrain['dist_to_type_1_mean'] = train['dist'] / train.groupby('type_1')['dist'].transform('mean')\ntest['dist_to_type_1_mean'] = test['dist'] / test.groupby('type_1')['dist'].transform('mean')\n\ntrain[f'molecule_type_dist_mean'] = train.groupby(['molecule_name', 'type'])['dist'].transform('mean')\ntest[f'molecule_type_dist_mean'] = test.groupby(['molecule_name', 'type'])['dist'].transform('mean')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Compute distance between mollecule barycentre coordinates and atoms coordinates"},{"metadata":{"trusted":true},"cell_type":"code","source":"train['dist_barycentre_x0']=np.abs(train['x_bar_x'] - train['x_0'])\ntrain['dist_barycentre_y0']=np.abs(train['y_bar_x'] - train['y_0'])\ntrain['dist_barycentre_z0']=np.abs(train['z_bar_x'] - train['z_0'])\ntrain['dist_barycentre_x1']=np.abs(train['x_bar_x'] - train['x_1'])\ntrain['dist_barycentre_y1']=np.abs(train['y_bar_x'] - train['y_1'])\ntrain['dist_barycentre_z1']=np.abs(train['z_bar_x'] - train['z_1'])\ntest['dist_barycentre_x0']=np.abs(test['x_bar_x'] - test['x_0'])\ntest['dist_barycentre_y0']=np.abs(test['y_bar_x'] - test['y_0'])\ntest['dist_barycentre_z0']=np.abs(test['z_bar_x'] - test['z_0'])\ntest['dist_barycentre_x1']=np.abs(test['x_bar_x'] - test['x_1'])\ntest['dist_barycentre_y1']=np.abs(test['y_bar_x'] - test['y_1'])\ntest['dist_barycentre_z1']=np.abs(test['z_bar_x'] - test['z_1'])\n\n\ntrain['squared_measure_0'] =  np.sqrt(\n    np.array(   \n    np.power(train['dist_barycentre_x0'],2) +\n    \n    np.power(train['dist_barycentre_y0'],2) +\n    \n    np.power(train['dist_barycentre_z0'],2) \n     ,dtype=float   \n    ))\n\ntrain['squared_measure_1'] =   np.sqrt(\n    np.array(\n    np.power(train['dist_barycentre_x1'],2) +\n    \n    np.power(train['dist_barycentre_y1'],2) +\n    \n    np.power(train['dist_barycentre_z1'],2) \n,dtype=float   \n    ))\n\ntest['squared_measure_0'] =  np.sqrt(\n    np.array(\n    np.power(test['dist_barycentre_x0'],2) +\n    \n    np.power(test['dist_barycentre_y0'],2) +\n    \n    np.power(test['dist_barycentre_z0'],2) \n,dtype=float   \n    ))\n\ntest['squared_measure_1'] =   np.sqrt(\n    np.array(\n    np.power(test['dist_barycentre_x1'],2) +\n    \n    np.power(test['dist_barycentre_y1'],2) +\n    \n    np.power(test['dist_barycentre_z1'],2) \n     ,dtype=float   \n    )  )","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.columns","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn import preprocessing\nfor f in [ 'type_0', 'type_1', 'type']:\n    lbl = preprocessing.LabelEncoder()\n    lbl.fit(list(train[f].values) + list(test[f].values))\n    train[f] = lbl.transform(list(train[f].values))\n    test[f] = lbl.transform(list(test[f].values))\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#collect residual garbage\ngc.collect()\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Training**"},{"metadata":{"trusted":true},"cell_type":"code","source":"import os\nimport time\nimport datetime\nimport json\nimport gc\nfrom numba import jit\n\nimport numpy as np\nimport pandas as pd\n\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom tqdm import tqdm_notebook\n\nimport lightgbm as lgb\nimport xgboost as xgb\nfrom catboost import CatBoostRegressor, CatBoostClassifier\nfrom sklearn import metrics\n\nfrom itertools import product\n\nimport altair as alt\nfrom altair.vega import v5\nfrom IPython.display import HTML\n\n# using ideas from this kernel: https://www.kaggle.com/notslush/altair-visualization-2018-stackoverflow-survey\ndef prepare_altair():\n    \"\"\"\n    Helper function to prepare altair for working.\n    \"\"\"\n\n    vega_url = 'https://cdn.jsdelivr.net/npm/vega@' + v5.SCHEMA_VERSION\n    vega_lib_url = 'https://cdn.jsdelivr.net/npm/vega-lib'\n    vega_lite_url = 'https://cdn.jsdelivr.net/npm/vega-lite@' + alt.SCHEMA_VERSION\n    vega_embed_url = 'https://cdn.jsdelivr.net/npm/vega-embed@3'\n    noext = \"?noext\"\n    \n    paths = {\n        'vega': vega_url + noext,\n        'vega-lib': vega_lib_url + noext,\n        'vega-lite': vega_lite_url + noext,\n        'vega-embed': vega_embed_url + noext\n    }\n    \n    workaround = f\"\"\"    requirejs.config({{\n        baseUrl: 'https://cdn.jsdelivr.net/npm/',\n        paths: {paths}\n    }});\n    \"\"\"\n    \n    return workaround\n    \n\ndef add_autoincrement(render_func):\n    # Keep track of unique <div/> IDs\n    cache = {}\n    def wrapped(chart, id=\"vega-chart\", autoincrement=True):\n        if autoincrement:\n            if id in cache:\n                counter = 1 + cache[id]\n                cache[id] = counter\n            else:\n                cache[id] = 0\n            actual_id = id if cache[id] == 0 else id + '-' + str(cache[id])\n        else:\n            if id not in cache:\n                cache[id] = 0\n            actual_id = id\n        return render_func(chart, id=actual_id)\n    # Cache will stay outside and \n    return wrapped\n           \n\n@add_autoincrement\ndef render(chart, id=\"vega-chart\"):\n    \"\"\"\n    Helper function to plot altair visualizations.\n    \"\"\"\n    chart_str = \"\"\"\n    <div id=\"{id}\"></div><script>\n    require([\"vega-embed\"], function(vg_embed) {{\n        const spec = {chart};     \n        vg_embed(\"#{id}\", spec, {{defaultStyle: true}}).catch(console.warn);\n        console.log(\"anything?\");\n    }});\n    console.log(\"really...anything?\");\n    </script>\n    \"\"\"\n    return HTML(\n        chart_str.format(\n            id=id,\n            chart=json.dumps(chart) if isinstance(chart, dict) else chart.to_json(indent=None)\n        )\n    )\n    \n\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        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.float64).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\n    \n\n@jit\ndef fast_auc(y_true, y_prob):\n    \"\"\"\n    fast roc_auc computation: https://www.kaggle.com/c/microsoft-malware-prediction/discussion/76013\n    \"\"\"\n    y_true = np.asarray(y_true)\n    y_true = y_true[np.argsort(y_prob)]\n    nfalse = 0\n    auc = 0\n    n = len(y_true)\n    for i in range(n):\n        y_i = y_true[i]\n        nfalse += (1 - y_i)\n        auc += y_i * nfalse\n    auc /= (nfalse * (n - nfalse))\n    return auc\n\n\ndef eval_auc(y_true, y_pred):\n    \"\"\"\n    Fast auc eval function for lgb.\n    \"\"\"\n    return 'auc', fast_auc(y_true, y_pred), True\n\n\ndef group_mean_log_mae(y_true, y_pred, types, floor=1e-9):\n    \"\"\"\n    Fast metric computation for this competition: https://www.kaggle.com/c/champs-scalar-coupling\n    Code is from this kernel: https://www.kaggle.com/uberkinder/efficient-metric\n    \"\"\"\n    maes = (y_true-y_pred).abs().groupby(types).mean()\n    return np.log(maes.map(lambda x: max(x, floor))).mean()\n    \n\ndef train_model_regression(X, X_test, y, params, folds, model_type='lgb', eval_metric='mae', columns=None, plot_feature_importance=False, model=None,\n                               verbose=10000, early_stopping_rounds=200, n_estimators=50000):\n    \"\"\"\n    A function to train a variety of regression models.\n    Returns dictionary with oof predictions, test predictions, scores and, if necessary, feature importances.\n    \n    :params: X - training data, can be pd.DataFrame or np.ndarray (after normalizing)\n    :params: X_test - test data, can be pd.DataFrame or np.ndarray (after normalizing)\n    :params: y - target\n    :params: folds - folds to split data\n    :params: model_type - type of model to use\n    :params: eval_metric - metric to use\n    :params: columns - columns to use. If None - use all columns\n    :params: plot_feature_importance - whether to plot feature importance of LGB\n    :params: model - sklearn model, works only for \"sklearn\" model type\n    :params: verbose - parameters for gradient boosting models\n    :params: early_stopping_rounds - parameters for gradient boosting models\n    :params: n_estimators - parameters for gradient boosting models\n    \n    \"\"\"\n    columns = X.columns if columns is None else columns\n    X_test = X_test[columns]\n    \n    # to set up scoring parameters\n    metrics_dict = {'mae': {'lgb_metric_name': 'mae',\n                        'catboost_metric_name': 'MAE',\n                        'sklearn_scoring_function': metrics.mean_absolute_error},\n                    'group_mae': {'lgb_metric_name': 'mae',\n                        'catboost_metric_name': 'MAE',\n                        'scoring_function': group_mean_log_mae},\n                    'mse': {'lgb_metric_name': 'mse',\n                        'catboost_metric_name': 'MSE',\n                        'sklearn_scoring_function': metrics.mean_squared_error}\n                    }\n\n    \n    result_dict = {}\n    \n    # out-of-fold predictions on train data\n    oof = np.zeros(len(X))\n    \n    # averaged predictions on train data\n    prediction = np.zeros(len(X_test))\n    \n    # list of scores on folds\n    scores = []\n    feature_importance = pd.DataFrame()\n    \n    # split and train on folds\n    for fold_n, (train_index, valid_index) in enumerate(folds.split(X)):\n        print(f'Fold {fold_n + 1} started at {time.ctime()}')\n        if type(X) == np.ndarray:\n            X_train, X_valid = X[columns][train_index], X[columns][valid_index]\n            y_train, y_valid = y[train_index], y[valid_index]\n        else:\n            X_train, X_valid = X[columns].iloc[train_index], X[columns].iloc[valid_index]\n            y_train, y_valid = y.iloc[train_index], y.iloc[valid_index]\n            \n        if model_type == 'lgb':\n            model = lgb.LGBMRegressor(**params, n_estimators = n_estimators, n_jobs = -1)\n            model.fit(X_train, y_train, \n                    eval_set=[(X_train, y_train), (X_valid, y_valid)], eval_metric=metrics_dict[eval_metric]['lgb_metric_name'],\n                    verbose=verbose, early_stopping_rounds=early_stopping_rounds)\n            \n            y_pred_valid = model.predict(X_valid)\n            y_pred = model.predict(X_test, num_iteration=model.best_iteration_)\n            \n        if model_type == 'xgb':\n            train_data = xgb.DMatrix(data=X_train, label=y_train, feature_names=X.columns)\n            valid_data = xgb.DMatrix(data=X_valid, label=y_valid, feature_names=X.columns)\n\n            watchlist = [(train_data, 'train'), (valid_data, 'valid_data')]\n            model = xgb.train(dtrain=train_data, num_boost_round=20000, evals=watchlist, early_stopping_rounds=200, verbose_eval=verbose, params=params)\n            y_pred_valid = model.predict(xgb.DMatrix(X_valid, feature_names=X.columns), ntree_limit=model.best_ntree_limit)\n            y_pred = model.predict(xgb.DMatrix(X_test, feature_names=X.columns), ntree_limit=model.best_ntree_limit)\n        \n        if model_type == 'sklearn':\n            model = model\n            model.fit(X_train, y_train)\n            \n            y_pred_valid = model.predict(X_valid).reshape(-1,)\n            score = metrics_dict[eval_metric]['sklearn_scoring_function'](y_valid, y_pred_valid)\n            print(f'Fold {fold_n}. {eval_metric}: {score:.4f}.')\n            print('')\n            \n            y_pred = model.predict(X_test).reshape(-1,)\n        \n        if model_type == 'cat':\n            model = CatBoostRegressor(iterations=20000,  eval_metric=metrics_dict[eval_metric]['catboost_metric_name'], **params,\n                                      loss_function=metrics_dict[eval_metric]['catboost_metric_name'])\n            model.fit(X_train, y_train, eval_set=(X_valid, y_valid), cat_features=[], use_best_model=True, verbose=False)\n\n            y_pred_valid = model.predict(X_valid)\n            y_pred = model.predict(X_test)\n        \n        oof[valid_index] = y_pred_valid.reshape(-1,)\n        if eval_metric != 'group_mae':\n            scores.append(metrics_dict[eval_metric]['sklearn_scoring_function'](y_valid, y_pred_valid))\n        else:\n            scores.append(metrics_dict[eval_metric]['scoring_function'](y_valid, y_pred_valid, X_valid['type']))\n\n        prediction += y_pred    \n        \n        if model_type == 'lgb' and plot_feature_importance:\n            # feature importance\n            fold_importance = pd.DataFrame()\n            fold_importance[\"feature\"] = columns\n            fold_importance[\"importance\"] = model.feature_importances_\n            fold_importance[\"fold\"] = fold_n + 1\n            feature_importance = pd.concat([feature_importance, fold_importance], axis=0)\n\n    prediction /= folds.n_splits\n    \n    print('CV mean score: {0:.4f}, std: {1:.4f}.'.format(np.mean(scores), np.std(scores)))\n    \n    result_dict['oof'] = oof\n    result_dict['prediction'] = prediction\n    result_dict['scores'] = scores\n    \n    if model_type == 'lgb':\n        if plot_feature_importance:\n            feature_importance[\"importance\"] /= folds.n_splits\n            cols = feature_importance[[\"feature\", \"importance\"]].groupby(\"feature\").mean().sort_values(\n                by=\"importance\", ascending=False)[:50].index\n\n            best_features = feature_importance.loc[feature_importance.feature.isin(cols)]\n\n            plt.figure(figsize=(16, 12));\n            sns.barplot(x=\"importance\", y=\"feature\", data=best_features.sort_values(by=\"importance\", ascending=False));\n            plt.title('LGB Features (avg over folds)');\n            \n            result_dict['feature_importance'] = feature_importance\n        \n    return result_dict\n    \n\n\ndef train_model_classification(X, X_test, y, params, folds, model_type='lgb', eval_metric='auc', columns=None, plot_feature_importance=False, model=None,\n                               verbose=10000, early_stopping_rounds=200, n_estimators=50000):\n    \"\"\"\n    A function to train a variety of classification models.\n    Returns dictionary with oof predictions, test predictions, scores and, if necessary, feature importances.\n    \n    :params: X - training data, can be pd.DataFrame or np.ndarray (after normalizing)\n    :params: X_test - test data, can be pd.DataFrame or np.ndarray (after normalizing)\n    :params: y - target\n    :params: folds - folds to split data\n    :params: model_type - type of model to use\n    :params: eval_metric - metric to use\n    :params: columns - columns to use. If None - use all columns\n    :params: plot_feature_importance - whether to plot feature importance of LGB\n    :params: model - sklearn model, works only for \"sklearn\" model type\n    :params: verbose - parameters for gradient boosting models\n    :params: early_stopping_rounds - parameters for gradient boosting models\n    :params: n_estimators - parameters for gradient boosting models\n    \n    \"\"\"\n    columns = X.columns if columns == None else columns\n    X_test = X_test[columns]\n    \n    # to set up scoring parameters\n    metrics_dict = {'auc': {'lgb_metric_name': eval_auc,\n                        'catboost_metric_name': 'AUC',\n                        'sklearn_scoring_function': metrics.roc_auc_score},\n                    }\n    \n    result_dict = {}\n    \n    # out-of-fold predictions on train data\n    oof = np.zeros((len(X), len(set(y.values))))\n    \n    # averaged predictions on train data\n    prediction = np.zeros((len(X_test), oof.shape[1]))\n    \n    # list of scores on folds\n    scores = []\n    feature_importance = pd.DataFrame()\n    \n    # split and train on folds\n    for fold_n, (train_index, valid_index) in enumerate(folds.split(X)):\n        print(f'Fold {fold_n + 1} started at {time.ctime()}')\n        if type(X) == np.ndarray:\n            X_train, X_valid = X[columns][train_index], X[columns][valid_index]\n            y_train, y_valid = y[train_index], y[valid_index]\n        else:\n            X_train, X_valid = X[columns].iloc[train_index], X[columns].iloc[valid_index]\n            y_train, y_valid = y.iloc[train_index], y.iloc[valid_index]\n            \n        if model_type == 'lgb':\n            model = lgb.LGBMClassifier(**params, n_estimators=n_estimators, n_jobs = -1)\n            model.fit(X_train, y_train, \n                    eval_set=[(X_train, y_train), (X_valid, y_valid)], eval_metric=metrics_dict[eval_metric]['lgb_metric_name'],\n                    verbose=verbose, early_stopping_rounds=early_stopping_rounds)\n            \n            y_pred_valid = model.predict_proba(X_valid)\n            y_pred = model.predict_proba(X_test, num_iteration=model.best_iteration_)\n            \n        if model_type == 'xgb':\n            train_data = xgb.DMatrix(data=X_train, label=y_train, feature_names=X.columns)\n            valid_data = xgb.DMatrix(data=X_valid, label=y_valid, feature_names=X.columns)\n\n            watchlist = [(train_data, 'train'), (valid_data, 'valid_data')]\n            model = xgb.train(dtrain=train_data, num_boost_round=n_estimators, evals=watchlist, early_stopping_rounds=early_stopping_rounds, verbose_eval=verbose, params=params)\n            y_pred_valid = model.predict(xgb.DMatrix(X_valid, feature_names=X.columns), ntree_limit=model.best_ntree_limit)\n            y_pred = model.predict(xgb.DMatrix(X_test, feature_names=X.columns), ntree_limit=model.best_ntree_limit)\n        \n        if model_type == 'sklearn':\n            model = model\n            model.fit(X_train, y_train)\n            \n            y_pred_valid = model.predict(X_valid).reshape(-1,)\n            score = metrics_dict[eval_metric]['sklearn_scoring_function'](y_valid, y_pred_valid)\n            print(f'Fold {fold_n}. {eval_metric}: {score:.4f}.')\n            print('')\n            \n            y_pred = model.predict_proba(X_test)\n        \n        if model_type == 'cat':\n            model = CatBoostClassifier(iterations=n_estimators, eval_metric=metrics_dict[eval_metric]['catboost_metric_name'], **params,\n                                      loss_function=metrics_dict[eval_metric]['catboost_metric_name'])\n            model.fit(X_train, y_train, eval_set=(X_valid, y_valid), cat_features=[], use_best_model=True, verbose=False)\n\n            y_pred_valid = model.predict(X_valid)\n            y_pred = model.predict(X_test)\n        \n        oof[valid_index] = y_pred_valid\n        scores.append(metrics_dict[eval_metric]['sklearn_scoring_function'](y_valid, y_pred_valid[:, 1]))\n\n        prediction += y_pred    \n        \n        if model_type == 'lgb' and plot_feature_importance:\n            # feature importance\n            fold_importance = pd.DataFrame()\n            fold_importance[\"feature\"] = columns\n            fold_importance[\"importance\"] = model.feature_importances_\n            fold_importance[\"fold\"] = fold_n + 1\n            feature_importance = pd.concat([feature_importance, fold_importance], axis=0)\n\n    prediction /= folds.n_splits\n    \n    print('CV mean score: {0:.4f}, std: {1:.4f}.'.format(np.mean(scores), np.std(scores)))\n    \n    result_dict['oof'] = oof\n    result_dict['prediction'] = prediction\n    result_dict['scores'] = scores\n    \n    if model_type == 'lgb':\n        if plot_feature_importance:\n            feature_importance[\"importance\"] /= folds.n_splits\n            cols = feature_importance[[\"feature\", \"importance\"]].groupby(\"feature\").mean().sort_values(\n                by=\"importance\", ascending=False)[:50].index\n\n            best_features = feature_importance.loc[feature_importance.feature.isin(cols)]\n\n            plt.figure(figsize=(16, 12));\n            sns.barplot(x=\"importance\", y=\"feature\", data=best_features.sort_values(by=\"importance\", ascending=False));\n            plt.title('LGB Features (avg over folds)');\n            \n            result_dict['feature_importance'] = feature_importance\n        \n    return result_dict\n\n# setting up altair\nworkaround = prepare_altair()\nHTML(\"\".join((\n    \"<script>\",\n    workaround,\n    \"</script>\",\n)))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Reformat variables"},{"metadata":{"trusted":true},"cell_type":"code","source":"train[['x_bar_x', 'y_bar_x', 'z_bar_x', 'dist_barycentre_x0', 'dist_barycentre_y0', 'dist_barycentre_z0','dist_barycentre_x1', 'dist_barycentre_y1', 'dist_barycentre_z1']] = train[['x_bar_x', \n'y_bar_x', 'z_bar_x', 'dist_barycentre_x0', 'dist_barycentre_y0', 'dist_barycentre_z0', 'dist_barycentre_x1', 'dist_barycentre_y1', 'dist_barycentre_z1']].apply(pd.to_numeric)\ntest[['x_bar_x', 'y_bar_x', 'z_bar_x', 'dist_barycentre_x0', 'dist_barycentre_y0', 'dist_barycentre_z0','dist_barycentre_x1', 'dist_barycentre_y1', 'dist_barycentre_z1']] = test[['x_bar_x', \n'y_bar_x', 'z_bar_x', 'dist_barycentre_x0', 'dist_barycentre_y0', 'dist_barycentre_z0', 'dist_barycentre_x1', 'dist_barycentre_y1', 'dist_barycentre_z1']].apply(pd.to_numeric)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.columns","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Prepare X and Y for training \n"},{"metadata":{"trusted":true},"cell_type":"code","source":"X = train.drop(['id', 'molecule_name', 'scalar_coupling_constant'], axis=1)\ny = train['scalar_coupling_constant']\nX_test = test.drop(['id', 'molecule_name'], axis=1)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X.columns","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"del train\ndel structures\n#del test\ndel structures_notreal\ndel structures_bary\ndel structures_aux_xbary\ndel structures_aux_xbary_\ndel train_p_0\ndel train_p_1\ndel test_p_0\ndel test_p_1\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"y.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X1=reduce_mem_usage(X)\n#y1=reduce_mem_usage(y)\nX_test1 =reduce_mem_usage(X_test)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X1.columns","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"del X\ndel X_test\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Split data based on 3 folds "},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.model_selection import KFold\nn_fold = 3\nfolds = KFold(n_splits=n_fold, shuffle=True, random_state=11)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"params = {'num_leaves': 128,\n          'min_child_samples': 79,\n          'objective': 'regression',\n          'max_depth': 15,\n          'learning_rate': 0.1,\n          \"boosting_type\": \"gbdt\",\n          \"subsample_freq\": 1,\n          \"subsample\": 0.9,\n          \"bagging_seed\": 11,\n          \"metric\": 'mae',\n          \"verbosity\": -1,\n          'reg_alpha': 0.1,\n          'reg_lambda': 0.3,\n          'colsample_bytree': 1.0,\n          \n         }\nresult_dict_lgb = train_model_regression(X=X1, X_test=X_test1, y=y, params=params, folds=folds, model_type='lgb', eval_metric='group_mae', plot_feature_importance=True,\n                                                      verbose=1000, early_stopping_rounds=200, n_estimators=10000)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission=pd.DataFrame({'id':test[\"id\"], 'scalar_coupling_constant':result_dict_lgb['prediction']})\nsubmission.to_csv('submission.csv',index=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","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}