Gobbledygook.

src/soap/fieldtensor.cpp

 #include "soap/fieldtensor.hpp"
 #include "soap/functions.hpp"
+#include "soap/linalg/numpy.hpp"
 #include <boost/math/special_functions/legendre.hpp>
 
 namespace soap {
@@ -26,7 +27,7 @@ AtomicSpectrumFT::AtomicSpectrumFT(Particle *center, int K, int L)
 
 AtomicSpectrumFT::~AtomicSpectrumFT() {
     GLOG() << "[~] Destruct " << this->getCenter()->getId() << std::endl;
-    // Deallocate body-order terms
+    // Deallocate body-order field terms
     for (int k=0; k <= _K; ++k) {
         field_map_t &fm = _body_map[k];
         GLOG() << "[~]   Deallocating k=" << k<< std::endl;
@@ -36,6 +37,27 @@ AtomicSpectrumFT::~AtomicSpectrumFT() {
             delete it->second;
         }
     }
+    // Deallocate contraction coefficients
+    for (auto it = _coeff_map.begin(); it != _coeff_map.end(); ++it) {
+        GLOG() << "[~]   Deallocating s1:s2 = " <<  it->first.first << ":" << it->first.second << std::endl;
+        delete it->second;
+    }
+}
+
+boost::python::list AtomicSpectrumFT::getTypes() {
+    boost::python::list types;
+    // The highest-body-order term should have the complete set
+    for (auto it = _body_map[_K].begin(); it != _body_map[_K].end(); ++it) {
+        std::string type = it->first;
+        types.append(type);
+    }
+    return types;
+}
+
+boost::python::object AtomicSpectrumFT::getCoefficientsNumpy(std::string s1, std::string s2) {
+    soap::linalg::numpy_converter npc(_numpy_t.c_str());
+    channel_t channel(s1, s2);
+    return npc.ublas_to_numpy< dtype_t >(*_coeff_map[channel]);
 }
 
 void AtomicSpectrumFT::addField(int k, std::string type, field_t &flm) {
@@ -143,24 +165,22 @@ void AtomicSpectrumFT::contract() {
                             }
                             int l1l2 = LambdaLambda_off_k+lambda1*Lambda1 + lambda2;
                             GLOG() << " Store " << lambda1 << ":" << lambda2 << ":" << l << " @ " << l1l2 << ":" << l << " = " << phi_l_s1s2_l1l2 <<  std::endl;
-                            coeffs(l1l2, l) = inv_alpha*phi_l_s1s2_l1l2;
-                        }
-                    }
-                }
-
-
-            }
-        }
+                            coeffs(l1l2, l) = inv_alpha*phi_l_s1s2_l1l2.real();
+                        } // l
+                    } // lambda 2
+                } // lambda 1
+            } // Channel type 2
+        } // Channel type 1
         GLOG() << std::endl;
-
-
-
     }
+    return;
 }
 
 void AtomicSpectrumFT::registerPython() {
     using namespace boost::python;
     class_<AtomicSpectrumFT, AtomicSpectrumFT*>("AtomicSpectrumFT", init<Particle*, int, int>())
+        .def("getTypes", &AtomicSpectrumFT::getTypes)
+        .def("getPower", &AtomicSpectrumFT::getCoefficientsNumpy)
         .def("getCenter", &AtomicSpectrumFT::getCenter, return_value_policy<reference_existing_object>());
 }
 
@@ -335,6 +355,8 @@ FTSpectrum::FTSpectrum(Structure &structure, Options &options)
     : _structure(&structure), _options(&options) {
     _cutoff = CutoffFunctionOutlet().create(_options->get<std::string>("radialcutoff.type"));
 	_cutoff->configure(*_options);
+    _L = _options->get<int>("fieldtensor.L");
+    _K = _options->get<int>("fieldtensor.K");
     return;
 }
 
@@ -429,8 +451,8 @@ void InteractUpdateSourceTarget(
 void FTSpectrum::compute() {
     GLOG() << "Computing FTSpectrum ..." << std::endl;
 
-    int K = 3; // TODO
-    int L = 3; // TODO
+    int K = _K;
+    int L = _L;
     Tlmlm T12(L);
 
     // CREATE ATOMIC SPECTRA

src/soap/fieldtensor.hpp

@@ -11,7 +11,7 @@ namespace soap {
 class AtomicSpectrumFT
 {
 public:
-
+    typedef double dtype_t;
     static const std::string _numpy_t;
 
     // FIELD MOMENTS
@@ -23,14 +23,16 @@ public:
 
     // CONTRACTIONS
     typedef std::pair<std::string, std::string> channel_t;
-    typedef ub::matrix<cmplx_t> coeff_t; // (k=1:0 k=1:l' k=2:l'l'', l)
-    typedef ub::zero_matrix<cmplx_t> coeff_zero_t;
+    typedef ub::matrix<dtype_t> coeff_t; // (k=1:0 k=1:l' k=2:l'l'', l)
+    typedef ub::zero_matrix<dtype_t> coeff_zero_t;
     typedef std::map<channel_t, coeff_t*> coeff_map_t;
 
     AtomicSpectrumFT(Particle *center, int K, int L);
    ~AtomicSpectrumFT();
     Particle *getCenter() { return _center; }
     std::string getType() { return _type; }
+    boost::python::list getTypes();
+    boost::python::object getCoefficientsNumpy(std::string s1, std::string s2);
     int getTypeIdx() { return _s; }
 
     static void registerPython();
@@ -52,6 +54,17 @@ private:
     int _L; // <- Angular momentum cutoff
     int _s; // <- Type index
     std::string _type; // <- Type string
+    std::list<std::string> _nb_types;
+};
+
+class FTBasis // TODO
+{
+public:
+    FTBasis(Options &options) : _K(-1), _L(-1) {};
+   ~FTBasis() {};
+private:
+    int _K;
+    int _L;
 };
 
 class FTSpectrum
@@ -74,6 +87,9 @@ private:
     Options *_options;
     CutoffFunction *_cutoff;
 
+    int _K;
+    int _L;
+
     atomic_array_t _atomic_array;
 };
 

src/soap/soapy/kernel.py

@@ -104,6 +104,71 @@ class TrajectoryLogger(object):
         self.ofs.close()
         return
 
+class KernelAdaptorFTD(object):
+    def __init__(self, options, types_global):
+        self.types = types_global
+        self.S = len(types_global)
+        return
+    def adapt(self, spectrum, return_pos_matrix=False):
+        IX = np.zeros((0,0), dtype='float64') # feature matrix
+        dimX = -1
+        IR = np.zeros((0,0), dtype='float64') # position matrix
+        types = []
+        for atomic_i in spectrum:
+            Xi_unnorm, Xi_norm = self.adaptScalar(atomic_i)
+            Ri = atomic_i.getCenter().pos
+            types.append(atomic_i.getCenter().type)
+            dimX = Xi_norm.shape[0]
+            if not IX.any():
+                IX = np.copy(Xi_norm) # TODO Is this necessary?
+                IX.resize((1, dimX))
+                IR = np.copy(Ri)
+                IR.resize((1, 3))
+            else:
+                i = IX.shape[0]
+                IX.resize((i+1, dimX))
+                IX[-1,:] = Xi_norm
+                IR.resize((i+1, 3))
+                IR[-1,:] = Ri
+        if return_pos_matrix:
+            return IX, IR, types
+        else:
+            return IX
+    def adaptScalar(self, atomic, epsilon=1e-20):
+        X = self.reduce(atomic, self.types)
+        X_mag = np.dot(X,X)**0.5
+        if X_mag < epsilon:
+            X_mag = 1.
+        X_norm = X/X_mag
+        return X, X_norm
+    def reduce(self, atomic, types_global, verbose=False):
+        types_atomic = atomic.getTypes()
+        S = len(types_global)
+        SS = S*S
+        S_atomic = len(types_atomic)
+        X = None
+        dim_ab = -1
+        dim_total = -1
+        # Channel - type a
+        for i in range(S_atomic):
+            a = types_atomic[i]
+            sa = types_global.index(a)
+            # Channel - type b
+            for j in range(S_atomic):
+                b = types_atomic[j]
+                sb = types_global.index(b)
+                x = atomic.getPower(a,b)
+                # Initialize aggregated array
+                if i == 0 and j == 0:
+                    dim_ab = x.shape[0]*x.shape[1]
+                    dim_total = SS*dim_ab
+                    X = np.zeros((dim_total,), dtype='float64')
+                # Locate in super-array X
+                i0 = (sa*S+sb)*dim_ab
+                i1 = i0 + dim_ab
+                X[i0:i1] = x.flatten()
+        return X
+
 class Xnklab(object):
     def __init__(self, atomic, types_global):
         self.types_global = types_global
@@ -683,7 +748,8 @@ KernelAdaptorFactory = {
 'specific-unique-dmap': KernelAdaptorSpecificUniqueDMap,
 'global-generic': KernelAdaptorGlobalGeneric,
 'global-specific': KernelAdaptorGlobalSpecific,
-'global-specific-energy': KernelAdaptorGlobalSpecificEnergy
+'global-specific-energy': KernelAdaptorGlobalSpecificEnergy,
+'ftd-specific': KernelAdaptorFTD
 }
 
 KernelFunctionFactory = {