Merge branch 'master' of https://github.com/capoe/soapxx

src/soap/SOAPRC.in

@@ -6,6 +6,7 @@ test "$is_csh" = 101 && goto CSH
 export PYTHONPATH="${PYTHONPATH}:@CMAKE_INSTALL_PREFIX@"
 export SOAP_ROOT="@CMAKE_INSTALL_PREFIX@"
 export LD_LIBRARY_PATH="{LD_LIBRARY_PATH}:@CMAKE_INSTALL_PREFIX@/soap"
+return
 
 # CSH/TCSH
 CSH:

src/soap/__init__.py

 from _soapxx import *
 from linalg import *
 from . import soapy
+from . import tools
 

src/soap/soapy/CMakeLists.txt

-install(FILES __init__.py elements.py util.py npthreads.py kernel.py simspace.py pca.py dimred.py lagraph.py lamatch.py math.py momo.py DESTINATION ${LOCAL_INSTALL_DIR}/soapy)
+install(FILES __init__.py learn.py elements.py util.py npthreads.py kernel.py simspace.py pca.py dimred.py lagraph.py lamatch.py math.py momo.py DESTINATION ${LOCAL_INSTALL_DIR}/soapy)
 

src/soap/soapy/__init__.py

@@ -5,4 +5,5 @@ from pca import *
 from util import *
 from math import *
 from elements import *
+import learn
 

src/soap/soapy/kernel.py

@@ -162,24 +162,37 @@ class KernelAdaptorSpecificUnique(object):
         self.types = types_global
         self.S = len(types_global)
         return
-    def adapt(self, spectrum):
-        IX = np.zeros((0.,0.), dtype='float64')
+    def adapt(self, spectrum, return_pos_matrix=False):
+        IX = np.zeros((0,0), dtype='float64') # feature matrix
         dimX = -1
+        IR = np.zeros((0,0), dtype='float64') # position matrix
+        types = []
         for atomic_i in spectrum:
             Xi_unnorm, Xi_norm = self.adaptScalar(atomic_i)
+            Ri = atomic_i.getCenter().pos
+            types.append(atomic_i.getCenter().type)
             dimX = Xi_norm.shape[0]
             if not IX.any():
                 IX = np.copy(Xi_norm) # TODO Is this necessary?
                 IX.resize((1, dimX))
+                IR = np.copy(Ri)
+                IR.resize((1, 3))
             else:
                 i = IX.shape[0]
                 IX.resize((i+1, dimX))
                 IX[-1,:] = Xi_norm
-            #print IX
-        return IX
-    def adaptScalar(self, atomic):
+                IR.resize((i+1, 3))
+                IR[-1,:] = Ri
+        if return_pos_matrix:
+            return IX, IR, types
+        else:
+            return IX
+    def adaptScalar(self, atomic, epsilon=1e-20):
         X = reduce_xnklab_atomic(atomic, self.types)
-        X_norm = X/np.dot(X,X)**0.5
+        X_mag = np.dot(X,X)**0.5
+        if X_mag < epsilon: 
+            X_mag = 1.
+        X_norm = X/X_mag
         return X, X_norm
 
 class KernelAdaptorSpecific(object):
@@ -203,7 +216,6 @@ class KernelAdaptorSpecific(object):
                 i = IX.shape[0]
                 IX.resize((i+1, dimX))
                 IX[-1,:] = Xi_norm
-            #print IX
         return IX
     def adaptScalar(self, atomic):
         xnklab_atomic = Xnklab(atomic, self.types)

src/soap/soapy/lagraph.py

@@ -680,10 +680,7 @@ class ParticleGraph(object):
                 #print z, z_idx
             P = np.array(ix)
             dim = P.shape[1]
-            if options_soap.get('kernel.adaptor') in [ 'global-generic', 'global-specific' ]:
-                pass
-            else:
-                assert P.shape[0] == n_atoms 
+            assert P.shape[0] == n_atoms 
             #print P.dot(P.T)
         elif descriptor_type == 'npy_load':
             folder = options_descriptor["folder"]

src/soap/soapy/learn.py

+import numpy as np
+
+def subsample_array(array, n_select, method='stride', stride_shift=0):
+    # Number of data items; selected, discarded items
+    if type(array) == list:
+        n_data = len(array)
+    else:
+        n_data = array.shape[0]
+    n_discard = n_data - n_select
+    assert n_discard >= 0 # Trying to select more data than available?
+    # Subsample
+    if method == 'random':
+        # Random access indices
+        idcs = np.arange(0, n_data)
+        np.random.shuffle(idcs)
+        # Create selections
+        if type(array) == list:
+            array_select = []
+            array_discard = []
+            for i in range(n_data):
+                if i < n_select:
+                    array_select.append(array[idcs[i]])
+                else:
+                    array_discard.append(array[idcs[i]])
+        else:
+            raise NotImplementedError("<subsample_data_array> Array type other than list")
+    elif method == 'stride':
+        # Construct index sets
+        idcs = np.arange(0, n_data)
+        idcs_sel = [ int(float(idcs.shape[0])/n_select*i) for i in range(n_select) ]
+        idcs_sel = np.array(idcs_sel)
+        # Take into account periodic shift
+        if stride_shift != 0:
+            idcs_sel = np.sort((idcs_sel+stride_shift) % idcs.shape[0])
+        # Index complement
+        mask = np.zeros(idcs.shape[0], dtype=bool)
+        mask[idcs_sel] = True
+        idcs_disc = idcs[~mask]
+        # Create selections
+        if type(array) == list:
+            array_select = []
+            array_discard = []
+            for i in idcs_sel:
+                array_select.append(array[i])
+            for i in idcs_disc:
+                array_discard.append(array[i])
+        else:
+            raise NotImplementedError("<subsample_data_array> Array type other than list")
+    return array_select, array_discard
+
+if __name__ == "__main__":
+    print "Data array"
+    a = range(20)
+    print a
+    print "Select with stride, shift=0"
+    a1, a2 = subsample_array(a, 5, 'stride', 0)
+    print a1, a2
+    print "Select with stride, shift=6"
+    a1, a2 = subsample_array(a, 5, 'stride', 6)
+    print a1, a2
+    print "Select random subsample"
+    a1, a2 = subsample_array(a, 5, 'random')
+    print a1, a2
+
+
+

src/soap/soapy/util.py

@@ -11,6 +11,30 @@ HARTREE_TO_KCALMOL = 627.509469
 
 MP_LOCK = mp.Lock()
 
+def mp_pool_compute_upper_triangle(
+        kfct,
+        g_list,
+        n_procs,
+        dtype='float64',
+        mplog=None,
+        **kwargs):
+    kfct_primed=fct.partial(kfct, **kwargs)
+    n_rows = len(g_list)
+    kmat = np.zeros((n_rows, n_rows), dtype=dtype)
+    for i in range(n_rows):
+        if mplog: mplog << mplog.back << "Computing row %d" % i << mplog.endl
+        g_pair_list = []
+        gi = g_list[i]
+        for j in range(i, n_rows):
+            gj = g_list[j]
+            g_pair_list.append([gi,gj])
+        pool = mp.Pool(processes=n_procs)
+        krow = pool.map(kfct_primed, g_pair_list)
+        pool.close()
+        pool.join()
+        kmat[i,i:] = krow
+    return kmat
+
 def mp_compute_vector(
         kfct,
         g_list,
@@ -42,12 +66,30 @@ def mp_compute_column_block(gi, gj_list, kfct):
         krow.append(k)
     return krow
 
+def compute_upper_triangle(
+        kfct,
+        g_list,
+        **kwargs):
+    dim = len(g_list)
+    kmat = np.zeros((dim,dim), dtype='float64')
+    kfct_primed = fct.partial(
+        kfct,
+        **kwargs)
+    for i in range(dim):
+        gi = g_list[i]
+        for j in range(i, dim):
+            gj = g_list[j]
+            kij = kfct(gi, gj)
+            kmat[i,j] = kij
+            kmat[j,i] = kij
+    return kmat
+
 def mp_compute_upper_triangle(
         kfct, 
         g_list, 
         n_procs, 
         n_blocks, 
-        log=None, 
+        mplog=None, 
         tstart_twall=(None,None), 
         backup=True,
         verbose=True,
@@ -63,7 +105,7 @@ def mp_compute_upper_triangle(
     n_blocks: number of column blocks onto which computation is split
     kwargs: keyword arguments supplied to kfct
     """
-    if not verbose: log=None
+    if not verbose: mplog=None
     t_start = tstart_twall[0]
     t_wall = tstart_twall[1]
     dim = len(g_list)
@@ -78,7 +120,7 @@ def mp_compute_upper_triangle(
         # Column start, column end
         c0 = col_div[0]
         c1 = col_div[-1]+1
-        if log: log << "Column block i[%d:%d] j[%d:%d]" % (0, c1, c0, c1) << log.endl
+        if mplog: mplog << "Column block i[%d:%d] j[%d:%d]" % (0, c1, c0, c1) << mplog.endl
         gj_list = g_list[c0:c1]
         gi_list = g_list[0:c1]
         # Prime kernel function
@@ -95,7 +137,7 @@ def mp_compute_upper_triangle(
         npyfile = 'out.block_i_%d_%d_j_%d_%d.npy' % (0, c1, c0, c1)
         # ... but first check for previous calculations of same slice
         if backup and npyfile in os.listdir('./'):
-            if log: log << "Load block from '%s'" % npyfile << log.endl
+            if mplog: mplog << "Load block from '%s'" % npyfile << mplog.endl
             kmat_column_block = np.load(npyfile)
         else:
             kmat_column_block = pool.map(mp_compute_column_block_primed, gi_list)
@@ -110,9 +152,9 @@ def mp_compute_upper_triangle(
         dt_block = t_out-t_in
         if t_start and t_wall:
             t_elapsed = t_out-t_start
-            if log: log << "Time elapsed =" << t_elapsed << " (wall time = %s) (maxmem = %d)" % (t_wall, resource.getrusage(resource.RUSAGE_SELF).ru_maxrss) << log.endl
+            if mplog: mplog << "Time elapsed =" << t_elapsed << " (wall time = %s) (maxmem = %d)" % (t_wall, resource.getrusage(resource.RUSAGE_SELF).ru_maxrss) << mplog.endl
             if col_div_idx+1 != len(col_div_list) and t_elapsed+dt_block > t_wall-dt_block:
-                log << "Wall time hit expected for next iteration, break ..." << log.endl
+                mplog << "Wall time hit expected for next iteration, break ..." << mplog.endl
                 break
             else: pass
     return kmat

src/soap/spectrum.cpp

@@ -47,6 +47,10 @@ void Spectrum::compute() {
     this->compute(_structure->particles(), _structure->particles());
 }
 
+void Spectrum::compute(Segment *center) {
+    this->compute(center->particles(), _structure->particles());
+}
+
 void Spectrum::compute(Segment *center, Segment *target) {
     this->compute(center->particles(), target->particles());
 }
@@ -139,7 +143,7 @@ AtomicSpectrum *Spectrum::computeAtomic(Particle *center, Structure::particle_ar
             weight0 *= _basis->getCutoff()->getCenterWeight();
         }
 
-        GLOG() << (*pit)->getType() << " " << dr.getX() << " " << dr.getY() << " " << dr.getZ() << std::endl;
+        GLOG() << (*pit)->getType() << " " << dr.getX() << " " << dr.getY() << " " << dr.getZ() << " " << (*pit)->getWeight() << std::endl;
 
         // COMPUTE EXPANSION & ADD TO SPECTRUM
         bool gradients = (is_image) ? false : _options->get<bool>("spectrum.gradients");
@@ -275,6 +279,7 @@ void Spectrum::load(std::string archfile) {
 void Spectrum::registerPython() {
     using namespace boost::python;
     void (Spectrum::*computeAll)() = &Spectrum::compute;
+    void (Spectrum::*computeSeg)(Segment*) = &Spectrum::compute;
     void (Spectrum::*computeSegPair)(Segment*, Segment*) = &Spectrum::compute;
     void (Spectrum::*computeCentersTargets)(Structure::particle_array_t&, Structure::particle_array_t&) = &Spectrum::compute;
 
@@ -283,6 +288,7 @@ void Spectrum::registerPython() {
     	.def(init<std::string>())
     	.def("__iter__", range<return_value_policy<reference_existing_object> >(&Spectrum::beginAtomic, &Spectrum::endAtomic))
 	    .def("compute", computeAll)
+        .def("compute", computeSeg)
 	    .def("compute", computeSegPair)
 	    .def("compute", computeCentersTargets)
 		.def("computePower", &Spectrum::computePower)

src/soap/spectrum.hpp

@@ -46,6 +46,7 @@ public:
 	void clean();
 
 	void compute();
+    void compute(Segment *centers);
 	void compute(Segment *centers, Segment *targets);
 	void compute(Structure::particle_array_t &sources, Structure::particle_array_t &targets);
 	AtomicSpectrum *computeAtomic(Particle *center);

src/soap/tools/CMakeLists.txt

-install(FILES __init__.py loadwrite.py inverse.py extract.py DESTINATION ${LOCAL_INSTALL_DIR}/tools)
+install(FILES __init__.py loadwrite.py inverse.py extract.py partition.py DESTINATION ${LOCAL_INSTALL_DIR}/tools)
 

src/soap/tools/__init__.py

 from loadwrite import *
 from inverse import *
 from extract import *
+import partition
 

src/soap/tools/loadwrite.py

 import os
 import numpy as np
+import itertools
+import partition
 from .. import _soapxx as soap
 from ..soapy import momo
-
+from ..soapy import elements
 
 try:
     import ase
@@ -10,6 +12,169 @@ try:
 except ImportError:
     print("Note: ase.io import failed. Install PYTHON-ASE to harvest full reader functionality.")
 
+def structures_from_xyz(xyz_file, **kwargs):
+    # Read xyz via ASE
+    ase_configs = ase.io.read(xyz_file, index=':')
+    return structures_from_ase(
+        ase_configs=ase_configs, 
+        **kwargs)
+
+def structures_from_ase(
+        configs, 
+        return_all=False, 
+        **kwargs):
+    configs_mod = []
+    structs = []
+    return_tuples = []
+    for config in configs:
+        return_tuple = \
+            structure_from_ase(
+                config,
+                **kwargs)
+        configs_mod.append(return_tuple[0])
+        structs.append(return_tuple[1])
+    if return_all:
+        return return_tuples
+    else:
+        return configs_mod, structs
+
+def structure_from_ase(
+        config, 
+        do_partition=False, 
+        add_fragment_com=False, 
+        use_center_of_geom=False, 
+        log=None):
+    # NOTE Center of mass is computed without considering PBC => Requires unwrapped coordinates
+    structure = None
+    frag_bond_matrix = None
+    atom_bond_matrix = None
+    frag_labels = []
+    atom_labels = []
+    # System properties
+    label = config.info['label']
+    positions = config.get_positions()
+    types = config.get_chemical_symbols()
+    if log: log << log.back << "Reading '%s'" % label << log.flush
+    # Simulation cell
+    if config.pbc.all(): 
+        box = np.array([config.cell[0], config.cell[1], config.cell[2]])
+    elif not config.pbc.any(): 
+        box = np.zeros((3,3))
+    else: 
+        raise NotImplementedError("<structures_from_xyz> Partial periodicity not implemented.")
+    # Partition 
+    if do_partition:
+        # Partition => Frags, top, reordered positions
+        frags, top = partition.PartitionStructure(types, positions, outfile_gro='%s.gro' % label)
+        atms = [ atm for atm in itertools.chain(*frags) ]
+        positions = [ atm.xyz for atm in itertools.chain(*frags) ]
+        types = [ atm.e for atm in itertools.chain(*frags) ]
+        # Generate fragment labels
+        for frag in frags:
+            label = '-'.join(atm.e for atm in frag)
+            frag_labels.append(label)
+        # Reorder map
+        id_reorder_map = {}
+        for atm in atms:
+            id_reorder_map[atm.id_initial] = atm.id
+        # Connectivity matrix: fragments
+        frag_bond_matrix = np.zeros((len(frags),len(frags)), dtype=bool)
+        for i in range(len(frags)):
+            frag_bond_matrix[i,i] = True
+            for j in range(i+1, len(frags)):
+                frags_are_bonded = False
+                for ai in frags[i]:
+                    for aj in frags[j]:
+                        if aj in ai.bonded:
+                            frags_are_bonded = True
+                            break
+                    if frags_are_bonded: break
+                frag_bond_matrix[i,j] = frags_are_bonded
+                frag_bond_matrix[j,i] = frags_are_bonded
+        # Connectivity matrix: atoms
+        atom_bond_matrix = np.zeros((len(positions),len(positions)), dtype=bool)
+        for i in range(len(atms)):
+            atom_bond_matrix[i,i] = True
+            ai = atms[i]
+            for j in range(i+1, len(atms)):
+                atoms_are_bonded = False
+                aj = atms[j]
+                if aj in ai.bonded:
+                    atoms_are_bonded = True
+                atom_bond_matrix[i,j] = atoms_are_bonded
+                atom_bond_matrix[j,i] = atoms_are_bonded
+        # Reorder ASE atoms
+        config = ase.Atoms(
+            sorted(config, key = lambda atm: id_reorder_map[atm.index+1]), 
+            info=config.info, 
+            cell=config.cell, 
+            pbc=config.pbc)
+    else:
+        top = [('SEG', 1, len(positions))]
+    # Check particle count consistent with top
+    atom_count = 0
+    for section in top:
+        atom_count += section[1]*section[2]
+    assert atom_count == len(positions) # Does topology match structure?
+    # Create segments, particles
+    structure = soap.Structure(label)
+    structure.box = box
+    atom_idx = 0
+    for section in top:
+        seg_type = section[0]
+        n_segs = section[1]
+        n_atoms = section[2]
+        for i in range(n_segs):
+            segment = structure.addSegment()
+            segment.name = seg_type
+            # Add particles, compute CoM
+            com = np.array([0.,0.,0.])
+            com_weight_total = 0
+            for j in range(n_atoms):
+                particle = structure.addParticle(segment)
+                particle.pos = positions[atom_idx]
+                particle.weight = 1.
+                particle.sigma = 0.5
+                particle.type = types[atom_idx]
+                atom_labels.append(particle.type)
+                # Compute CoMass/CoGeom
+                if use_center_of_geom:
+                    com_weight = 1.
+                else:
+                    com_weight = elements.periodic_table[particle.type].mass
+                com_weight_total += com_weight
+                com = com + com_weight*positions[atom_idx]
+                atom_idx += 1
+            com = com/com_weight_total
+            # Add CoM particle if requested
+            if add_fragment_com:
+                segment = structure.addSegment()
+                segment.name = "%s.COM" % seg_type
+                particle = structure.addParticle(segment)
+                particle.pos = com
+                particle.weight = 0.
+                particle.sigma = 0.5
+                particle.type = "COM"
+    if log: log << log.endl
+    return config, structure, top, frag_bond_matrix, atom_bond_matrix, frag_labels, atom_labels
+
+def join_structures_as_segments(structs):
+            # NOTE Segments of the individual structures will be lost
+            # NOTE Box will be adopted from first structure in list
+            label = ':'.join([s.label for s in structs])
+            struct = soap.Structure(label)
+            struct.box = structs[0].box
+            for struct_i in structs:
+                seg = struct.addSegment()
+                seg.name = struct_i.label
+                for p_i in struct_i.particles:
+                    p = struct.addParticle(seg)
+                    p.pos = p_i.pos
+                    p.weight = p_i.weight
+                    p.sigma = p_i.sigma
+                    p.type = p_i.type
+            return struct
+
 def write_xyz(xyz_file, structure):
     ofs = open(xyz_file, 'w')
     ofs.write('%d\n%s\n' % (len(list(structure.particles)), structure.label))