Systep exploration.

test/test_kernelpot/kernel.py

@@ -150,9 +150,41 @@ class KernelFunctionDot(object):
             return (abs(1.-K))**0.5
         else:
             return K
-   
+    def computeBlockDot(self, IX, return_distance=False):
+        return self.computeBlock(IX, return_distance)
+
+class KernelFunctionDotHarmonic(object):
+    """
+    C_i = ( d^2 [ IX.X ]^xi - mu_i )^2
+    """
+    def __init__(self, options):
+        self.delta = float(options.get('kernel.delta'))
+        self.xi = float(options.get('kernel.xi'))
+        self.mu = float(options.get('kernel.mu'))
+        self.kfctdot = KernelFunctionDot(options)
+        return
+    def compute(self, IX, X):
+        if type(self.mu) == float:
+            mu = np.zeros((IX.shape[0]))
+            np.fill(mu, self.mu)
+            self.mu = mu
+        C = self.kfctdot.compute(IX, X)
+        return (C-self.mu)**2
+    def computeDerivativeOuter(self, IX, X):
+        if type(self.mu) == float:
+            mu = np.zeros((IX.shape[0]))
+            mu.fill(self.mu)
+            self.mu = mu
+        C = self.kfctdot.compute(IX, X)
+        IC = self.kfctdot.computeDerivativeOuter(IX, X)
+        return 2*(C-self.mu)*IC
+    def computeBlockDot(self, IX, return_distance=False):
+        return self.kfctdot.computeBlock(IX, return_distance)
+        
+        
+
 KernelAdaptorFactory = { 'generic': KernelAdaptorGeneric, 'global-generic': KernelAdaptorGlobalGeneric }     
-KernelFunctionFactory = { 'dot':KernelFunctionDot }
+KernelFunctionFactory = { 'dot':KernelFunctionDot, 'dot-harmonic': KernelFunctionDotHarmonic }
 
 class KernelPotential(object):
     def __init__(self, options):
@@ -174,6 +206,14 @@ class KernelPotential(object):
         # INCLUSIONS / EXCLUSIONS
         # -> Already be enforced when computing spectra
         return
+    def computeKernelMatrix(self, return_distance=False):
+        return self.kernelfct.computeBlock(self.IX, return_distance)
+    def computeDotKernelMatrix(self, return_distance=False):
+        # Even if the, e.g., dot-harmonic kernel is used for optimization, 
+        # the dot-product based kernel matrix may still be relevant
+        return self.kernelfct.computeBlockDot(self.IX, return_distance)
+    def nConfigs(self):
+        return self.IX.shape[0]
     def importAcquire(self, IX_acqu, alpha):
         n_acqu = IX_acqu.shape[0]
         dim_acqu = IX_acqu.shape[1]
@@ -285,7 +325,7 @@ class KernelPotential(object):
             logging.info("  Center %d" % (pid))
             # neighbour-pid-independent kernel "prevector" (outer derivative)
             X_unnorm, X_norm = self.adaptor.adaptScalar(atomic)
-            dIC = self.kernelfct.computeDerivativeOuter(self.IX, X_norm) # TODO This must be X_norm!
+            dIC = self.kernelfct.computeDerivativeOuter(self.IX, X_norm)
             alpha_dIC = self.alpha.dot(dIC)
             for nb_pid in nb_pids:
                 # Force on neighbour
@@ -364,7 +404,7 @@ def apply_force_norm_step(structure, forces, scale, constrain_particles=[]):
         part.pos = part.pos + scale*forces[idx]
     return [ part.pos for part in structure ]
 
-def evaluate_energy(positions, structure, kernelpot, opt_pids, verbose=False, ofs=None):
+def evaluate_energy(positions, structure, kernelpot, opt_pids, verbose=False, ofs=None, average=False):
     if verbose: print "Energy"
     # Impose positions
     pid_pos = positions.reshape((opt_pids.shape[0],3))
@@ -376,13 +416,14 @@ def evaluate_energy(positions, structure, kernelpot, opt_pids, verbose=False, of
         if verbose: print part.id, part.type, part.pos
     # Evaluate energy function
     energy = kernelpot.computeEnergy(structure)
+    if average: energy /= kernelpot.nConfigs()
     # Log
     if ofs: ofs.logFrame(structure)
     if verbose: print energy
     print energy
     return energy
 
-def evaluate_energy_gradient(positions, structure, kernelpot, opt_pids, verbose=False, ofs=None):
+def evaluate_energy_gradient(positions, structure, kernelpot, opt_pids, verbose=False, ofs=None, average=False):
     if verbose: print "Forces"
     # Adjust positions
     pid_pos = positions.reshape((opt_pids.shape[0],3))
@@ -398,6 +439,8 @@ def evaluate_energy_gradient(positions, structure, kernelpot, opt_pids, verbose=
     opt_pidcs = opt_pids-1
     gradients_short = gradients[opt_pidcs]
     gradients_short = gradients_short.flatten()
+    if average:
+        gradients_short /= kernelpot.nConfigs()
     if verbose: print gradients_short
     #forces[2] = 0. # CONSTRAIN TO Z=0 PLANE
     return gradients_short

test/test_sampling/0_sample.py

+#! /usr/bin/env python
+import soap
+import soap.tools
+soap.silence()
+
+import os
+import sys
+import numpy as np
+import logging
+import io
+
+from momo import osio, endl, flush
+from kernel import KernelPotential, TrajectoryLogger
+from kernel import perturb_positions, random_positions
+from kernel import evaluate_energy, evaluate_energy_gradient
+from pca import PCA, IPCA
+from dimred import dimred
+
+logging.basicConfig(
+    format='[%(asctime)s] %(message)s',
+    datefmt='%I:%M:%S',
+    level=logging.ERROR)
+
+# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+select_struct_label = 'config.xyz'
+exclude_centers = []
+adaptor_type = 'global-generic'
+kernel_type = 'dot-harmonic'
+alpha = +1.
+mu = 0.9
+sc = 0.1 # 1.5 # 0.1
+average_energy = True
+lj_sigma = 0.02 # 0.00001 # 0.02
+# <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+# NOTE adaptor_type = 'generic' => exclude center if not constrained in optimization
+
+def sample_random_x0(struct, exclude_pids_rnd, opt_pids, b0, b1):
+    # Sample
+    x0 = [ part.pos for part in struct ]
+    x0 = random_positions(struct, exclude_pid=exclude_pids_rnd, b0=b0, b1=b1)
+    x0 = np.array(x0)    
+    # Short-list
+    opt_pids = np.array(opt_pids)
+    opt_pidcs = opt_pids-1
+    x0 = x0[opt_pidcs]
+    x0 = x0.flatten()    
+    return x0
+
+class LJRepulsive(object):
+    def __init__(self, sigma=lj_sigma):
+        self.sigma = sigma            
+    def computeEnergy(self, struct):
+        coords = np.zeros((struct.n_particles, 3))
+        for idx, part in enumerate(struct):
+            coords[idx,:] = part.pos            
+        E = 0.
+        for i in range(struct.n_particles):
+            for j in range(i):
+                r = np.sqrt(np.sum((coords[i, :] - coords[j, :]) ** 2))
+                E += (self.sigma/r)**12
+        return E
+    def computeGradient(self, struct):
+        coords = np.zeros((struct.n_particles, 3))
+        for idx, part in enumerate(struct):
+            coords[idx,:] = part.pos            
+        grad = np.zeros(coords.shape)
+        for i in range(struct.n_particles):
+            for j in range(i):
+                dr = coords[i, :] - coords[j, :]
+                r = np.sqrt(np.sum(dr**2))
+                g = 12./self.sigma*(self.sigma/r)**13
+                grad[i, :] += -g * dr / r
+                grad[j, :] += g * dr / r
+        return grad
+
+def evaluate_energy_local(positions, structure, kernelpot, opt_pids, verbose=False, ofs=None, average=False):
+    if verbose: print "Energy"
+    # Impose positions
+    pid_pos = positions.reshape((opt_pids.shape[0],3))
+    for pidx, pid in enumerate(opt_pids):
+        pos = pid_pos[pidx,:]
+        particle = structure.getParticle(pid)        
+        particle.pos = pos        
+    for part in structure:
+        if verbose: print part.id, part.type, part.pos
+    # Evaluate energy function
+    energy = kernelpot.computeEnergy(structure)
+    if average: energy /= kernelpot.nConfigs()
+    # Evaluate additional potentials
+    lj = LJRepulsive()
+    energy_add = lj.computeEnergy(struct)
+    energy += energy_add
+    # Log
+    if ofs: ofs.logFrame(structure)
+    if verbose: print energy
+    print energy
+    return energy
+
+def evaluate_energy_gradient_local(positions, structure, kernelpot, opt_pids, verbose=False, ofs=None, average=False):
+    if verbose: print "Forces"
+    # Adjust positions
+    pid_pos = positions.reshape((opt_pids.shape[0],3))
+    for pidx, pid in enumerate(opt_pids):
+        pos = pid_pos[pidx,:]
+        particle = structure.getParticle(pid)        
+        particle.pos = pos        
+    for part in structure:
+        if verbose: print part.id, part.type, part.pos
+    # Evaluate forces
+    forces = kernelpot.computeForces(structure)
+    gradients = -1.*np.array(forces)
+    if average:
+        gradients = gradients/kernelpot.nConfigs()
+    # Evaluate additional potentials
+    lj = LJRepulsive()
+    gradients_add = lj.computeGradient(struct)
+    gradients = gradients + gradients_add
+    # Short-list
+    opt_pidcs = opt_pids-1
+    gradients_short = gradients[opt_pidcs]
+    gradients_short = gradients_short.flatten()
+    if verbose: print gradients_short
+    #forces[2] = 0. # CONSTRAIN TO Z=0 PLANE
+    return gradients_short
+
+# OPTIONS
+options = soap.Options()
+options.excludeCenters(['H'])
+options.excludeTargets(['H'])
+options.excludeCenterIds(exclude_centers)
+options.excludeTargetIds([])
+# Spectrum
+options.set('spectrum.gradients', True)
+options.set('spectrum.2l1_norm', False)                     # <- pull here (True/False)
+# Basis
+options.set('radialbasis.type', 'gaussian')
+options.set('radialbasis.mode', 'adaptive')
+options.set('radialbasis.N', 9)
+options.set('radialbasis.sigma', 0.5)
+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.5)
+options.set('angularbasis.type', 'spherical-harmonic')
+options.set('angularbasis.L', 6)
+# Kernel
+options.set('kernel.adaptor', adaptor_type)                 # <- pull here (generic/global-generic)
+options.set('kernel.type', kernel_type)
+options.set('kernel.delta', 1.)
+options.set('kernel.xi', 2.)                                # <- pull here (1., ..., 4., ...)
+options.set('kernel.mu', mu)
+# Cube files
+options.set('densitygrid.N', 20)
+options.set('densitygrid.dx', 0.15)
+
+
+# LOAD STRUCTURES
+structures = [ 
+    soap.tools.setup_structure_ase(c.config_file, c.atoms)
+    for c in soap.tools.ase_load_all('in.configs')
+]
+struct_dict = {}
+for struct in structures:
+    print struct.label
+    struct_dict[struct.label] = struct
+
+# SELECT STRUCTURE
+struct = struct_dict[select_struct_label]
+options.excludeCenterIds(exclude_centers)
+
+# SETUP KERNEL (=> BASIS, KERNELPOT WITH ADAPTOR)
+kernelpot_buffer = KernelPotential(options)
+kernelpot_buffer.acquire(struct, alpha)
+kernelpot = KernelPotential(options)
+kernelpot.acquire(struct, alpha) # <- repulsive
+np.savetxt('out.x_ref.txt', kernelpot.IX)
+
+# MANAGE EXCLUSIONS (RANDOMIZATION, OPTIMIZATION)
+exclude_pids_rnd = [1]
+opt_pids = [ i for i in range(2, struct.n_particles+1) ]
+opt_pids = np.array(opt_pids)
+#if adaptor_type == 'global-generic':
+#    exclude_pids_rnd = []
+#    opt_pids = [ i for i in range(1,struct.n_particles+1) ]
+#elif adaptor_type == 'generic':
+#    ...
+
+# TRAJECTORY LOGGER
+ofs_opt = TrajectoryLogger('out.opt.xyz')
+ofs_opt.logFrame(struct)
+ofs_acquired = TrajectoryLogger('out.acq.xyz')
+ofs_acquired.logFrame(struct)
+ofs_warned = TrajectoryLogger('out.warn.xyz')
+
+# OPTIMIZER INTERFACE
+import scipy.optimize
+f = evaluate_energy_local
+fprime = evaluate_energy_gradient_local
+args = (struct, kernelpot, opt_pids, False, ofs_opt, average_energy)
+
+print "Select label:    ", select_struct_label
+print "Exclude centers: ", exclude_centers
+print "Adaptor type:    ", adaptor_type
+
+pos_opt_list = []
+e_opt_list = []
+
+
+
+
+for mu in [ 0.9, 0.7, 0.5, 0.3 ]:
+    kernelpot.kernelfct.mu = mu   
+    
+    
+    n_iterations = 10
+    skip_warnings = False
+
+    iter_count = 0
+    while iter_count < n_iterations:
+        # RANDOMIZE INITIAL CONFIG & SETUP X0
+        x0 = sample_random_x0(struct, exclude_pids_rnd, opt_pids, -1.*sc, +1.*sc)    
+        print "Positions after initial randomization"
+        for part in struct:
+            print part.id, part.pos
+        print "Optimization arguments (initial)"
+        for pos in x0.reshape((len(opt_pids),3)):
+            print pos
+        print type(x0)
+        
+
+        x0_opt, f_opt, n_calls, n_grad_calls, warnflag = scipy.optimize.fmin_cg(
+            f=f, 
+            x0=x0, 
+            fprime=fprime, 
+            args=args, 
+            gtol=1e-6,
+            full_output=True)
+        if warnflag > 0:
+            print "Warnflag =", warnflag, " => try again."
+            ofs_warned.logFrame(struct)
+            if skip_warnings: continue
+        
+        print "Acquire structure ..."
+        ofs_acquired.logFrame(struct)
+        #kernelpot.acquire(struct, alpha)
+        kernelpot_buffer.acquire(struct, alpha)
+        iter_count += 1
+        print "Kernel size = %d" % kernelpot.nConfigs()
+        
+        e, prj = kernelpot.computeEnergy(struct, return_prj_mat=True)
+        print prj
+        
+        pos = np.array([ part.pos for part in struct ])
+        print "Positions"
+        print pos
+        pos_opt_list.append(pos)
+        e_opt_list.append(f_opt)
+        
+        #raw_input('...')
+
+# CLOSE LOGGER
+ofs_opt.close()
+ofs_acquired.close()
+ofs_warned.close()
+
+np.savetxt('out.e_opt_list.txt', np.array(e_opt_list))
+np.savetxt('out.kernelpot.kernel.txt', kernelpot.computeDotKernelMatrix())
+np.savetxt('out.kernelpot.buffer.kernel.txt', kernelpot_buffer.computeDotKernelMatrix())
+np.savetxt('out.kernelpot.buffer.ix.txt', kernelpot_buffer.IX)
+np.savetxt('out.kernelpot.ix.txt', kernelpot.IX)
+
+if kernelpot.nConfigs() > 1: dimred(kernelpot, 'mds', 'out.dimred.txt', symmetrize=True)
+if kernelpot_buffer.nConfigs() > 1: dimred(kernelpot_buffer, 'mds', 'out.dimred.buffer.txt', symmetrize=True)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

test/test_sampling/dimred.py

+#! /usr/bin/env python
+import numpy as np
+import sklearn.manifold
+import sklearn.decomposition
+
+def dimred_matrix(method, kmat=None, distmat=None, outfile=None, ix=None, symmetrize=False):
+    if symmetrize:
+        kmat = 0.5*(kmat+kmat.T)
+        dmat = 0.5*(dmat+dmat.T)
+    if method == 'mds':
+        # MDS
+        # http://scikit-learn.org/stable/modules/generated/sklearn.manifold.MDS.html#sklearn.manifold.MDS
+        mds = sklearn.manifold.MDS(
+            n_components=2,
+            metric=False,
+            verbose=1,
+            n_init=10,
+            dissimilarity='precomputed')
+        positions = mds.fit_transform(distmat)
+    elif method == 'isomap':
+        isomap = sklearn.manifold.Isomap(
+            n_neighbors=5,
+            n_components=3,
+            eigen_solver='auto',
+            path_method='auto',
+            neighbors_algorithm='auto')
+        positions = isomap.fit_transform(ix)
+    elif method == 'kernelpca':
+        kernelpca = sklearn.decomposition.KernelPCA(
+            n_components=None,
+            kernel='precomputed',
+            eigen_solver='auto',
+            max_iter=None,
+            remove_zero_eig=True)
+        positions = kernelpca.fit_transform(kmat)
+    else: raise NotImplementedError(method)
+    if outfile: np.savetxt(outfile, positions)
+    return positions
+
+def dimred(kernelpot, method, outfile, symmetrize=False):
+    ix = kernelpot.IX
+    kmat = kernelpot.kernelfct.computeBlockDot(ix, return_distance=False)
+    distmat = kernelpot.kernelfct.computeBlockDot(ix, return_distance=True)
+    positions = dimred_matrix(method, kmat, distmat, outfile, ix, symmetrize=symmetrize)
+    return positions
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

test/test_sampling/in.configs/config.xyz

+3
+
+C 0.0 0.0 0.0
+C 1.5 0.0 0.0
+C 0.0 1.5 0.0

test/test_sampling/symmetrize.py

+import numpy as np
+
+IX = np.loadtxt('out.kernelpot.buffer.ix.txt')
+K = np.loadtxt('out.kernelpot.buffer.kernel.txt')
+D = (abs(1.-K))**0.5
+
+Dsym = 0.5*(D+D.T)
+Ksym = 0.5*(K+K.T)
+
+from dimred import dimred_matrix
+dimred_matrix('mds', kmat=K, distmat=Dsym, outfile='tmp.txt')
+#dimred_matrix('kernelpca', kmat=K, distmat=Dsym, outfile='tmp.txt')
+#dimred_matrix('isomap', kmat=Ksym, distmat=Dsym, outfile='tmp.txt', ix=IX, symmetrize=False)