Bob bug (#33)

* Fixed ordering in Bob, added integration test

* Tightened threshold in Bob integration test to 2.8 kcal/mol

Author: andersx

qml/compound.py

@@ -214,7 +214,7 @@ def generate_atomic_coulomb_matrix(self, size = 23, sorting = "row-norm",
             sorting = sorting, central_cutoff = central_cutoff, central_decay = central_decay,
             interaction_cutoff = interaction_cutoff, interaction_decay = interaction_decay)
 
-    def generate_bob(self, asize = {"O":3, "C":7, "N":3, "H":16, "S":1}):
+    def generate_bob(self, size=23, asize = {"O":3, "C":7, "N":3, "H":16, "S":1}):
         """ Creates a Bag of Bonds (BOB) representation of a molecule.
             The representation expands on the coulomb matrix representation.
             For each element a bag (vector) is constructed for self interactions

qml/representations.py

@@ -30,14 +30,41 @@
 from .frepresentations import fgenerate_local_coulomb_matrix
 from .frepresentations import fgenerate_atomic_coulomb_matrix
 from .frepresentations import fgenerate_eigenvalue_coulomb_matrix
-from .frepresentations import fgenerate_bob
 
 from .data import NUCLEAR_CHARGE
 
 from .slatm import get_boa
 from .slatm import get_sbop
 from .slatm import get_sbot
 
+def vector_to_matrix(v):
+    """ Converts a representation from 1D vector to 2D square matrix.
+    :param v: 1D input representation.
+    :type v: numpy array 
+    :return: Square matrix representation.
+    :rtype: numpy array 
+    """
+
+    if not (np.sqrt(8*v.shape[0]+1) == int(np.sqrt(8*v.shape[0]+1))):
+        print("ERROR: Can not make a square matrix.")
+        exit(1)
+
+    n = v.shape[0]
+    l = (-1 + int(np.sqrt(8*n+1)))//2
+    M = np.empty((l,l))
+
+    index = 0
+    for i in range(l):
+        for j in range(l):
+            if j > i:
+                continue
+
+            M[i,j] = v[index]
+            M[j,i] = M[i,j]
+
+            index += 1
+    return M
+
 def generate_coulomb_matrix(nuclear_charges, coordinates, size = 23, sorting = "row-norm"):
     """ Creates a Coulomb Matrix representation of a molecule.
         Sorting of the elements can either be done by ``sorting="row-norm"`` or ``sorting="unsorted"``.
@@ -234,7 +261,7 @@ def generate_eigenvalue_coulomb_matrix(nuclear_charges, coordinates, size = 23):
     return fgenerate_eigenvalue_coulomb_matrix(nuclear_charges,
         coordinates, size)
 
-def generate_bob(nuclear_charges, coordinates, atomtypes, asize = {"O":3, "C":7, "N":3, "H":16, "S":1}):
+def generate_bob(nuclear_charges, coordinates, atomtypes, size=23, asize = {"O":3, "C":7, "N":3, "H":16, "S":1}):
     """ Creates a Bag of Bonds (BOB) representation of a molecule.
         The representation expands on the coulomb matrix representation.
         For each element a bag (vector) is constructed for self interactions
@@ -259,12 +286,45 @@ def generate_bob(nuclear_charges, coordinates, atomtypes, asize = {"O":3, "C":7,
         :type nuclear_charges: numpy array
         :param coordinates: 3D Coordinates of the atoms in the molecule
         :type coordinates: numpy array
+        :param size: The maximum number of atoms in the representation
+        :type size: integer
         :param asize: The maximum number of atoms of each element type supported by the representation
-        :type size: dictionary
+        :type asize: dictionary
 
         :return: 1D representation
         :rtype: numpy array
     """
+    natoms = len(nuclear_charges)
+
+    coulomb_matrix = fgenerate_unsorted_coulomb_matrix(nuclear_charges, coordinates, size)
+
+    coulomb_matrix = vector_to_matrix(coulomb_matrix)
+    descriptor = []
+    atomtypes = np.asarray(atomtypes)
+    for atom1, size1 in sorted(asize.items()):
+        pos1 = np.where(atomtypes == atom1)[0]
+        feature_vector = np.zeros(size1)
+        feature_vector[:pos1.size] = np.diag(coulomb_matrix)[pos1]
+        feature_vector.sort()
+        descriptor.append(feature_vector[:])
+        for atom2, size2 in sorted(asize.items()):
+            if atom1 > atom2:
+                continue
+            if atom1 == atom2:
+                size = size1*(size1-1)//2
+                feature_vector = np.zeros(size)
+                sub_matrix = coulomb_matrix[np.ix_(pos1,pos1)]
+                feature_vector[:pos1.size*(pos1.size-1)//2] = sub_matrix[np.triu_indices(pos1.size, 1)]
+                feature_vector.sort()
+                descriptor.append(feature_vector[:])
+            else:
+                pos2 = np.where(atomtypes == atom2)[0]
+                feature_vector = np.zeros(size1*size2)
+                feature_vector[:pos1.size*pos2.size] = coulomb_matrix[np.ix_(pos1,pos2)].ravel()
+                feature_vector.sort()
+                descriptor.append(feature_vector[:])
+
+    return np.concatenate(descriptor)
 
     n = 0
     atoms = sorted(asize, key=asize.get)