Skip to content
GitLab
Commit 5e1b95c2 authored by Carl Poelking's avatar Carl Poelking
Browse files

Configuration kernel PCA.

parent 3d8e7a67
Hide whitespace changes
Inline Side-by-side
CMakeLists.txt
View file @ 5e1b95c2
......@@ -10,7 +10,7 @@ set(CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/cmake)
# BUILD OPTIONS
enable_language(CXX)
#set(CMAKE_CXX_COMPILER "/usr/local/shared/intel/compilers_and_libraries_2016.0.109/linux/bin/intel64/icc")
set(CMAKE_CXX_COMPILER "/usr/local/shared/intel/compilers_and_libraries_2016.2.181/linux/bin/intel64/icc")
message("C++ compiler: " ${CMAKE_CXX_COMPILER} " " ${CMAKE_CXX_COMPILER_ID})
option(BUILD_SHARED_LIBS "Build shared libs" ON)
if(${CMAKE_VERSION} VERSION_GREATER 3.1)
......
build.sh
View file @ 5e1b95c2
#! /bin/bash
mkdir -p build
cd build
cmake .. && make && make install
cmake .. -DCMAKE_CXX_COMPILER_ID=Intel && make && make install
cd ..
test/test_kernelpot/kernel.py
View file @ 5e1b95c2
......@@ -166,6 +166,28 @@ class KernelPotential(object):
# INCLUSIONS / EXCLUSIONS
# -> Already be enforced when computing spectra
return
def importAcquire(self, IX_acqu, alpha):
n_acqu = IX_acqu.shape[0]
dim_acqu = IX_acqu.shape[1]
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("Imported %d environments." % n_acqu)
return
def acquire(self, structure, alpha):
logging.info("Acquire ...")
spectrum = soap.Spectrum(structure, self.options, self.basis)
......
test/test_pca/config_pca.py 0 → 100755
View file @ 5e1b95c2
#! /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.DEBUG)
# 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', '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')
]
# SETUP KERNEL (=> BASIS, KERNELPOT WITH ADAPTOR)
basis = soap.Basis(options)
kernelpot = KernelPotential(options)
# FILL KERNEL
generate = False
if generate:
for struct in structures:
print struct.label
kernelpot.acquire(struct, 1.)
print kernelpot.IX.shape
np.savetxt('out.kernelpot.ix.txt', kernelpot.IX)
else:
IX = np.loadtxt('out.kernelpot.ix.txt')
kernelpot.importAcquire(IX, 1.)
# KERNEL PCA
pca = PCA()
pca.compute(IX, normalize_mean=False, normalize_std=False)
#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)
np.savetxt('out.pca.eigenvec_0.txt', X_eig_0)
np.savetxt('out.pca.eigenvec_1.txt', X_eig_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)
test/test_pca/pca.py
View file @ 5e1b95c2
#! /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
# 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', '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')
]
# SETUP KERNEL (=> BASIS, KERNELPOT WITH ADAPTOR)
basis = soap.Basis(options)
kernelpot = KernelPotential(options)
for struct in structures:
print struct.label
kernelpot.acquire(struct, 1.)
print kernelpot.IX.shape
np.savetxt('kernelpot.ix.txt', kernelpot.IX)
# Load XX from text file
# Run PCA (normalize, diagonalize)
# Project structures onto PCA => plot
# Reconstruct eigenstructures (unnormalize, load into kernel, vary centers, optimize)
def PCAtrain(Xtrain):
#### IMPLEMENT HERE
#### Note that for the eigen decomposition you can use the numpy function 'np.linalg.eigh'
mean = X_train.mean(0)
std = X_train.std(0)
print(mean.shape)
print(X_train.shape)
Xtrain = (Xtrain - mean) / std
Sigma = np.dot(Xtrain.T,Xtrain)/Xtrain.shape[0]
singular_values, pc_vectors = np.linalg.eigh(Sigma)
idx = singular_values.argsort()[::-1]
singular_values = singular_values[idx]
pc_vectors = pc_vectors[:,idx]
return pc_vectors, singular_values, mean, std
class PCA(object):
def __init__(self):
self.IX = None
self.IX_norm = None
self.X_mean = None
self.X_std = None
return
def compute(self, IX, normalize_mean=True, normalize_std=True):
# Normalize
self.IX = IX
self.X_mean = IX.mean(0)
self.X_std= IX.std(0)
if not normalize_mean:
self.X_mean.fill(0.)
if not normalize_std:
self.X_std.fill(1.)
self.IX_norm = (IX - self.X_mean)/self.X_std
# Correlation matrix / moment matrix
Sigma = np.dot(self.IX_norm.T,self.IX_norm)/self.IX_norm.shape[0]
eigvals, eigvecs = np.linalg.eigh(Sigma)
# Sort largest -> smallest eigenvalue
idcs = eigvals.argsort()[::-1]
eigvals = eigvals[idcs]
eigvecs = eigvecs[:,idcs]
eigvecs.transpose()
# Store eigenvalues, eigenvectors
self.eigvals = eigvals # eigvals[i] => i-th eigenvalue, where
self.eigvecs = eigvecs # eigvecs[i] => i-th eigenvector
eigvals_cumulative = np.cumsum(eigvals)
np.savetxt('out.pca.sigma.txt', Sigma)
np.savetxt('out.pca.eigvals_cum.txt', eigvals_cumulative/eigvals_cumulative[-1])
return
def expand(self, X, idx_cutoff=None):
assert False
eigvecs_sel = self.eigvecs[0:idx_cutoff] if idx_cutoff else self.eigvecs
eigvals_sel = self.eigvals[0:idx_cutoff] if idx_cutoff else self.eigvals
# Normalize
X_norm = X-self.X_mean
X_norm = X_norm/self.X_std
# Expand
eigencoeffs = eigvecs_sel.dot(X_norm)
# Reconstruct
X_recon = eigencoeffs.dot(eigvecs_sel)
X_recon = self.unnorm(X_recon)
return X_recon
def unnorm(self, X):
return X*self.X_std+self.X_mean
def expandBlock(self, IX, idx_cutoff=None):
eigvecs_sel = self.eigvecs[0:idx_cutoff] if idx_cutoff else self.eigvecs
eigvals_sel = self.eigvals[0:idx_cutoff] if idx_cutoff else self.eigvals
# Normalize
IX_norm = IX-self.X_mean
IX_norm = IX_norm/self.X_std
# Expand
eigencoeffs = IX_norm.dot(eigvecs_sel.T)
return eigencoeffs
def reconstructEigenvec(self, idx):
X_recon = self.eigvecs[idx]
X_recon = self.unnorm(X_recon)
return X_recon
def reconstructBlock(self, eigencoeffs):
# Order such that: eigencoeffs[i] => eigencoefficients of sample i
n_samples = eigencoeffs.shape[0]
n_components = eigencoeffs.shape[1]
idx_cutoff = n_components
eigvecs_sel = self.eigvecs[0:idx_cutoff]
eigvals_sel = self.eigvals[0:idx_cutoff]
# Reconstruct
X_recon = eigencoeffs.dot(eigvecs_sel)
X_recon = self.unnormBlock(X_recon)
return X_recon
def unnormBlock(self, IX):
return IX*self.X_std+self.X_mean
def dot(self, e1, e2):
# e1, e1 => vectors with eigencoefficients of two samples
assert e1.shape == e2.shape
n_components = e1.shape[0]
assert False and not "SENSE"
class IPCA(PCA):
def compute(self, IX, normalize_mean=True, normalize_std=True):
print "Compute invariant"
# Normalize
self.IX = IX
self.X_mean = IX.mean(0)
self.X_std= IX.std(0)
if not normalize_mean:
self.X_mean.fill(0.)
if not normalize_std:
self.X_std.fill(1.)
self.IX_norm = (IX - self.X_mean)/self.X_std
# Correlation matrix / moment matrix
Sigma = np.dot(self.IX_norm.T,self.IX_norm)/self.IX_norm.shape[0]
# Invariance step
d = (1./np.diag(Sigma))**0.5
D = np.zeros(Sigma.shape)
np.fill_diagonal(D, d)
Sigma = D.dot(Sigma).dot(D)
eigvals, eigvecs = np.linalg.eigh(Sigma)
# Sort largest -> smallest eigenvalue
idcs = eigvals.argsort()[::-1]
eigvals = eigvals[idcs]
eigvecs = eigvecs[:,idcs]
eigvecs.transpose()
# Store eigenvalues, eigenvectors
self.eigvals = eigvals # eigvals[i] => i-th eigenvalue, where
self.eigvecs = eigvecs # eigvecs[i] => i-th eigenvector
eigvals_cumulative = np.cumsum(eigvals)
np.savetxt('out.pca.sigma.txt', Sigma)
np.savetxt('out.pca.eigvals_cum.txt', eigvals_cumulative/eigvals_cumulative[-1])
return
......
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