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

More kernel exploríng, added build file.

parent 91f6f42d
apply
build
soap
./soap
Makefile
test/old
test/back
......
# HEADERS
file(GLOB HEADERS *.hpp)
install(FILES ${HEADERS} DESTINATION ${LOCAL_INSTALL_DIR}/include/soap/base)
......@@ -43,6 +43,7 @@ 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;
_exclude_center_list.append(type);
}
return;
}
......@@ -51,20 +52,49 @@ 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;
_exclude_target_list.append(type);
}
return;
}
bool Options::doExcludeCenter(std::string &type) {
map_exclude_t::iterator it = _exclude_center.find(type);
auto 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);
auto it = _exclude_target.find(type);
return (it == _exclude_target.end()) ? false : true;
}
void Options::excludeCenterIds(boost::python::list &types) {
for (int i = 0; i < boost::python::len(types); ++i) {
int pid = boost::python::extract<int>(types[i]);
_exclude_center_id[pid] = true;
_exclude_center_id_list.append(pid);
}
return;
}
void Options::excludeTargetIds(boost::python::list &types) {
for (int i = 0; i < boost::python::len(types); ++i) {
int pid = boost::python::extract<int>(types[i]);
_exclude_target_id[pid] = true;
_exclude_target_id_list.append(pid);
}
return;
}
bool Options::doExcludeCenterId(int pid) {
auto it = _exclude_center_id.find(pid);
return (it == _exclude_center_id.end()) ? false : true;
}
bool Options::doExcludeTargetId(int pid) {
auto it = _exclude_target_id.find(pid);
return (it == _exclude_target_id.end()) ? false : true;
}
void Options::registerPython() {
using namespace boost::python;
void (Options::*set_int)(std::string, int) = &Options::set;
......@@ -74,10 +104,16 @@ void Options::registerPython() {
//void (Options::*set_bool)(std::string, bool) = &Options::set;
class_<Options, Options*>("Options")
.add_property("exclude_center_list", &Options::getExcludeCenterList)
.add_property("exclude_target_list", &Options::getExcludeTargetList)
.add_property("exclude_center_id_list", &Options::getExcludeCenterIdList)
.add_property("exclude_target_id_list", &Options::getExcludeTargetIdList)
.def("__str__", &Options::summarizeOptions)
.def("summarizeOptions", &Options::summarizeOptions)
.def("excludeCenters", &Options::excludeCenters)
.def("excludeTargets", &Options::excludeTargets)
.def("excludeCenterIds", &Options::excludeCenterIds)
.def("excludeTargetIds", &Options::excludeTargetIds)
.def("set", set_int)
.def("set", set_double)
.def("set", set_string)
......
......@@ -13,6 +13,7 @@ class Options
public:
typedef std::map<std::string, std::string> map_options_t;
typedef std::map<std::string, bool> map_exclude_t;
typedef std::map<int, bool> map_exclude_id_t;
// CONSTRUCT & DESTRUCT
Options();
......@@ -28,7 +29,7 @@ public:
void set(std::string key, double value) { this->set(key, boost::lexical_cast<std::string>(value)); }
bool hasKey(std::string key) { return (_key_value_map.find(key) == _key_value_map.end()) ? false : true; }
//void set(std::string key, bool value) { this->set(key, boost::lexical_cast<std::string>(value)); }
void configureCenters(boost::python::list center_excludes) { _center_excludes = center_excludes; }
//void configureCenters(boost::python::list center_excludes) { _center_excludes = center_excludes; }
std::string summarizeOptions();
// EXCLUSIONS TODO Move this to a distinct class
......@@ -36,6 +37,15 @@ public:
void excludeTargets(boost::python::list &types);
bool doExcludeCenter(std::string &type);
bool doExcludeTarget(std::string &type);
boost::python::list getExcludeCenterList() { return _exclude_center_list; }
boost::python::list getExcludeTargetList() { return _exclude_target_list; }
void excludeCenterIds(boost::python::list &pids);
void excludeTargetIds(boost::python::list &pids);
bool doExcludeCenterId(int pid);
bool doExcludeTargetId(int pid);
boost::python::list getExcludeCenterIdList() { return _exclude_center_id_list; }
boost::python::list getExcludeTargetIdList() { return _exclude_target_id_list; }
// PYTHON
static void registerPython();
......@@ -46,6 +56,9 @@ public:
arch & _key_value_map;
arch & _exclude_center;
arch & _exclude_target;
arch & _exclude_center_id;
arch & _exclude_target_id;
return;
}
......@@ -55,6 +68,13 @@ private:
map_exclude_t _exclude_center;
map_exclude_t _exclude_target;
boost::python::list _exclude_center_list;
boost::python::list _exclude_target_list;
map_exclude_id_t _exclude_center_id;
map_exclude_id_t _exclude_target_id;
boost::python::list _exclude_center_id_list;
boost::python::list _exclude_target_id_list;
};
}
......
......@@ -59,7 +59,8 @@ void Spectrum::compute(Structure::particle_array_t &centers, Structure::particle
Structure::particle_it_t pit;
for (pit = centers.begin(); pit != centers.end(); ++pit) {
// Continue if exclusion defined ...
if (_options->doExcludeCenter((*pit)->getType())) continue;
if (_options->doExcludeCenter((*pit)->getType()) ||
_options->doExcludeCenterId((*pit)->getId())) continue;
// Compute ...
AtomicSpectrum *atomic_spectrum = this->computeAtomic(*pit, targets);
this->addAtomic(atomic_spectrum);
......@@ -94,7 +95,8 @@ AtomicSpectrum *Spectrum::computeAtomic(Particle *center, Structure::particle_ar
for (pit = targets.begin(); pit != targets.end(); ++pit) {
// CHECK FOR EXCLUSIONS
if (_options->doExcludeTarget((*pit)->getType())) continue;
if (_options->doExcludeTarget((*pit)->getType()) ||
_options->doExcludeTargetId((*pit)->getId())) continue;
// FIND DISTANCE & DIRECTION, CHECK CUTOFF
vec dr = _structure->connect(center->getPos(), (*pit)->getPos());
......
......@@ -147,7 +147,7 @@ class KernelPotential:
logging.info("Compute energy from %d atomic environments ..." % n_acqu)
for n in range(n_acqu):
X = IX_acqu[n]
print "X_MAG", np.dot(X,X)
#print "X_MAG", np.dot(X,X)
ic = self.kernelfct.compute(self.IX, X)
energy += self.alpha.dot(ic)
#print "Projection", ic
......@@ -207,6 +207,15 @@ def perturb_positions(structure, exclude_pid=[]):
dz = np.random.uniform(-1.,1.)
part.pos = part.pos + 0.1*np.array([dx,dy,dz])
return [ part.pos for part in structure ]
def random_positions(structure, exclude_pid=[]):
for part in structure:
if part.id in exclude_pid: continue
dx = np.random.uniform(-1.,1.)
dy = np.random.uniform(-1.,1.)
dz = np.random.uniform(-1.,1.)
part.pos = np.array([dx,dy,dz])
return [ part.pos for part in structure ]
def apply_force_step(structure, forces, scale, constrain_particles=[]):
max_step = 0.05
......@@ -216,7 +225,24 @@ def apply_force_step(structure, forces, scale, constrain_particles=[]):
if df > max_f: max_f = df
if max_f > max_step: scale = scale*max_step/max_f
else: pass
print "Scale =", scale
#print "Scale =", scale
idx = -1
for part in structure:
idx += 1
if part.id in constrain_particles:
print "Skip force step, pid =", part.id
continue
#if np.random.uniform(0.,1.) > 0.5: continue
part.pos = part.pos + scale*forces[idx]
return [ part.pos for part in structure ]
def apply_force_norm_step(structure, forces, scale, constrain_particles=[]):
min_step = 0.5
max_f = 0.0
for f in forces:
df = np.dot(f,f)**0.5
if df > max_f: max_f = df
scale = min_step/max_f
idx = -1
for part in structure:
idx += 1
......
......@@ -5,7 +5,15 @@ import soap.tools
import os
import numpy as np
import logging
from momo import osio, endl, flush
from kernel import KernelPotential, restore_positions, apply_force_step
logging.basicConfig(
format='[%(asctime)s] %(message)s',
datefmt='%I:%M:%S',
level=logging.ERROR)
verbose = False
options = soap.Options()
options.excludeCenters([])
......@@ -18,7 +26,7 @@ options.set('radialbasis.integration_steps', 15)
options.set('radialcutoff.Rc', 6.8)
options.set('radialcutoff.Rc_width', 0.5)
options.set('radialcutoff.type', 'heaviside')
options.set('radialcutoff.center_weight', 0.)
options.set('radialcutoff.center_weight', 1.)
options.set('angularbasis.type', 'spherical-harmonic')
options.set('angularbasis.L', 6) # 6
options.set('densitygrid.N', 20)
......@@ -26,219 +34,8 @@ options.set('densitygrid.dx', 0.15)
options.set('spectrum.gradients', True)
options.set('kernel.adaptor', 'generic')
options.set('kernel.type', 'dot')
options.set('kernel.delta', 0.1)
options.set('kernel.xi', 1.)
class KernelAdaptorGeneric:
def __init__(self, options):
return
def adapt(self, spectrum, pid_list=None):
# IX = [
# X1 ->,
# X2 ->,
# ...
# ]
IX = np.zeros((0,0), dtype='complex128')
dimX = -1
# TODO More adaptors:
# TODO Particle ID mask, center-type mask
# TODO For type-resolved adaptors, increase dimensionality accordingly
for atomic_i in spectrum:
#print atomic_i.getCenter().pos
if pid_list:
pid = atomic_i.getCenter().id
if not pid in pid_list:
logging.debug("Skip, adapt atomic, pid = %d" % pid)
continue
Xi = np.real(atomic_i.getPower("","").array)
dimX = Xi.shape[0]*Xi.shape[1]
Xi = Xi.reshape((dimX))
#print "dim(X) =", dimX
if not IX.any():
IX = np.copy(Xi)
IX.resize((1, dimX))
else:
i = IX.shape[0]
IX.resize((i+1, dimX))
IX[-1,:] = Xi
#print IX
return IX
def adaptScalar(self, atomic):
X = np.real(atomic.getPower("","").array)
dimX = X.shape[0]*X.shape[1]
X = np.real(X.reshape((dimX)))
return X
def adaptGradients(self, atomic, nb_pid):
dxnkl_pid = atomic.getPowerGradGeneric(nb_pid)
dX_dx = dxnkl_pid.getArrayGradX()
dX_dy = dxnkl_pid.getArrayGradY()
dX_dz = dxnkl_pid.getArrayGradZ()
dimX = dX_dx.shape[0]*dX_dx.shape[1]
dX_dx = np.real(dX_dx.reshape((dimX)))
dX_dy = np.real(dX_dy.reshape((dimX)))
dX_dz = np.real(dX_dz.reshape((dimX)))
return dX_dx, dX_dy, dX_dz
class KernelFunctionDot(object):
def __init__(self, options):
self.delta = float(options.get('kernel.delta'))
self.xi = float(options.get('kernel.xi'))
return
def computeDot(self, IX, X, xi, delta):
return delta**2 * np.dot(IX,X)**xi
def compute(self, IX, X):
return self.computeDot(IX, X, self.xi, self.delta)
def computeDerivativeOuter(self, IX, X):
c = self.computeDot(IX, X, self.xi-1, self.delta)
return self.xi*np.diag(c).dot(IX)
KernelAdaptorFactory = { 'generic': KernelAdaptorGeneric }
KernelFunctionFactory = { 'dot':KernelFunctionDot }
class KernelPotential:
def __init__(self, options):
logging.info("Construct kernel potential ...")
self.basis = soap.Basis(options)
self.options = options
# CORE DATA
self.IX = None
self.alpha = None
self.dimX = None # <- Also used as flag set by first ::acquire call
# KERNEL
logging.info("Choose kernel function ...")
self.kernelfct = KernelFunctionFactory[options.get('kernel.type')](options)
# ADAPTOR
logging.info("Choose adaptor ...")
self.adaptor = KernelAdaptorFactory[options.get('kernel.adaptor')](options)
# INCLUSIONS / EXCLUSIONS
self.pid_list_acquire = [1]
self.pid_list_force = [1]
return
def acquire(self, structure, alpha):
logging.info("Acquire ...")
spectrum = soap.Spectrum(structure, self.options, self.basis)
spectrum.compute()
spectrum.computePower()
spectrum.computePowerGradients()
# New X's
logging.info("Adapt spectrum ...")
IX_acqu = self.adaptor.adapt(spectrum, self.pid_list_acquire)
n_acqu = IX_acqu.shape[0]
dim_acqu = IX_acqu.shape[1]
# New alpha's
alpha_acqu = np.zeros((n_acqu))
alpha_acqu.fill(alpha)
if not self.dimX:
# First time ...
self.dimX = dim_acqu
self.IX = IX_acqu
self.alpha = alpha_acqu
else:
# Check and extend ...
assert self.dimX == dim_acqu # Acquired descr. should match linear dim. of previous descr.'s
I = self.IX.shape[0]
self.IX.resize((I+n_acqu, self.dimX))
self.IX[I:I+n_acqu,:] = IX_acqu
self.alpha.resize((I+n_acqu))
self.alpha[I:I+n_acqu] = alpha_acqu
#print self.alpha
#print self.IX
logging.info("Acquired %d environments." % n_acqu)
return
def computeEnergy(self, structure):
logging.info("Start energy ...")
energy = 0.0
spectrum = soap.Spectrum(structure, self.options, self.basis)
spectrum.compute()
spectrum.computePower()
spectrum.computePowerGradients()
IX_acqu = self.adaptor.adapt(spectrum, self.pid_list_acquire)
n_acqu = IX_acqu.shape[0]
dim_acqu = IX_acqu.shape[1]
logging.info("Compute energy from %d atomic environments ..." % n_acqu)
for n in range(n_acqu):
X = IX_acqu[n]
print "X_MAG", np.dot(X,X)
ic = self.kernelfct.compute(self.IX, X)
energy += self.alpha.dot(ic)
print "Projection", ic
return energy
def computeForces(self, structure, verbose=False):
logging.info("Start forces ...")
if verbose:
for p in structure: print p.pos
forces = [ np.zeros((3)) for i in range(structure.n_particles) ]
logging.info("Compute forces on %d particles ..." % structure.n_particles)
# Compute X's, dX's
spectrum = soap.Spectrum(structure, self.options, self.basis)
spectrum.compute()
spectrum.computePower()
spectrum.computePowerGradients()
# Extract & compute force
for atomic in spectrum:
pid = atomic.getCenter().id
if not pid in self.pid_list_force:
logging.debug("Skip forces derived from environement with pid = %d" % pid)
continue
#if pid != 1: continue
nb_pids = atomic.getNeighbourPids()
logging.info(" Center %d" % (pid))
# neighbour-pid-independent kernel "prevector" (outer derivative)
X = self.adaptor.adaptScalar(atomic)
dIC = self.kernelfct.computeDerivativeOuter(self.IX, X)
alpha_dIC = self.alpha.dot(dIC)
for nb_pid in nb_pids:
# Force on neighbour
logging.info(" -> Nb %d" % (nb_pid))
dX_dx, dX_dy, dX_dz = self.adaptor.adaptGradients(atomic, nb_pid)
force_x = -alpha_dIC.dot(dX_dx)
force_y = -alpha_dIC.dot(dX_dy)
force_z = -alpha_dIC.dot(dX_dz)
forces[nb_pid-1][0] += force_x
forces[nb_pid-1][1] += force_y
forces[nb_pid-1][2] += force_z
#print forces
#raw_input('...')
return forces
def restore_positions(structure, positions):
idx = 0
for part in structure:
part.pos = positions[idx]
idx += 1
return [ part.pos for part in structure ]
def apply_force_step(structure, forces):
max_step = 0.05
max_f = 0.0
for f in forces:
df = np.dot(f,f)**0.5
if df > max_f: max_f = df
if max_f > max_step: scale = max_step/max_f
else: scale = 1.
idx = 0
for part in structure:
#if np.random.uniform(0.,1.) > 0.5: continue
part.pos = part.pos + scale*forces[idx]
idx += 1
return [ part.pos for part in structure ]
logging.basicConfig(
format='[%(asctime)s] %(message)s',
datefmt='%I:%M:%S',
level=logging.ERROR)
verbose = False
options.set('kernel.delta', 1.)
options.set('kernel.xi', 4.)
# STRUCTURE
xyzfile = 'config.xyz'
......@@ -255,6 +52,13 @@ soap.silence()
kernelpot = KernelPotential(options)
kernelpot.acquire(structure, 1.)
positions = [
np.array([0., 0., 0.]),
np.array([1.75, 0., 0.]),
np.array([0., 1.75, 0.]) ]
restore_positions(structure, positions)
#kernelpot.acquire(structure, 1.)
print "O'", kernelpot.computeEnergy(structure)
positions = [
np.array([0., 0., 0.]),
......@@ -290,15 +94,15 @@ print "D", kernelpot.computeEnergy(structure)
# MAP
Nx = 20
Ny = 20
dx = 0.2
dy = 0.2
dx = 0.25
dy = 0.25
ofs = open('out', 'w')
for i in range(20):
for j in range(20):
for i in range(Nx+1):
for j in range(Ny+1):
x = i*dx - 0.5*Nx*dx
y = j*dy - 0.5*Ny*dy
print x, y
print i, j
positions = [
np.array([0., 0., 0.]),
np.array([1.5, 0., 0.]),
......@@ -309,3 +113,4 @@ for i in range(20):
ofs.write('\n')
ofs.close()
......@@ -7,18 +7,20 @@ import numpy as np
import logging
from momo import osio, endl, flush
from kernel import KernelPotential, restore_positions, apply_force_step, perturb_positions
from kernel import KernelPotential, restore_positions, apply_force_step, apply_force_norm_step, perturb_positions, random_positions
logging.basicConfig(
format='[%(asctime)s] %(message)s',
datefmt='%I:%M:%S',
level=logging.DEBUG)
verbose = True
level=logging.ERROR)
verbose = False
exclude_center_pids = [2,3]
exclude_perturb_pids = [1]
constrain_pids = []
perturb_initial = False
random_initial = True
......@@ -28,11 +30,11 @@ options.excludeTargets([])
options.excludeCenterIds(exclude_center_pids)
options.excludeTargetIds([])
options.set('radialbasis.type', 'gaussian')
options.set('radialbasis.mode', 'adaptive')
options.set('radialbasis.mode', 'adaptive') # equispaced or adaptive
options.set('radialbasis.N', 9) # 9
options.set('radialbasis.sigma', 0.5) # 0.9
options.set('radialbasis.sigma', 1.5) # 0.9
options.set('radialbasis.integration_steps', 15)
options.set('radialcutoff.Rc', 6.8)
options.set('radialcutoff.Rc', 4.)
options.set('radialcutoff.Rc_width', 0.5)
options.set('radialcutoff.type', 'heaviside')
options.set('radialcutoff.center_weight', 1.)
......@@ -61,7 +63,7 @@ for p in positions_0: print p
# SETUP KERNEL
soap.silence()
kernelpot = KernelPotential(options)
kernelpot.acquire(structure, 1.)
kernelpot.acquire(structure, -1.)
#positions = [
# np.array([0., 0., 0.]),
......@@ -72,9 +74,15 @@ kernelpot.acquire(structure, 1.)
#restore_positions(structure, positions_0)
# PERTURB POSITIONS
positions_0_pert = perturb_positions(structure, exclude_pid=exclude_perturb_pids)
print "Positions after initial perturbation"
for p in positions_0_pert: print p
if perturb_initial:
positions_0_pert = perturb_positions(structure, exclude_pid=exclude_perturb_pids)
print "Positions after initial perturbation"
for p in positions_0_pert: print p
if random_initial:
positions_0_pert = random_positions(structure)
print "Positions after initial randomization"
for p in positions_0_pert: print p
raw_input('...')
......@@ -82,6 +90,12 @@ acquire_restore_every_nth = -1
ofs = open('log.xyz', 'w')
# Write first frame
ofs.write('%d\n\n' % structure.n_particles)
for p in structure.particles:
r = p.pos
ofs.write('%s %+1.7f %+1.7f %+1.7f\n' % (p.type, r[0], r[1], r[2]))
# Energy (int)
e_in = kernelpot.computeEnergy(structure)
osio << osio.my << "Energy(in)" << e_in << osio.endl
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment