Repository: aspuru-guzik-group/selfies

Issue:
Feature request: Retrieve mapping betwewen SMILES and SELFIES tokens
Is it possible to get a map which SMILES tokens were used to generate which SELFIES tokens (or v.v.)? 

I am looking for a feature like this:
```py
>>>smiles = 'CCO'
>>>encoder(smiles, get_mapping=True)
([C][C][O], [0,1,2])
```
In this simple example `[0,1,2]` would imply that the first SMILES token (`C`) is mapped to the first selfies token (`[C]`) and so on.

**Motivation**:
I think this feature could be very useful to close the gap between RDKit and SELFIES. One example are scaffolds. Say we have a molecule, want to retrieve its scaffold and decorate it with a generative model. With SMILES it's easy (see example below) but with SELFIES it's not possible (as far as I understand).


**My questions:**
1) Is it, in principle, possible to obtain such a mapping?
2) If yes, is there already a way to obtain it with the current package?
3) If no, is this a feature that seems worth implementing?

**Discussion:**
Such a mapping would imply a standardized way of splitting the strings into tokens. Fortunately, we have `split_selfies` already, but regarding SMILES, I think that the tokenizer from the _Found in Translation_  [paper](https://pubs.rsc.org/lv/content/articlehtml/2018/sc/c8sc02339e)!) could be a good choice since it's used widely. (I'm using that tokenizer in the example below.)



**==== EXAMPLE ===**
This is just the appendix to the post. It's an example for how to retrieve which SMILES tokens constitute the scaffold of a given molecule. As it appears to me, this is currently not possible with SELFIES.

First, some boring setup:
```py
from rdkit import Chem
from selfies import encoder, decoder, split_selfies
from rdkit.Chem.Scaffolds.MurckoScaffold import GetScaffoldForMol
from pytoda.smiles.processing import tokenize_smiles
import re 

# Setup tokenizer
NON_ATOM_CHARS = set(list(map(str, range(1, 10))) + ['/', '\\', '(', ')', '#', '=', '.', ':', '-'])
regexp = re.compile(
    r'(\[[^\]]+]|Br?|Cl?|N|O|S|P|F|I|b|c|n|o|s|p|\(|\)|\.|=|#|'
    r'-|\+|\\\\|\/|:|~|@|\?|>|\*|\$|\%[0-9]{2}|[0-9])'
)
smiles_tokenizer = lambda smi: [token for token in regexp.split(smi) if token]
```

Example molecule (left) and RDKit-extracted scaffold (right):
![Screenshot 2021-05-06 at 11 59 07](https://user-images.githubusercontent.com/15703818/117286474-891b5100-ae69-11eb-9bf0-3e5026a66f83.png)

```
smiles = 'CCOc1[nH]c(N=Cc2ccco2)c(C#N)c1C#N'
mol = Chem.MolFromSmiles(smiles)
atom_symbols = [atom.GetSymbol() for atom in mol.GetAtoms()]
scaffold = GetScaffoldForMol(mol)
# List of ints pointing to scaffold atoms as they occur in SMILES
scaffold_atoms = mol.GetSubstructMatches(scaffold)[0]
smiles_tokens =  smiles_tokenizer(smiles)


atom_id = -1
for token in smiles_tokens:
    if token not in NON_ATOM_CHARS:
        # Found atom
        atom_id += 1
        if atom_id in scaffold_atoms:
            print(token, '--> on scaffold')
        else:
            print(token, '--> not on scaffold')
    else:
        # Non-Atom-Chars
        if (atom_id in scaffold_atoms and atom_id+1 in scaffold_atoms) or atom_id==scaffold_atoms[-1]:
            print(token, '--> on scaffold')
        else:
            print(token, '--> not on scaffold')
    
```
Output will be:

```
C --> not on scaffold
C --> not on scaffold
O --> not on scaffold
c --> on scaffold
1 --> on scaffold
[nH] --> on scaffold
c --> on scaffold
( --> on scaffold
N --> on scaffold
= --> on scaffold
C --> on scaffold
c --> on scaffold
2 --> on scaffold
c --> on scaffold
c --> on scaffold
c --> on scaffold
o --> on scaffold
2 --> on scaffold
) --> on scaffold
c --> on scaffold
( --> not on scaffold
C --> not on scaffold
# --> not on scaffold
N --> not on scaffold
) --> not on scaffold
c --> on scaffold
1 --> on scaffold
C --> not on scaffold
# --> not on scaffold
N --> not on scaffold
```

Trying to achieve the same with SELFIES does not seem to work. This is because `selfies.encoder` does not fully preserve the order of the tokens passed. It preserves it to large extents (which is great) but around ring symbols it usually breaks. I feel like I would need to reverse-engineer the context free grammar to solve this. 
Here would be the tokens in SMILES and SELFIES respectively:
```
C [C]
C [C]
O [O]
c [C]
1 [NHexpl]
[nH] [C]
c [Branch1_1]
( [Branch2_3]
N [N]
= [=C]
C [C]
c [=C]
2 [C]
c [=C]
c [O]
c [Ring1]
o [Branch1_1]
2 [=C]
) [Branch1_1]
c [Ring1]
( [C]
C [#N]
# [C]
N [Expl=Ring1]
) [=C]
c [C]
1 [#N]
C
#
N
```

Discussion:

Patch:
diff --git a/selfies/utils/smiles_utils.py b/selfies/utils/smiles_utils.py
index a93bfeb..c3d6aed 100644
--- a/selfies/utils/smiles_utils.py
+++ b/selfies/utils/smiles_utils.py
@@ -227,7 +227,6 @@ def _derive_mol_from_tokens(mol, smiles, tokens, i):
     prev_stack.append(tok)
     while tokens:
         tok = tokens.popleft()
-        i += 1
         bond_char = tok.extract_bond_char(smiles)
         symbol, symbol_type = tok.extract_symbol(smiles), tok.token_type
         prev_atom = prev_stack[-1]
@@ -241,7 +240,7 @@ def _derive_mol_from_tokens(mol, smiles, tokens, i):
                 err_msg = "invalid atom symbol '{}'".format(symbol)
                 raise SMILESParserError(smiles, err_msg, tok.start_idx)
 
-            curr = _attach_atom(mol, bond_char, curr, prev_atom, i, tok)
+            curr, i = _attach_atom(mol, bond_char, curr, prev_atom, i, tok)
             prev_stack.pop()
             prev_stack.append(curr)
             chain_start = False
@@ -277,6 +276,7 @@ def _derive_mol_from_tokens(mol, smiles, tokens, i):
         else:
             # should not happen
             raise Exception("invalid symbol type")
+        i += 1
 
     if len(mol) == 0:
         err_idx = (len(smiles) if (tok is None) else tok.start_idx) - 1
@@ -295,6 +295,8 @@ def _derive_mol_from_tokens(mol, smiles, tokens, i):
 
 def _attach_atom(mol, bond_char, atom, prev_atom, i, tok):
     is_root = (prev_atom is None)
+    if bond_char:
+        i += 1
     o = mol.add_atom(atom, mark_root=is_root)
     mol.add_attribution(o, [(i, str(tok))])
     if not is_root:
@@ -304,7 +306,7 @@ def _attach_atom(mol, bond_char, atom, prev_atom, i, tok):
             order = 1.5  # handle implicit aromatic bonds, e.g. cc
         o = mol.add_bond(src=src, dst=dst, order=order, stereo=stereo)
         mol.add_attribution(o, [(i, str(tok))])
-    return atom
+    return atom, i
 
 
 def _make_ring_bonds(mol, smiles, ltoken, latom, lpos, rtoken, ratom):


Test patch:
diff --git a/tests/test_selfies.py b/tests/test_selfies.py
index 29980a2..a1c7bfd 100644
--- a/tests/test_selfies.py
+++ b/tests/test_selfies.py
@@ -148,10 +148,14 @@ def test_decoder_attribution():
 
 
 def test_encoder_attribution():
-    s, am = sf.encoder(
-        "C1([O-])C=CC=C1Cl", attribute=True)
+    smiles = "C1([O-])C=CC=C1Cl"
+    indices = [0, 3, 3, 3, 5, 7, 8, 10, None, None, 12]
+    s, am = sf.encoder(smiles, attribute=True)
     # check that Cl lined up
-    for ta in am:
+    for i, ta in enumerate(am):
+        if ta[1]:
+            assert indices[i] == ta[1][0][0], \
+                f'found {ta[1]}; should be {indices[i]}'
         if ta[0] == '[Cl]':
             for i, v in ta[1]:
                 if v == 'Cl':