Skip to content
GitLab
Commit 4bc8634d authored by Carl Poelking's avatar Carl Poelking
Browse files

PowerSpectrum wrapper, supporting descriptor maps.

parent 81ec4f49
Hide whitespace changes
Inline Side-by-side
src/soap/atomicspectrum.cpp
View file @ 4bc8634d
......@@ -47,6 +47,22 @@ AtomicSpectrum::~AtomicSpectrum() {
delete _xnkl_generic_incoherent;
_xnkl_generic_incoherent = NULL;
}
// CLEAN PID-RESOLVED GRADIENTS
this->prunePidData();
}
void AtomicSpectrum::null() {
_center = NULL;
_center_id = -1;
_center_pos = vec(0,0,0);
_center_type = "?";
_basis = NULL;
_qnlm_generic = NULL;
_xnkl_generic_coherent = NULL;
_xnkl_generic_incoherent = NULL;
}
void AtomicSpectrum::prunePidData() {
// CLEAN PID-RESOLVED GRADIENTS
// ... Neighbour density expansions with gradients
for (auto it = _map_pid_qnlm.begin(); it != _map_pid_qnlm.end(); ++it) delete it->second.second;
......@@ -65,17 +81,6 @@ AtomicSpectrum::~AtomicSpectrum() {
}
}
void AtomicSpectrum::null() {
_center = NULL;
_center_id = -1;
_center_pos = vec(0,0,0);
_center_type = "?";
_basis = NULL;
_qnlm_generic = NULL;
_xnkl_generic_coherent = NULL;
_xnkl_generic_incoherent = NULL;
}
void AtomicSpectrum::invert(map_xnkl_t &map_xnkl, xnkl_t *xnkl_generic_coherent, std::string type1, std::string type2) {
int N = _basis->getRadBasis()->N();
......
src/soap/atomicspectrum.hpp
View file @ 4bc8634d
......@@ -36,12 +36,12 @@ public:
typedef std::map<int, map_xnkl_t> map_pid_xnkl_t; // <- id=>type=>xnkl
typedef std::map<int, xnkl_t*> map_pid_xnkl_gc_t; // <- id=>xnkl_generic_coherent
AtomicSpectrum(Particle *center, Basis *basis);
AtomicSpectrum(Basis *basis);
AtomicSpectrum() { this->null(); }
~AtomicSpectrum();
void null();
void prunePidData();
void write(std::ostream &ofs);
void invert(map_xnkl_t &map_xnkl, xnkl_t *xnkl_generic_coherent, std::string type1, std::string type2);
void invert(xnkl_t *xnkl, std::string type1, std::string type2);
......
src/soap/soapy/graph.py
View file @ 4bc8634d
......@@ -52,13 +52,14 @@ class Graph(object):
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)
if type(self.P_avg) != type(None):
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(
......@@ -103,9 +104,10 @@ class Graph(object):
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
if 'feature_dmap_avg' in h5f:
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
if 'feature_mat_avg' in h5f:
......
src/soap/soapy/kernel.py
View file @ 4bc8634d
......@@ -22,7 +22,10 @@ class DescriptorMap(object):
"""
def __init__(self):
self.dmap = {}
self.T = self # Transpose
#self.T = self # Transpose
@property
def T(self):
return self
def dot(self, other):
dot = 0.0
for key in self.dmap:
......@@ -294,7 +297,7 @@ def reduce_xnklab_atomic(atomic, types_global, verbose=False):
return X
class KernelAdaptorSpecificUniqueDMap(object):
def __init__(self, options, types_global):
def __init__(self, options, types_global=None):
return
def adapt(self, spectrum, return_pos_matrix=False):
IDMap = DescriptorMapMatrix() # List of <DescriptorMaps>
......@@ -469,6 +472,53 @@ class KernelAdaptorGlobalSpecific(object):
X_norm = X/np.dot(X,X)**0.5
return X, X_norm
class KernelAdaptorGlobalSpecificDMap(object):
def __init__(self, options, types_global=None):
return
def adapt(self, spectrum, return_pos_matrix=False):
IDMap = DescriptorMapMatrix() # List of <DescriptorMaps>
dimX = -1
IR = np.zeros((1,3), dtype='float64') # position matrix
types = []
atomic_i = spectrum.getGlobal()
IDMap.append(self.adaptScalar(atomic_i))
Ri = np.array([0., 0., 0.])
types.append("global")
IR[0,:] = 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 KernelAdaptorGeneric(object):
def __init__(self, options, types_global=None):
return
......@@ -742,25 +792,26 @@ class KernelFunctionDot3HarmonicDist(object):
return 2*(D_A - D_B)*(0.5/D_A*(-1.)*IC_A - 0.5/D_B*(-1.)*IC_B) + 2*(D_A - D_AB)*0.5/D_A*(-1.)*IC_A + 2*(D_B - D_AB)*0.5/D_B*(-1.)*IC_B + 2*(D_A-self.d0)*1./D_A*(-1.)*IC_A + 2*(D_B-self.d0)*1./D_B*(-1.)*IC_B
KernelAdaptorFactory = {
'generic': KernelAdaptorGeneric,
'specific': KernelAdaptorSpecific,
'specific-unique': KernelAdaptorSpecificUnique,
'specific-unique-dmap': KernelAdaptorSpecificUniqueDMap,
'global-generic': KernelAdaptorGlobalGeneric,
'global-specific': KernelAdaptorGlobalSpecific,
'global-specific-energy': KernelAdaptorGlobalSpecificEnergy,
'ftd-specific': KernelAdaptorFTD
'generic': KernelAdaptorGeneric,
'specific': KernelAdaptorSpecific,
'specific-unique': KernelAdaptorSpecificUnique,
'specific-unique-dmap': KernelAdaptorSpecificUniqueDMap,
'global-generic': KernelAdaptorGlobalGeneric,
'global-specific': KernelAdaptorGlobalSpecific,
'global-specific-dmap': KernelAdaptorGlobalSpecificDMap,
'global-specific-energy': KernelAdaptorGlobalSpecificEnergy,
'ftd-specific': KernelAdaptorFTD
}
KernelFunctionFactory = {
'gaussian-diff': KernelFunctionGaussianDiff,
'dot': KernelFunctionDot,
'dot-shifted': KernelFunctionDotShifted,
'dot-harmonic': KernelFunctionDotHarmonic,
'dot-harmonic-dist': KernelFunctionDotHarmonicDist,
'dot-lj': KernelFunctionDotLj,
'dot-3-harmonic': KernelFunctionDot3Harmonic,
'dot-3-harmonic-dist': KernelFunctionDot3HarmonicDist
'gaussian-diff': KernelFunctionGaussianDiff,
'dot': KernelFunctionDot,
'dot-shifted': KernelFunctionDotShifted,
'dot-harmonic': KernelFunctionDotHarmonic,
'dot-harmonic-dist': KernelFunctionDotHarmonicDist,
'dot-lj': KernelFunctionDotLj,
'dot-3-harmonic': KernelFunctionDot3Harmonic,
'dot-3-harmonic-dist': KernelFunctionDot3HarmonicDist
}
class KernelPotential(object):
......
src/soap/spectrum.cpp
View file @ 4bc8634d
......@@ -24,6 +24,11 @@ Spectrum::Spectrum(std::string archfile) :
this->load(archfile);
}
Spectrum::Spectrum() :
_log(NULL), _options(NULL), _structure(NULL), _basis(NULL), _own_basis(true), _global_atomic(NULL) {
;
}
Spectrum::~Spectrum() {
delete _log;
_log = NULL;
......@@ -281,6 +286,25 @@ void Spectrum::save(std::string archfile) {
return;
}
std::string Spectrum::saves(bool prune) {
if (prune) {
// Delete all pid-resolved data from atomic spectra.
// Note that this data is required to compute gradients.
// Pruning will no longer allow computing gradients
// for serialised objects.
for (auto it = _atomspec_array.begin(); it != _atomspec_array.end(); ++it) {
(*it)->prunePidData();
}
if (_global_atomic) {
_global_atomic->prunePidData();
}
}
std::stringstream bstream;
boost::archive::binary_oarchive arch(bstream);
arch << (*this);
return bstream.str();
}
void Spectrum::load(std::string archfile) {
std::ifstream ifs(archfile.c_str());
boost::archive::binary_iarchive arch(ifs);
......@@ -288,6 +312,14 @@ void Spectrum::load(std::string archfile) {
return;
}
Spectrum &Spectrum::loads(std::string bstr) {
std::stringstream bstream;
bstream << bstr;
boost::archive::binary_iarchive arch(bstream);
arch >> (*this);
return (*this);
}
void Spectrum::registerPython() {
using namespace boost::python;
void (Spectrum::*computeAll)() = &Spectrum::compute;
......@@ -298,6 +330,7 @@ void Spectrum::registerPython() {
class_<Spectrum>("Spectrum", init<Structure &, Options &>())
.def(init<Structure &, Options &, Basis &>())
.def(init<std::string>())
.def(init<>())
.def("__iter__", range<return_value_policy<reference_existing_object> >(&Spectrum::beginAtomic, &Spectrum::endAtomic))
.def("__len__", &Spectrum::length)
.def("compute", computeAll)
......@@ -306,6 +339,7 @@ void Spectrum::registerPython() {
.def("compute", computeCentersTargets)
.def("computePower", &Spectrum::computePower)
.def("computePowerGradients", &Spectrum::computePowerGradients)
.def("deleteGlobal", &Spectrum::deleteGlobal)
.def("computeGlobal", &Spectrum::computeGlobal, return_value_policy<reference_existing_object>())
.def("addAtomic", &Spectrum::addAtomic)
.def("getAtomic", &Spectrum::getAtomic, return_value_policy<reference_existing_object>())
......@@ -313,6 +347,8 @@ void Spectrum::registerPython() {
.def("saveAndClean", &Spectrum::saveAndClean)
.def("save", &Spectrum::save)
.def("load", &Spectrum::load)
.def("saves", &Spectrum::saves)
.def("loads", &Spectrum::loads, ref_existing())
.def("writeDensityOnGrid", &Spectrum::writeDensityOnGrid)
.def("writeDensityCubeFile", &Spectrum::writeDensityCubeFile)
.def("writeDensity", &Spectrum::writeDensity)
......
src/soap/spectrum.hpp
View file @ 4bc8634d
......@@ -28,6 +28,7 @@ public:
typedef std::vector<AtomicSpectrum*>::iterator atomic_it_t;
typedef std::map<std::string, atomspec_array_t> map_atomspec_array_t;
Spectrum();
Spectrum(std::string archfile);
Spectrum(Structure &structure, Options &options);
Spectrum(Structure &structure, Options &options, Basis &basis);
......@@ -42,7 +43,9 @@ public:
void saveAndClean() { std::cout << "spectrum::save&clean" << std::endl; }
void save(std::string archfile);
std::string saves(bool prune = true);
void load(std::string archfile);
Spectrum &loads(std::string bstr);
void clean();
int length() { return _atomspec_array.size(); }
......@@ -52,7 +55,8 @@ public:
void compute(Structure::particle_array_t &sources, Structure::particle_array_t &targets);
AtomicSpectrum *computeAtomic(Particle *center);
AtomicSpectrum *computeAtomic(Particle *center, Structure::particle_array_t &targets);
AtomicSpectrum *computeGlobal();
AtomicSpectrum *computeGlobal();
void deleteGlobal() { if (_global_atomic) { delete _global_atomic; _global_atomic = NULL; } }
void addAtomic(AtomicSpectrum *atomspec);
AtomicSpectrum *getAtomic(int slot_idx, std::string center_type);
AtomicSpectrum *getGlobal() { assert(_global_atomic && "Compute first"); return _global_atomic; }
......
src/soap/tools/loadwrite.py
View file @ 4bc8634d
......@@ -51,7 +51,8 @@ def structure_from_ase(
frag_labels = []
atom_labels = []
# System properties
label = str(config.info['label'])
if 'label' in config.info: label = str(config.info['label'])
else: label = "?"
positions = config.get_positions()
types = config.get_chemical_symbols()
if log: log << log.back << "Reading '%s'" % label << log.flush
......
test/test_nomadsim/data/structures.xyz 0 → 100644
View file @ 4bc8634d
83
label=amk
O -0.2481000 -1.5435000 -0.0318000
O -1.0143000 +1.2262000 -0.0189000
O +1.6187000 +0.4791000 +2.3746000
O -0.7154000 -3.6608000 +0.9232000
O -1.9669000 +3.0045000 +1.2262000
O -0.3871000 -2.6479000 -2.5716000
O -3.5554000 +0.0746000 -0.3225000
O -4.3071000 -3.3663000 +0.1595000
O -5.0786000 +2.2906000 -1.3926000
O -3.4770000 +4.6550000 -1.7036000
O +4.4438000 -1.5561000 -1.7036000
O -1.6527000 -3.7443000 +3.4833000
O +6.8243000 -1.3551000 -0.5031000
N +3.3036000 -0.6804000 +0.1187000
N +2.5248000 +2.2729000 +0.1922000
N -3.1638000 -3.6565000 -2.4905000
N -0.5372000 +5.4195000 +1.3645000
N +6.2986000 +2.8934000 -0.1960000
C +0.9202000 -1.0504000 +0.6256000
C +0.2427000 +1.4238000 +0.6315000
C +1.9786000 -0.7591000 -0.4492000
C +0.5548000 +0.1783000 +1.4760000
C +1.3141000 +1.7359000 -0.4302000
C +1.6412000 +0.4919000 -1.2745000
C -0.1857000 -2.9586000 -0.2015000
C -2.0964000 +1.6431000 +0.8116000
C -0.9272000 -3.3750000 -1.4736000
C -2.4343000 -3.1227000 -1.3342000
C -3.4606000 +1.4101000 +0.1602000
C -2.9594000 -3.7821000 -0.0481000
C -3.7099000 +2.3830000 -0.9905000
C -2.1001000 -3.3866000 +1.1635000
C -3.4360000 +3.8217000 -0.5476000
C -2.0719000 +3.9387000 +0.1457000
C +4.4357000 -1.0908000 -0.5659000
C -2.4938000 -4.1564000 +2.4216000
C -1.8458000 +5.3284000 +0.7348000
C +5.7132000 -0.8806000 +0.2389000
C +5.9422000 +0.5901000 +0.6082000
C +6.0680000 +1.5124000 -0.6101000
H +1.3060000 -1.8138000 +1.3174000
H +0.2011000 +2.2954000 +1.2966000
H +2.0009000 -1.5892000 -1.1650000
H -0.3098000 -0.0930000 +2.0906000
H +0.9222000 +2.5018000 -1.1112000
H +0.7665000 +0.2667000 -1.8997000
H +2.4566000 +0.7186000 -1.9739000
H +0.8563000 -3.2821000 -0.3155000
H -2.0762000 +1.0334000 +1.7214000
H +3.3970000 -0.2830000 +1.0483000
H -0.7599000 -4.4451000 -1.6543000
H -2.6251000 -2.0454000 -1.3213000
H -4.2503000 +1.5445000 +0.9105000
H +2.2996000 +3.1249000 +0.7040000
H +3.1901000 +2.5414000 -0.5312000
H -2.9521000 -4.8736000 -0.1628000
H -3.1116000 +2.1137000 -1.8688000
H -2.2084000 -2.3168000 +1.3844000
H +1.3290000 +1.2160000 +2.9388000
H -4.2341000 +4.1580000 +0.1255000
H -1.2718000 +3.7447000 -0.5813000
H -3.5346000 -3.9635000 +2.6986000
H -2.3500000 -5.2328000 +2.2804000
H -1.9200000 +6.0939000 -0.0444000
H -2.6075000 +5.5441000 +1.4931000
H -4.1545000 -3.4350000 -2.3998000
H -2.8396000 -3.1950000 -3.3394000
H +5.6351000 -1.4860000 +1.1497000
H +0.5562000 -2.8739000 -2.6384000
H -3.3704000 -0.5141000 +0.4283000
H -4.3113000 -2.4018000 +0.2774000
H -5.6353000 +2.5131000 -0.6270000
H -2.7520000 +4.3864000 -2.2934000
H +6.8726000 +0.6593000 +1.1877000
H +5.1400000 +0.9504000 +1.2624000
H -0.4631000 +4.7238000 +2.1057000
H -0.4383000 +6.3289000 +1.8135000
H -1.9182000 -4.2477000 +4.2718000
H +7.5985000 -1.3396000 +0.0851000
H +5.1634000 +1.4714000 -1.2245000
H +6.9062000 +1.1938000 -1.2386000
H +6.3280000 +3.4922000 -1.0202000
H +5.5117000 +3.2180000 +0.3648000
45
label=cef
S -0.5586000 +0.8986000 -2.0988000
S +6.0436000 -1.5788000 -0.4902000
O +1.3241000 +3.3320000 +0.4702000
O +0.9713000 +0.7594000 +2.2896000
O -4.8654000 -0.3485000 -0.7008000
O -2.0761000 +0.1273000 +2.9040000
O +3.5679000 +1.1153000 +0.9287000
O -1.7889000 -1.9676000 +2.0326000
O -6.2321000 -2.1827000 -0.4408000
N -0.6238000 +0.8095000 +0.5539000
N +2.0028000 +1.3479000 -0.7677000
N -7.0943000 -0.1389000 -1.0458000
C +0.9670000 +2.0023000 +0.0430000
C -0.4526000 +1.7713000 -0.5197000
C +0.5039000 +1.1137000 +1.2326000
C -1.6782000 -0.0412000 +0.6443000
C -2.5005000 -0.2626000 -0.4089000
C -2.2661000 +0.3298000 -1.7968000
C -3.7371000 -1.1278000 -0.3645000
C +3.2287000 +0.9577000 -0.2430000
C -1.8484000 -0.7600000 +1.9156000
C +1.3280000 +4.3199000 -0.5493000
C +4.1325000 +0.2877000 -1.2584000
C +4.4254000 -1.1154000 -0.8049000
C +3.5405000 -2.1493000 -0.5909000
C -6.0766000 -0.9953000 -0.7055000
C +4.1998000 -3.3375000 -0.1627000
C +5.5633000 -3.1667000 -0.0658000
H -1.1285000 +2.6293000 -0.4236000
H +1.8295000 +1.1804000 -1.7556000
H -2.4760000 -0.4239000 -2.5632000
H -2.9494000 +1.1708000 -1.9576000
H -3.6261000 -1.9520000 -1.0801000
H -3.8992000 -1.5652000 +0.6253000
H +1.6717000 +5.2618000 -0.1136000
H +0.3200000 +4.4684000 -0.9454000
H +2.0110000 +4.0403000 -1.3564000
H +3.6652000 +0.2557000 -2.2488000
H +5.0394000 +0.8964000 -1.3484000
H -2.1959000 -0.3111000 +3.7733000
H +2.4708000 -2.0626000 -0.7340000
H +3.6994000 -4.2703000 +0.0627000
H +6.3002000 -3.8984000 +0.2338000
H -8.0567000 -0.4555000 -1.0970000
H -6.9150000 +0.8370000 -1.2568000
test/test_nomadsim/run.py 0 → 100755
View file @ 4bc8634d
#! /usr/bin/env python
import soap
import ase.io
import json
import numpy as np
import pickle
import h5py
import time
import datetime
import os
import psutil
log = soap.soapy.momo.osio
options_default = {
"spectrum.gradients": False,
"spectrum.global": False,
"spectrum.2l1_norm": False, # TODO Tune this
"radialbasis.type" : "gaussian",
"radialbasis.mode" : "adaptive", # TODO 'equispaced' may be better suited
"radialbasis.N" : 9, # TODO Tune this
"radialbasis.sigma": 0.5, # TODO Tune this
"radialbasis.integration_steps": 15,
"radialcutoff.Rc": 3.5, # TODO Tune this
"radialcutoff.Rc_width": 0.5,
"radialcutoff.type": "shifted-cosine",
"radialcutoff.center_weight": 1.0,
"angularbasis.type": "spherical-harmonic",
"angularbasis.L": 6, # TODO Tune this
"kernel.adaptor": "specific-unique-dmap",
"exclude_centers": [],
"exclude_targets": [],
"type_list": ["C", "H", "N", "O", "S"]
}
class PowerSpectrum(object):
# TODO Need option for normalisation of gsdmap, sdmap
# ----
# TODO Certain header settings will throw errors:
# TODO e.g., store_gcmap = False and store_gsd, store_gsdmap = True
# ----
# TODO Global descriptor computed twice if options["spectrum.global"]
# TODO and store_gcmap both true
# ----
# NOTE Output hdf5 becomes very large with options["spectrum.gradients"] true
settings = {
'cxx_compute_power' : True,
'store_cxx_serial' : False,
'store_cmap' : True,
'store_gcmap' : True,
'store_sd' : False,
'store_gsd' : False,
'store_sdmap' : True,
'store_gsdmap' : True,
'dtype': 'float32' # TODO implement
}
verbose = False
log = log
def __init__(self, config=None, options=None, label="?"):
if type(config) != type(None):
if self.verbose: self.log << self.log.mg << "Initialising power spectrum '%s'" % label << self.log.endl
self.label = label
self.config = config
self.options = options
# Cxx structure and spectrum
self.has_structure = False
self.structure = None
self.has_spectrum = False
self.spectrum = None
self.has_spectrum_global = False
self.spectrum_global = None
# Density expansion coefficients
self.has_cnlm = False
self.cmap = None
self.gcmap = None
# Power expansion coefficients
self.has_pnkl = False
self.sd = None
self.gsd = None
self.gsdmap = None
self.sdmap = None
if type(config) != type(None):
self.compute(config, options)
return
def compute(self, config, options):
if self.has_cnlm: pass
else:
# CONVERT ASE ATOMS TO SOAPXX TOPOLOGY
self.config, self.structure, top, frag_bond_mat, atom_bond_mat, frag_labels, atom_labels = \
soap.tools.structure_from_ase(
config,
do_partition=False,
add_fragment_com=False,
log=None)
self.has_structure = True
# COMPUTE SPECTRUM
options_cxx = self.getCxxOptions(options)
self.spectrum = soap.Spectrum(self.structure, options_cxx)
if self.verbose: self.log << "[cc] Computing density expansion" << self.log.endl
self.spectrum.compute()
self.has_spectrum = True
# EXTRACT CNLM ATOMIC
if PowerSpectrum.settings["store_cmap"]:
if self.verbose: self.log << "[py] Storing coefficient map (-> cmap)" << self.log.endl
self.cmap = []
for atomic in self.spectrum:
cmap = {}
types = atomic.getTypes()
for t in types:
cnlm = atomic.getLinear(t).array
cmap[t] = cnlm
self.cmap.append(cmap)
# EXTRACT CNLM GLOBAL
if PowerSpectrum.settings["store_gcmap"]:
if self.verbose: self.log << "[cc] Computing global coefficient map" << self.log.endl
self.gcmap = []
self.spectrum_global = self.spectrum.computeGlobal()
if self.verbose: self.log << "[py] Storing global coefficient map (-> gcmap)" << self.log.endl
types = atomic.getTypes()
for t in types:
self.gcmap.append({ t: atomic.getLinear(t).array })
self.has_spectrum_global = True
self.has_cnlm = True
# COMPUTE POWER SPECTRUM?
if PowerSpectrum.settings["cxx_compute_power"]:
self.computePower()
return
def computePower(self):
if not self.has_cnlm:
raise RuntimeError("Cannot compute power spectrum without computing spectrum first")
if self.has_pnkl: pass
else:
if self.verbose: self.log << "[cc] Computing power spectrum" << self.log.endl
self.spectrum.computePower()
if self.options['spectrum.gradients']:
if self.verbose: self.log << "[cc] ... gradients ..." << self.log.endl
self.spectrum.computePowerGradients()
if self.options['spectrum.global']:
if self.verbose: self.log << "[cc] ... global avg ..." << self.log.endl
self.spectrum.deleteGlobal()
self.spectrum.computeGlobal()
if PowerSpectrum.settings["store_sd"]:
if self.verbose: self.log << "[py] Storing specific descriptor (-> sd)" << self.log.endl
adaptor = soap.soapy.kernel.KernelAdaptorFactory['specific-unique'](
self.spectrum.options,
types_global=self.options['type_list'])
IX, IR, types = adaptor.adapt(self.spectrum, return_pos_matrix=True)
self.sd = IX
if PowerSpectrum.settings["store_gsd"]:
if self.verbose: self.log << "[py] Storing global specific descriptor (-> gsd)" << self.log.endl
adaptor = soap.soapy.kernel.KernelAdaptorFactory['global-specific'](
self.spectrum.options,