From ea6e13c77d6b4c486cd802e6b61da47098dca800 Mon Sep 17 00:00:00 2001
From: capoe <cp605@cam.ac.uk>
Date: Thu, 14 Sep 2017 15:04:21 +0100
Subject: [PATCH] Graph dMap support, DEBUG option.
---
src/soap/fieldtensor.cpp | 42 +++++++++++++-------
src/soap/soapy/graph.py | 85 ++++++++++++++++++++++++++++++++++++----
2 files changed, 106 insertions(+), 21 deletions(-)
diff --git a/src/soap/fieldtensor.cpp b/src/soap/fieldtensor.cpp
index 56fe120..83e18d2 100644
--- a/src/soap/fieldtensor.cpp
+++ b/src/soap/fieldtensor.cpp
@@ -202,8 +202,10 @@ void AtomicSpectrumFT::contractDeep() {
int lm = l*l+l+m;
phi_l_s1s2 += f1(0, lm)*std::conj(f2(0, lm));
}
- GLOG() << l << ":" << l << std::endl;
- GLOG() << "Save @ " << offset << std::endl;
+ #ifdef DBG
+ GLOG() << l << ":" << l << std::endl;
+ GLOG() << "Save @ " << offset << std::endl;
+ #endif
coeffs(offset, 0) = pow(this->getCenter()->getSigma(), 2*l)*inv_alpha_l*phi_l_s1s2.real();
offset += 1;
} // l
@@ -224,8 +226,10 @@ void AtomicSpectrumFT::contractDeep() {
phi_l1_l1l1 += w111 * f1(0,lm1) * f2(0,lm2) * f2(0,lm3);
}}}
// Store l1l1l2_s1s1s2
- GLOG() << l1 << "::" << l1 << ":" << l1 << " " << w111_sq << std::endl;
- GLOG() << "Save @ " << offset << std::endl;
+ #ifdef DBG
+ GLOG() << l1 << "::" << l1 << ":" << l1 << " " << w111_sq << std::endl;
+ GLOG() << "Save @ " << offset << std::endl;
+ #endif
coeffs(offset, 0) =
sqrt(inv_alpha_l1)*sqrt(inv_alpha_l1)*sqrt(inv_alpha_l1)
* pow(this->getCenter()->getSigma(), l1+l1+l1+0.5) *
@@ -260,8 +264,10 @@ void AtomicSpectrumFT::contractDeep() {
phi_l1_l1l2 += w112 * f1(0,lm1) * f2(0,lm2) * f2(0,lm3);
}}}
// Store l1l1l2_s1s1s2
- GLOG() << l1 << ":" << l1 << "::" << l2 << " " << w112_sq << std::endl;
- GLOG() << "Save @ " << offset << std::endl;
+ #ifdef DBG
+ GLOG() << l1 << ":" << l1 << "::" << l2 << " " << w112_sq << std::endl;
+ GLOG() << "Save @ " << offset << std::endl;
+ #endif
coeffs(offset, 0) =
sqrt(inv_alpha_l1)*sqrt(inv_alpha_l1)*sqrt(inv_alpha_l2)
* pow(this->getCenter()->getSigma(), l1+l1+l2+0.5) *
@@ -269,8 +275,10 @@ void AtomicSpectrumFT::contractDeep() {
* phi_l1l1_l2.real();
offset += 1;
// Store l1l1l2_s1s2s2
- GLOG() << l1 << "::" << l1 << ":" << l2 << " " << w112_sq << std::endl;
- GLOG() << "Save @ " << offset << std::endl;
+ #ifdef DBG
+ GLOG() << l1 << "::" << l1 << ":" << l2 << " " << w112_sq << std::endl;
+ GLOG() << "Save @ " << offset << std::endl;
+ #endif
coeffs(offset, 0) =
sqrt(inv_alpha_l1)*sqrt(inv_alpha_l1)*sqrt(inv_alpha_l2)
* pow(this->getCenter()->getSigma(), l1+l1+l2+0.5) *
@@ -303,8 +311,10 @@ void AtomicSpectrumFT::contractDeep() {
phi_l3_l1l2 += w123 * f1(0,lm3) * f2(0,lm1) * f2(0,lm2);
}}}
// Store (l1)(l2l3)
- GLOG() << l1 << "::" << l2 << ":" << l3 << " " << w123_sq << std::endl;
- GLOG() << "Save @ " << offset << std::endl;
+ #ifdef DBG
+ GLOG() << l1 << "::" << l2 << ":" << l3 << " " << w123_sq << std::endl;
+ GLOG() << "Save @ " << offset << std::endl;
+ #endif
coeffs(offset, 0) =
sqrt(inv_alpha_l1)*sqrt(inv_alpha_l2)*sqrt(inv_alpha_l3)
* pow(this->getCenter()->getSigma(), l1+l2+l3+0.5) *
@@ -312,8 +322,10 @@ void AtomicSpectrumFT::contractDeep() {
* phi_l1_l2l3.real();
offset += 1;
// Store (l2)(l3l1)
- GLOG() << l2 << "::" << l3 << ":" << l1 << " " << w123_sq << std::endl;
- GLOG() << "Save @ " << offset << std::endl;
+ #ifdef DBG
+ GLOG() << l2 << "::" << l3 << ":" << l1 << " " << w123_sq << std::endl;
+ GLOG() << "Save @ " << offset << std::endl;
+ #endif
coeffs(offset, 0) =
sqrt(inv_alpha_l1)*sqrt(inv_alpha_l2)*sqrt(inv_alpha_l3)
* pow(this->getCenter()->getSigma(), l1+l2+l3+0.5) *
@@ -321,8 +333,10 @@ void AtomicSpectrumFT::contractDeep() {
* phi_l2_l3l1.real();
offset += 1;
// Store (l3)(l1l2)
- GLOG() << l3 << "::" << l1 << ":" << l2 << " " << w123_sq << std::endl;
- GLOG() << "Save @ " << offset << std::endl;
+ #ifdef DBG
+ GLOG() << l3 << "::" << l1 << ":" << l2 << " " << w123_sq << std::endl;
+ GLOG() << "Save @ " << offset << std::endl;
+ #endif
coeffs(offset, 0) =
sqrt(inv_alpha_l1)*sqrt(inv_alpha_l2)*sqrt(inv_alpha_l3)
* pow(this->getCenter()->getSigma(), l1+l2+l3+0.5) *
diff --git a/src/soap/soapy/graph.py b/src/soap/soapy/graph.py
index 9f0743a..68b717b 100644
--- a/src/soap/soapy/graph.py
+++ b/src/soap/soapy/graph.py
@@ -19,38 +19,109 @@ class Graph(object):
idx=-1,
label='',
feature_mat=None,
+ feature_mat_avg=None,
position_mat=None,
connectivity_mat=None,
vertex_info=[],
graph_info={}):
+ # Labels
self.idx = idx
self.label = label
+ self.graph_info = graph_info
+ # Vertex data: descriptors
self.P = feature_mat
+ self.P_avg = feature_mat_avg
+ self.P_type_str = str(type(self.P))
+ # Vertex data: positions, labels
self.R = position_mat
- self.C = connectivity_mat
self.vertex_info = vertex_info
- self.graph_info = graph_info
+ # Edge data
+ self.C = connectivity_mat
return
def save_to_h5(self, h5f, dtype='float32'):
group = h5f.create_group('%06d' % self.idx)
- group.create_dataset('feature_mat', data=self.P, compression='gzip', dtype=dtype)
+ if self.P_type_str == "<class 'soap.soapy.kernel.DescriptorMapMatrix'>":
+ # Save list of descriptor maps
+ g0 = group.create_group('feature_dmap')
+ for dmap_idx, dmap in enumerate(self.P):
+ g1 = g0.create_group('%d' % dmap_idx)
+ for key in dmap:
+ g1.create_dataset(
+ key,
+ data=dmap[key],
+ compression='gzip',
+ dtype=dtype)
+ # Save averaged descriptor map
+ g0_avg = group.create_group('feature_dmap_avg')
+ for key in self.P_avg:
+ g0_avg.create_dataset(
+ key,
+ data=self.P_avg[key],
+ compression='gzip',
+ dtype=dtype)
+ elif self.P_type_str == "<type 'numpy.ndarray'>":
+ # Save numpy arrays
+ group.create_dataset(
+ 'feature_mat',
+ data=self.P,
+ compression='gzip',
+ dtype=dtype)
+ group.create_dataset(
+ 'feature_mat_avg',
+ data=self.P_avg,
+ compression='gzip',
+ dtype=dtype)
+ else: raise NotImplementedError(self.P_type_str)
group.create_dataset('position_mat', data=self.R)
group.create_dataset('connectivity_mat', data=self.C)
group.attrs['idx'] = self.idx
group.attrs['label'] = self.label
group.attrs['vertex_info'] = json.dumps(self.vertex_info)
group.attrs['graph_info'] = json.dumps(self.graph_info)
+ group.attrs['P_type_str'] = self.P_type_str
return
def load_from_h5(self, h5f):
self.idx = h5f.attrs['idx']
self.label = h5f.attrs['label']
self.vertex_info = json.loads(h5f.attrs['vertex_info'])
self.graph_info = json.loads(h5f.attrs['graph_info'])
- self.P = h5f['feature_mat'].value
+ # Determine feature matrix type
+ if 'P_type_str'in h5f.attrs:
+ self.P_type_str = h5f.attrs['P_type_str']
+ else:
+ self.P_type_str = "<type 'numpy.ndarray'>"
+ if self.P_type_str == "<class 'soap.soapy.kernel.DescriptorMapMatrix'>":
+ # Load list of descriptor maps
+ self.P = soap.soapy.kernel.DescriptorMapMatrix()
+ g0 = h5f['feature_dmap']
+ for i in range(len(g0)):
+ Pi = soap.soapy.kernel.DescriptorMap()
+ g1 = g0['%d' % i]
+ for key in g1:
+ Pi[key] = g1[key].value
+ self.P.append(Pi)
+ # Load averaged descriptor map
+ self.P_avg = soap.soapy.kernel.DescriptorMap()
+ g0_avg = h5f['feature_dmap_avg']
+ for key in g0_avg:
+ self.P_avg[key] = g0_avg[key].value
+ elif self.P_type_str == "<type 'numpy.ndarray'>":
+ self.P = h5f['feature_mat'].value
+ self.P_avg = h5f['feature_mat_avg'].value
+ else: raise NotImplementedError(self.P_type_str)
self.R = h5f['position_mat'].value
self.C = h5f['connectivity_mat'].value
return self
+def load_graphs(hdf5, n=-1):
+ h = h5py.File(hdf5, 'r')
+ gsec = h['graphs']
+ if n == None or n < 0: n = len(gsec)
+ graphs = [ Graph().load_from_h5(
+ gsec['%06d' % i]) for i in range(len(gsec)) if i < n ]
+ h.close()
+ return graphs
+
def calculate_graph_connectivity(graph, zero_diagonal):
"""
Suited for periodic-table elements:
@@ -297,7 +368,7 @@ def mp_compute_graph(
vertex_info = frag_labels
else:
connectivity_mat = atom_bond_mat
- vertex_info = atom_labels
+ vertex_info = type_vec
graph = Graph(
idx = config.info['idx'],
label = str(config.info['label']),
@@ -332,7 +403,7 @@ def mp_compute_kernel_block_hdf5(block, kernel, log, dtype_result, h5_file):
#soap.soapy.util.MP_LOCK.release()
return mp_compute_kernel_block([g_rows, g_cols], kernel, log, dtype_result)
-def mp_compute_kernel_block(block, kernel, log, dtype_result):
+def mp_compute_kernel_block(block, kernel, log, dtype_result, symmetric=True):
gi_block = block[0]
gj_block = block[1]
if log:
@@ -342,7 +413,7 @@ def mp_compute_kernel_block(block, kernel, log, dtype_result):
kmat = np.zeros((len(gi_block),len(gj_block)), dtype=dtype_result)
for i,gi in enumerate(gi_block):
for j,gj in enumerate(gj_block):
- if gi.idx > gj.idx: pass
+ if symmetric and gi.idx > gj.idx: pass
else: kmat[i,j] = kernel.compute(gi,gj)
return kmat
--
GitLab