Commit ae9521a1 authored by Carl Poelking's avatar Carl Poelking
Browse files

Real Ylm in inversion, Wigner normalization.

parent 1820246b
......@@ -8,43 +8,43 @@
namespace soap {
AtomicSpectrum::AtomicSpectrum(Particle *center, Basis *basis) {
this->null();
_center = center;
_center_pos = center->getPos();
_center_type = center->getType();
_basis = basis;
_qnlm_generic = new BasisExpansion(_basis);
this->null();
_center = center;
_center_pos = center->getPos();
_center_type = center->getType();
_basis = basis;
_qnlm_generic = new BasisExpansion(_basis);
}
AtomicSpectrum::~AtomicSpectrum() {
// Clean Qnlm's
for (map_qnlm_t::iterator it = _map_qnlm.begin(); it != _map_qnlm.end(); ++it) delete it->second;
_map_qnlm.clear();
if (_qnlm_generic) {
delete _qnlm_generic;
_qnlm_generic = NULL;
}
// Clean Xnkl's
for (map_xnkl_t::iterator it = _map_xnkl.begin(); it != _map_xnkl.end(); ++it) delete it->second;
_map_xnkl.clear();
if (_xnkl_generic_coherent) {
delete _xnkl_generic_coherent;
_xnkl_generic_coherent = NULL;
}
if (_xnkl_generic_incoherent) {
delete _xnkl_generic_incoherent;
_xnkl_generic_incoherent = NULL;
}
// Clean Qnlm's
for (map_qnlm_t::iterator it = _map_qnlm.begin(); it != _map_qnlm.end(); ++it) delete it->second;
_map_qnlm.clear();
if (_qnlm_generic) {
delete _qnlm_generic;
_qnlm_generic = NULL;
}
// Clean Xnkl's
for (map_xnkl_t::iterator it = _map_xnkl.begin(); it != _map_xnkl.end(); ++it) delete it->second;
_map_xnkl.clear();
if (_xnkl_generic_coherent) {
delete _xnkl_generic_coherent;
_xnkl_generic_coherent = NULL;
}
if (_xnkl_generic_incoherent) {
delete _xnkl_generic_incoherent;
_xnkl_generic_incoherent = NULL;
}
}
void AtomicSpectrum::null() {
_center = NULL;
_center_pos = vec(0,0,0);
_center_type = "?";
_basis = NULL;
_qnlm_generic = NULL;
_xnkl_generic_coherent = NULL;
_xnkl_generic_incoherent = NULL;
_center = NULL;
_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) {
......@@ -53,7 +53,7 @@ void AtomicSpectrum::invert(map_xnkl_t &map_xnkl, xnkl_t *xnkl_generic_coherent,
int L = _basis->getAngBasis()->L();
assert(xnkl_generic_coherent->getBasis() == _basis &&
"Trying to invert spectrum linked against foreign basis.");
"Trying to invert spectrum linked against foreign basis.");
PowerExpansion::coeff_t &xnkl = xnkl_generic_coherent->getCoefficients();
BasisExpansion::coeff_t &qnlm = _qnlm_generic->getCoefficients();
......@@ -65,11 +65,11 @@ void AtomicSpectrum::invert(map_xnkl_t &map_xnkl, xnkl_t *xnkl_generic_coherent,
// 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.);
}
}
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;
......@@ -79,19 +79,19 @@ void AtomicSpectrum::invert(map_xnkl_t &map_xnkl, xnkl_t *xnkl_generic_coherent,
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.);
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.);
// 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.);
// double xnnl = sqrt(xnkl(n*N+n,l).real());
// qnlm(n, l*l+l) = std::complex<double>(xnnl, 0.);
// }
// }
......@@ -99,96 +99,96 @@ void AtomicSpectrum::invert(map_xnkl_t &map_xnkl, xnkl_t *xnkl_generic_coherent,
}
void AtomicSpectrum::addQnlm(std::string type, qnlm_t &nb_expansion) {
assert(nb_expansion.getBasis() == _basis &&
"Should not sum expansions linked against different bases.");
map_qnlm_t::iterator it = _map_qnlm.find(type);
if (it == _map_qnlm.end()) {
_map_qnlm[type] = new BasisExpansion(_basis);
it = _map_qnlm.find(type);
}
it->second->add(nb_expansion);
_qnlm_generic->add(nb_expansion);
return;
assert(nb_expansion.getBasis() == _basis &&
"Should not sum expansions linked against different bases.");
map_qnlm_t::iterator it = _map_qnlm.find(type);
if (it == _map_qnlm.end()) {
_map_qnlm[type] = new BasisExpansion(_basis);
it = _map_qnlm.find(type);
}
it->second->add(nb_expansion);
_qnlm_generic->add(nb_expansion);
return;
}
AtomicSpectrum::qnlm_t *AtomicSpectrum::getQnlm(std::string type) {
if (type == "") {
return _qnlm_generic;
}
map_qnlm_t::iterator it = _map_qnlm.find(type);
if (it == _map_qnlm.end()) {
throw soap::base::OutOfRange("AtomicSpectrum: No such type '" + type + "'");
}
else {
return it->second;
}
if (type == "") {
return _qnlm_generic;
}
map_qnlm_t::iterator it = _map_qnlm.find(type);
if (it == _map_qnlm.end()) {
throw soap::base::OutOfRange("AtomicSpectrum: No such type '" + type + "'");
}
else {
return it->second;
}
}
void AtomicSpectrum::computePower() {
// Specific (i.e., type-dependent)
map_qnlm_t::iterator it1;
map_qnlm_t::iterator it2;
for (it1 = _map_qnlm.begin(); it1 != _map_qnlm.end(); ++it1) {
for (it2 = _map_qnlm.begin(); it2 != _map_qnlm.end(); ++it2) {
type_pair_t types(it1->first, it2->first);
GLOG() << " " << types.first << ":" << types.second << std::flush;
PowerExpansion *powex = new PowerExpansion(_basis);
powex->computeCoefficients(it1->second, it2->second);
_map_xnkl[types] = powex;
}
}
// Generic coherent
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;
_xnkl_generic_incoherent = new PowerExpansion(_basis);
map_xnkl_t::iterator it;
for (it = _map_xnkl.begin(); it != _map_xnkl.end(); ++it) {
_xnkl_generic_incoherent->add(it->second);
}
GLOG() << std::endl;
// Specific (i.e., type-dependent)
map_qnlm_t::iterator it1;
map_qnlm_t::iterator it2;
for (it1 = _map_qnlm.begin(); it1 != _map_qnlm.end(); ++it1) {
for (it2 = _map_qnlm.begin(); it2 != _map_qnlm.end(); ++it2) {
type_pair_t types(it1->first, it2->first);
GLOG() << " " << types.first << ":" << types.second << std::flush;
PowerExpansion *powex = new PowerExpansion(_basis);
powex->computeCoefficients(it1->second, it2->second);
_map_xnkl[types] = powex;
}
}
// Generic coherent
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;
_xnkl_generic_incoherent = new PowerExpansion(_basis);
map_xnkl_t::iterator it;
for (it = _map_xnkl.begin(); it != _map_xnkl.end(); ++it) {
_xnkl_generic_incoherent->add(it->second);
}
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;
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()) {
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;
map_xnkl_t::iterator it = _map_xnkl.find(types);
if (it == _map_xnkl.end()) {
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);
type_pair_t types(type1, type2);
return this->getXnkl(types);
}
void AtomicSpectrum::registerPython() {
......@@ -199,13 +199,13 @@ void AtomicSpectrum::registerPython() {
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>());
.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>());
}
}
......
......@@ -113,7 +113,7 @@ void BasisExpansion::writeDensity(
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") %
ofs << (boost::format("n %1$2d l %2$2d m %3$+2d %4$+1.7e %5$+1.7e") %
n % l %m % c_nlm_real % c_nlm_imag) << std::endl;
}
}
......@@ -169,7 +169,7 @@ void BasisExpansion::writeDensityOnGrid(
% (*pit)->getTypeId() % (dr.x()*conv) % (dr.y()*conv) % (dr.z()*conv);
}
GLOG() << "Fill grid " << std::flush;
GLOG() << "'" << filename << "': " << "Fill grid " << std::flush;
double int_density_dr = 0.0;
int ijk = 0;
for (int i = -I; i <= I; ++i) {
......@@ -255,7 +255,9 @@ void BasisExpansion::registerPython() {
class_<BasisExpansion, BasisExpansion*>("BasisExpansion", init<Basis*>())
.add_property("array", &BasisExpansion::getCoefficientsNumpy, &BasisExpansion::setCoefficientsNumpy)
.def("getArray", &BasisExpansion::getCoefficientsNumpy);
.def("getArray", &BasisExpansion::getCoefficientsNumpy)
.def("writeDensity", &BasisExpansion::writeDensity)
.def("writeDensityOnGrid", &BasisExpansion::writeDensityOnGrid);
}
}
......
......@@ -4,48 +4,48 @@
namespace soap {
PowerExpansion::PowerExpansion(Basis *basis) :
_basis(basis),
_L(basis->getAngBasis()->L()),
_N(basis->getRadBasis()->N()) {
_coeff = coeff_zero_t(_N*_N, _L+1);
_basis(basis),
_L(basis->getAngBasis()->L()),
_N(basis->getRadBasis()->N()) {
_coeff = coeff_zero_t(_N*_N, _L+1);
}
void PowerExpansion::computeCoefficients(BasisExpansion *basex1, BasisExpansion *basex2) {
if (!_basis) throw soap::base::APIError("PowerExpansion::computeCoefficients, basis not initialised.");
BasisExpansion::coeff_t &coeff1 = basex1->getCoefficients();
BasisExpansion::coeff_t &coeff2 = basex2->getCoefficients();
for (int n = 0; n < _N; ++n) {
for (int k = 0; k < _N; ++k) {
for (int l = 0; l < (_L+1); ++l) {
//std::cout << n << " " << k << " " << l << " : " << std::flush;
std::complex<double> c_nkl = 0.0;
for (int m = -l; m <= l; ++m) {
//std::cout << m << " " << std::flush;
c_nkl += coeff1(n, l*l+l+m)*std::conj(coeff2(k, l*l+l+m));
}
_coeff(n*_N+k, l) = c_nkl;
//std::cout << std::endl;
}
}
}
//throw soap::base::APIError("");
return;
if (!_basis) throw soap::base::APIError("PowerExpansion::computeCoefficients, basis not initialised.");
BasisExpansion::coeff_t &coeff1 = basex1->getCoefficients();
BasisExpansion::coeff_t &coeff2 = basex2->getCoefficients();
for (int n = 0; n < _N; ++n) {
for (int k = 0; k < _N; ++k) {
for (int l = 0; l < (_L+1); ++l) {
//std::cout << n << " " << k << " " << l << " : " << std::flush;
std::complex<double> c_nkl = 0.0;
for (int m = -l; m <= l; ++m) {
//std::cout << m << " " << std::flush;
c_nkl += coeff1(n, l*l+l+m)*std::conj(coeff2(k, l*l+l+m));
}
_coeff(n*_N+k, l) = 2.*sqrt(2.)*M_PI/sqrt(2.*l+1)*c_nkl; // Normalization = sqrt(8\pi^2/(2l+1))
//std::cout << std::endl;
}
}
}
//throw soap::base::APIError("");
return;
}
void PowerExpansion::add(PowerExpansion *other) {
assert(other->_basis == _basis &&
"Should not sum expansions linked against different bases.");
_coeff = _coeff + other->_coeff;
return;
assert(other->_basis == _basis &&
"Should not sum expansions linked against different bases.");
_coeff = _coeff + other->_coeff;
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.");
}
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() {
......@@ -54,42 +54,42 @@ boost::python::object PowerExpansion::getCoefficientsNumpy() {
}
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);
}
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;
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;
}
}
}
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();
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);
.add_property("array", &PowerExpansion::getCoefficientsNumpy, &PowerExpansion::setCoefficientsNumpy)
.def("getArray", &PowerExpansion::getCoefficientsNumpy);
}
}
......
......@@ -4,10 +4,9 @@ import numpy.linalg
from momo import osio, endl, flush
def compute_X(Q, P, M):
# Q is N x (L+1) matrix
# M = diag(1,2,2,...)
return Q.dot(M.dot(Q.T)) + P.dot(M.dot(P.T))
def compute_X(Qr):
# Q is N x (2*L+1) real-valued matrix
return Qr.dot(Qr.T)
def compute_usv(Q):
U, s, V = np.linalg.svd(Q, full_matrices=True)
......@@ -25,15 +24,15 @@ def compute_usv(Q):
return U, S, S_inv, V
def compute_rms_dist(X, Y):
D = X-Y
return (np.sum(D**2)/(D.shape[0]*D.shape[1]))**0.5
def compute_rms_dist(X, Y, eps=1e-10):
D = (X-Y)**2/(Y*Y+eps*eps)
return (np.sum(D)/(D.shape[0]*D.shape[1]))**0.5
def create_random(N, M):
def create_random(N, M, lower=-0.1, upper=0.1):
R = np.zeros((N,M))
for i in range(N):
for j in range(M):
R[i][j] = np.random.uniform(0.,1.)
R[i][j] = np.random.uniform(lower, upper)
return R
def load_X(tabfile):
......@@ -64,82 +63,85 @@ def setup_M(L1):
return M, M_inv
def invert_xnkl_aa_fixed_l(X, Y, N, l, Q_return, P_return):
assert X.shape == Y.shape == (N,N)
assert X.shape == Y.shape == (N,N)
debug = False
l1 = l+1
l21 = 2*l+1
# Real and imaginary l-expansions
Q_n = create_random(N,l1)
P_n = create_random(N,l1)
P_n[:,0] = 0. # Yl0 is real for all l
# Metric to reflect symmetries Y
M, M_inv = setup_M(l1)
with_P = False
if not with_P: P_n[:,:] = 0.
print "Starting l = %d, Q0+iP0 =" % (l1-1)
for i in range(N):
print Q_n[i], "+ i*", P_n[i]
# Real-valued coefficients (associated with real-valued Ylm)
Qr_n = create_random(N,l21,-0.1,0.1)
if l == 0:
pass
else:
for m in np.arange(-l,l+1):
if m != -1:
Qr_n[:,l+m] = 0.
else:
pass
#Qr_n[:,l+m] = 1.
print "Starting l = %d, Qr_0 =" % l
print Qr_n
i = 0
max_i = 512
max_i = 4096
converged = False
omega_sor = 0.5
tolerance = 1e-12
while True:
i += 1
# Compute Q_n