Work on integral / inversion examples.

src/atomicspectrum.cpp

@@ -176,7 +176,7 @@ void AtomicSpectrum::write(std::ostream &ofs) {
 AtomicSpectrum::xnkl_t *AtomicSpectrum::getXnkl(type_pair_t &types) {
     map_xnkl_t::iterator it = _map_xnkl.find(types);
     if (it == _map_xnkl.end()) {
-        if (types.first == "g" and types.second == "c") {
+        if (types.first == "" and types.second == "") {
             return _xnkl_generic_coherent;
         }
         else {

src/base/logger.hpp

@@ -27,10 +27,12 @@ else (*plog)()
 class LogBuffer : public std::stringbuf {
 
 public:
-	LogBuffer() : std::stringbuf(),
+	LogBuffer() : std::stringbuf(), _silent(false),
                 _errorPreface("ERROR "), _warnPreface("WARNING "),
                 _infoPreface(""), _dbgPreface("DEBUG "),
                 _writePreface(true) {}
+
+    void silence() { _silent = true; }
         
         // sets the log level (needed for output)
 	void setLogLevel(TLogLevel LogLevel) { _LogLevel = LogLevel; }
@@ -82,6 +84,7 @@ private:
   // Multithreading
   bool _maverick;
   bool _writePreface;
+  bool _silent;
   
   std::string _timePreface;
   std::string _errorPreface;
@@ -92,6 +95,7 @@ private:
 
 protected:
 	virtual int sync() {
+            if (_silent) return 0;
             
             std::ostringstream _message;
 
@@ -177,6 +181,10 @@ public:
 	}
         
         void setReportLevel( TLogLevel ReportLevel ) { _ReportLevel = ReportLevel; }
+        void silence() {
+            _silent = true;
+            dynamic_cast<LogBuffer*>(rdbuf())->silence();
+        }
         void setMultithreading( bool maverick ) { 
             _maverick = maverick;
             dynamic_cast<LogBuffer *>( rdbuf() )->setMultithreading( _maverick );
@@ -203,6 +211,7 @@ private:
     
     // if true, only a single processor job is executed
     bool      _maverick;
+    bool      _silent;
     
     std::string Messages() {
         return dynamic_cast<LogBuffer *>( rdbuf() )->Messages();

src/bindings.cpp

@@ -22,5 +22,8 @@ BOOST_PYTHON_MODULE(_soapxx)
     soap::RadialBasisFactory::registerAll();
     soap::AngularBasisFactory::registerAll();
     soap::CutoffFunctionFactory::registerAll();
+
+
+    boost::python::def("silence", &soap::GLOG_SILENCE);
 }
 

src/globals.cpp

@@ -4,5 +4,11 @@ namespace soap {
 
 Logger GLOG;
 
+void GLOG_SILENCE() {
+    GLOG() << "Silencing logger ..." << std::endl;
+    GLOG.silence();
+    return;
+}
+
 }
 

src/globals.hpp

@@ -7,6 +7,8 @@ namespace soap {
 
 extern Logger GLOG;
 
+void GLOG_SILENCE();
+
 namespace constants {
 
 const double ANGSTROM_TO_BOHR = 1./0.52917721067;

src/options.cpp

@@ -38,6 +38,32 @@ std::string Options::summarizeOptions() {
 	return info;
 }
 
+void Options::excludeCenters(boost::python::list &types) {
+    for (int i = 0; i < boost::python::len(types); ++i) {
+        std::string type = boost::python::extract<std::string>(types[i]);
+        _exclude_center[type] = true;
+    }
+    return;
+}
+
+void Options::excludeTargets(boost::python::list &types) {
+    for (int i = 0; i < boost::python::len(types); ++i) {
+        std::string type = boost::python::extract<std::string>(types[i]);
+        _exclude_target[type] = true;
+    }
+    return;
+}
+
+bool Options::doExcludeCenter(std::string &type) {
+    map_exclude_t::iterator it = _exclude_center.find(type);
+    return (it == _exclude_center.end()) ? false : true;
+}
+
+bool Options::doExcludeTarget(std::string &type) {
+    map_exclude_t::iterator it = _exclude_target.find(type);
+    return (it == _exclude_target.end()) ? false : true;
+}
+
 void Options::registerPython() {
 	using namespace boost::python;
 	void (Options::*set_int)(std::string, int) = &Options::set;
@@ -47,7 +73,8 @@ void Options::registerPython() {
 	class_<Options, Options*>("Options")
         .def("__str__", &Options::summarizeOptions)
 	    .def("summarizeOptions", &Options::summarizeOptions)
-		.def("configureCenters", &Options::configureCenters)
+		.def("excludeCenters", &Options::excludeCenters)
+		.def("excludeTargets", &Options::excludeTargets)
 		.def("set", set_int)
 		.def("set", set_double)
 		.def("set", set_string);

src/options.hpp

@@ -12,6 +12,7 @@ class Options
 {
 public:
 	typedef std::map<std::string, std::string> map_options_t;
+	typedef std::map<std::string, bool> map_exclude_t;
 
 	// CONSTRUCT & DESTRUCT
 	Options();
@@ -28,6 +29,12 @@ public:
 	void configureCenters(boost::python::list center_excludes) { _center_excludes = center_excludes; }
 	std::string summarizeOptions();
 
+	// EXCLUSIONS TODO Move this to a distinct class
+	void excludeCenters(boost::python::list &types);
+	void excludeTargets(boost::python::list &types);
+	bool doExcludeCenter(std::string &type);
+	bool doExcludeTarget(std::string &type);
+
 	// PYTHON
 	static void registerPython();
 
@@ -35,12 +42,17 @@ public:
 	template<class Archive>
 	void serialize(Archive &arch, const unsigned int version) {
 		arch & _key_value_map;
+		arch & _exclude_center;
+		arch & _exclude_target;
 		return;
 	}
 
 private:
 	boost::python::list _center_excludes;
 	map_options_t _key_value_map;
+
+	map_exclude_t _exclude_center;
+	map_exclude_t _exclude_target;
 };
 
 }

src/spectrum.cpp

@@ -49,6 +49,9 @@ void Spectrum::compute() {
 
     Structure::particle_it_t pit;
     for (pit = _structure->beginParticles(); pit != _structure->endParticles(); ++pit) {
+        // Continue if exclusion defined ...
+        if (_options->doExcludeCenter((*pit)->getType())) continue;
+        // Compute ...
     	AtomicSpectrum *atomic_spectrum = this->computeAtomic(*pit);
     	//atomic_spectrum->getReduced()->writeDensityOnGrid("density.expanded.cube", _options, _structure, *pit, true);
     	//atomic_spectrum->getReduced()->writeDensityOnGrid("density.explicit.cube", _options, _structure, *pit, false);
@@ -156,6 +159,9 @@ AtomicSpectrum *Spectrum::computeAtomic(Particle *center) {
     Structure::particle_it_t pit;
     for (pit = _structure->beginParticles(); pit != _structure->endParticles(); ++pit) {
 
+        // CHECK FOR EXCLUSIONS
+        if (_options->doExcludeTarget((*pit)->getType())) continue;
+
     	// FIND DISTANCE & DIRECTION, APPLY CUTOFF (= WEIGHT REDUCTION)
     	vec dr = _structure->connect(center->getPos(), (*pit)->getPos());
     	double r = soap::linalg::abs(dr);
@@ -249,6 +255,7 @@ void Spectrum::registerPython() {
     class_<Spectrum>("Spectrum", init<Structure &, Options &>())
     	.def(init<Structure &, Options &, Basis &>())
     	.def(init<std::string>())
+    	.def("__iter__", range<return_value_policy<reference_existing_object> >(&Spectrum::beginAtomic, &Spectrum::endAtomic))
 	    .def("compute", &Spectrum::compute)
 		.def("computePower", &Spectrum::computePower)
 		.def("addAtomic", &Spectrum::addAtomic)
@@ -263,6 +270,8 @@ void Spectrum::registerPython() {
 		.add_property("options", make_function(&Spectrum::getOptions, ref_existing()))
 		.add_property("basis", make_function(&Spectrum::getBasis, ref_existing()))
 		.add_property("structure", make_function(&Spectrum::getStructure, ref_existing()));
+    class_<atomspec_array_t>("AtomicSpectrumContainer")
+           .def(vector_indexing_suite<atomspec_array_t>());
 }
 
 /* STORAGE, BASIS, COMPUTATION, PARALLELIZATION */

src/spectrum.hpp

@@ -25,6 +25,7 @@ class Spectrum
 {
 public:
 	typedef std::vector<AtomicSpectrum*> atomspec_array_t;
+	typedef std::vector<AtomicSpectrum*>::iterator atomic_it_t;
 	typedef std::map<std::string, atomspec_array_t> map_atomspec_array_t;
 
 	Spectrum(std::string archfile);
@@ -36,6 +37,9 @@ public:
     Options *getOptions() { return _options; }
     Basis *getBasis() { return _basis; }
 
+    atomic_it_t beginAtomic() { return _atomspec_array.begin(); }
+    atomic_it_t endAtomic() { return _atomspec_array.end(); }
+
 	void saveAndClean() { std::cout << "spectrum::save&clean" << std::endl; }
 	void save(std::string archfile);
 	void load(std::string archfile);

src/tools/inverse.py

@@ -75,7 +75,7 @@ def invert_xnkl_aa_fixed_l(X, Y, N, l, Q_return, P_return):
         pass
     else:
         for m in np.arange(-l,l+1):
-            if m != -1:
+            if m != 0:
                 Qr_n[:,l+m] = 0.
             else:
                 pass
@@ -162,7 +162,7 @@ def invert_xnkl_aa(xnkl_cmplx, basis, l_smaller=0, l_larger=-1):
     
     for l in range(L1):
         if l < l_smaller: continue
-        if l == l_larger+1: break
+        if l and l == l_larger+1: break
 
         conv = (2.*l+1.)**0.5/(2.*2.**0.5*np.pi)
         

src/tools/loadwrite.py

@@ -17,7 +17,7 @@ def write_xyz(xyz_file, structure):
     ofs.close()
     return    
 
-def ase_load_all(folder):
+def ase_load_all(folder, log=None):
     cwd = os.getcwd()
     os.chdir(folder)
     config_listdir = sorted(os.listdir('./'))
@@ -25,10 +25,14 @@ def ase_load_all(folder):
     for config_file in config_listdir:
         if os.path.isfile(config_file):
             config_files.append(config_file)
-    ase_config_list = AseConfigList(config_files)
+    ase_config_list = AseConfigList(config_files, log=log)
     os.chdir(cwd)
     return ase_config_list
 
+def ase_load_single(config_file, log=None):
+    ase_config_list = AseConfigList([config_file], log=log)
+    return ase_config_list[0]
+
 def setup_structure_ase(label, ase_config):
     # DEFINE SYSTEM
     structure = soap.Structure(label)    
@@ -102,6 +106,8 @@ class AseConfigList(object):
         return self.configs[idx]
     def __iter__(self):
         return self.configs.__iter__()
+    def __len__(self):
+        return len(self.configs)
 
 class AseConfig(object):
     def __init__(self,