{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Summary\n\nThis kernel uses `DScribe`'s descriptor generation to get features for Champs. I refered [@educatedguessing's kernel here](https://www.kaggle.com/educatedguessing/diy-atom-centered-symmetry-functions) to generate the labels for the features. In his/her kernel, the feature generation is querying the `structures` dataframe, here I use `xyz` file reading, thus much faster.\n\nThe features generated are: **Atom-centered Symmetry Functions (ACSF)**, the columns inovling Fluorine are removed from get-go to avoid taking too much memory.\n\nACSFs are atom-based, in the actual features for the two atoms in a coupling pair, I used:\n\n- Simple mapping the features to each atom to `acsf_feat_x` for the Hydrogen atom, `acsf_feat_y` for the Carbon/Nitrogen/a second Hydrogen atom.\n- Charged-weighted harmonic mean $ab/(a+b)$, geometric mean $\\sqrt{ab}$ of the two. The charge weighted was suggested in [WACSF - Weighted Atom-Centered Symmetry Functions as Descriptors in Machine Learning Potentials](https://arxiv.org/abs/1712.05861). Since ACSFs all have exponential weights, it is not sensible to use simple arithmetic mean.\n- Various distance of $G^2$ and $G^4$ columns for the same atom combinations.\n\nThese features alone, together with feature selection techniques, are able to achieve -2.05 in the private LB using LGBM.\n\nReference:\n- [DIY: Atom-Centered Symmetry Functions](https://www.kaggle.com/educatedguessing/diy-atom-centered-symmetry-functions)"},{"metadata":{"trusted":true},"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nfrom numpy.linalg import eig\nfrom numpy.linalg import svd\nfrom sklearn.decomposition import PCA\nfrom tqdm import tqdm_notebook, tqdm\nimport gc\nimport os\n\nimport time, copy\n\nimport matplotlib.pyplot as plt\n%matplotlib inline\nimport seaborn as sns\nsns.set()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import ase\nfrom ase import Atoms\nimport ase.visualize\nfrom ase.build import molecule","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from dscribe.descriptors import ACSF, LMBTR, SOAP","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def view(molecule):\n    # Select a molecule\n    mol = structures[structures['molecule_name'] == molecule]\n    \n    # Get atomic coordinates\n    xcart = mol.iloc[:, 3:].values\n    \n    # Get atomic symbols\n    symbols = mol.iloc[:, 2].values\n    \n    # Display molecule\n    system = Atoms(positions=xcart, symbols=symbols)\n    print('Molecule Name: %s.' %molecule)\n    return ase.visualize.view(system, viewer=\"x3d\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"structures = pd.read_csv('../input/structures.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"all_molecule_names = structures['molecule_name'].unique()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"_DEBUG = True\nfolder = '../input/structures/'","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Testing I/O for one file"},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nf = open(f\"{folder}/dsgdb9nsd_133885.xyz\")\npositions = []\nsymbols = []\nfor row, line in enumerate(f):\n    print(row, line.replace('\\n', ''))\n    \n    fields = line.split(' ')\n    if row < 2:\n        continue\n    # Then rows of atomic positions and chemical symbols.\n    else:\n        positions.append(fields[1:4])\n        print(f\"{fields[0]} -> {symbols}\\n\")\n        symbols.append(fields[0])       \n\nprint(len(symbols))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# ACSF generation examples"},{"metadata":{},"cell_type":"markdown","source":"## Feature count\n\nLet $N := (\\# \\,\\text{unique atoms})$ which is the minimum among the number of `species` argument in `ACSF` and the actual unique atoms in a molecule, then the total number of features for *one* atom is:\n\n$$\nN\\times \\Big(\\underbrace{1}_{G^1}\n+ \\underbrace{(\\#\\, \\eta, R_S \\text{ couplets})}_{G^2}\\Big)\n+ \\underbrace{N(N+1)/2\\times (\\#\\, \\eta, \\zeta, \\lambda \\text{ triplets})}_{G^4}\n$$\n\nFor the $i$-th atom, let $\\mathcal{S}$ be set of atoms of the same symbol other than this atom of interest:\n$$G_{i}^{1}=\\sum_{j\\in\\mathcal{S} }\\; f_{\\mathrm{c}}\\left(R_{i j}\\right),$$\nand similarly,\n$$G_{i}^{2}=\\sum_{j\\in \\mathcal{S}}\\; e^{-\\eta\\left(R_{i j}-R_{s}\\right)^{2}} \\cdot f_{c}\\left(R_{i j}\\right).$$\n$$\n G_{i}^{4}= 2^{1-\\zeta} \\sum_{j, k \\in \\mathcal{S}\\backslash\\{ i\\} }^{\\text { all }}\\left(1+\\lambda \\cos \\theta_{i j k}\\right)^{\\zeta} \\cdot e^{-\\eta\\left(R_{i j}^{2}+R_{i k}^{2}+R_{j k}^{2}\\right)} \\cdot f_{\\mathrm{c}}\\left(R_{i j}\\right) \\cdot f_{\\mathrm{c}}\\left(R_{i k}\\right) \\cdot f_{\\mathrm{c}}\\left(R_{j k}\\right) \n$$\n\nThe arrangement of the features is by the feature being radial or angular, then by the `species` parameter. First $G^1$ with $G^2$ for `H`, `C`, `N`, `O`, `F`, then the $G^4$ for every two combinations of atom in `H`, `C`, `N`, `O`, `F`.\n\n"},{"metadata":{},"cell_type":"markdown","source":"### Testing G^1\n\nExample 1: water"},{"metadata":{"trusted":true},"cell_type":"code","source":"#For all ACSF functions R_c\nrcut = 10.0\n\nspecies = [\"H\", \"C\", \"N\", \"O\", \"F\"]\n\ng2_params = None\ng4_params = None\n\ng3_params = None\ng5_params = None\n\nacsf = ACSF(\n    species=species, # order of atomic number\n    rcut=rcut,\n    g2_params=g2_params,\n    g4_params=g4_params,\n)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"water = molecule(\"H2O\")\n\n# Create ACSF output\nacsf_water = acsf.create(water, positions=[0,1,2])\n\nprint(acsf_water.shape, \"\\n\")\nprint(water.get_all_distances(), \"\\n\")\nprint(f\"Atom O: {acsf_water[0]} \\n\")\nprint(f\"Atom H: {acsf_water[1]} \\n\")\nprint(f\"Atom H: {acsf_water[2]} \\n\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Testing G^2"},{"metadata":{"trusted":true},"cell_type":"code","source":"#G2 - eta/R_s couples:\ng2_params = [[0.5, 2], [0.1, 2]]\n\n#G4 - eta/zeta/lambda triplets:\n# g4_params = [[1, 4,  1], [0.1, 4,  1], [0.01, 4,  1]]\n\nacsf = ACSF(\n    species=[\"H\", \"C\", \"N\", \"O\", \"F\"], # order of atomic number\n    rcut=rcut,\n    g2_params=g2_params,\n    g4_params=None,\n)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"acsf_water = acsf.create(water, positions=[0,1,2])\n\nprint(acsf_water.shape, \"\\n\")\nprint(water.get_all_distances(), \"\\n\")\nprint(f\"Atom O\\n {acsf_water[0]} \\n\")\nprint(f\"Atom H\\n {acsf_water[1]} \\n\")\nprint(f\"Atom H\\n {acsf_water[2]} \\n\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Testing G^4"},{"metadata":{"trusted":true},"cell_type":"code","source":"#G2 - eta/R_s couples:\ng2_params = [[2, 1], [0.5, 1]]\n\n#G4 - eta/zeta/lambda triplets:\ng4_params = [[1, 4,  1], [0.1, 4,  1], [0.01, 4,  1]]\n\nacsf = ACSF(\n    species=species, # order of atomic number\n    rcut=rcut,\n    g2_params=g2_params,\n    g4_params=g4_params,\n)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"acsf_water = acsf.create(water, positions=[0,1,2])\n\nprint(acsf_water.shape, \"\\n\")\nprint(f\"Atom O\\n {acsf_water[0]} \\n\")\nprint(f\"Atom H\\n {acsf_water[1]} \\n\")\nprint(f\"Atom H\\n {acsf_water[2]} \\n\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Example 2: formic acid"},{"metadata":{"trusted":true},"cell_type":"code","source":"acid = molecule(\"HCOOH\")\nase.visualize.view(acid, viewer=\"x3d\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Create ACSF output\nacsf_acid = acsf.create(acid)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"k=1\nprint(acsf_acid.shape, \"\\n\")\nprint(acid.get_atomic_numbers(),\"\\n\")\nprint(acid.get_all_distances(), \"\\n\")\nprint(f\"All atoms {acid.get_chemical_symbols()}\")\nprint(f\"Atom {acid.get_chemical_symbols()[k]}:\\n {acsf_acid[k]} \\n\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Labeling features\n\nThe unsung hero: [DIY: Atom-Centered Symmetry Functions](https://www.kaggle.com/educatedguessing/diy-atom-centered-symmetry-functions)."},{"metadata":{"trusted":true},"cell_type":"code","source":"tmp_system = Atoms(species, [[0,0,0]]*len(species))\n\nnr_to_symbol = {number:symbol for symbol, number in\n                    zip(tmp_system.get_chemical_symbols(),tmp_system.get_atomic_numbers())\n                    }\natomic_numbers = sorted(tmp_system.get_atomic_numbers())\n\nlabel_feats = []\n\ng_params={\n        'g1': [rcut],\n        'g2': g2_params,\n        'g3': g3_params,\n        'g4': g4_params,\n        'g5': g5_params\n    }\n\nfor atom_nr in atomic_numbers:\n    atom_id = nr_to_symbol[atom_nr]\n    \n    for g in [\"g1\", \"g2\", \"g3\"]:\n        params = g_params[g]\n        if params is not None:\n            for para in params:\n                label_feats.append(f'acsf_{g}_{atom_id}_{para}')\n                \nfor atom_nr in atomic_numbers:\n    atom_id = nr_to_symbol[atom_nr]\n    for i in range(0, atom_nr+1): # this is worth noting\n        if i in atomic_numbers:\n            atom_id_2 = nr_to_symbol[i]\n            for g in [\"g4\",\"g5\"]:\n                params = g_params[g]\n                if params is not None:\n                    for para in params:\n                        label_feats.append(f'acsf_{g}_{atom_id}_{atom_id_2}_{para}')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def create_feature_labels(species, rcut,\n                          g2_params=None,\n                          g3_params=None,\n                          g4_params=None,\n                          g5_params=None,\n                          transform_to_symbols=True):\n    #sub function to transform from atom numbers to chemical symbols\n    def get_atom_id(atom_nr, transform_to_symbols):\n\n        if transform_to_symbols == True:\n            atom_id = nr_to_symbol[atom_nr]\n        else:\n            atom_id = atom_nr\n        return atom_id\n\n    feature_label = []\n\n    g_params={\n        'g1': [rcut],\n        'g2': g2_params,\n        'g3': g3_params,\n        'g4': g4_params,\n        'g5': g5_params\n    }\n\n    # create_atom_numbers -> symbol dict\n    tmp_system = ase.Atoms(species, [[0,0,0]]*len(species))\n\n    nr_to_symbol = {number:symbol for symbol, number in\n                    zip(tmp_system.get_chemical_symbols(),tmp_system.get_atomic_numbers())\n                    }\n\n    atomic_numbers = sorted(tmp_system.get_atomic_numbers())\n\n    for atom_nr in atomic_numbers:\n        atom_id = get_atom_id(atom_nr, transform_to_symbols)\n        for g in [\"g1\", \"g2\", \"g3\"]:\n            params = g_params[g]\n            if params is not None:\n                for para in params:\n                    feature_label.append(f'acsf_{g}_{atom_id}_{para}')\n\n    for atom_nr in atomic_numbers:\n        atom_id = get_atom_id(atom_nr, transform_to_symbols)\n        for i in range(0, atom_nr+1):\n            if i in atomic_numbers:\n                atom_id_2 = get_atom_id(i, transform_to_symbols)\n                for g in [\"g4\",\"g5\"]:\n                    params = g_params[g]\n                    if params is not None:\n                        for para in params:\n                            feature_label.append(f'acsf_{g}_{atom_id}_{atom_id_2}_{para}')\n\n    return feature_label","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"labels_tmp = create_feature_labels([\"H\", \"C\", \"O\", \"N\", \"F\"], rcut,\n                          g2_params=g2_params,\n                          g3_params=None,\n                          g4_params=g4_params,\n                          g5_params=None,\n                          transform_to_symbols=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"## sanity check\nlabels_tmp == label_feats","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"acid_feat = pd.DataFrame(acsf_acid, columns=label_feats)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"acid_feat.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Benchmark using few molecules"},{"metadata":{"trusted":true},"cell_type":"code","source":"filenames = [folder + i + '.xyz' for i in all_molecule_names[:100]]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Final parameters for CHAMPS competition\n\nBased on the original paper: \n\n- Jörg Behler. Atom-centered symmetry functions for constructing high-dimensional neural network potentials. J. Chem. Phys., 134(7):074106, 2011.\n\nI decided to use the following paramters to cover all three coupling types."},{"metadata":{"trusted":true},"cell_type":"code","source":"rcut = 9.0\n\nspecies = [\"H\", \"C\", \"N\", \"O\", \"F\"]\n\n#G2 - eta/R_s couples:\ng2_params = [[5, 1], [2, 1], [0.5, 1],\n             [5, 3], [2, 3], [0.5, 3]]\n\n#G4 - eta/zeta/lambda triplets:\ng4_params = [[0.5, 2,  1], [0.5, 6,  1], [0.5, 16,  1], \n             [0.05, 2,  1],[0.05, 6,  1], [0.05, 16,  1], \n             [0.5, 2, -1], [0.5, 6, -1], [0.5, 16, -1], \n             [0.05, 2, -1], [0.05, 6, -1], [0.05, 16, -1]]\n\nacsf = ACSF(\n    species=species, # order of atomic number\n    rcut=rcut,\n    g2_params=g2_params,\n    g4_params=g4_params,\n)\n\nlabel_feats = create_feature_labels([\"H\", \"C\", \"O\", \"N\", \"F\"], rcut,\n                          g2_params=g2_params,\n                          g3_params=None,\n                          g4_params=g4_params,\n                          g5_params=None,\n                          transform_to_symbols=True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Iterating in xyz file, numpy array then dataframe"},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nimages = []\n\n#Open and parse file.\nfor filename in tqdm_notebook(filenames):\n    # Define lists for output.\n    positions = []\n    symbols = []\n    with open(filename) as f:\n        for row, line in enumerate(f):\n            fields = line.split(' ')\n            # Each file contains a 3 line header.\n            if row < 2:\n                continue\n            # Then rows of atomic positions and chemical symbols.\n            else:\n                positions.append(fields[1:4])\n                symbols.append(fields[0])\n    # Make an atoms object from each file.\n    atom = Atoms(positions=np.array(positions, dtype=float),\n                  symbols=symbols)\n    features = acsf.create(atom, n_jobs=2) \n    # structure of return is [[#acsf features] for each position in molecule_system]\n    images.append(features)\n\nstructures_acsf = pd.DataFrame(np.concatenate(images),columns=label_feats)\nprint(structures_acsf.shape,'\\n')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Query the dataframe to generate ACSF\n\nSlow. Slow. Slow.\n\nhttps://www.kaggle.com/educatedguessing/diy-atom-centered-symmetry-functions"},{"metadata":{"trusted":true},"cell_type":"code","source":"def calculate_symmetric_functions(df_structure, rcut, g2_params=None,\n                                  g3_params=None,\n                                  g4_params=None,\n                                  g5_params=None):\n\n    species = [\"H\", \"C\", \"O\", \"N\", \"F\"]\n\n    acsf = ACSF(\n        species=species,\n        rcut=rcut,\n        g2_params=g2_params,\n        g3_params=g3_params,\n        g4_params=g4_params,\n        g5_params=g5_params,\n    )\n\n    structure_molecules = df_structure.molecule_name.unique()\n\n    acsf_feature_labels = create_feature_labels(species=species,\n                                                rcut=rcut,\n                                                g2_params=g2_params,\n                                                g3_params=g3_params,\n                                                g4_params=g4_params,\n                                                g5_params=g5_params,\n                                                )\n\n    df_structure= df_structure.reindex(columns = df_structure.columns.tolist() + acsf_feature_labels)\n\n    df_structure = df_structure.sort_values(['molecule_name','atom_index'])\n\n    acsf_structure_chunks = calculate_acsf_in_chunks(structure_molecules, df_structure, acsf, acsf_feature_labels)\n\n    acsf_structure = pd.DataFrame().append(acsf_structure_chunks)\n\n    return acsf_structure\n\ndef calculate_acsf_in_chunks(structure_molecules, \n                             df_structure, acsf, acsf_feature_labels, \n                             step_size=2000):\n\n    mol_counter = 0\n    max_counter = len(structure_molecules)\n    all_chunks = []\n    \n    while mol_counter*step_size < max_counter:\n\n        tmp_molecules = structure_molecules[mol_counter*step_size:(mol_counter+1)*step_size]\n\n        tmp_structure = df_structure.loc[df_structure.molecule_name.isin(tmp_molecules),:].copy()\n\n        tmp_results = calculate_acsf_multiple_molecules(tmp_molecules, tmp_structure, acsf, acsf_feature_labels)\n\n        all_chunks.append(tmp_results.copy())\n\n        mol_counter += 1\n\n    return all_chunks\n\n\ndef calculate_acsf_multiple_molecules(molecule_names, df_structure, acsf, acsf_feature_labels):\n\n    #acsf_feature_labels = [f'acsf_{nr}' for nr in range(0, acsf.get_number_of_features())]\n    #df_molecules = df_structure.loc[df_structure.molecule_name.isin(molecule_names),:].copy()\n    counter = 0\n    for molecule_name in molecule_names:\n\n        df_molecule = df_structure.loc[df_structure.molecule_name == molecule_name,:]\n        acsf_values = calculate_acsf_single_molecule(df_molecule, acsf)\n\n\n        df_structure.loc[df_structure.molecule_name==molecule_name, acsf_feature_labels] \\\n        = copy.copy(acsf_values)\n\n        counter += 1\n\n    return df_structure\n\ndef calculate_acsf_single_molecule(df_molecule, acsf):\n\n    molecule_atoms = df_molecule.loc[:, 'atom']\n    molecule_positions = df_molecule.loc[:, ['x','y','z']]\n\n    molecule_system = Atoms(symbols=molecule_atoms, positions=molecule_positions)\n\n    return acsf.create(molecule_system, n_jobs=2)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"struct_small = structures.loc[structures.molecule_name.isin(all_molecule_names[:100])].copy()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nstruct_small = calculate_symmetric_functions(struct_small, \n                              rcut, \n                              g2_params=g2_params, \n                              g4_params=g4_params)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"struct_small.head(5)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"## sanity check\nnp.allclose(struct_small[label_feats], structures_acsf)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Checking correlation of ACSF features"},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_correlated_cols(df,threshold=0.98):\n    '''\n    threshold: threshold to remove correlated variables\n    '''\n    \n    # Absolute value correlation matrix\n    corr_matrix = df.corr().abs()\n    \n    # Getting the upper triangle of correlations\n    upper = corr_matrix.where(np.triu(np.ones(corr_matrix.shape), k=1).astype(np.bool))\n    \n    # Select columns with correlations above threshold\n    cols_to_drop = [column for column in upper.columns if any(upper[column] > threshold)]\n    \n    print('There are {} columns to remove.'.format(len(cols_to_drop)))\n    return cols_to_drop","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"corr_cols = get_correlated_cols(structures_acsf, threshold=0.97)\nprint(corr_cols)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"np.abs(structures_acsf.corrwith(structures_acsf[corr_cols[1]])).sort_values(ascending=False)[:10]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"np.abs(structures_acsf.corrwith(structures_acsf.iloc[:,0])).sort_values(ascending=False)[:10]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Generation for all molecules"},{"metadata":{"trusted":true},"cell_type":"code","source":"filenames = [folder + i + '.xyz' for i in all_molecule_names]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nimages = []\n\n#Open and parse file.\nfor filename in tqdm_notebook(filenames):\n    # Define lists for output.\n    positions = []\n    symbols = []\n    with open(filename) as f:\n        for row, line in enumerate(f):\n            fields = line.split(' ')\n            # Each file contains a 3 line header.\n            if row < 2:\n                continue\n            # Then rows of atomic positions and chemical symbols.\n            else:\n                positions.append(fields[1:4])\n                symbols.append(fields[0])\n    # Make an atoms object from each file.\n    atom = Atoms(positions=np.array(positions, dtype=float),\n                  symbols=symbols)\n    features = acsf.create(atom, n_jobs=2) \n    # structure of return is [[#acsf features] for each position in molecule_system]\n    images.append(features)\n\nstructures_acsf = pd.DataFrame(np.concatenate(images),columns=label_feats)\nprint('\\n', structures_acsf.shape,'\\n')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"## drop Fluorine columns\nF_columns = [col for col in structures_acsf.columns if '_F_' in col]\nstructures_acsf = structures_acsf.drop(columns=F_columns)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"structures_acsf.to_csv('structures_acsf.csv', index=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Addendum: feature generation functions\n\nAfter mapping the `structures_acsf` to the `train` and `test`, the following the feature generation functions."},{"metadata":{"trusted":true},"cell_type":"code","source":"TOL = 1e-12\n\ndef get_chi2_distance(v1, v2):\n    '''\n    all columns must be non-negative\n    compute the weighted Chi-square distance\n    '''  \n    diff = ((v1 - v2)**2)/(v1+v2+TOL)\n    \n    return diff.sum(axis=1)\n\ndef get_angular_distance(v1, v2):\n    '''\n    Compute the cosine distance along axis 1\n    inputs: 2 n by m array\n    '''\n    \n    cosine = (v1*v2).sum(axis=1)/(norm(v1,axis=1)*norm(v2,axis=1)+TOL)\n    \n    return cosine\n\ndef get_tanimoto_distance(v1, v2):\n    '''\n    Compute the Tanimoto similarity\n    '''\n    a = (v1*v2).sum(axis=1)\n    b = (v1*v1).sum(axis=1)\n    c = (v2*v2).sum(axis=1)\n    \n    return a/(b + c - a + TOL)\n    \n\ndef add_acsf_features(df):\n    \n    acsf_cols = []\n    for col in df.columns:\n        if 'acsf' in col:\n            acsf_cols.append(col)\n            \n    #### G1 difference features\n    g1_cols = [col for col in acsf_cols if 'g1' in col]\n    g1_cols_atom0 = [col for col in g1_cols if 'x' in col]\n    g1_cols_atom1 = [col for col in g1_cols if 'y' in col]\n    \n    v1 = df[g1_cols_atom0].values\n    v2 = df[g1_cols_atom1].values\n    \n    df['acsf_g1_diff'] = get_chi2_distance(v1, v2)\n    df['acsf_g1_cos'] = get_angular_distance(v1, v2)\n    df['acsf_g1_tanimoto'] = get_tanimoto_distance(v1, v2)\n    \n    #### G2 difference features\n    g2_cols = [col for col in acsf_cols if 'g2' in col]\n    for symbol in ['H', 'C', 'N', 'O', 'F']:\n        \n        g2_cols_atom0 = [col for col in g2_cols if 'x' in col if symbol in col]\n        g2_cols_atom1 = [col for col in g2_cols if 'y' in col if symbol in col]\n        \n        v1 = df[g2_cols_atom0].values\n        v2 = df[g2_cols_atom1].values\n        \n        df['acsf_g2_diff_'+str(symbol)] = get_chi2_distance(v1, v2)\n        df['acsf_g2_cos_'+str(symbol)] = get_angular_distance(v1, v2)\n        df['acsf_g2_tanimoto_'+str(symbol)] = get_tanimoto_distance(v1, v2)\n        \n        \n    #### G4 difference features\n    g4_cols = [col for col in acsf_cols if 'g4' in col]\n    \n    g4_pairs = []\n    all_symbol = ['H', 'C', 'N', 'O' ]\n    for i, s in enumerate(all_symbol):\n        for j in range(i+1):\n            g4_pairs.append(str(s)+'_'+str(all_symbol[j]))\n            \n    for pair in g4_pairs:\n        \n        g4_cols_atom0 = [col for col in g4_cols if 'x' in col if symbol in col]\n        g4_cols_atom1 = [col for col in g4_cols if 'y' in col if symbol in col]\n        \n        v1 = df[g4_cols_atom0].values\n        v2 = df[g4_cols_atom1].values\n        \n        df['acsf_g4_diff_'+str(pair)] = get_chi2_distance(v1, v2)\n        df['acsf_g4_cos_'+str(pair)] = get_angular_distance(v1, v2)\n        df['acsf_g4_tanimoto_'+str(pair)] = get_tanimoto_distance(v1, v2)\n    \n    return df\n\ndef add_prod_features(df, cols=None, weights=None):\n    if weights is not None and isinstance(weights, pd.DataFrame):\n        weights = weights.values\n    if cols is not None:\n        for col in cols:\n            if col+'_x' in df.columns and col+'_y' in df.columns:\n                df[col+'_prod'] = np.sqrt(weights[:,0]*weights[:,1]*df[col+'_x']*df[col+'_y'])\n    return df\n\ndef add_mean_features(df, cols=None, weights=None):\n    if weights is not None and isinstance(weights, pd.DataFrame):\n        weights = weights.values\n    if cols is not None:\n        for col in cols:\n            if col+'_x' in df.columns and col+'_y' in df.columns:\n                val_atom_0 = weights[:,0]*df[col+'_x']\n                val_atom_1 = weights[:,1]*df[col+'_y']\n                val_atom_0 = np.abs(val_atom_0)\n                val_atom_1 = np.abs(val_atom_1)\n                val_atom_0[val_atom_0<1e-13] = 1e-13\n                val_atom_1[val_atom_1<1e-13] = 1e-13\n                df[col+'_hmean'] = hmean(np.c_[val_atom_0,val_atom_1], axis=1)\n    return df","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}