Commit 1820246b authored by Carl Poelking's avatar Carl Poelking
Browse files

Inverse framework.

parent 07b25c69
......@@ -47,6 +47,57 @@ void AtomicSpectrum::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();
int L = _basis->getAngBasis()->L();
assert(xnkl_generic_coherent->getBasis() == _basis &&
"Trying to invert spectrum linked against foreign basis.");
PowerExpansion::coeff_t &xnkl = xnkl_generic_coherent->getCoefficients();
BasisExpansion::coeff_t &qnlm = _qnlm_generic->getCoefficients();
type_pair_t types(type1, type2);
if (type1 == "g" && type2 == "c") xnkl = xnkl_generic_coherent->getCoefficients();
else if (type1 == "g" && type2 == "i") throw soap::base::NotImplemented("::invert g/i");
else xnkl = map_xnkl[types]->getCoefficients();
// ZERO QNLM TO BE SAFE
for (int n= 0; n < N; ++n) {
for (int l = 0; l <= L; ++l) {
for (int m = -l; m <= l; ++m) {
qnlm(n, l*l+l+m) = std::complex<double>(0.,0.);
}
}
}
// Q000 from X000, then Qk00 = X0k0/Q000
typedef std::complex<double> cmplx;
int n = 0;
int k = 0;
int l = 0;
int m = 0;
qnlm(n, l*l+l+m) = cmplx(sqrt(xnkl(n*N+k, l).real()), 0.);
for (k = 1; k < N; ++k) {
qnlm(k, l*l+l+m) = cmplx(xnkl(n*N+k, l).real()/qnlm(0, l*l+l+m).real(), 0.);
}
// // FILL Qn00's USING Xnn0's
// for (int n = 0; n < N; ++n) {
// double xnn0 = sqrt(xnkl(n*N+n,0).real());
// std::cout << xnn0 << std::endl;
// qnlm(n, 0) = std::complex<double>(xnn0, 0.);
// }
// for (int n = 0; n < N; ++n) {
// for (int l = 0; l <= L; ++l) {
// double xnnl = sqrt(xnkl(n*N+n,l).real());
// qnlm(n, l*l+l) = std::complex<double>(xnnl, 0.);
// }
// }
return;
}
void AtomicSpectrum::addQnlm(std::string type, qnlm_t &nb_expansion) {
assert(nb_expansion.getBasis() == _basis &&
"Should not sum expansions linked against different bases.");
......@@ -90,6 +141,7 @@ void AtomicSpectrum::computePower() {
GLOG() << " g/c" << std::flush;
if (_xnkl_generic_coherent) delete _xnkl_generic_coherent;
_xnkl_generic_coherent = new PowerExpansion(_basis);
_xnkl_generic_coherent->computeCoefficients(_qnlm_generic, _qnlm_generic);
// Generic incoherent
GLOG() << " g/i" << std::flush;
if (_xnkl_generic_incoherent) delete _xnkl_generic_incoherent;
......@@ -101,13 +153,60 @@ void AtomicSpectrum::computePower() {
GLOG() << std::endl;
}
void AtomicSpectrum::write(std::ostream &ofs) {
throw soap::base::NotImplemented("AtomicSpectrum::write");
map_qnlm_t::iterator it;
for (it = _map_qnlm.begin(); it != _map_qnlm.end(); ++it) {
std::cout << "TYPE" << it->first << std::endl;
qnlm_t *qnlm = it->second;
qnlm_t::coeff_t &coeff = qnlm->getCoefficients();
int L = qnlm->getBasis()->getAngBasis()->L();
int N = qnlm->getBasis()->getRadBasis()->N();
for (int n = 0; n < N; ++n) {
for (int l = 0; l < (L+1); ++l) {
for (int m = -l; m <= l; ++l) {
std::cout << n << " " << l << " " << m << " " << coeff(n,l*l+l+m) << std::endl;
}
}
}
}
return;
}
AtomicSpectrum::xnkl_t *AtomicSpectrum::getXnkl(type_pair_t &types) {
map_xnkl_t::iterator it = _map_xnkl.find(types);
if (it == _map_xnkl.end()) {
throw soap::base::OutOfRange("AtomicSpectrum: No such type pair '" + types.first + ":" + types.second + "'");
if (types.first == "g" and types.second == "c") {
return _xnkl_generic_coherent;
}
else {
throw soap::base::OutOfRange("AtomicSpectrum: No such type pair '" + types.first + ":" + types.second + "'");
}
}
return it->second;
}
AtomicSpectrum::xnkl_t *AtomicSpectrum::getPower(std::string type1, std::string type2) {
type_pair_t types(type1, type2);
return this->getXnkl(types);
}
void AtomicSpectrum::registerPython() {
using namespace boost::python;
// Does not work with boost::python ??
//xnkl_t *(AtomicSpectrum::*getXnkl_string_string)(std::string, std::string) = &AtomicSpectrum::getPower;
class_<AtomicSpectrum, AtomicSpectrum*>("AtomicSpectrum", init<>())
.def(init<Particle*, Basis*>())
.def("addLinear", &AtomicSpectrum::addQnlm)
.def("getLinear", &AtomicSpectrum::getQnlm, return_value_policy<reference_existing_object>())
.def("computePower", &AtomicSpectrum::computePower)
.def("getPower", &AtomicSpectrum::getPower, return_value_policy<reference_existing_object>())
.def("getCenter", &AtomicSpectrum::getCenter, return_value_policy<reference_existing_object>())
.def("getCenterType", &AtomicSpectrum::getCenterType, return_value_policy<reference_existing_object>())
.def("getCenterPos", &AtomicSpectrum::getCenterPos, return_value_policy<reference_existing_object>());
}
}
......@@ -32,6 +32,9 @@ public:
AtomicSpectrum() { this->null(); }
~AtomicSpectrum();
void null();
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);
// CENTER & BASIS METHODS
Particle *getCenter() { return _center; }
......@@ -42,12 +45,16 @@ public:
void addQnlm(std::string type, qnlm_t &nb_expansion);
qnlm_t *getQnlm(std::string type);
qnlm_t *getQnlmGeneric() { return _qnlm_generic; }
// XNLM METHODS
// XNKL METHODS
void computePower();
xnkl_t *getPower(std::string type1, std::string type2);
xnkl_t *getXnkl(type_pair_t &types);
map_xnkl_t &getXnklMap() { return _map_xnkl; }
xnkl_t *getXnklGenericCoherent() { return _xnkl_generic_coherent; }
xnkl_t *getXnklGenericIncoherent() { return _xnkl_generic_incoherent; }
static void registerPython();
template<class Archive>
void serialize(Archive &arch, const unsigned int version) {
// Center & basis
......
#include "basis.hpp"
#include <math.h>
#include "linalg/numpy.hpp"
#include "basis.hpp"
namespace soap {
// =====
......@@ -30,6 +33,13 @@ Basis::~Basis() {
_cutoff = NULL;
}
void Basis::registerPython() {
using namespace boost::python;
class_<Basis, Basis*>("Basis", init<>())
.add_property("N", make_function(&Basis::N, copy_const()))
.add_property("L", make_function(&Basis::L, copy_const()));
}
// ==============
// BasisExpansion
// ==============
......@@ -77,6 +87,42 @@ void BasisExpansion::conjugate() {
}
}
void BasisExpansion::writeDensity(
std::string filename,
Options *options,
Structure *structure,
Particle *center) {
if (_angbasis == NULL || _radbasis == NULL) {
throw soap::base::APIError("<BasisExpansion::writeDensityOnGrid> "
"Object not linked against basis.");
}
std::ofstream ofs;
ofs.open(filename.c_str(), std::ofstream::out);
if (!ofs.is_open()) {
throw soap::base::IOError("Bad file handle: " + filename);
}
double sum_intensity = 0.0;
BasisExpansion::coeff_t &coeff = this->getCoefficients();
for (int n = 0; n < _radbasis->N(); ++n) {
for (int l = 0; l <= _angbasis->L(); ++l) {
for (int m = -l; m <= l; ++m) {
std::complex<double> c_nlm = coeff(n, l*l+l+m);
double c_nlm_real = c_nlm.real();
double c_nlm_imag = c_nlm.imag();
sum_intensity += c_nlm_real*c_nlm_real + c_nlm_imag*c_nlm_imag;
ofs << (boost::format("%1$2d %2$2d %3$+2d %4$+1.7e %5$+1.7e") %
n % l %m % c_nlm_real % c_nlm_imag) << std::endl;
}
}
}
GLOG() << "<BasisExpansion::writeDensity> Summed intensity = " << sum_intensity << std::endl;
ofs.close();
}
void BasisExpansion::writeDensityOnGrid(
std::string filename,
Options *options,
......@@ -188,5 +234,29 @@ void BasisExpansion::writeDensityOnGrid(
return;
}
void BasisExpansion::setCoefficientsNumpy(boost::python::object &np_array) {
soap::linalg::numpy_converter npc("complex128");
npc.numpy_to_ublas< std::complex<double> >(np_array, _coeff);
int N = _basis->getRadBasis()->N();
int L = _basis->getAngBasis()->L();
if (_coeff.size1() != N ||
_coeff.size2() != (L+1)*(L+1)) {
throw soap::base::APIError("<BasisExpansion::setCoefficientsNumpy> Matrix size not consistent with basis.");
}
}
boost::python::object BasisExpansion::getCoefficientsNumpy() {
soap::linalg::numpy_converter npc("complex128");
return npc.ublas_to_numpy< std::complex<double> >(_coeff);
}
void BasisExpansion::registerPython() {
using namespace boost::python;
class_<BasisExpansion, BasisExpansion*>("BasisExpansion", init<Basis*>())
.add_property("array", &BasisExpansion::getCoefficientsNumpy, &BasisExpansion::setCoefficientsNumpy)
.def("getArray", &BasisExpansion::getCoefficientsNumpy);
}
}
......@@ -28,6 +28,10 @@ public:
RadialBasis *getRadBasis() { return _radbasis; }
AngularBasis *getAngBasis() { return _angbasis; }
CutoffFunction *getCutoff() { return _cutoff; }
const int &N() { return _radbasis->N(); }
const int &L() { return _angbasis->L(); }
static void registerPython();
template<class Archive>
void serialize(Archive &arch, const unsigned int version) {
......@@ -63,9 +67,15 @@ public:
void computeCoefficients(double r, vec d, double weight, double sigma);
void add(BasisExpansion &other) { _coeff = _coeff + other._coeff; }
void conjugate();
void writeDensity(std::string filename, Options *options,
Structure *structure, Particle *center);
void writeDensityOnGrid(std::string filename, Options *options,
Structure *structure, Particle *center, bool fromExpansion);
void setCoefficientsNumpy(boost::python::object &array);
boost::python::object getCoefficientsNumpy();
static void registerPython();
template<class Archive>
void serialize(Archive &arch, const unsigned int version) {
arch & _basis;
......
......@@ -7,11 +7,17 @@ namespace soap {
BOOST_PYTHON_MODULE(_soapxx)
{
using namespace boost::python;
soap::Particle::registerPython();
soap::Segment::registerPython();
soap::Structure::registerPython();
soap::Segment::registerPython();
soap::Particle::registerPython();
soap::Options::registerPython();
soap::Spectrum::registerPython();
soap::Basis::registerPython();
soap::AtomicSpectrum::registerPython();
soap::BasisExpansion::registerPython();
soap::PowerExpansion::registerPython();
soap::RadialBasisFactory::registerAll();
soap::AngularBasisFactory::registerAll();
......
#ifndef _SOAP_LINALG_NUMPY
#define _SOAP_LINALG_NUMPY
// See https://mail.python.org/pipermail/cplusplus-sig/2006-December/011409.html
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/python.hpp>
namespace soap { namespace linalg {
/** Converts between numpy arrays and boost ublas matrices. */
struct numpy_converter
{
//using namespace boost::python;
//typedef boost::python::object object;
//typedef boost::python::tuple tuple;
boost::python::object numpy;
boost::python::object array_type;
boost::python::object array_function;
boost::python::object dtype;
/** Constructor. dtype_name determines type of matrices created. */
numpy_converter( const char * dtype_name = "float64" )
{
PyObject* module = ::PyImport_Import( boost::python::object( "numpy" ).ptr() );
if( ! module )
{
throw std::logic_error( "Could not import numpy" );
}
numpy = boost::python::object( boost::python::handle<>( module ) );
array_type = numpy.attr( "ndarray" );
array_function = numpy.attr( "array" );
set_dtype( dtype_name );
}
/** Set which dtype the created numpy matrices have. */
void set_dtype( const char * dtype_name = "float64" )
{
dtype = numpy.attr( dtype_name );
}
/** Convert a numpy matrix to a ublas one. */
template< typename T >
boost::numeric::ublas::matrix< T > &
numpy_to_ublas(
boost::python::object a,
boost::numeric::ublas::matrix< T > & m )
{
boost::python::tuple shape( a.attr("shape") );
if( boost::python::len( shape ) != 2 )
{
throw std::logic_error( "numeric::array must have 2 dimensions" );
}
m.resize(
boost::python::extract< unsigned >( shape[0] ),
boost::python::extract< unsigned >( shape[1] ) );
for( unsigned i = 0; i < m.size1(); ++i )
{
for( unsigned j = 0; j < m.size2(); ++j )
{
m( i, j ) = boost::python::extract< T >( a[
boost::python::make_tuple( i, j ) ] );
}
}
return m;
}
/** Convert a ublas matrix to a numpy matrix. */
template< typename T >
boost::python::object
ublas_to_numpy(
const boost::numeric::ublas::matrix< T > & m )
{
//create a numpy array to put it in
boost::python::object result(
array_type(
boost::python::make_tuple( m.size1(), m.size2() ),
dtype ) );
//copy the elements
for( unsigned i = 0; m.size1() != i; ++i )
{
for( unsigned j = 0; m.size2() != j; ++j )
{
result[ boost::python::make_tuple( i, j ) ] = m( i, j );
}
}
return result;
}
};
}} /* CLOSE NAMESPACE */
#endif
......@@ -44,7 +44,8 @@ void Options::registerPython() {
void (Options::*set_double)(std::string, double) = &Options::set;
void (Options::*set_string)(std::string, std::string) = &Options::set;
class_<Options>("Options")
class_<Options, Options*>("Options")
.def("__str__", &Options::summarizeOptions)
.def("summarizeOptions", &Options::summarizeOptions)
.def("configureCenters", &Options::configureCenters)
.def("set", set_int)
......
#include "power.hpp"
#include "linalg/numpy.hpp"
namespace soap {
......@@ -30,6 +31,7 @@ void PowerExpansion::computeCoefficients(BasisExpansion *basex1, BasisExpansion
//throw soap::base::APIError("");
return;
}
void PowerExpansion::add(PowerExpansion *other) {
assert(other->_basis == _basis &&
"Should not sum expansions linked against different bases.");
......@@ -37,5 +39,58 @@ void PowerExpansion::add(PowerExpansion *other) {
return;
}
void PowerExpansion::setCoefficientsNumpy(boost::python::object &np_array) {
soap::linalg::numpy_converter npc("complex128");
npc.numpy_to_ublas< std::complex<double> >(np_array, _coeff);
if (_coeff.size1() != _N*_N ||
_coeff.size2() != _L+1) {
throw soap::base::APIError("<PowerExpansion::setCoefficientsNumpy> Matrix size not consistent with basis.");
}
}
boost::python::object PowerExpansion::getCoefficientsNumpy() {
soap::linalg::numpy_converter npc("complex128");
return npc.ublas_to_numpy< std::complex<double> >(_coeff);
}
void PowerExpansion::writeDensity(
std::string filename,
Options *options,
Structure *structure,
Particle *center) {
std::ofstream ofs;
ofs.open(filename.c_str(), std::ofstream::out);
if (!ofs.is_open()) {
throw soap::base::IOError("Bad file handle: " + filename);
}
double sum_intensity = 0.0;
PowerExpansion::coeff_t &coeff = this->getCoefficients();
for (int n = 0; n < _N; ++n) {
for (int k = 0; k < _N; ++k) {
for (int l = 0; l <= _L; ++l) {
std::complex<double> c_nkl = coeff(n*_N+k, l);
double c_nkl_real = c_nkl.real();
double c_nkl_imag = c_nkl.imag();
sum_intensity += c_nkl_real*c_nkl_real + c_nkl_imag*c_nkl_imag;
ofs << (boost::format("%1$2d %2$2d %3$+2d %4$+1.7e %5$+1.7e") %
n % k % l % c_nkl_real % c_nkl_imag) << std::endl;
}
}
}
GLOG() << "<PowerExpansion::writeDensity> Summed intensity = " << sum_intensity << std::endl;
ofs.close();
}
void PowerExpansion::registerPython() {
using namespace boost::python;
class_<PowerExpansion, PowerExpansion*>("PowerExpansion", init<Basis*>())
.add_property("array", &PowerExpansion::getCoefficientsNumpy, &PowerExpansion::setCoefficientsNumpy)
.def("getArray", &PowerExpansion::getCoefficientsNumpy);
}
}
......@@ -23,8 +23,15 @@ public:
PowerExpansion() : _basis(NULL), _L(-1), _N(-1) {;}
PowerExpansion(Basis *basis);
Basis *getBasis() { return _basis; }
coeff_t &getCoefficients() { return _coeff; }
void computeCoefficients(BasisExpansion *basex1, BasisExpansion *basex2);
void add(PowerExpansion *other);
void writeDensity(std::string filename, Options *options, Structure *structure, Particle *center);
void setCoefficientsNumpy(boost::python::object &np_array);
boost::python::object getCoefficientsNumpy();
static void registerPython();
template<class Archive>
void serialize(Archive &arch, const unsigned int version) {
......
......@@ -7,23 +7,30 @@
namespace soap {
Spectrum::Spectrum(Structure &structure, Options &options)
: _log(NULL), _options(&options), _structure(&structure) {
Spectrum::Spectrum(Structure &structure, Options &options) :
_log(NULL), _options(&options), _structure(&structure), _own_basis(true) {
GLOG() << "Configuring spectrum ..." << std::endl;
// CREATE & CONFIGURE BASIS
_basis = new Basis(&options);
}
Spectrum::Spectrum(Structure &structure, Options &options, Basis &basis) :
_log(NULL), _options(&options), _structure(&structure), _basis(&basis), _own_basis(false) {
;
}
Spectrum::Spectrum(std::string archfile) :
_log(NULL), _options(NULL), _structure(NULL), _basis(NULL) {
_log(NULL), _options(NULL), _structure(NULL), _basis(NULL), _own_basis(true) {
this->load(archfile);
}
Spectrum::~Spectrum() {
delete _log;
_log = NULL;
delete _basis;
_basis = NULL;
if (_own_basis) {
delete _basis;
_basis = NULL;
}
atomspec_array_t::iterator it;
for (it = _atomspec_array.begin(); it != _atomspec_array.end(); ++it) {
delete *it;
......@@ -45,11 +52,11 @@ void Spectrum::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);
this->add(atomic_spectrum);
this->addAtomic(atomic_spectrum);
}
}
void Spectrum::writeDensityOnGrid(int slot_idx, std::string center_type, std::string density_type) {
AtomicSpectrum *Spectrum::getAtomic(int slot_idx, std::string center_type) {
AtomicSpectrum *atomic_spectrum = NULL;
// FIND SPECTRUM
if (center_type == "") {
......@@ -79,6 +86,11 @@ void Spectrum::writeDensityOnGrid(int slot_idx, std::string center_type, std::st
}
}
}
return atomic_spectrum;
}
void Spectrum::writeDensityOnGrid(int slot_idx, std::string center_type, std::string density_type) {
AtomicSpectrum *atomic_spectrum = this->getAtomic(slot_idx, center_type);
// WRITE CUBE FILES
if (atomic_spectrum) {
atomic_spectrum->getQnlm(density_type)->writeDensityOnGrid(
......@@ -89,6 +101,37 @@ void Spectrum::writeDensityOnGrid(int slot_idx, std::string center_type, std::st
return;
}
void Spectrum::writeDensityOnGridInverse(int slot_idx, std::string center_type, std::string type1, std::string type2) {
AtomicSpectrum *atomic_spectrum = this->getAtomic(slot_idx, center_type);
AtomicSpectrum inverse_atomic_spectrum(atomic_spectrum->getCenter(), atomic_spectrum->getBasis());
inverse_atomic_spectrum.invert(atomic_spectrum->getXnklMap(), atomic_spectrum->getXnklGenericCoherent(), type1, type2);
inverse_atomic_spectrum.getQnlm("")->writeDensityOnGrid(
"density.inverse.cube", _options, _structure, inverse_atomic_spectrum.getCenter(), true);
inverse_atomic_spectrum.getQnlm("")->writeDensity(
"density.inverse.coeff", _options, _structure, inverse_atomic_spectrum.getCenter());
return;
}
void Spectrum::writeDensity(int slot_idx, std::string center_type, std::string density_type) {
AtomicSpectrum *atomic_spectrum = this->getAtomic(slot_idx, center_type);
// WRITE COEFF FILE
if (atomic_spectrum) {
atomic_spectrum->getQnlm(density_type)->writeDensity(
"density.expanded.coeff", _options, _structure, atomic_spectrum->getCenter());
}
return;
}
void Spectrum::writePowerDensity(int slot_idx, std::string center_type, std::string type1, std::string type2) {
AtomicSpectrum *atomic_spectrum = this->getAtomic(slot_idx, center_type);
if (atomic_spectrum) {
AtomicSpectrum::type_pair_t types(type1, type2);
atomic_spectrum->getXnkl(types)->writeDensity(
"density.power.coeff", _options, _structure, atomic_spectrum->getCenter());
}
return;
}
void Spectrum::computePower() {
atomspec_array_t::iterator it;
for (it = _atomspec_array.begin(); it != _atomspec_array.end(); ++it) {
......@@ -173,7 +216,7 @@ AtomicSpectrum *Spectrum::computeAtomic(Particle *center) {
}
void Spectrum::add(AtomicSpectrum *atomspec) {
void Spectrum::addAtomic(AtomicSpectrum *atomspec) {
assert(atomspec->getBasis() == _basis &&
"Should not append atomic spectrum linked against different basis.");
std::string atomspec_type = atomspec->getCenterType();
......@@ -204,13 +247,22 @@ void Spectrum::load(std::string archfile) {
void Spectrum::registerPython() {