Support for key-value descriptors.

src/soap/linalg/bindings.cpp

 #include "soap/linalg/bindings.hpp"
+#include "soap/linalg/kernel.hpp"
 
 namespace soap { namespace linalg {
 
@@ -10,5 +11,6 @@ BOOST_PYTHON_MODULE(_linalg) {
     boost::python::numeric::array::set_module_and_type("numpy", "ndarray");
     soap::linalg::vec::registerPython();
     soap::linalg::matrix::registerPython();
+    soap::linalg::KernelModule::registerPython();
 }
 

src/soap/linalg/permanent/CMakeLists.txt

@@ -9,7 +9,7 @@ endforeach()
 # ADD & LINK LIBRARY
 add_library(_permanent ${LOCAL_SOURCES})
 target_compile_options(_permanent PRIVATE "-O3" PRIVATE "-Ofast" PRIVATE "-march=native")
-target_link_libraries(_permanent ${PYTHON_LIBRARIES})
+target_link_libraries(_permanent ${Boost_LIBRARIES} ${PYTHON_LIBRARIES})
 set_target_properties(_permanent PROPERTIES PREFIX "" SUFFIX ".so" LIBRARY_OUTPUT_DIRECTORY .)
 
 # INSTALL

src/soap/linalg/permanent/permanent.cpp

@@ -12,6 +12,7 @@
 #include <random>
 #include <cmath>
 #include <iostream>
+#include "../numpy.hpp"
 // Array access macros.
 
 #define SM(x0, x1) (*(npy_double*) (( (char*) PyArray_DATA(matrix) + \
@@ -19,6 +20,8 @@
                     (x1) * PyArray_STRIDES(matrix)[1])))
 #define SM_shape(x0) (int) PyArray_DIM(matrix, x0)
 
+namespace ub = boost::numeric::ublas;
+
 template <typename T>
 class Matrix
 {
@@ -244,3 +247,14 @@ static PyObject *regmatch(PyObject *self, PyObject *args) {
   npy_double p = _shmatch(matrix,gamma,eps);
   return PyFloat_FromDouble(p);
 }
+
+char const* greet()
+{
+   return "hello, world";
+}
+
+BOOST_PYTHON_MODULE(_reboost)
+{
+    using namespace boost::python;
+    def("greet", greet);
+}

src/soap/soapy/dimred.py

@@ -4,9 +4,12 @@ import sklearn.manifold
 import sklearn.decomposition
 
 def dimred_matrix(method, kmat=None, distmat=None, outfile=None, ix=None, symmetrize=False, prj_dimension=2):
+    dmat = distmat
     if symmetrize:
-        kmat = 0.5*(kmat+kmat.T)
-        dmat = 0.5*(dmat+dmat.T)
+        if type(kmat) != type(None):
+            kmat = 0.5*(kmat+kmat.T)
+        if type(dmat) != type(None):
+            dmat = 0.5*(dmat+dmat.T)
     if method == 'mds':
         # MDS
         # http://scikit-learn.org/stable/modules/generated/sklearn.manifold.MDS.html#sklearn.manifold.MDS

src/soap/soapy/kernel.py

@@ -14,6 +14,82 @@ from momo import osio, endl, flush
 # TODO PCA + identification of eigenstructures
 # TODO Sample Bethe tree
 
+class DescriptorMap(object):
+    """
+    Structured data type that stores 
+    descriptors as key-value pairs
+    and emulates basic algebraic properties
+    """
+    def __init__(self):
+        self.dmap = {}
+        self.T = self # Transpose
+    def dot(self, other):
+        dot = 0.0
+        for key in self.dmap:
+            if key in other.dmap:
+                dot += self[key].dot(other[key])
+        return dot
+    def add(self, other):
+        for key in other:
+            if key in self:
+                self.dmap[key] = self.dmap[key] + other[key]
+            else:
+                self.dmap[key] = other[key]
+        return
+    def normalise(self):
+        norm = self.dot(self)**0.5
+        for key in self.dmap:
+            self.dmap[key] = self.dmap[key] / norm
+        return
+    def __div__(self, scalar):
+        for k in self:
+            self[k] = self[k]/scalar
+        return self
+    def __mul__(self, scalar):
+        for k in self:
+            self[k] = self[k]*scalar
+        return self
+    def __getitem__(self, key):
+        return self.dmap[key]
+    def __setitem__(self, key, value):
+        self.dmap[key] = value
+    def __iter__(self):
+        return iter(self.dmap)
+    def __len__(self):
+        return len(self.dmap)
+    def __contains__(self, key):
+        return key in self.dmap
+
+class DescriptorMapMatrix(object):
+    """
+    List of <DescriptorMap>s
+    """
+    def __init__(self):
+        self.mat = []
+        self.T = self # Transpose
+    def append(self, other):
+        self.mat.append(other)
+    def dot(self, other, dtype='float64'):
+        n = len(self)
+        m = len(other)
+        k = np.zeros((n,m), dtype=dtype)
+        for i in range(n):
+            for j in range(m):
+                k[i,j] = self[i].dot(other[j])
+        return k
+    def sum(self, axis=0):
+        assert axis == 0
+        dmap_summed = soap.soapy.DescriptorMap()
+        for dmap in self:
+            dmap_summed.add(dmap)
+        return dmap_summed
+    def __len__(self):
+        return len(self.mat)
+    def __getitem__(self, i):
+        return self.mat[i]
+    def __iter__(self):
+        return iter(self.mat)
+
 class TrajectoryLogger(object):
     def __init__(self, outfile, mode='w'):
         self.ofs = open(outfile, mode)
@@ -152,6 +228,53 @@ def reduce_xnklab_atomic(atomic, types_global, verbose=False):
             #print a,b,sa,sb,i0,i1
     return X
 
+class KernelAdaptorSpecificUniqueDMap(object):
+    def __init__(self, options, types_global):
+        return
+    def adapt(self, spectrum, return_pos_matrix=False):
+        IDMap = DescriptorMapMatrix() # List of <DescriptorMaps>
+        dimX = -1
+        IR = np.zeros((len(spectrum),3), dtype='float64') # position matrix
+        types = []
+        for idx, atomic_i in enumerate(spectrum):
+            IDMap.append(self.adaptScalar(atomic_i))
+            Ri = atomic_i.getCenter().pos
+            types.append(atomic_i.getCenter().type)
+            IR[idx,:] = Ri
+        if return_pos_matrix:
+            return IDMap, IR, types
+        else:
+            return IDMap
+    def adaptScalar(self, atomic, epsilon=1e-20):
+        dmap = DescriptorMap()
+        types_atomic = atomic.getTypes()
+        S = len(types_atomic)
+        N = atomic.basis.N
+        L = atomic.basis.L
+        # SPECIES a = b
+        for i in range(S):
+            a = types_atomic[i]
+            xnklaa = atomic.getPower(a,a).array.real
+            xnklaa_red = np.zeros(((N*N+N)/2, L+1))
+            # Select where n <= k
+            for i in range(N):
+                for j in range(i, N):
+                    ij_red = i*N - (i*i-i)/2 + j-i
+                    xnklaa_red[ij_red] = xnklaa[i*N+j]
+            dmap['%s:%s' % (a,a)] = xnklaa_red.flatten()
+        # SPECIES a != b
+        for i in range(S):
+            for j in range(i+1, S):
+                # Select all
+                a = types_atomic[i]
+                b = types_atomic[j]
+                xnklab = atomic.getPower(a,b).array.real
+                xnklab_red = xnklab
+                pair_ab = '%s:%s' % ((a,b) if a < b else (b,a))
+                dmap[pair_ab] = xnklab_red.flatten()
+        dmap.normalise()
+        return dmap
+
 class KernelAdaptorSpecificUnique(object):
     """
     Extracts Xnklab from spectrum, while removing duplicates
@@ -557,6 +680,7 @@ KernelAdaptorFactory = {
 'generic': KernelAdaptorGeneric,
 'specific': KernelAdaptorSpecific,
 'specific-unique': KernelAdaptorSpecificUnique,
+'specific-unique-dmap': KernelAdaptorSpecificUniqueDMap,
 'global-generic': KernelAdaptorGlobalGeneric,
 'global-specific': KernelAdaptorGlobalSpecific,
 'global-specific-energy': KernelAdaptorGlobalSpecificEnergy

src/soap/soapy/lamatch.py

@@ -158,4 +158,35 @@ GraphKernelFactory = {
 'global': compare_graphs_global,
 'laplacian': lagraph.compare_graphs_laplacian_kernelized
 }
+
+# =====================
+# NOMAD Structure Tools
+# =====================
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
 

src/soap/soapy/pca.py

@@ -17,6 +17,7 @@ def pca_compute(IX, eps=0., log=None, norm_div_std=True, norm_sub_mean=True):
     """
     # Normalize: mean, std
     if log: log << "PCA: Normalize ..." << log.endl
+    IX_norm = IX
     if norm_sub_mean:
         IX_norm = IX - IX.mean(0)
     if norm_div_std:

src/soap/spectrum.cpp

@@ -107,9 +107,9 @@ AtomicSpectrum *Spectrum::computeAtomic(Particle *center, Structure::particle_ar
     int nb_max = na_nb_nc[1];
     int nc_max = na_nb_nc[2];
 
-    GLOG() << box_a << " " << box_b << " " << box_c << std::endl;
-    GLOG() << rc << std::endl;
-    GLOG() << na_max << " " << nb_max << " " << nc_max << std::endl;
+    //GLOG() << box_a << " " << box_b << " " << box_c << std::endl;
+    //GLOG() << rc << std::endl;
+    //GLOG() << na_max << " " << nb_max << " " << nc_max << std::endl;
 
     // CREATE BLANK
     AtomicSpectrum *atomic_spectrum = new AtomicSpectrum(center, this->_basis);
@@ -299,6 +299,7 @@ void Spectrum::registerPython() {
     	.def(init<Structure &, Options &, Basis &>())
     	.def(init<std::string>())
     	.def("__iter__", range<return_value_policy<reference_existing_object> >(&Spectrum::beginAtomic, &Spectrum::endAtomic))
+        .def("__len__", &Spectrum::length)
 	    .def("compute", computeAll)
         .def("compute", computeSeg)
 	    .def("compute", computeSegPair)

src/soap/spectrum.hpp

@@ -44,6 +44,7 @@ public:
 	void save(std::string archfile);
 	void load(std::string archfile);
 	void clean();
+    int length() { return _atomspec_array.size(); }
 
 	void compute();
     void compute(Segment *centers);