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

Eigstruct and NLDR tests.

parent 5e1b95c2
......@@ -142,7 +142,14 @@ class KernelFunctionDot(object):
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)
return self.xi*np.diag(c).dot(IX)
def computeBlock(self, IX, return_distance=False):
# Order such that: IX[i] => fingerprint observation i
K = self.delta**2 * IX.dot(IX.T)**self.xi
if return_distance:
return (abs(1.-K))**0.5
else:
return K
KernelAdaptorFactory = { 'generic': KernelAdaptorGeneric, 'global-generic': KernelAdaptorGlobalGeneric }
KernelFunctionFactory = { 'dot':KernelFunctionDot }
......@@ -156,6 +163,7 @@ class KernelPotential(object):
self.IX = None
self.alpha = None
self.dimX = None # <- Also used as flag set by first ::acquire call
self.labels = []
# KERNEL
logging.info("Choose kernel function ...")
self.kernelfct = KernelFunctionFactory[options.get('kernel.type')](options)
......@@ -188,7 +196,7 @@ class KernelPotential(object):
#print self.IX
logging.info("Imported %d environments." % n_acqu)
return
def acquire(self, structure, alpha):
def acquire(self, structure, alpha, label=None):
logging.info("Acquire ...")
spectrum = soap.Spectrum(structure, self.options, self.basis)
spectrum.compute()
......@@ -204,6 +212,9 @@ class KernelPotential(object):
# New alpha's
alpha_acqu = np.zeros((n_acqu))
alpha_acqu.fill(alpha)
# Labels
labels = [ label for i in range(n_acqu) ]
self.labels = self.labels + labels
if not self.dimX:
# First time ...
self.dimX = dim_acqu
......@@ -367,7 +378,8 @@ def evaluate_energy(positions, structure, kernelpot, opt_pids, verbose=False, of
energy = kernelpot.computeEnergy(structure)
# Log
if ofs: ofs.logFrame(structure)
if verbose: print energy
if verbose: print energy
print energy
return energy
def evaluate_energy_gradient(positions, structure, kernelpot, opt_pids, verbose=False, ofs=None):
......
This diff is collapsed.
#! /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, KernelFunctionFactory, TrajectoryLogger
from kernel import perturb_positions, random_positions
from kernel import evaluate_energy, evaluate_energy_gradient
# KERNEL OPTIONS
options = soap.Options()
options.set('kernel.type', 'dot')
options.set('kernel.delta', 1.)
options.set('kernel.xi', 4.)
method = 'kernelpca'
# KERNEL FUNCTION
kernelfct = KernelFunctionFactory[options.get('kernel.type')](options)
# DISTANCE MATRIX
ix = np.loadtxt('in.kernelpot.ix.txt')
kmat = kernelfct.computeBlock(ix, return_distance=False)
distmat = kernelfct.computeBlock(ix, return_distance=True)
# PROJECTION
import sklearn.manifold
import sklearn.decomposition
if method == 'mds':
# MDS
# http://scikit-learn.org/stable/modules/generated/sklearn.manifold.MDS.html#sklearn.manifold.MDS
mds = sklearn.manifold.MDS(
n_components=3,
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)
# SAVE POSITIONS
np.savetxt('out.positions.txt', positions)
......@@ -70,7 +70,7 @@ generate = False
if generate:
for struct in structures:
print struct.label
kernelpot.acquire(struct, 1.)
kernelpot.acquire(struct, 1., label=struct.label)
print kernelpot.IX.shape
np.savetxt('out.kernelpot.ix.txt', kernelpot.IX)
else:
......@@ -137,7 +137,8 @@ print np.dot(X_eig_0, X_eig_0)
print np.dot(X_eig_0, X_eig_1)
print np.dot(X_eig_1, X_eig_1)
X_eig_vecs = pca.reconstructEigenvecs()
np.savetxt('out.pca.xeigenvecs.txt', X_eig_vecs)
......
#! /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
# Load XX from text file
# Run PCA (normalize, diagonalize)
# Project structures onto PCA => plot
# TODO Reconstruct eigenstructures (unnormalize, load into kernel, vary centers, optimize)
logging.basicConfig(
format='[%(asctime)s] %(message)s',
datefmt='%I:%M:%S',
level=logging.ERROR)
# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
idx_eig_target = 199
N = 5vmd
select_struct_label = 'config_C_%02d.xyz' % N
exclude_pids_rnd = [1]
exclude_centers = [ i for i in range(2, N+1) ]
opt_pids = [ i for i in range(2, N+1) ]
x_array_file = 'out.pca.xeigenvecs.txt'
x_array_file = 'out.kernelpot.ix.txt'
# <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
# 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', 'generic') # <- pull here (generic/global-generic)
options.set('kernel.type', 'dot')
options.set('kernel.delta', 1.)
options.set('kernel.xi', 4.) # <- pull here (1., ..., 4., ...)
# 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('configs_n_carbons')
]
struct_dict = {}
for struct in structures:
print struct.label
struct_dict[struct.label] = struct
# LOAD EIGENSTRUCTURES
print "Load fingerprints from '%s'" % x_array_file
IX_eig = np.loadtxt(x_array_file)
x1 = IX_eig[0]
x2 = IX_eig[4]
print np.dot(x1,x1), np.dot(x1,x2), np.dot(x2,x2)
# SELECT STRUCTURE, X_TARGET
print "Target eigenstructure with index", idx_eig_target
print "Exclude centers", exclude_centers
print "Select label", select_struct_label
X_target = IX_eig[idx_eig_target]
struct = struct_dict[select_struct_label]
options.excludeCenterIds(exclude_centers)
# SETUP KERNEL (=> BASIS, KERNELPOT WITH ADAPTOR)
basis = soap.Basis(options)
kernelpot = KernelPotential(options)
kernelpot.importAcquire(np.array([X_target]), -1.) # <- attractive
sc = 2.
# RANDOMIZE INITIAL CONFIG
positions_0 = [ part.pos for part in struct ]
positions_0_rnd = random_positions(struct, exclude_pid=exclude_pids_rnd, b0=-1.*sc, b1=1.*sc)
print "Positions after initial randomization"
for part in struct:
print part.id, part.pos
# EXTRACT POSITIONS (OPTIMIZATION TARGETS)
opt_pids = np.array(opt_pids)
opt_pidcs = opt_pids-1
positions = np.array(positions_0_rnd)
positions_short = positions[opt_pidcs]
positions_short = positions_short.flatten()
print "Optimization arguments (initial)"
for pos in positions_short.reshape((len(opt_pids),3)):
print pos
raw_input('...')
# TRAJECTORY LOGGER
ofs = TrajectoryLogger('opt.xyz')
ofs.logFrame(struct)
# OPTIMIZER INTERFACE
import scipy.optimize
f = evaluate_energy
x0 = positions_short
fprime = evaluate_energy_gradient
args = (struct, kernelpot, opt_pids, False, ofs)
positions_opt = scipy.optimize.fmin_cg(f, x0, fprime=fprime, args=args, gtol=1e-6)
#positions_opt = scipy.optimize.minimize(f, x0, args=args, method='CG', jac=fprime, tol=1e-5)
#positions_opt = scipy.optimize.fmin_cg(f, x0, args=args, gtol=1e-5)
#scipy.optimize.minimize(f, x0, args=args, method='CG', jac=fprime, tol=1e-99)
#scipy.optimize.minimize(f, x0, args=args, method='CG', tol=1e-5)
ofs.close()
print "Optimized positions:", positions_opt
kernelpot.acquire(struct, -1.)
np.savetxt('opt_X.txt', kernelpot.IX[-1])
......@@ -65,6 +65,9 @@ class PCA(object):
X_recon = self.eigvecs[idx]
X_recon = self.unnorm(X_recon)
return X_recon
def reconstructEigenvecs(self):
X_eig_recon = self.unnormBlock(self.eigvecs)
return X_eig_recon
def reconstructBlock(self, eigencoeffs):
# Order such that: eigencoeffs[i] => eigencoefficients of sample i
n_samples = eigencoeffs.shape[0]
......
#! /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
# Load XX from text file
# Run PCA (normalize, diagonalize)
# Project structures onto PCA => plot
# TODO Reconstruct eigenstructures (unnormalize, load into kernel, vary centers, optimize)
logging.basicConfig(
format='[%(asctime)s] %(message)s',
datefmt='%I:%M:%S',
level=logging.ERROR)
# >>>>>>>>>>>>>>>>
norm_mean = True
norm_std = False
adaptor_type = 'global-generic'
# <<<<<<<<<<<<<<<<
# OPTIONS
options = soap.Options()
options.excludeCenters(['H'])
options.excludeTargets(['H'])
options.excludeCenterIds([])
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', 'dot')
options.set('kernel.delta', 1.)
options.set('kernel.xi', 4.) # <- pull here (1., ..., 4., ...)
# 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('configs')
]
# SETUP KERNEL (=> BASIS, KERNELPOT WITH ADAPTOR)
basis = soap.Basis(options)
kernelpot = KernelPotential(options)
# FILL KERNEL
generate = True
if generate:
for struct in structures:
print struct.label
kernelpot.acquire(struct, 1., label=struct.label)
print kernelpot.IX.shape
np.savetxt('out.kernelpot.ix.txt', kernelpot.IX)
ofs = open('out.kernelpot.labels.txt', 'w')
for idx,label in enumerate(kernelpot.labels):
ofs.write('%3d %s\n' % (idx,label))
ofs.close()
else:
IX = np.loadtxt('out.kernelpot.ix.txt')
kernelpot.importAcquire(IX, 1.)
# KERNEL PCA
pca = PCA()
pca.compute(kernelpot.IX, normalize_mean=norm_mean, normalize_std=norm_std)
#pca = IPCA()
#pca.compute(IX, normalize_mean=False, normalize_std=False)
# =============================
# CHECK COMPONENT NORMALIZATION
# =============================
ones_vec = np.zeros(567)
ones_vec.fill(1.)
np.savetxt('out.pca.unnorm.txt', pca.unnormBlock(ones_vec))
# =================
# INDEX CUTOFF SCAN
# =================
K_full = None
idx_ref = 567
sel = 7
for idx_cutoff in [idx_ref, 500, 400, 300, 200, 100, 50, 40, 30, 20, 10, 5, 4, 3, 2, 1]:
print idx_cutoff
I_eigencoeffs = pca.expandBlock(kernelpot.IX, idx_cutoff)
IX_recon = pca.reconstructBlock(I_eigencoeffs)
K = IX_recon.dot(IX_recon.T)
eigen_K = I_eigencoeffs.dot(I_eigencoeffs.T)
if idx_cutoff == idx_ref:
K_full = K
K_slice = K[sel,:]
K_slice_full = K_full[sel,:]
def stretch_K(K):
k_min = np.min(K)
k_max = np.max(K)
return (K-k_min)/(k_max-k_min)
K_slice = stretch_K(K_slice)
K_slice_full = stretch_K(K_slice_full)
KK = np.array([K_slice,K_slice_full]).T
np.savetxt('out.pca.kernel_recon_cutoff_%03d.txt' % (idx_cutoff), KK)
# ========================
# EXPAND KERNEL STRUCTURES
# ========================
I_eigencoeffs = pca.expandBlock(kernelpot.IX)
np.savetxt('out.pca.eigencoeffs.txt', I_eigencoeffs)
# =================
# SAVE EIGENVECTORS
# =================
X_eig_0 = pca.reconstructEigenvec(0)
X_eig_1 = pca.reconstructEigenvec(1)
print np.dot(X_eig_0, X_eig_0)
print np.dot(X_eig_0, X_eig_1)
print np.dot(X_eig_1, X_eig_1)
X_eig_vecs = pca.reconstructEigenvecs()
np.savetxt('out.pca.xeigenvecs.txt', X_eig_vecs)
#! /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
# Load XX from text file
# Run PCA (normalize, diagonalize)
# Project structures onto PCA => plot
# TODO Reconstruct eigenstructures (unnormalize, load into kernel, vary centers, optimize)
logging.basicConfig(
format='[%(asctime)s] %(message)s',
datefmt='%I:%M:%S',
level=logging.ERROR)
# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
idx_eig_target = int(sys.argv[1])
N = int(sys.argv[2])
select_struct_label = 'config_C_%02d.xyz' % N
exclude_pids_rnd = [1]
exclude_centers = [ i for i in range(2, N+1) ]
opt_pids = [ i for i in range(2, N+1) ]
x_array_file = 'out.pca.xeigenvecs.txt'
#x_array_file = 'out.kernelpot.ix.txt'
adaptor_type = 'global-generic'
sc = 0.1 # 1.5 # 0.1 # 1.5
# <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
if adaptor_type == 'global-generic':
exclude_pids_rnd = []
exclude_centers = []
opt_pids = [ i for i in range(1,N+1) ]
# 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', 'dot')
options.set('kernel.delta', 1.)
options.set('kernel.xi', 2.) # <- pull here (1., ..., 4., ...)
# 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('configs_n_carbons')
]
struct_dict = {}
for struct in structures:
print struct.label
struct_dict[struct.label] = struct
# LOAD EIGENSTRUCTURES
print "Load fingerprints from '%s'" % x_array_file
IX_eig = np.loadtxt(x_array_file)
x1 = IX_eig[0]
x2 = IX_eig[4]
print np.dot(x1,x1), np.dot(x1,x2), np.dot(x2,x2)
# SELECT STRUCTURE, X_TARGET
print "Target eigenstructure with index", idx_eig_target
print "Exclude centers", exclude_centers
print "Select label", select_struct_label
X_target = IX_eig[idx_eig_target]
struct = struct_dict[select_struct_label]
options.excludeCenterIds(exclude_centers)
# SETUP KERNEL (=> BASIS, KERNELPOT WITH ADAPTOR)
basis = soap.Basis(options)
kernelpot = KernelPotential(options)
kernelpot.importAcquire(np.array([X_target]), -1.) # <- attractive
# RANDOMIZE INITIAL CONFIG
positions_0 = [ part.pos for part in struct ]
positions_0_rnd = random_positions(struct, exclude_pid=exclude_pids_rnd, b0=-1.*sc, b1=1.*sc)
print "Positions after initial randomization"
for part in struct:
print part.id, part.pos
# EXTRACT POSITIONS (OPTIMIZATION TARGETS)
opt_pids = np.array(opt_pids)
opt_pidcs = opt_pids-1
positions = np.array(positions_0_rnd)
positions_short = positions[opt_pidcs]
positions_short = positions_short.flatten()
print "Optimization arguments (initial)"
for pos in positions_short.reshape((len(opt_pids),3)):
print pos
#raw_input('...')