From af028a7e08fa3faf4dbe0599d06c9de42af85427 Mon Sep 17 00:00:00 2001
From: Carl Poelking <cp605@cam.ac.uk>
Date: Mon, 27 Mar 2017 18:50:12 +0100
Subject: [PATCH] Gobbledygook.
---
src/soap/fieldtensor.cpp | 50 ++++++++++++++++++++--------
src/soap/fieldtensor.hpp | 22 +++++++++++--
src/soap/soapy/kernel.py | 70 ++++++++++++++++++++++++++++++++++++++--
3 files changed, 123 insertions(+), 19 deletions(-)
diff --git a/src/soap/fieldtensor.cpp b/src/soap/fieldtensor.cpp
index 21e228b..fa9fa4a 100644
--- a/src/soap/fieldtensor.cpp
+++ b/src/soap/fieldtensor.cpp
@@ -1,5 +1,6 @@
#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
diff --git a/src/soap/fieldtensor.hpp b/src/soap/fieldtensor.hpp
index b78f572..87f9121 100644
--- a/src/soap/fieldtensor.hpp
+++ b/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;
};
diff --git a/src/soap/soapy/kernel.py b/src/soap/soapy/kernel.py
index 6c2b317..62c63dd 100644
--- a/src/soap/soapy/kernel.py
+++ b/src/soap/soapy/kernel.py
@@ -16,7 +16,7 @@ from momo import osio, endl, flush
class DescriptorMap(object):
"""
- Structured data type that stores
+ Structured data type that stores
descriptors as key-value pairs
and emulates basic algebraic properties
"""
@@ -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 = {
--
GitLab