From 76441bafdb1cc3dda580c411db7af4626051b57b Mon Sep 17 00:00:00 2001 From: Jonathan Vandermause Date: Sat, 22 Feb 2020 12:19:53 -0500 Subject: [PATCH 01/24] update y vector when adding training structures --- src/kernels.cpp | 6 +++--- src/sparse_gp.cpp | 19 ++++++++++++++++++- src/sparse_gp.h | 2 +- src/structure.cpp | 8 ++++---- src/structure.h | 8 ++++---- tests/test_kernels.cpp | 6 +++--- tests/test_sparse_gp.cpp | 40 ++++++++++++++++++++++------------------ tests/test_structure.cpp | 4 ++-- 8 files changed, 57 insertions(+), 36 deletions(-) diff --git a/src/kernels.cpp b/src/kernels.cpp index 50751fe51..2ab379eb6 100644 --- a/src/kernels.cpp +++ b/src/kernels.cpp @@ -86,7 +86,7 @@ Eigen::VectorXd TwoBodyKernel :: env_struc(const LocalEnvironment & env1, LocalEnvironment env_curr; for (int i = 0; i < struc1.noa; i ++){ - env_curr = struc1.environment_descriptors[i]; + env_curr = struc1.local_environments[i]; cent2 = env_curr.central_species; for (int m = 0; m < env1.rs.size(); m ++){ @@ -304,7 +304,7 @@ Eigen::VectorXd ThreeBodyKernel :: env_struc(const LocalEnvironment & env1, int c1 = env1.central_species; for (int i = 0; i < struc1.noa; i ++){ - env2 = struc1.environment_descriptors[i]; + env2 = struc1.local_environments[i]; c2 = env2.central_species; for (int m = 0; m < env1.three_body_indices.size(); m ++){ @@ -522,7 +522,7 @@ Eigen::VectorXd DotProductKernel LocalEnvironment env_curr; for (int i = 0; i < struc1.noa; i ++){ - env_curr = struc1.environment_descriptors[i]; + env_curr = struc1.local_environments[i]; // Check that the environments have the same central species. if (env1.central_species != env_curr.central_species) continue; diff --git a/src/sparse_gp.cpp b/src/sparse_gp.cpp index b6c344099..2a39348e5 100644 --- a/src/sparse_gp.cpp +++ b/src/sparse_gp.cpp @@ -124,5 +124,22 @@ void SparseGP :: add_training_structure(StructureDescriptor training_structure){ // Store training structure. training_structures.push_back(training_structure); - // TODO: update y vector. + // Update y vector. + Eigen::VectorXd labels = Eigen::VectorXd::Zero(n_labels); + + if (training_structure.energy.size() != 0){ + labels.head(1) = training_structure.energy; + } + + if (training_structure.forces.size() != 0){ + labels.segment(1, n_atoms * 3) = training_structure.forces; + } + + if (training_structure.stresses.size() != 0){ + stresses.conservativeResize(stresses.size() + 6); + labels.tail(6) = training_structure.stresses; + } + + y.conservativeResize(y.size() + n_labels); + y.tail(n_labels) = labels; } diff --git a/src/sparse_gp.h b/src/sparse_gp.h index 333da6258..d860db89e 100644 --- a/src/sparse_gp.h +++ b/src/sparse_gp.h @@ -12,11 +12,11 @@ class SparseGP{ public: Eigen::MatrixXd Kuu, Kuf, Sigma; Eigen::VectorXd y, alpha, kernel_hyperparameters; + Eigen::VectorXd energies, forces, stresses; std::vector kernels; std::vector sparse_environments; std::vector training_structures; - std::vector energies, forces, stresses; double energy_norm, forces_norm, stresses_norm, energy_offset, forces_offset, stresses_offset; diff --git a/src/structure.cpp b/src/structure.cpp index 3ebd19ed5..b1074f4ba 100644 --- a/src/structure.cpp +++ b/src/structure.cpp @@ -124,7 +124,7 @@ void StructureDescriptor :: compute_environments(){ for (int i = 0; i < noa; i ++){ env = LocalEnvironment(*this, i, cutoff); - environment_descriptors.push_back(env); + local_environments.push_back(env); } } @@ -134,7 +134,7 @@ void StructureDescriptor :: compute_nested_environments(){ for (int i = 0; i < noa; i ++){ env = LocalEnvironment(*this, i, cutoff, nested_cutoffs); - environment_descriptors.push_back(env); + local_environments.push_back(env); } } @@ -144,7 +144,7 @@ void StructureDescriptor :: compute_descriptors(){ for (int i = 0; i < noa; i ++){ env = LocalEnvironment(*this, i, cutoff, descriptor_calculator); - environment_descriptors.push_back(env); + local_environments.push_back(env); } } @@ -155,6 +155,6 @@ void StructureDescriptor :: nested_descriptors(){ for (int i = 0; i < noa; i ++){ env = LocalEnvironment(*this, i, cutoff, nested_cutoffs, descriptor_calculator); - environment_descriptors.push_back(env); + local_environments.push_back(env); } } diff --git a/src/structure.h b/src/structure.h index e6969ee5f..0ea448ab9 100644 --- a/src/structure.h +++ b/src/structure.h @@ -31,14 +31,14 @@ class Structure{ class StructureDescriptor : public Structure{ public: DescriptorCalculator * descriptor_calculator; - std::vector environment_descriptors; + std::vector local_environments; double cutoff; std::vector nested_cutoffs; // Make structure labels empty by default. - std::vector energy {}; - std::vector forces {}; - std::vector stresses {}; + Eigen::VectorXd energy; + Eigen::VectorXd forces; + Eigen::VectorXd stresses; StructureDescriptor(); diff --git a/tests/test_kernels.cpp b/tests/test_kernels.cpp index 44ac222f4..2e59c87d4 100644 --- a/tests/test_kernels.cpp +++ b/tests/test_kernels.cpp @@ -65,7 +65,7 @@ class KernelTest : public ::testing::Test{ test_struc_2 = StructureDescriptor(cell, species, positions_2, cutoff, nested_cutoffs, &desc1); - test_env = test_struc_2.environment_descriptors[0]; + test_env = test_struc_2.local_environments[0]; kernel = DotProductKernel(kernel_hyperparameters); two_body_kernel = TwoBodyKernel(length_scale, cutoff_string, @@ -80,8 +80,8 @@ class KernelTest : public ::testing::Test{ }; TEST_F(KernelTest, NormTest){ - LocalEnvironment env1 = test_struc.environment_descriptors[0]; - LocalEnvironment env2 = test_struc.environment_descriptors[1]; + LocalEnvironment env1 = test_struc.local_environments[0]; + LocalEnvironment env2 = test_struc.local_environments[1]; double kern_val = kernel.env_env(env1, env1); EXPECT_NEAR(kern_val, 1, THRESHOLD); } diff --git a/tests/test_sparse_gp.cpp b/tests/test_sparse_gp.cpp index 5327f6ea5..c9096aeb7 100644 --- a/tests/test_sparse_gp.cpp +++ b/tests/test_sparse_gp.cpp @@ -5,10 +5,15 @@ class SparseTest : public ::testing::Test{ public: // structure + int n_atoms = 2; Eigen::MatrixXd cell{3, 3}, cell_2{3, 3}; - std::vector species {0, 1, 0, 1, 0}; - Eigen::MatrixXd positions{5, 3}, positions_2{5, 3}, positions_3{5, 3}; + std::vector species {0, 1}; + Eigen::MatrixXd positions{n_atoms, 3}, positions_2{n_atoms, 3}, + positions_3{n_atoms, 3}; StructureDescriptor test_struc, test_struc_2, test_struc_3; + + // labels + Eigen::VectorXd energy{1}, forces{n_atoms * 3}, stresses{6}; // descriptor std::string radial_string = "chebyshev"; @@ -37,22 +42,23 @@ class SparseTest : public ::testing::Test{ 0, 0, 10; positions << -0.68402216, 0.54343671, -0.52961224, - 0.33045915, 0.40010388, 0.59849816, - -0.92832825, 0.06239221, 0.51601996, - 0.75120489, -0.39606988, -0.34000017, - -0.8242705, -0.73860995, 0.92679555; + 0.33045915, 0.40010388, 0.59849816; positions_2 << 0.69955637, -0.41619112, -0.51725003, - 0.43189622, 0.88548458, -0.74495343, - 0.31395126, -0.32179606, -0.35013419, - 0.08793497, -0.70567732, -0.3811633, - 0.35585787, -0.87190223, 0.06770428; + 0.43189622, 0.88548458, -0.74495343; + + energy = Eigen::VectorXd::Random(1); + forces = Eigen::VectorXd::Random(n_atoms * 3); + stresses = Eigen::VectorXd::Random(6); desc1 = B2_Calculator(radial_string, cutoff_string, radial_hyps, cutoff_hyps, descriptor_settings); test_struc = StructureDescriptor(cell, species, positions, cutoff, nested_cutoffs, &desc1); + test_struc.energy = energy; + test_struc.forces = forces; + test_struc.stresses = stresses; two_body_kernel = TwoBodyKernel(length_scale, cutoff_string, cutoff_hyps); @@ -68,18 +74,16 @@ class SparseTest : public ::testing::Test{ TEST_F(SparseTest, AddSparseTest){ SparseGP sparse_gp = SparseGP(kernels); - LocalEnvironment env1 = test_struc.environment_descriptors[0]; - LocalEnvironment env2 = test_struc.environment_descriptors[1]; - LocalEnvironment env3 = test_struc.environment_descriptors[2]; + LocalEnvironment env1 = test_struc.local_environments[0]; + LocalEnvironment env2 = test_struc.local_environments[1]; sparse_gp.add_sparse_environment(env1); sparse_gp.add_sparse_environment(env2); - sparse_gp.add_sparse_environment(env3); - - test_struc.energy = std::vector {10}; - test_struc.forces = std::vector (15, 10); - test_struc.stresses = std::vector {1, 2, 3, 4, 5, 6}; sparse_gp.add_training_structure(test_struc); + sparse_gp.add_training_structure(test_struc); + + std::cout << sparse_gp.y; + // sparse_gp.add_training_structure(test_struc); // std::cout << sparse_gp.Kuf << std::endl; // std::cout << sparse_gp.training_structures.size() << std::endl; diff --git a/tests/test_structure.cpp b/tests/test_structure.cpp index af046dfa1..6f6ed6c3a 100644 --- a/tests/test_structure.cpp +++ b/tests/test_structure.cpp @@ -70,11 +70,11 @@ TEST_F(StructureTest, StructureDescriptor){ for (int j = 0; j < desc1.descriptor_vals.size(); j ++){ EXPECT_EQ(desc1.descriptor_vals(j), - test_struc.environment_descriptors[i] + test_struc.local_environments[i] .descriptor_vals(j)); for (int k = 0; k < test_struc.species.size(); k ++){ EXPECT_EQ(desc1.descriptor_force_dervs(k, j), - test_struc.environment_descriptors[i] + test_struc.local_environments[i] .descriptor_force_dervs(k, j)); } } From 7289b9bb54318d4574065613bdffaac91f214507 Mon Sep 17 00:00:00 2001 From: Jonathan Vandermause Date: Sat, 22 Feb 2020 12:34:19 -0500 Subject: [PATCH 02/24] add test of covariance updates --- tests/test_sparse_gp.cpp | 20 +++++++++++++++----- 1 file changed, 15 insertions(+), 5 deletions(-) diff --git a/tests/test_sparse_gp.cpp b/tests/test_sparse_gp.cpp index c9096aeb7..bfb9bd85c 100644 --- a/tests/test_sparse_gp.cpp +++ b/tests/test_sparse_gp.cpp @@ -72,7 +72,7 @@ class SparseTest : public ::testing::Test{ } }; -TEST_F(SparseTest, AddSparseTest){ +TEST_F(SparseTest, UpdateK){ SparseGP sparse_gp = SparseGP(kernels); LocalEnvironment env1 = test_struc.local_environments[0]; LocalEnvironment env2 = test_struc.local_environments[1]; @@ -80,11 +80,21 @@ TEST_F(SparseTest, AddSparseTest){ sparse_gp.add_sparse_environment(env2); sparse_gp.add_training_structure(test_struc); + + test_struc.stresses = Eigen::VectorXd {}; sparse_gp.add_training_structure(test_struc); - std::cout << sparse_gp.y; + EXPECT_EQ(sparse_gp.Kuu.rows(), 2); + EXPECT_EQ(sparse_gp.Kuf.rows(), 2); + EXPECT_EQ(sparse_gp.Kuf.cols(), 2 + 6 * n_atoms + 6); + EXPECT_EQ(sparse_gp.y.size(), 2 + 6 * n_atoms + 6); + EXPECT_EQ(test_struc.forces, sparse_gp.y.segment(1, 3 * n_atoms)); - // sparse_gp.add_training_structure(test_struc); - // std::cout << sparse_gp.Kuf << std::endl; - // std::cout << sparse_gp.training_structures.size() << std::endl; + double kern_val = 0; + double kern_curr; + for (int i = 0; i < kernels.size(); i ++){ + kern_curr = kernels[i] -> env_env(env1, env1); + kern_val += kern_curr; + } + EXPECT_EQ(kern_val, sparse_gp.Kuu(0,0)); } From 7a771fed8bd408e14b65da871ba4138b73e24245 Mon Sep 17 00:00:00 2001 From: Jonathan Vandermause Date: Sat, 22 Feb 2020 14:45:37 -0500 Subject: [PATCH 03/24] set kernel hyperparameters attribute --- src/kernels.cpp | 12 +++++++----- src/kernels.h | 5 ++--- tests/test_kernels.cpp | 3 +-- tests/test_sparse_gp.cpp | 3 +-- 4 files changed, 11 insertions(+), 12 deletions(-) diff --git a/src/kernels.cpp b/src/kernels.cpp index 2ab379eb6..7cae2fb6f 100644 --- a/src/kernels.cpp +++ b/src/kernels.cpp @@ -19,6 +19,8 @@ TwoBodyKernel :: TwoBodyKernel(double ls, const std::string & cutoff_function, ls2 = 1 / (ls * ls); this->cutoff_hyps = cutoff_hyps; + kernel_hyperparameters.push_back(ls); + if (cutoff_function == "quadratic"){ this->cutoff_pointer = quadratic_cutoff; } @@ -159,6 +161,8 @@ ThreeBodyKernel :: ThreeBodyKernel(double ls, ls2 = 1 / (ls * ls); this->cutoff_hyps = cutoff_hyps; + kernel_hyperparameters.push_back(ls); + if (cutoff_function == "quadratic"){ this->cutoff_pointer = quadratic_cutoff; } @@ -481,12 +485,10 @@ void ThreeBodyKernel :: update_kernel_vector(Eigen::VectorXd & kernel_vector, DotProductKernel :: DotProductKernel() {}; -DotProductKernel :: DotProductKernel(std::vector - kernel_hyperparameters) : Kernel(kernel_hyperparameters){ +DotProductKernel :: DotProductKernel(double power){ - signal_variance = kernel_hyperparameters[0]; - sig2 = signal_variance * signal_variance; - power = kernel_hyperparameters[1]; + this->power = power; + kernel_hyperparameters.push_back(power); } double DotProductKernel :: env_env(const LocalEnvironment & env1, diff --git a/src/kernels.h b/src/kernels.h index 7126f932a..7e63f400b 100644 --- a/src/kernels.h +++ b/src/kernels.h @@ -24,12 +24,11 @@ class Kernel{ class DotProductKernel : public Kernel{ public: - double signal_variance, power, sig2; + double power; DotProductKernel(); - // Hyperparameters: [signal_variance, power] - DotProductKernel(std::vector kernel_hyperparameters); + DotProductKernel(double power); double env_env(const LocalEnvironment & env1, const LocalEnvironment & env2); diff --git a/tests/test_kernels.cpp b/tests/test_kernels.cpp index 2e59c87d4..7792a5354 100644 --- a/tests/test_kernels.cpp +++ b/tests/test_kernels.cpp @@ -35,7 +35,6 @@ class KernelTest : public ::testing::Test{ double signal_variance = 2; double length_scale = 1; double power = 2; - std::vector kernel_hyperparameters {signal_variance, power}; DotProductKernel kernel; TwoBodyKernel two_body_kernel; ThreeBodyKernel three_body_kernel; @@ -67,7 +66,7 @@ class KernelTest : public ::testing::Test{ nested_cutoffs, &desc1); test_env = test_struc_2.local_environments[0]; - kernel = DotProductKernel(kernel_hyperparameters); + kernel = DotProductKernel(power); two_body_kernel = TwoBodyKernel(length_scale, cutoff_string, cutoff_hyps); three_body_kernel = ThreeBodyKernel(length_scale, cutoff_string, diff --git a/tests/test_sparse_gp.cpp b/tests/test_sparse_gp.cpp index bfb9bd85c..23e06e89a 100644 --- a/tests/test_sparse_gp.cpp +++ b/tests/test_sparse_gp.cpp @@ -29,7 +29,6 @@ class SparseTest : public ::testing::Test{ double signal_variance = 2; double length_scale = 1; double power = 2; - std::vector kernel_hyperparameters {signal_variance, power}; DotProductKernel kernel; TwoBodyKernel two_body_kernel; ThreeBodyKernel three_body_kernel; @@ -64,7 +63,7 @@ class SparseTest : public ::testing::Test{ cutoff_hyps); three_body_kernel = ThreeBodyKernel(length_scale, cutoff_string, cutoff_hyps); - many_body_kernel = DotProductKernel(kernel_hyperparameters); + many_body_kernel = DotProductKernel(power); kernels = std::vector {&two_body_kernel, &three_body_kernel, From fe6e990d61b15636af5fc8a59a4c35e85cb94eb7 Mon Sep 17 00:00:00 2001 From: Jonathan Vandermause Date: Sat, 22 Feb 2020 15:17:26 -0500 Subject: [PATCH 04/24] add sigma to kernels --- src/kernels.cpp | 82 ++++++++++++++++++++++++---------------- src/kernels.h | 16 ++++---- tests/test_kernels.cpp | 12 +++--- tests/test_sparse_gp.cpp | 10 ++--- 4 files changed, 69 insertions(+), 51 deletions(-) diff --git a/src/kernels.cpp b/src/kernels.cpp index 7cae2fb6f..b3364df25 100644 --- a/src/kernels.cpp +++ b/src/kernels.cpp @@ -11,15 +11,19 @@ Kernel :: Kernel(std::vector kernel_hyperparameters){ TwoBodyKernel :: TwoBodyKernel() {}; -TwoBodyKernel :: TwoBodyKernel(double ls, const std::string & cutoff_function, +TwoBodyKernel :: TwoBodyKernel(double sigma, double ls, + const std::string & cutoff_function, std::vector cutoff_hyps){ + this->sigma = sigma; + sig2 = sigma * sigma; + this->ls = ls; ls1 = 1 / (2 * ls * ls); ls2 = 1 / (ls * ls); this->cutoff_hyps = cutoff_hyps; - kernel_hyperparameters.push_back(ls); + kernel_hyperparameters = std::vector {sigma, ls}; if (cutoff_function == "quadratic"){ this->cutoff_pointer = quadratic_cutoff; @@ -64,7 +68,7 @@ double TwoBodyKernel :: env_env(const LocalEnvironment & env1, } } } - return kern; + return sig2 * kern; } Eigen::VectorXd TwoBodyKernel :: env_struc(const LocalEnvironment & env1, @@ -120,7 +124,7 @@ Eigen::VectorXd TwoBodyKernel :: env_struc(const LocalEnvironment & env1, // energy kernel c1 = rdiff * rdiff; c2 = exp(-c1 * ls1); - kernel_vector(0) += fi * fj * c2 / 2; + kernel_vector(0) += sig2 * fi * fj * c2 / 2; // helper constants c3 = c2 * ls2 * fi * fj * rdiff; @@ -128,21 +132,27 @@ Eigen::VectorXd TwoBodyKernel :: env_struc(const LocalEnvironment & env1, // fx + exx, exy, exz stress components fx = xrel * c3 + c4 * xrel; - kernel_vector(1 + 3 * i) += fx; - kernel_vector(no_elements - 6) -= fx * xval * vol_inv / 2; - kernel_vector(no_elements - 5) -= fx * yval * vol_inv / 2; - kernel_vector(no_elements - 4) -= fx * zval * vol_inv / 2; - + kernel_vector(1 + 3 * i) += sig2 * fx; + kernel_vector(no_elements - 6) -= + sig2 * fx * xval * vol_inv / 2; + kernel_vector(no_elements - 5) -= + sig2 * fx * yval * vol_inv / 2; + kernel_vector(no_elements - 4) -= + sig2 * fx * zval * vol_inv / 2; + // fy + eyy, eyz stress components fy = yrel * c3 + c4 * yrel; - kernel_vector(2 + 3 * i) += fy; - kernel_vector(no_elements - 3) -= fy * yval * vol_inv / 2; - kernel_vector(no_elements - 2) -= fy * zval * vol_inv / 2; + kernel_vector(2 + 3 * i) += sig2 * fy; + kernel_vector(no_elements - 3) -= + sig2 * fy * yval * vol_inv / 2; + kernel_vector(no_elements - 2) -= + sig2 * fy * zval * vol_inv / 2; // fz + ezz stress component fz = zrel * c3 + c4 * zrel; - kernel_vector(3 + 3 * i) += fz; - kernel_vector(no_elements - 1) -= fz * zval * vol_inv / 2; + kernel_vector(3 + 3 * i) += sig2 * fz; + kernel_vector(no_elements - 1) -= + sig2 * fz * zval * vol_inv / 2; } } } @@ -153,15 +163,18 @@ Eigen::VectorXd TwoBodyKernel :: env_struc(const LocalEnvironment & env1, ThreeBodyKernel :: ThreeBodyKernel() {}; -ThreeBodyKernel :: ThreeBodyKernel(double ls, +ThreeBodyKernel :: ThreeBodyKernel(double sigma, double ls, const std::string & cutoff_function, std::vector cutoff_hyps){ + this->sigma = sigma; + sig2 = sigma * sigma; + this->ls = ls; ls1 = 1 / (2 * ls * ls); ls2 = 1 / (ls * ls); this->cutoff_hyps = cutoff_hyps; - kernel_hyperparameters.push_back(ls); + kernel_hyperparameters = std::vector {sigma, ls}; if (cutoff_function == "quadratic"){ this->cutoff_pointer = quadratic_cutoff; @@ -281,7 +294,7 @@ double ThreeBodyKernel :: env_env(const LocalEnvironment & env1, } } - return kern; + return sig2 * kern; } Eigen::VectorXd ThreeBodyKernel :: env_struc(const LocalEnvironment & env1, @@ -453,42 +466,45 @@ void ThreeBodyKernel :: update_kernel_vector(Eigen::VectorXd & kernel_vector, double p1 = r11 * r11 + r22 * r22 + r33 * r33; double p2 = exp(-p1 * ls1); - kernel_vector(0) += p2 * fi * fj / 3; + kernel_vector(0) += sig2 * p2 * fi * fj / 3; double p3 = p2 * ls2 * fi * fj; double p4 = p2 * fi; double fx1 = p3 * r11 * xrel1 + p4 * fdjx1; double fx2 = p3 * r22 * xrel2 + p4 * fdjx2; - kernel_vector(1 + 3 * i) += fx1 + fx2; + kernel_vector(1 + 3 * i) += sig2 * (fx1 + fx2); kernel_vector(no_elements - 6) -= - (fx1 * xval1 + fx2 * xval2) * vol_inv / 2; + sig2 * (fx1 * xval1 + fx2 * xval2) * vol_inv / 2; kernel_vector(no_elements - 5) -= - (fx1 * yval1 + fx2 * yval2) * vol_inv / 2; + sig2 * (fx1 * yval1 + fx2 * yval2) * vol_inv / 2; kernel_vector(no_elements - 4) -= - (fx1 * zval1 + fx2 * zval2) * vol_inv / 2; + sig2 * (fx1 * zval1 + fx2 * zval2) * vol_inv / 2; double fy1 = p3 * r11 * yrel1 + p4 * fdjy1; double fy2 = p3 * r22 * yrel2 + p4 * fdjy2; - kernel_vector(2 + 3 * i) += fy1 + fy2; + kernel_vector(2 + 3 * i) += sig2 * (fy1 + fy2); kernel_vector(no_elements - 3) -= - (fy1 * yval1 + fy2 * yval2) * vol_inv / 2; + sig2 * (fy1 * yval1 + fy2 * yval2) * vol_inv / 2; kernel_vector(no_elements - 2) -= - (fy1 * zval1 + fy2 * zval2) * vol_inv / 2; + sig2 * (fy1 * zval1 + fy2 * zval2) * vol_inv / 2; double fz1 = p3 * r11 * zrel1 + p4 * fdjz1; double fz2 = p3 * r22 * zrel2 + p4 * fdjz2; - kernel_vector(3 + 3 * i) += fz1 + fz2; + kernel_vector(3 + 3 * i) += sig2 * (fz1 + fz2); kernel_vector(no_elements - 1) -= - (fz1 * zval1 + fz2 * zval2) * vol_inv / 2; + sig2 * (fz1 * zval1 + fz2 * zval2) * vol_inv / 2; } DotProductKernel :: DotProductKernel() {}; -DotProductKernel :: DotProductKernel(double power){ +DotProductKernel :: DotProductKernel(double sigma, double power){ + this->sigma = sigma; + sig2 = sigma * sigma; this->power = power; - kernel_hyperparameters.push_back(power); + + kernel_hyperparameters = std::vector {sigma, power}; } double DotProductKernel :: env_env(const LocalEnvironment & env1, @@ -500,7 +516,7 @@ double DotProductKernel :: env_env(const LocalEnvironment & env1, double d1 = env1.descriptor_norm; double d2 = env2.descriptor_norm; - return pow(dot / (d1 * d2), power); + return sig2 * pow(dot / (d1 * d2), power); } Eigen::VectorXd DotProductKernel @@ -553,8 +569,8 @@ Eigen::VectorXd DotProductKernel stress_kern += dval * s1; } - kern_vec(0) = en_kern; - kern_vec.segment(1, struc1.noa * 3) = -force_kern; - kern_vec.tail(6) = -stress_kern * vol_inv; + kern_vec(0) = sig2 * en_kern; + kern_vec.segment(1, struc1.noa * 3) = -sig2 * force_kern; + kern_vec.tail(6) = -sig2 * stress_kern * vol_inv; return kern_vec; } diff --git a/src/kernels.h b/src/kernels.h index 7e63f400b..ab632a6f0 100644 --- a/src/kernels.h +++ b/src/kernels.h @@ -24,11 +24,11 @@ class Kernel{ class DotProductKernel : public Kernel{ public: - double power; + double sigma, sig2, power; DotProductKernel(); - DotProductKernel(double power); + DotProductKernel(double sigma, double power); double env_env(const LocalEnvironment & env1, const LocalEnvironment & env2); @@ -40,13 +40,14 @@ class DotProductKernel : public Kernel{ class TwoBodyKernel : public Kernel{ public: - double ls, ls1, ls2; + double sigma, sig2, ls, ls1, ls2; void (*cutoff_pointer)(double *, double, double, std::vector); std::vector cutoff_hyps; TwoBodyKernel(); - TwoBodyKernel(double ls, const std::string & cutoff_function, + TwoBodyKernel(double sigma, double ls, + const std::string & cutoff_function, std::vector cutoff_hyps); double env_env(const LocalEnvironment & env1, @@ -58,13 +59,14 @@ class TwoBodyKernel : public Kernel{ class ThreeBodyKernel : public Kernel{ public: - double ls, ls1, ls2; + double sigma, sig2, ls, ls1, ls2; void (*cutoff_pointer)(double *, double, double, std::vector); std::vector cutoff_hyps; ThreeBodyKernel(); - - ThreeBodyKernel(double ls, const std::string & cutoff_function, + + ThreeBodyKernel(double sigma, double ls, + const std::string & cutoff_function, std::vector cutoff_hyps); double env_env(const LocalEnvironment & env1, diff --git a/tests/test_kernels.cpp b/tests/test_kernels.cpp index 7792a5354..0c931fc0a 100644 --- a/tests/test_kernels.cpp +++ b/tests/test_kernels.cpp @@ -66,11 +66,11 @@ class KernelTest : public ::testing::Test{ nested_cutoffs, &desc1); test_env = test_struc_2.local_environments[0]; - kernel = DotProductKernel(power); - two_body_kernel = TwoBodyKernel(length_scale, cutoff_string, - cutoff_hyps); - three_body_kernel = ThreeBodyKernel(length_scale, cutoff_string, - cutoff_hyps); + kernel = DotProductKernel(signal_variance, power); + two_body_kernel = TwoBodyKernel(signal_variance, length_scale, + cutoff_string, cutoff_hyps); + three_body_kernel = ThreeBodyKernel(signal_variance, length_scale, + cutoff_string, cutoff_hyps); kern_vec = kernel.env_struc(test_env, test_struc); kern_vec_2 = two_body_kernel.env_struc(test_env, test_struc); @@ -82,7 +82,7 @@ TEST_F(KernelTest, NormTest){ LocalEnvironment env1 = test_struc.local_environments[0]; LocalEnvironment env2 = test_struc.local_environments[1]; double kern_val = kernel.env_env(env1, env1); - EXPECT_NEAR(kern_val, 1, THRESHOLD); + EXPECT_NEAR(kern_val, signal_variance*signal_variance, THRESHOLD); } TEST_F(KernelTest, ForceTest){ diff --git a/tests/test_sparse_gp.cpp b/tests/test_sparse_gp.cpp index 23e06e89a..c152a7d15 100644 --- a/tests/test_sparse_gp.cpp +++ b/tests/test_sparse_gp.cpp @@ -59,11 +59,11 @@ class SparseTest : public ::testing::Test{ test_struc.forces = forces; test_struc.stresses = stresses; - two_body_kernel = TwoBodyKernel(length_scale, cutoff_string, - cutoff_hyps); - three_body_kernel = ThreeBodyKernel(length_scale, cutoff_string, - cutoff_hyps); - many_body_kernel = DotProductKernel(power); + two_body_kernel = TwoBodyKernel(signal_variance, length_scale, + cutoff_string, cutoff_hyps); + three_body_kernel = ThreeBodyKernel(signal_variance, length_scale, + cutoff_string, cutoff_hyps); + many_body_kernel = DotProductKernel(signal_variance, power); kernels = std::vector {&two_body_kernel, &three_body_kernel, From 034e2038c4c11670f4330d6cf29673f605749af2 Mon Sep 17 00:00:00 2001 From: Jonathan Vandermause Date: Sat, 22 Feb 2020 15:37:35 -0500 Subject: [PATCH 05/24] set hyperparameters in sparse gp constructor --- src/sparse_gp.cpp | 32 +++++++++++++++++++++++++++++++- src/sparse_gp.h | 7 +++++-- tests/test_sparse_gp.cpp | 6 +++++- 3 files changed, 41 insertions(+), 4 deletions(-) diff --git a/src/sparse_gp.cpp b/src/sparse_gp.cpp index 2a39348e5..4bb2edeee 100644 --- a/src/sparse_gp.cpp +++ b/src/sparse_gp.cpp @@ -4,8 +4,38 @@ SparseGP :: SparseGP(){} -SparseGP :: SparseGP(std::vector kernels){ +SparseGP :: SparseGP(std::vector kernels, double sigma_e, + double sigma_f, double sigma_s){ + this->kernels = kernels; + + // Count hyperparameters. + int n_hyps = 0; + for (int i = 0; i < kernels.size(); i ++){ + n_hyps += kernels[i]->kernel_hyperparameters.size(); + } + + // Set the kernel hyperparameters. + hyperparameters = Eigen::VectorXd::Zero(n_hyps + 3); + std::vector hyps_curr; + int hyp_counter = 0; + for (int i = 0; i < kernels.size(); i ++){ + hyps_curr = kernels[i]->kernel_hyperparameters; + + for (int j = 0; j < hyps_curr.size(); j ++){ + hyperparameters(hyp_counter) = hyps_curr[j]; + hyp_counter ++; + } + } + + // Set the noise hyperparameters. + hyperparameters(n_hyps) = sigma_e; + hyperparameters(n_hyps+1) = sigma_f; + hyperparameters(n_hyps+2) = sigma_s; + + this->sigma_e = sigma_e; + this->sigma_f = sigma_f; + this->sigma_s = sigma_s; } void SparseGP :: add_sparse_environment(LocalEnvironment env){ diff --git a/src/sparse_gp.h b/src/sparse_gp.h index d860db89e..115f76714 100644 --- a/src/sparse_gp.h +++ b/src/sparse_gp.h @@ -11,7 +11,7 @@ class SparseGP{ public: Eigen::MatrixXd Kuu, Kuf, Sigma; - Eigen::VectorXd y, alpha, kernel_hyperparameters; + Eigen::VectorXd y, alpha, hyperparameters; Eigen::VectorXd energies, forces, stresses; std::vector kernels; @@ -21,9 +21,12 @@ class SparseGP{ double energy_norm, forces_norm, stresses_norm, energy_offset, forces_offset, stresses_offset; + double sigma_e, sigma_f, sigma_s; + SparseGP(); - SparseGP(std::vector); + SparseGP(std::vector, double sigma_e, double sigma_f, + double sigma_s); void add_sparse_environment(LocalEnvironment env); void add_training_structure(StructureDescriptor training_structure); diff --git a/tests/test_sparse_gp.cpp b/tests/test_sparse_gp.cpp index c152a7d15..4bc03425c 100644 --- a/tests/test_sparse_gp.cpp +++ b/tests/test_sparse_gp.cpp @@ -72,7 +72,11 @@ class SparseTest : public ::testing::Test{ }; TEST_F(SparseTest, UpdateK){ - SparseGP sparse_gp = SparseGP(kernels); + double sigma_e = 0; + double sigma_f = 0; + double sigma_s = 0; + + SparseGP sparse_gp = SparseGP(kernels, sigma_e, sigma_f, sigma_s); LocalEnvironment env1 = test_struc.local_environments[0]; LocalEnvironment env2 = test_struc.local_environments[1]; sparse_gp.add_sparse_environment(env1); From 7607214e95334d25711b61caddc7808c90950103 Mon Sep 17 00:00:00 2001 From: Jonathan Vandermause Date: Sat, 22 Feb 2020 16:54:51 -0500 Subject: [PATCH 06/24] update noise matrix when adding training structures --- src/sparse_gp.cpp | 13 +++++++++++-- src/sparse_gp.h | 5 ++--- tests/test_sparse_gp.cpp | 6 +++--- 3 files changed, 16 insertions(+), 8 deletions(-) diff --git a/src/sparse_gp.cpp b/src/sparse_gp.cpp index 4bb2edeee..fdd3489de 100644 --- a/src/sparse_gp.cpp +++ b/src/sparse_gp.cpp @@ -154,22 +154,31 @@ void SparseGP :: add_training_structure(StructureDescriptor training_structure){ // Store training structure. training_structures.push_back(training_structure); - // Update y vector. + // Update y vector and noise matrix. Eigen::VectorXd labels = Eigen::VectorXd::Zero(n_labels); + Eigen::VectorXd noise_vector = Eigen::VectorXd::Zero(n_labels); if (training_structure.energy.size() != 0){ labels.head(1) = training_structure.energy; + noise_vector(0) = 1 / (sigma_e * sigma_e); } if (training_structure.forces.size() != 0){ labels.segment(1, n_atoms * 3) = training_structure.forces; + noise_vector.segment(1, n_atoms * 3) = + Eigen::VectorXd::Constant(n_atoms * 3, 1 / (sigma_f * sigma_f)); } if (training_structure.stresses.size() != 0){ - stresses.conservativeResize(stresses.size() + 6); labels.tail(6) = training_structure.stresses; + noise_vector.tail(6) = + Eigen::VectorXd::Constant(6, 1 / (sigma_s * sigma_s)); } y.conservativeResize(y.size() + n_labels); y.tail(n_labels) = labels; + + noise.conservativeResize(prev_cols + n_labels); + noise.tail(n_labels) = noise_vector; + noise_matrix = noise.asDiagonal(); } diff --git a/src/sparse_gp.h b/src/sparse_gp.h index 115f76714..cb6bd692e 100644 --- a/src/sparse_gp.h +++ b/src/sparse_gp.h @@ -10,9 +10,8 @@ class SparseGP{ public: - Eigen::MatrixXd Kuu, Kuf, Sigma; - Eigen::VectorXd y, alpha, hyperparameters; - Eigen::VectorXd energies, forces, stresses; + Eigen::MatrixXd Kuu, Kuf, Sigma, noise_matrix; + Eigen::VectorXd y, alpha, hyperparameters, noise; std::vector kernels; std::vector sparse_environments; diff --git a/tests/test_sparse_gp.cpp b/tests/test_sparse_gp.cpp index 4bc03425c..7c1c300c0 100644 --- a/tests/test_sparse_gp.cpp +++ b/tests/test_sparse_gp.cpp @@ -72,9 +72,9 @@ class SparseTest : public ::testing::Test{ }; TEST_F(SparseTest, UpdateK){ - double sigma_e = 0; - double sigma_f = 0; - double sigma_s = 0; + double sigma_e = 1; + double sigma_f = 2; + double sigma_s = 3; SparseGP sparse_gp = SparseGP(kernels, sigma_e, sigma_f, sigma_s); LocalEnvironment env1 = test_struc.local_environments[0]; From 025b93a5bba18e975f46acb35137ce2c677f96bd Mon Sep 17 00:00:00 2001 From: Jonathan Vandermause Date: Sat, 22 Feb 2020 19:11:54 -0500 Subject: [PATCH 07/24] implement mean prediction --- src/sparse_gp.cpp | 34 ++++++++++++++++++++-- src/sparse_gp.h | 2 ++ tests/test_sparse_gp.cpp | 62 +++++++++++++++++++++++++++++----------- 3 files changed, 79 insertions(+), 19 deletions(-) diff --git a/src/sparse_gp.cpp b/src/sparse_gp.cpp index fdd3489de..566256fa3 100644 --- a/src/sparse_gp.cpp +++ b/src/sparse_gp.cpp @@ -158,14 +158,16 @@ void SparseGP :: add_training_structure(StructureDescriptor training_structure){ Eigen::VectorXd labels = Eigen::VectorXd::Zero(n_labels); Eigen::VectorXd noise_vector = Eigen::VectorXd::Zero(n_labels); + label_count = 0; if (training_structure.energy.size() != 0){ labels.head(1) = training_structure.energy; noise_vector(0) = 1 / (sigma_e * sigma_e); + label_count ++; } if (training_structure.forces.size() != 0){ - labels.segment(1, n_atoms * 3) = training_structure.forces; - noise_vector.segment(1, n_atoms * 3) = + labels.segment(label_count, n_atoms * 3) = training_structure.forces; + noise_vector.segment(label_count, n_atoms * 3) = Eigen::VectorXd::Constant(n_atoms * 3, 1 / (sigma_f * sigma_f)); } @@ -182,3 +184,31 @@ void SparseGP :: add_training_structure(StructureDescriptor training_structure){ noise.tail(n_labels) = noise_vector; noise_matrix = noise.asDiagonal(); } + +void SparseGP::update_alpha(){ + Eigen::MatrixXd sigma_inv = Kuu + Kuf * noise_matrix * Kuf.transpose(); + Sigma = sigma_inv.inverse(); + alpha = Sigma * Kuf * noise_matrix * y; +} + +Eigen::VectorXd SparseGP::predict(StructureDescriptor test_structure){ + int n_atoms = test_structure.noa; + int n_out = 1 + 3 * n_atoms + 6; + int n_sparse = sparse_environments.size(); + int n_kernels = kernels.size(); + Eigen::MatrixXd kern_mat = Eigen::MatrixXd::Zero(n_out, n_sparse); + + LocalEnvironment sparse_env; + Eigen::VectorXd kernel_vector; + for (int i = 0; i < n_sparse; i ++){ + sparse_env = sparse_environments[i]; + kernel_vector = Eigen::VectorXd::Zero(n_out); + for (int j = 0; j < n_kernels; j ++){ + kernel_vector += + kernels[j] -> env_struc(sparse_env, test_structure); + } + kern_mat.col(i) = kernel_vector; + } + + return kern_mat * alpha; +} diff --git a/src/sparse_gp.h b/src/sparse_gp.h index cb6bd692e..fc258db94 100644 --- a/src/sparse_gp.h +++ b/src/sparse_gp.h @@ -29,6 +29,8 @@ class SparseGP{ void add_sparse_environment(LocalEnvironment env); void add_training_structure(StructureDescriptor training_structure); + void update_alpha(); + Eigen::VectorXd predict(StructureDescriptor test_structure); }; #endif \ No newline at end of file diff --git a/tests/test_sparse_gp.cpp b/tests/test_sparse_gp.cpp index 7c1c300c0..a2c374e1a 100644 --- a/tests/test_sparse_gp.cpp +++ b/tests/test_sparse_gp.cpp @@ -8,9 +8,8 @@ class SparseTest : public ::testing::Test{ int n_atoms = 2; Eigen::MatrixXd cell{3, 3}, cell_2{3, 3}; std::vector species {0, 1}; - Eigen::MatrixXd positions{n_atoms, 3}, positions_2{n_atoms, 3}, - positions_3{n_atoms, 3}; - StructureDescriptor test_struc, test_struc_2, test_struc_3; + Eigen::MatrixXd positions{n_atoms, 3}; + StructureDescriptor test_struc; // labels Eigen::VectorXd energy{1}, forces{n_atoms * 3}, stresses{6}; @@ -26,9 +25,9 @@ class SparseTest : public ::testing::Test{ B2_Calculator desc1; // kernel - double signal_variance = 2; + double signal_variance = 1; double length_scale = 1; - double power = 2; + double power = 1; DotProductKernel kernel; TwoBodyKernel two_body_kernel; ThreeBodyKernel three_body_kernel; @@ -36,18 +35,15 @@ class SparseTest : public ::testing::Test{ std::vector kernels; SparseTest(){ - cell << 10, 0, 0, - 0, 10, 0, - 0, 0, 10; + cell << 100, 0, 0, + 0, 100, 0, + 0, 0, 100; - positions << -0.68402216, 0.54343671, -0.52961224, - 0.33045915, 0.40010388, 0.59849816; + positions << 0, 0, 0, + 1, 0, 0; - positions_2 << 0.69955637, -0.41619112, -0.51725003, - 0.43189622, 0.88548458, -0.74495343; - - energy = Eigen::VectorXd::Random(1); - forces = Eigen::VectorXd::Random(n_atoms * 3); + energy << 0.01; + forces << -1, 0, 0, 1, 0, 0; stresses = Eigen::VectorXd::Random(6); desc1 = B2_Calculator(radial_string, cutoff_string, @@ -65,9 +61,12 @@ class SparseTest : public ::testing::Test{ cutoff_string, cutoff_hyps); many_body_kernel = DotProductKernel(signal_variance, power); + // kernels = + // std::vector {&two_body_kernel, &three_body_kernel, + // &many_body_kernel}; + kernels = - std::vector {&two_body_kernel, &three_body_kernel, - &many_body_kernel}; + std::vector {&two_body_kernel}; } }; @@ -101,3 +100,32 @@ TEST_F(SparseTest, UpdateK){ } EXPECT_EQ(kern_val, sparse_gp.Kuu(0,0)); } + +TEST_F(SparseTest, Predict){ + double sigma_e = 0.1; + double sigma_f = 0.01; + double sigma_s = 1; + + SparseGP sparse_gp = SparseGP(kernels, sigma_e, sigma_f, sigma_s); + LocalEnvironment env1 = test_struc.local_environments[0]; + LocalEnvironment env2 = test_struc.local_environments[1]; + sparse_gp.add_sparse_environment(env1); + + // test_struc.energy = Eigen::VectorXd {}; + // test_struc.forces = Eigen::VectorXd {}; + test_struc.stresses = Eigen::VectorXd {}; + + sparse_gp.add_training_structure(test_struc); + + sparse_gp.update_alpha(); + + Eigen::VectorXd pred_vals = sparse_gp.predict(test_struc); + + std::cout << "predicted values:" << std::endl; + std::cout << pred_vals << std::endl; + + std::cout << sparse_gp.Kuu << std::endl; + std::cout << sparse_gp.Kuf << std::endl; + std::cout << sparse_gp.Sigma << std::endl; + std::cout << sparse_gp.alpha << std::endl; +} \ No newline at end of file From 926436e1a03b2089efd557481040b86fd729fef7 Mon Sep 17 00:00:00 2001 From: Jonathan Vandermause Date: Sun, 23 Feb 2020 15:40:31 -0500 Subject: [PATCH 08/24] start kernel and sparse gp bindings --- src/ace_binding.cpp | 69 ++++++++++++++++------------------------ src/kernels.h | 2 ++ tests/test_sparse_gp.cpp | 12 +++---- 3 files changed, 35 insertions(+), 48 deletions(-) diff --git a/src/ace_binding.cpp b/src/ace_binding.cpp index e7b45ab3e..f4b5825eb 100644 --- a/src/ace_binding.cpp +++ b/src/ace_binding.cpp @@ -2,6 +2,8 @@ #include "structure.h" #include "local_environment.h" #include "descriptor.h" +#include "sparse_gp.h" +#include "kernels.h" #include #include @@ -12,49 +14,8 @@ namespace py = pybind11; -// Define spherical harmonics class. -class SphericalHarmonics{ - public: - double x, y, z; - int lmax; - std::vector Y, Yx, Yy, Yz; - - SphericalHarmonics(double x, double y, double z, int lmax); -}; - -SphericalHarmonics :: SphericalHarmonics(double x, double y, double z, - int lmax){ - int number_of_harmonics = (lmax + 1) * (lmax + 1); - - this->x = x; - this->y = y; - this->z = z; - this->lmax = lmax; - - // Initialize spherical harmonic vectors. - Y = std::vector(number_of_harmonics, 0); - Yx = std::vector(number_of_harmonics, 0); - Yy = std::vector(number_of_harmonics, 0); - Yz = std::vector(number_of_harmonics, 0); - - get_Y(Y, Yx, Yy, Yz, x, y, z, lmax); -}; - PYBIND11_MODULE(ace, m){ - // Bind the spherical harmonics class. - py::class_(m, "SphericalHarmonics") - .def(py::init()) - // Make attributes accessible. - .def_readwrite("x", &SphericalHarmonics::x) - .def_readwrite("y", &SphericalHarmonics::y) - .def_readwrite("z", &SphericalHarmonics::z) - .def_readwrite("lmax", &SphericalHarmonics::lmax) - .def_readwrite("Y", &SphericalHarmonics::Y) - .def_readwrite("Yx", &SphericalHarmonics::Yx) - .def_readwrite("Yy", &SphericalHarmonics::Yy) - .def_readwrite("Yz", &SphericalHarmonics::Yz); - - // Bind the structure class. + // Structure py::class_(m, "Structure") .def(py::init &, @@ -67,6 +28,7 @@ PYBIND11_MODULE(ace, m){ .def_readwrite("volume", &Structure::volume) .def("wrap_positions", &Structure::wrap_positions); + // Local environment py::class_(m, "LocalEnvironment") .def(py::init()) .def_readwrite("sweep", &LocalEnvironment::sweep) @@ -84,6 +46,7 @@ PYBIND11_MODULE(ace, m){ .def_readwrite("zs", &LocalEnvironment::zs) .def_readwrite("cutoff", &LocalEnvironment::cutoff); + // Descriptor calculators py::class_(m, "B2_Calculator") .def(py::init &, const std::vector &, @@ -96,4 +59,26 @@ PYBIND11_MODULE(ace, m){ .def_readwrite("descriptor_stress_dervs", &B2_Calculator::descriptor_stress_dervs) .def("compute", &B2_Calculator::compute); + + // Kernel functions + py::class_(m, "Kernel") + .def("env_env", &Kernel::env_env) + .def("env_struc", &Kernel::env_struc); + + py::class_(m, "TwoBodyKernel") + .def(py::init>()); + + py::class_(m, "ThreeBodyKernel") + .def(py::init>()); + + py::class_(m, "DotProductKernel") + .def(py::init()); + + // Sparse GP + py::class_(m, "SparseGP") + .def(py::init, double, double, double>()) + .def("add_sparse_environment", &SparseGP::add_sparse_environment) + .def("add_training_structure", &SparseGP::add_training_structure) + .def("update_alpha", &SparseGP::update_alpha) + .def("predict", &SparseGP::predict); } diff --git a/src/kernels.h b/src/kernels.h index ab632a6f0..caf0bb506 100644 --- a/src/kernels.h +++ b/src/kernels.h @@ -20,6 +20,8 @@ class Kernel{ virtual Eigen::VectorXd env_struc(const LocalEnvironment & env1, const StructureDescriptor & struc1) = 0; + + virtual ~Kernel() = default; }; class DotProductKernel : public Kernel{ diff --git a/tests/test_sparse_gp.cpp b/tests/test_sparse_gp.cpp index a2c374e1a..4d00befdd 100644 --- a/tests/test_sparse_gp.cpp +++ b/tests/test_sparse_gp.cpp @@ -121,11 +121,11 @@ TEST_F(SparseTest, Predict){ Eigen::VectorXd pred_vals = sparse_gp.predict(test_struc); - std::cout << "predicted values:" << std::endl; - std::cout << pred_vals << std::endl; + // std::cout << "predicted values:" << std::endl; + // std::cout << pred_vals << std::endl; - std::cout << sparse_gp.Kuu << std::endl; - std::cout << sparse_gp.Kuf << std::endl; - std::cout << sparse_gp.Sigma << std::endl; - std::cout << sparse_gp.alpha << std::endl; + // std::cout << sparse_gp.Kuu << std::endl; + // std::cout << sparse_gp.Kuf << std::endl; + // std::cout << sparse_gp.Sigma << std::endl; + // std::cout << sparse_gp.alpha << std::endl; } \ No newline at end of file From d8a612426dbb009c47b0da5abd936739b1a7480d Mon Sep 17 00:00:00 2001 From: Jonathan Vandermause Date: Sun, 23 Feb 2020 19:29:21 -0500 Subject: [PATCH 09/24] more python bindings --- src/ace_binding.cpp | 63 +++++++++++++++++++++++++++++++++++++-------- src/descriptor.h | 2 ++ 2 files changed, 54 insertions(+), 11 deletions(-) diff --git a/src/ace_binding.cpp b/src/ace_binding.cpp index f4b5825eb..647dfe518 100644 --- a/src/ace_binding.cpp +++ b/src/ace_binding.cpp @@ -28,6 +28,20 @@ PYBIND11_MODULE(ace, m){ .def_readwrite("volume", &Structure::volume) .def("wrap_positions", &Structure::wrap_positions); + py::class_(m, "StructureDescriptor") + .def(py::init &, + const Eigen::MatrixXd &, double>()) + .def(py::init &, + const Eigen::MatrixXd &, double, + DescriptorCalculator *>()) + .def_readwrite("local_environments", + &StructureDescriptor::local_environments) + .def_readwrite("energy", &StructureDescriptor::energy) + .def_readwrite("forces", &StructureDescriptor::forces) + .def_readwrite("stresses", &StructureDescriptor::stresses) + .def_readwrite("cutoff", &StructureDescriptor::cutoff) + .def_readwrite("nested_cutoffs", &StructureDescriptor::nested_cutoffs); + // Local environment py::class_(m, "LocalEnvironment") .def(py::init()) @@ -44,21 +58,39 @@ PYBIND11_MODULE(ace, m){ .def_readwrite("xs", &LocalEnvironment::xs) .def_readwrite("ys", &LocalEnvironment::ys) .def_readwrite("zs", &LocalEnvironment::zs) - .def_readwrite("cutoff", &LocalEnvironment::cutoff); + .def_readwrite("cutoff", &LocalEnvironment::cutoff) + .def_readwrite("descriptor_vals", &LocalEnvironment::descriptor_vals) + .def_readwrite("descriptor_force_dervs", + &LocalEnvironment::descriptor_force_dervs) + .def_readwrite("descriptor_stress_dervs", + &LocalEnvironment::descriptor_stress_dervs); + Eigen::VectorXd descriptor_vals; + Eigen::MatrixXd descriptor_force_dervs, descriptor_stress_dervs, + force_dot, stress_dot; // Descriptor calculators - py::class_(m, "B2_Calculator") - .def(py::init &, const std::vector &, - const std::vector &>()) - .def_readwrite("radial_basis", &B2_Calculator::radial_basis) - .def_readwrite("cutoff_function", &B2_Calculator::cutoff_function) - .def_readwrite("descriptor_vals", &B2_Calculator::descriptor_vals) + py::class_(m, "DescriptorCalculator") + .def("compute", &DescriptorCalculator::compute) + .def_readwrite("radial_basis", &DescriptorCalculator::radial_basis) + .def_readwrite("cutoff_function", + &DescriptorCalculator::cutoff_function) + .def_readwrite("descriptor_vals", + &DescriptorCalculator::descriptor_vals) .def_readwrite("descriptor_force_dervs", - &B2_Calculator::descriptor_force_dervs) + &DescriptorCalculator::descriptor_force_dervs) .def_readwrite("descriptor_stress_dervs", - &B2_Calculator::descriptor_stress_dervs) - .def("compute", &B2_Calculator::compute); + &DescriptorCalculator::descriptor_stress_dervs) + .def_readwrite("radial_hyps", + &DescriptorCalculator::radial_hyps) + .def_readwrite("cutoff_hyps", + &DescriptorCalculator::cutoff_hyps) + .def_readwrite("descriptor_settings", + &DescriptorCalculator::descriptor_settings); + + py::class_(m, "B2_Calculator") + .def(py::init &, const std::vector &, + const std::vector &>()); // Kernel functions py::class_(m, "Kernel") @@ -77,6 +109,15 @@ PYBIND11_MODULE(ace, m){ // Sparse GP py::class_(m, "SparseGP") .def(py::init, double, double, double>()) + .def_readwrite("Kuu", &SparseGP::Kuu) + .def_readwrite("Kuf", &SparseGP::Kuf) + .def_readwrite("Sigma", &SparseGP::Sigma) + .def_readwrite("noise", &SparseGP::noise) + .def_readwrite("noise_matrix", &SparseGP::noise_matrix) + .def_readwrite("kernels", &SparseGP::kernels) + .def_readwrite("alpha", &SparseGP::alpha) + .def_readwrite("y", &SparseGP::y) + .def_readwrite("hyperparameters", &SparseGP::hyperparameters) .def("add_sparse_environment", &SparseGP::add_sparse_environment) .def("add_training_structure", &SparseGP::add_training_structure) .def("update_alpha", &SparseGP::update_alpha) diff --git a/src/descriptor.h b/src/descriptor.h index b00f67a48..85b4fb7d3 100644 --- a/src/descriptor.h +++ b/src/descriptor.h @@ -29,6 +29,8 @@ class DescriptorCalculator{ const std::vector & descriptor_settings); virtual void compute(const LocalEnvironment & env) = 0; + + virtual ~DescriptorCalculator() = default; }; void B2_descriptor( From 58f3ba8f2d6cba3852f7a1684eab0c80371149f6 Mon Sep 17 00:00:00 2001 From: Jonathan Vandermause Date: Mon, 24 Feb 2020 16:25:27 -0500 Subject: [PATCH 10/24] add learning curve script --- src/ace_binding.cpp | 3 + timing/learning_curve.py | 122 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 125 insertions(+) create mode 100644 timing/learning_curve.py diff --git a/src/ace_binding.cpp b/src/ace_binding.cpp index 647dfe518..b00f9f232 100644 --- a/src/ace_binding.cpp +++ b/src/ace_binding.cpp @@ -34,6 +34,9 @@ PYBIND11_MODULE(ace, m){ .def(py::init &, const Eigen::MatrixXd &, double, DescriptorCalculator *>()) + .def(py::init &, + const Eigen::MatrixXd &, double, std::vector, + DescriptorCalculator *>()) .def_readwrite("local_environments", &StructureDescriptor::local_environments) .def_readwrite("energy", &StructureDescriptor::energy) diff --git a/timing/learning_curve.py b/timing/learning_curve.py new file mode 100644 index 000000000..0dfa088df --- /dev/null +++ b/timing/learning_curve.py @@ -0,0 +1,122 @@ +import numpy as np +import ace +import pickle + +# create two and three body kernels +kernel_2b = ace.TwoBodyKernel(0.02, 0.60, "quadratic", []) +kernel_3b = ace.ThreeBodyKernel(0.013, 0.2, "quadratic", []) + +# create many body kernel +sigma = 0.1 +power = 4 +kernel_mb = ace.DotProductKernel(sigma, power) + +# create descriptor calculator +radial_basis = "chebyshev" +cutoff_function = "quadratic" +cutoff = 4. +nested_cutoffs = [4., 4., 4.] +radial_hyps = [0, 4.] +cutoff_hyps = [] +descriptor_settings = [5, 1, 1] +b2_calculator = \ + ace.B2_Calculator(radial_basis, cutoff_function, radial_hyps, + cutoff_hyps, descriptor_settings) + +# create sparse gp +sigma_e = 0.1 +sigma_f = 0.1 +sigma_s = 0.1 +sparse_gp = ace.SparseGP([kernel_3b, kernel_mb], + sigma_e, sigma_f, sigma_s) + +# define species (based on poscar) +# C, F, H, O, S +species = [0] * 16 + [1] * 36 + [2] * 164 + [3] * 96 + [4] * 4 + +# prepare data +file_base = '/Users/jonpvandermause/Research/Bosch/10-4-19/Sulfonate_AIMD/'\ + 'extract_run' +file_tail = '.pkl' +file_nos = [8, 1, 7, 2, 6, 3, 5] +snapshots = [100, 200, 300, 400, 500] +test_snapshot = 300 + +# create test structure +test_no = 4 +data_file = '{}{}{}'.format(file_base, test_no, file_tail) +open_file = open(data_file, 'rb') +vasp_data = pickle.load(open_file) +cell = vasp_data['lattice'][test_snapshot] +positions = vasp_data['atoms'][test_snapshot] +test_forces = vasp_data['forces'][test_snapshot] +test_struc = \ + ace.StructureDescriptor(cell, species, positions, cutoff, + nested_cutoffs, b2_calculator) +test_struc.forces = test_forces.reshape(-1) + + +# train the model +for snapshot in snapshots: + for file_no in file_nos: + data_file = '{}{}{}'.format(file_base, file_no, file_tail) + open_file = open(data_file, 'rb') + vasp_data = pickle.load(open_file) + cell = vasp_data['lattice'][snapshot] + positions = vasp_data['atoms'][snapshot] + forces = vasp_data['forces'][snapshot] + train_struc = \ + ace.StructureDescriptor(cell, species, positions, cutoff, + nested_cutoffs, b2_calculator) + train_struc.forces = forces.reshape(-1) + + # choose sparse points randomly (equal number for each species) + sparse_inds = np.array([], dtype=int) + sparse_inds = \ + np.append(sparse_inds, + np.random.choice(np.arange(0, 16, 1), + size=4, replace=False)) + sparse_inds = \ + np.append(sparse_inds, + np.random.choice(np.arange(16, 52, 1), + size=4, replace=False)) + sparse_inds = \ + np.append(sparse_inds, + np.random.choice(np.arange(52, 216, 1), + size=4, replace=False)) + sparse_inds = \ + np.append(sparse_inds, + np.random.choice(np.arange(216, 312, 1), + size=4, replace=False)) + sparse_inds = \ + np.append(sparse_inds, + np.random.choice(np.arange(312, 316, 1), + size=4, replace=False)) + + # add training structure and sparse environments + print('adding training structure...') + sparse_gp.add_training_structure(train_struc) + print('adding sparse environments...') + for ind in sparse_inds: + print(ind) + sparse_gp.add_sparse_environment(train_struc.local_environments[ind]) + + # compute test error + print('updating alpha...') + sparse_gp.update_alpha() + + print('predicting on test structure...') + force_pred = sparse_gp.predict(test_struc) + mae = np.mean(np.abs(test_forces.reshape(-1) - force_pred[1:-6])) + print(mae) + +# # save maes and stds +# maes = np.array(maes) +# stds = np.array(stds) +# np.save('maes', maes) +# np.save('stds', stds) + +# # pickle the gp model +# gp_name = 'Sulfonate_3p5.gp' +# gp_file = open(gp_name, 'wb') +# pickle.dump(gp_model, gp_file) From b5e63aa6a0d8b53519ed233b1ffc5264aa53cae5 Mon Sep 17 00:00:00 2001 From: Jonathan Vandermause Date: Mon, 24 Feb 2020 21:31:28 -0500 Subject: [PATCH 11/24] parallelize predict --- CMakeLists.txt | 5 +++++ src/sparse_gp.cpp | 2 ++ 2 files changed, 7 insertions(+) diff --git a/CMakeLists.txt b/CMakeLists.txt index db1981926..71a359760 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -16,6 +16,11 @@ endif() set(CMAKE_CXX_FLAGS_RELEASE "-O3") +find_package(OpenMP) +if(OpenMP_CXX_FOUND) + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}") +endif() + find_package(Eigen3 REQUIRED) INCLUDE_DIRECTORIES ("${EIGEN3_INCLUDE_DIR}") diff --git a/src/sparse_gp.cpp b/src/sparse_gp.cpp index 566256fa3..8eb7a63a3 100644 --- a/src/sparse_gp.cpp +++ b/src/sparse_gp.cpp @@ -200,6 +200,8 @@ Eigen::VectorXd SparseGP::predict(StructureDescriptor test_structure){ LocalEnvironment sparse_env; Eigen::VectorXd kernel_vector; + + #pragma omp parallel for for (int i = 0; i < n_sparse; i ++){ sparse_env = sparse_environments[i]; kernel_vector = Eigen::VectorXd::Zero(n_out); From 13eb30e259457482f9df20782eefe355095efdf9 Mon Sep 17 00:00:00 2001 From: Jonathan Vandermause Date: Mon, 24 Feb 2020 21:38:19 -0500 Subject: [PATCH 12/24] parallelize add_training_structure --- src/sparse_gp.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/sparse_gp.cpp b/src/sparse_gp.cpp index 8eb7a63a3..518d03657 100644 --- a/src/sparse_gp.cpp +++ b/src/sparse_gp.cpp @@ -121,6 +121,8 @@ void SparseGP :: add_training_structure(StructureDescriptor training_structure){ LocalEnvironment sparse_env; Eigen::VectorXd kernel_vector; int label_count; + + #pragma omp parallel for for (int i = 0; i < n_sparse; i ++){ sparse_env = sparse_environments[i]; kernel_vector = Eigen::VectorXd::Zero(1 + 3 * n_atoms + 6); From e3912d159c284b95b576e627d2fef84d248a563a Mon Sep 17 00:00:00 2001 From: Jonathan Vandermause Date: Tue, 3 Mar 2020 22:37:18 -0500 Subject: [PATCH 13/24] fix predict seg fault --- CMakeLists.txt | 6 +---- parallel_tests/CMakeLists.txt | 6 +++++ parallel_tests/par_test.cpp | 27 ++++++++++++++++++++ src/CMakeLists.txt | 3 +++ src/sparse_gp.cpp | 22 ++++++++++++++++- src/sparse_gp.h | 1 + tests/test_sparse_gp.cpp | 46 +++++++++++++++++++++++++---------- 7 files changed, 92 insertions(+), 19 deletions(-) create mode 100644 parallel_tests/CMakeLists.txt create mode 100644 parallel_tests/par_test.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 71a359760..c3943d04b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -16,11 +16,6 @@ endif() set(CMAKE_CXX_FLAGS_RELEASE "-O3") -find_package(OpenMP) -if(OpenMP_CXX_FOUND) - set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}") -endif() - find_package(Eigen3 REQUIRED) INCLUDE_DIRECTORIES ("${EIGEN3_INCLUDE_DIR}") @@ -36,6 +31,7 @@ add_subdirectory(tests) add_subdirectory(timing) add_subdirectory(pybind11) add_subdirectory(googletest) +add_subdirectory(parallel_tests) # make python bindings pybind11_add_module(ace src/ace_binding.cpp src/y_grad.cpp src/structure.cpp diff --git a/parallel_tests/CMakeLists.txt b/parallel_tests/CMakeLists.txt new file mode 100644 index 000000000..8f3bb90c0 --- /dev/null +++ b/parallel_tests/CMakeLists.txt @@ -0,0 +1,6 @@ +# make library +add_executable(par_test par_test.cpp) + +find_package(OpenMP REQUIRED) +target_link_libraries(par_test PRIVATE OpenMP::OpenMP_CXX) + diff --git a/parallel_tests/par_test.cpp b/parallel_tests/par_test.cpp new file mode 100644 index 000000000..e076a1023 --- /dev/null +++ b/parallel_tests/par_test.cpp @@ -0,0 +1,27 @@ +#include +#include +#include +#include +#include + +int main(){ + auto t1 = std::chrono::high_resolution_clock::now(); + + int mat_size = 5; + Eigen::MatrixXd kern_mat = Eigen::MatrixXd::Zero(mat_size, mat_size); + Eigen::VectorXd vec_curr(5); + + #pragma omp parallel for + for (int i = 0; i < mat_size; i ++){ + std::this_thread::sleep_for(std::chrono::seconds(1)); + vec_curr << i, i+1, i+2, i+3, i+4; + kern_mat.col(i) = vec_curr; + } + + auto t2 = std::chrono::high_resolution_clock::now(); + auto tot_time = + std::chrono::duration_cast(t2-t1).count(); + std::cout << tot_time << std::endl; + std::cout << kern_mat << std::endl; + +} diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index b862536ae..7aaadba59 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -3,4 +3,7 @@ add_library(Ace y_grad.cpp radial.cpp cutoffs.cpp single_bond.cpp structure.cpp local_environment.cpp descriptor.cpp kernels.cpp sparse_gp.cpp) +find_package(OpenMP REQUIRED) +target_link_libraries(Ace PRIVATE OpenMP::OpenMP_CXX) + set(ACE_INCLUDE_DIR "${ACE_ROOT_DIR}/src" PARENT_SCOPE) diff --git a/src/sparse_gp.cpp b/src/sparse_gp.cpp index 518d03657..b2147f26f 100644 --- a/src/sparse_gp.cpp +++ b/src/sparse_gp.cpp @@ -204,6 +204,26 @@ Eigen::VectorXd SparseGP::predict(StructureDescriptor test_structure){ Eigen::VectorXd kernel_vector; #pragma omp parallel for + for (int i = 0; i < n_sparse; i ++){ + for (int j = 0; j < n_kernels; j ++){ + kern_mat.col(i) += + kernels[j] -> env_struc(sparse_environments[i], test_structure); + } + } + + return kern_mat * alpha; +} + +Eigen::VectorXd SparseGP::predict_serial(StructureDescriptor test_structure){ + int n_atoms = test_structure.noa; + int n_out = 1 + 3 * n_atoms + 6; + int n_sparse = sparse_environments.size(); + int n_kernels = kernels.size(); + Eigen::MatrixXd kern_mat = Eigen::MatrixXd::Zero(n_out, n_sparse); + + LocalEnvironment sparse_env; + Eigen::VectorXd kernel_vector; + for (int i = 0; i < n_sparse; i ++){ sparse_env = sparse_environments[i]; kernel_vector = Eigen::VectorXd::Zero(n_out); @@ -215,4 +235,4 @@ Eigen::VectorXd SparseGP::predict(StructureDescriptor test_structure){ } return kern_mat * alpha; -} +} \ No newline at end of file diff --git a/src/sparse_gp.h b/src/sparse_gp.h index fc258db94..c8d46019c 100644 --- a/src/sparse_gp.h +++ b/src/sparse_gp.h @@ -31,6 +31,7 @@ class SparseGP{ void add_training_structure(StructureDescriptor training_structure); void update_alpha(); Eigen::VectorXd predict(StructureDescriptor test_structure); + Eigen::VectorXd predict_serial(StructureDescriptor test_structure); }; #endif \ No newline at end of file diff --git a/tests/test_sparse_gp.cpp b/tests/test_sparse_gp.cpp index 4d00befdd..d5a349c0f 100644 --- a/tests/test_sparse_gp.cpp +++ b/tests/test_sparse_gp.cpp @@ -1,6 +1,7 @@ #include "gtest/gtest.h" #include "sparse_gp.h" #include "descriptor.h" +#include class SparseTest : public ::testing::Test{ public: @@ -19,7 +20,7 @@ class SparseTest : public ::testing::Test{ std::string cutoff_string = "cosine"; std::vector radial_hyps {0, 5}; std::vector cutoff_hyps; - std::vector descriptor_settings {2, 5, 5}; + std::vector descriptor_settings {2, 10, 10}; double cutoff = 5; std::vector nested_cutoffs {5, 5}; B2_Calculator desc1; @@ -61,12 +62,12 @@ class SparseTest : public ::testing::Test{ cutoff_string, cutoff_hyps); many_body_kernel = DotProductKernel(signal_variance, power); - // kernels = - // std::vector {&two_body_kernel, &three_body_kernel, - // &many_body_kernel}; - kernels = - std::vector {&two_body_kernel}; + std::vector {&two_body_kernel, &three_body_kernel, + &many_body_kernel}; + + // kernels = + // std::vector {&two_body_kernel}; } }; @@ -102,6 +103,8 @@ TEST_F(SparseTest, UpdateK){ } TEST_F(SparseTest, Predict){ + // Eigen::initParallel(); + double sigma_e = 0.1; double sigma_f = 0.01; double sigma_s = 1; @@ -110,19 +113,36 @@ TEST_F(SparseTest, Predict){ LocalEnvironment env1 = test_struc.local_environments[0]; LocalEnvironment env2 = test_struc.local_environments[1]; sparse_gp.add_sparse_environment(env1); + sparse_gp.add_sparse_environment(env2); - // test_struc.energy = Eigen::VectorXd {}; - // test_struc.forces = Eigen::VectorXd {}; test_struc.stresses = Eigen::VectorXd {}; sparse_gp.add_training_structure(test_struc); - sparse_gp.update_alpha(); + Eigen::VectorXd pred_vals; - Eigen::VectorXd pred_vals = sparse_gp.predict(test_struc); - - // std::cout << "predicted values:" << std::endl; - // std::cout << pred_vals << std::endl; + // predict in parallel + auto t1 = std::chrono::high_resolution_clock::now(); + for (int i = 0; i < 1000; i ++){ + pred_vals = sparse_gp.predict(test_struc); + } + auto t2 = std::chrono::high_resolution_clock::now(); + auto tot_time = + std::chrono::duration_cast(t2-t1).count(); + std::cout << tot_time << std::endl; + + // // predict in serial + // t1 = std::chrono::high_resolution_clock::now(); + // for (int i = 0; i < 1000; i ++){ + // Eigen::VectorXd pred_vals = sparse_gp.predict_serial(test_struc); + // } + // t2 = std::chrono::high_resolution_clock::now(); + // tot_time = + // std::chrono::duration_cast(t2-t1).count(); + // std::cout << tot_time << std::endl; + + std::cout << "predicted values:" << std::endl; + std::cout << pred_vals << std::endl; // std::cout << sparse_gp.Kuu << std::endl; // std::cout << sparse_gp.Kuf << std::endl; From 1fdd7666d1ef9986671dc2f253f7f8ce90e56777 Mon Sep 17 00:00:00 2001 From: Jonathan Vandermause Date: Fri, 6 Mar 2020 12:56:06 -0500 Subject: [PATCH 14/24] parallelize add_sparse_environment --- src/sparse_gp.cpp | 167 ++++++++++++++++++++++++++++++++++++++++++---- src/sparse_gp.h | 6 ++ 2 files changed, 161 insertions(+), 12 deletions(-) diff --git a/src/sparse_gp.cpp b/src/sparse_gp.cpp index b2147f26f..31c5a0a19 100644 --- a/src/sparse_gp.cpp +++ b/src/sparse_gp.cpp @@ -44,12 +44,11 @@ void SparseGP :: add_sparse_environment(LocalEnvironment env){ int n_sparse = sparse_environments.size(); int n_kernels = kernels.size(); - LocalEnvironment sparse_env; Eigen::VectorXd uu_vector = Eigen::VectorXd::Zero(n_sparse + 1); + #pragma omp parallel for for (int i = 0; i < n_sparse; i ++){ - sparse_env = sparse_environments[i]; for (int j = 0; j < n_kernels; j ++){ - uu_vector(i) += kernels[j] -> env_env(sparse_env, env); + uu_vector(i) += kernels[j] -> env_env(sparse_environments[i], env); } } @@ -69,31 +68,96 @@ void SparseGP :: add_sparse_environment(LocalEnvironment env){ int n_labels = Kuf.cols(); int n_strucs = training_structures.size(); Eigen::VectorXd uf_vector = Eigen::VectorXd::Zero(n_labels); - StructureDescriptor training_structure; Eigen::VectorXd kernel_vector; int label_count = 0; int n_atoms; + #pragma omp parallel for for (int i = 0; i < n_strucs; i ++){ - training_structure = training_structures[i]; - n_atoms = training_structure.noa; + n_atoms = training_structures[i].noa; kernel_vector = Eigen::VectorXd::Zero(1 + 3 * n_atoms + 6); for (int j = 0; j < n_kernels; j ++){ kernel_vector += kernels[j] -> env_struc(env, - training_structure); + training_structures[i]); } - if (training_structure.energy.size() != 0){ + if (training_structures[i].energy.size() != 0){ uf_vector(label_count) = kernel_vector(0); label_count += 1; } - if (training_structure.forces.size() != 0){ + if (training_structures[i].forces.size() != 0){ uf_vector.segment(label_count, n_atoms * 3) = kernel_vector.segment(1, n_atoms * 3); label_count += n_atoms * 3; } - if (training_structure.stresses.size() != 0){ + if (training_structures[i].stresses.size() != 0){ + uf_vector.segment(label_count, 6) = + kernel_vector.tail(6); + label_count += 6; + } + } + + // Update Kuf matrix. + Kuf.conservativeResize(Kuf.rows()+1, Kuf.cols()); + Kuf.row(Kuf.rows()-1) = uf_vector; + + // Store sparse environment. + sparse_environments.push_back(env); +} + +void SparseGP :: add_sparse_environment_serial(LocalEnvironment env){ + // Compute kernels between new environment and previous sparse + // environments. + int n_sparse = sparse_environments.size(); + int n_kernels = kernels.size(); + + Eigen::VectorXd uu_vector = Eigen::VectorXd::Zero(n_sparse + 1); + for (int i = 0; i < n_sparse; i ++){ + for (int j = 0; j < n_kernels; j ++){ + uu_vector(i) += kernels[j] -> env_env(sparse_environments[i], env); + } + } + + // Compute self kernel. + double self_kernel = 0; + for (int j = 0; j < n_kernels; j ++){ + self_kernel += kernels[j] -> env_env(env, env); + } + + // Update Kuu matrix. + Kuu.conservativeResize(Kuu.rows()+1, Kuu.cols()+1); + Kuu.col(Kuu.cols()-1) = uu_vector; + Kuu.row(Kuu.rows()-1) = uu_vector; + Kuu(n_sparse, n_sparse) = self_kernel; + + // Compute kernels between new environment and training structures. + int n_labels = Kuf.cols(); + int n_strucs = training_structures.size(); + Eigen::VectorXd uf_vector = Eigen::VectorXd::Zero(n_labels); + Eigen::VectorXd kernel_vector; + int label_count = 0; + int n_atoms; + for (int i = 0; i < n_strucs; i ++){ + n_atoms = training_structures[i].noa; + kernel_vector = Eigen::VectorXd::Zero(1 + 3 * n_atoms + 6); + for (int j = 0; j < n_kernels; j ++){ + kernel_vector += kernels[j] -> env_struc(env, + training_structures[i]); + } + + if (training_structures[i].energy.size() != 0){ + uf_vector(label_count) = kernel_vector(0); + label_count += 1; + } + + if (training_structures[i].forces.size() != 0){ + uf_vector.segment(label_count, n_atoms * 3) = + kernel_vector.segment(1, n_atoms * 3); + label_count += n_atoms * 3; + } + + if (training_structures[i].stresses.size() != 0){ uf_vector.segment(label_count, 6) = kernel_vector.tail(6); label_count += 6; @@ -108,6 +172,84 @@ void SparseGP :: add_sparse_environment(LocalEnvironment env){ sparse_environments.push_back(env); } +void SparseGP :: add_training_structure_serial(StructureDescriptor + training_structure){ + + int n_labels = training_structure.energy.size() + + training_structure.forces.size() + training_structure.stresses.size(); + int n_atoms = training_structure.noa; + int n_sparse = sparse_environments.size(); + int n_kernels = kernels.size(); + + // Calculate kernels between sparse environments and training structure. + Eigen::MatrixXd kernel_block = Eigen::MatrixXd::Zero(n_sparse, n_labels); + LocalEnvironment sparse_env; + Eigen::VectorXd kernel_vector; + int label_count; + + for (int i = 0; i < n_sparse; i ++){ + kernel_vector = Eigen::VectorXd::Zero(1 + 3 * n_atoms + 6); + for (int j = 0; j < n_kernels; j ++){ + kernel_vector += kernels[j] -> env_struc(sparse_environments[i], + training_structure); + } + + // Update kernel block. + label_count = 0; + if (training_structure.energy.size() != 0){ + kernel_block(i, 0) = kernel_vector(0); + label_count += 1; + } + + if (training_structure.forces.size() != 0){ + kernel_block.row(i).segment(label_count, n_atoms * 3) = + kernel_vector.segment(1, n_atoms * 3); + } + + if (training_structure.stresses.size() != 0){ + kernel_block.row(i).tail(6) = kernel_vector.tail(6); + } + } + + // Add kernel block to Kuf. + int prev_cols = Kuf.cols(); + Kuf.conservativeResize(n_sparse, prev_cols + n_labels); + Kuf.block(0, prev_cols, n_sparse, n_labels) = kernel_block; + + // Store training structure. + training_structures.push_back(training_structure); + + // Update y vector and noise matrix. + Eigen::VectorXd labels = Eigen::VectorXd::Zero(n_labels); + Eigen::VectorXd noise_vector = Eigen::VectorXd::Zero(n_labels); + + label_count = 0; + if (training_structure.energy.size() != 0){ + labels.head(1) = training_structure.energy; + noise_vector(0) = 1 / (sigma_e * sigma_e); + label_count ++; + } + + if (training_structure.forces.size() != 0){ + labels.segment(label_count, n_atoms * 3) = training_structure.forces; + noise_vector.segment(label_count, n_atoms * 3) = + Eigen::VectorXd::Constant(n_atoms * 3, 1 / (sigma_f * sigma_f)); + } + + if (training_structure.stresses.size() != 0){ + labels.tail(6) = training_structure.stresses; + noise_vector.tail(6) = + Eigen::VectorXd::Constant(6, 1 / (sigma_s * sigma_s)); + } + + y.conservativeResize(y.size() + n_labels); + y.tail(n_labels) = labels; + + noise.conservativeResize(prev_cols + n_labels); + noise.tail(n_labels) = noise_vector; + noise_matrix = noise.asDiagonal(); +} + void SparseGP :: add_training_structure(StructureDescriptor training_structure){ int n_labels = training_structure.energy.size() + @@ -124,10 +266,9 @@ void SparseGP :: add_training_structure(StructureDescriptor training_structure){ #pragma omp parallel for for (int i = 0; i < n_sparse; i ++){ - sparse_env = sparse_environments[i]; kernel_vector = Eigen::VectorXd::Zero(1 + 3 * n_atoms + 6); for (int j = 0; j < n_kernels; j ++){ - kernel_vector += kernels[j] -> env_struc(sparse_env, + kernel_vector += kernels[j] -> env_struc(sparse_environments[i], training_structure); } @@ -203,6 +344,8 @@ Eigen::VectorXd SparseGP::predict(StructureDescriptor test_structure){ LocalEnvironment sparse_env; Eigen::VectorXd kernel_vector; + // Compute the kernel between the test structure and each sparse + // environment, parallelizing over environments. #pragma omp parallel for for (int i = 0; i < n_sparse; i ++){ for (int j = 0; j < n_kernels; j ++){ diff --git a/src/sparse_gp.h b/src/sparse_gp.h index c8d46019c..79a100965 100644 --- a/src/sparse_gp.h +++ b/src/sparse_gp.h @@ -28,8 +28,14 @@ class SparseGP{ double sigma_s); void add_sparse_environment(LocalEnvironment env); + void add_sparse_environment_serial(LocalEnvironment env); + void add_training_structure(StructureDescriptor training_structure); + void add_training_structure_serial(StructureDescriptor + training_structure); + void update_alpha(); + Eigen::VectorXd predict(StructureDescriptor test_structure); Eigen::VectorXd predict_serial(StructureDescriptor test_structure); }; From 20161ed2ac6003b6067c696de8111ed88dd97e54 Mon Sep 17 00:00:00 2001 From: Jonathan Vandermause Date: Fri, 6 Mar 2020 14:13:21 -0500 Subject: [PATCH 15/24] fix structure constructor bug --- src/ace_binding.cpp | 9 +++- src/structure.cpp | 3 +- tests/test_sparse_gp.cpp | 106 +++++++++++++++++++-------------------- 3 files changed, 63 insertions(+), 55 deletions(-) diff --git a/src/ace_binding.cpp b/src/ace_binding.cpp index b00f9f232..c2e83eeec 100644 --- a/src/ace_binding.cpp +++ b/src/ace_binding.cpp @@ -34,6 +34,8 @@ PYBIND11_MODULE(ace, m){ .def(py::init &, const Eigen::MatrixXd &, double, DescriptorCalculator *>()) + .def(py::init &, + const Eigen::MatrixXd &, double, std::vector>()) .def(py::init &, const Eigen::MatrixXd &, double, std::vector, DescriptorCalculator *>()) @@ -122,7 +124,12 @@ PYBIND11_MODULE(ace, m){ .def_readwrite("y", &SparseGP::y) .def_readwrite("hyperparameters", &SparseGP::hyperparameters) .def("add_sparse_environment", &SparseGP::add_sparse_environment) + .def("add_sparse_environment_serial", + &SparseGP::add_sparse_environment_serial) .def("add_training_structure", &SparseGP::add_training_structure) + .def("add_training_structure_serial", + &SparseGP::add_training_structure_serial) .def("update_alpha", &SparseGP::update_alpha) - .def("predict", &SparseGP::predict); + .def("predict", &SparseGP::predict) + .def("predict_serial", &SparseGP::predict_serial); } diff --git a/src/structure.cpp b/src/structure.cpp index b1074f4ba..7fd118c47 100644 --- a/src/structure.cpp +++ b/src/structure.cpp @@ -87,7 +87,8 @@ StructureDescriptor :: StructureDescriptor(const Eigen::MatrixXd & cell, StructureDescriptor :: StructureDescriptor(const Eigen::MatrixXd & cell, const std::vector & species, const Eigen::MatrixXd & positions, - double cutoff, std::vector nested_cutoffs){ + double cutoff, std::vector nested_cutoffs) + : Structure(cell, species, positions){ this->cutoff = cutoff; this->nested_cutoffs = nested_cutoffs; diff --git a/tests/test_sparse_gp.cpp b/tests/test_sparse_gp.cpp index d5a349c0f..ddfd36188 100644 --- a/tests/test_sparse_gp.cpp +++ b/tests/test_sparse_gp.cpp @@ -51,7 +51,7 @@ class SparseTest : public ::testing::Test{ radial_hyps, cutoff_hyps, descriptor_settings); test_struc = StructureDescriptor(cell, species, positions, cutoff, - nested_cutoffs, &desc1); + nested_cutoffs); test_struc.energy = energy; test_struc.forces = forces; test_struc.stresses = stresses; @@ -62,12 +62,12 @@ class SparseTest : public ::testing::Test{ cutoff_string, cutoff_hyps); many_body_kernel = DotProductKernel(signal_variance, power); - kernels = - std::vector {&two_body_kernel, &three_body_kernel, - &many_body_kernel}; - // kernels = - // std::vector {&two_body_kernel}; + // std::vector {&two_body_kernel, &three_body_kernel, + // &many_body_kernel}; + + kernels = + std::vector {&two_body_kernel, &three_body_kernel}; } }; @@ -102,50 +102,50 @@ TEST_F(SparseTest, UpdateK){ EXPECT_EQ(kern_val, sparse_gp.Kuu(0,0)); } -TEST_F(SparseTest, Predict){ - // Eigen::initParallel(); - - double sigma_e = 0.1; - double sigma_f = 0.01; - double sigma_s = 1; - - SparseGP sparse_gp = SparseGP(kernels, sigma_e, sigma_f, sigma_s); - LocalEnvironment env1 = test_struc.local_environments[0]; - LocalEnvironment env2 = test_struc.local_environments[1]; - sparse_gp.add_sparse_environment(env1); - sparse_gp.add_sparse_environment(env2); - - test_struc.stresses = Eigen::VectorXd {}; - - sparse_gp.add_training_structure(test_struc); - sparse_gp.update_alpha(); - Eigen::VectorXd pred_vals; - - // predict in parallel - auto t1 = std::chrono::high_resolution_clock::now(); - for (int i = 0; i < 1000; i ++){ - pred_vals = sparse_gp.predict(test_struc); - } - auto t2 = std::chrono::high_resolution_clock::now(); - auto tot_time = - std::chrono::duration_cast(t2-t1).count(); - std::cout << tot_time << std::endl; - - // // predict in serial - // t1 = std::chrono::high_resolution_clock::now(); - // for (int i = 0; i < 1000; i ++){ - // Eigen::VectorXd pred_vals = sparse_gp.predict_serial(test_struc); - // } - // t2 = std::chrono::high_resolution_clock::now(); - // tot_time = - // std::chrono::duration_cast(t2-t1).count(); - // std::cout << tot_time << std::endl; - - std::cout << "predicted values:" << std::endl; - std::cout << pred_vals << std::endl; - - // std::cout << sparse_gp.Kuu << std::endl; - // std::cout << sparse_gp.Kuf << std::endl; - // std::cout << sparse_gp.Sigma << std::endl; - // std::cout << sparse_gp.alpha << std::endl; -} \ No newline at end of file +// TEST_F(SparseTest, Predict){ +// // Eigen::initParallel(); + +// double sigma_e = 0.1; +// double sigma_f = 0.01; +// double sigma_s = 1; + +// SparseGP sparse_gp = SparseGP(kernels, sigma_e, sigma_f, sigma_s); +// LocalEnvironment env1 = test_struc.local_environments[0]; +// LocalEnvironment env2 = test_struc.local_environments[1]; +// sparse_gp.add_sparse_environment(env1); +// sparse_gp.add_sparse_environment(env2); + +// test_struc.stresses = Eigen::VectorXd {}; + +// sparse_gp.add_training_structure(test_struc); +// sparse_gp.update_alpha(); +// Eigen::VectorXd pred_vals; + +// // predict in parallel +// auto t1 = std::chrono::high_resolution_clock::now(); +// for (int i = 0; i < 1000; i ++){ +// pred_vals = sparse_gp.predict(test_struc); +// } +// auto t2 = std::chrono::high_resolution_clock::now(); +// auto tot_time = +// std::chrono::duration_cast(t2-t1).count(); +// // std::cout << tot_time << std::endl; + +// // // predict in serial +// // t1 = std::chrono::high_resolution_clock::now(); +// // for (int i = 0; i < 1000; i ++){ +// // Eigen::VectorXd pred_vals = sparse_gp.predict_serial(test_struc); +// // } +// // t2 = std::chrono::high_resolution_clock::now(); +// // tot_time = +// // std::chrono::duration_cast(t2-t1).count(); +// // std::cout << tot_time << std::endl; + +// // std::cout << "predicted values:" << std::endl; +// // std::cout << pred_vals << std::endl; + +// // std::cout << sparse_gp.Kuu << std::endl; +// // std::cout << sparse_gp.Kuf << std::endl; +// // std::cout << sparse_gp.Sigma << std::endl; +// // std::cout << sparse_gp.alpha << std::endl; +// } \ No newline at end of file From c954c7093fda10cf6fb9f6b778287b19d1bee8d8 Mon Sep 17 00:00:00 2001 From: Jonathan Vandermause Date: Fri, 6 Mar 2020 16:03:49 -0500 Subject: [PATCH 16/24] link openmp and python library --- CMakeLists.txt | 2 ++ src/ace_binding.cpp | 5 +++-- 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index c3943d04b..0e480101b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -38,3 +38,5 @@ pybind11_add_module(ace src/ace_binding.cpp src/y_grad.cpp src/structure.cpp src/local_environment.cpp src/descriptor.cpp src/radial.cpp src/cutoffs.cpp src/single_bond.cpp src/kernels.cpp src/sparse_gp.cpp) +find_package(OpenMP REQUIRED) +target_link_libraries(ace PRIVATE OpenMP::OpenMP_CXX) diff --git a/src/ace_binding.cpp b/src/ace_binding.cpp index c2e83eeec..967b30c8f 100644 --- a/src/ace_binding.cpp +++ b/src/ace_binding.cpp @@ -123,9 +123,10 @@ PYBIND11_MODULE(ace, m){ .def_readwrite("alpha", &SparseGP::alpha) .def_readwrite("y", &SparseGP::y) .def_readwrite("hyperparameters", &SparseGP::hyperparameters) - .def("add_sparse_environment", &SparseGP::add_sparse_environment) + .def("add_sparse_environment", + &SparseGP::add_sparse_environment) .def("add_sparse_environment_serial", - &SparseGP::add_sparse_environment_serial) + &SparseGP::add_sparse_environment_serial) .def("add_training_structure", &SparseGP::add_training_structure) .def("add_training_structure_serial", &SparseGP::add_training_structure_serial) From 154a2ee6f8307ecdb31366d363413c8122159f7e Mon Sep 17 00:00:00 2001 From: Jonathan Vandermause Date: Sat, 7 Mar 2020 16:39:13 -0500 Subject: [PATCH 17/24] allow for multiple many body descriptors --- src/ace_binding.cpp | 24 +++++--- src/descriptor.cpp | 20 ++++--- src/descriptor.h | 11 ++-- src/kernels.cpp | 86 +++++++++++++++----------- src/kernels.h | 3 +- src/local_environment.cpp | 99 ++++++++++++++++++++---------- src/local_environment.h | 48 +++++++++------ src/structure.cpp | 36 +++++++---- src/structure.h | 23 ++++--- tests/test_descriptor.cpp | 15 ++++- tests/test_kernels.cpp | 37 ++++++++---- tests/test_local_environment.cpp | 16 +++-- tests/test_sparse_gp.cpp | 100 ++++++++++++++++--------------- tests/test_structure.cpp | 13 ++-- 14 files changed, 332 insertions(+), 199 deletions(-) diff --git a/src/ace_binding.cpp b/src/ace_binding.cpp index 967b30c8f..c93fde854 100644 --- a/src/ace_binding.cpp +++ b/src/ace_binding.cpp @@ -31,21 +31,27 @@ PYBIND11_MODULE(ace, m){ py::class_(m, "StructureDescriptor") .def(py::init &, const Eigen::MatrixXd &, double>()) - .def(py::init &, - const Eigen::MatrixXd &, double, - DescriptorCalculator *>()) + // n-body .def(py::init &, const Eigen::MatrixXd &, double, std::vector>()) + // many-body + .def(py::init &, + const Eigen::MatrixXd &, double, std::vector, + std::vector>()) + // n-body + many-body .def(py::init &, const Eigen::MatrixXd &, double, std::vector, - DescriptorCalculator *>()) + std::vector, + std::vector>()) .def_readwrite("local_environments", &StructureDescriptor::local_environments) .def_readwrite("energy", &StructureDescriptor::energy) .def_readwrite("forces", &StructureDescriptor::forces) .def_readwrite("stresses", &StructureDescriptor::stresses) .def_readwrite("cutoff", &StructureDescriptor::cutoff) - .def_readwrite("nested_cutoffs", &StructureDescriptor::nested_cutoffs); + .def_readwrite("n_body_cutoffs", &StructureDescriptor::n_body_cutoffs) + .def_readwrite("many_body_cutoffs", + &StructureDescriptor::many_body_cutoffs); // Local environment py::class_(m, "LocalEnvironment") @@ -70,9 +76,9 @@ PYBIND11_MODULE(ace, m){ .def_readwrite("descriptor_stress_dervs", &LocalEnvironment::descriptor_stress_dervs); - Eigen::VectorXd descriptor_vals; - Eigen::MatrixXd descriptor_force_dervs, descriptor_stress_dervs, - force_dot, stress_dot; + // Eigen::VectorXd descriptor_vals; + // Eigen::MatrixXd descriptor_force_dervs, descriptor_stress_dervs, + // force_dot, stress_dot; // Descriptor calculators py::class_(m, "DescriptorCalculator") .def("compute", &DescriptorCalculator::compute) @@ -109,7 +115,7 @@ PYBIND11_MODULE(ace, m){ .def(py::init>()); py::class_(m, "DotProductKernel") - .def(py::init()); + .def(py::init()); // Sparse GP py::class_(m, "SparseGP") diff --git a/src/descriptor.cpp b/src/descriptor.cpp index a32c3abf4..b7b4c5366 100644 --- a/src/descriptor.cpp +++ b/src/descriptor.cpp @@ -12,13 +12,15 @@ DescriptorCalculator::DescriptorCalculator( const std::string & radial_basis, const std::string & cutoff_function, const std::vector & radial_hyps, const std::vector & cutoff_hyps, - const std::vector & descriptor_settings){ + const std::vector & descriptor_settings, + int descriptor_index){ this->radial_basis = radial_basis; this->cutoff_function = cutoff_function; this->radial_hyps = radial_hyps; this->cutoff_hyps = cutoff_hyps; this->descriptor_settings = descriptor_settings; + this->descriptor_index = descriptor_index; if (radial_basis == "chebyshev"){ this->radial_pointer = chebyshev; @@ -102,9 +104,10 @@ B1_Calculator :: B1_Calculator(const std::string & radial_basis, const std::string & cutoff_function, const std::vector & radial_hyps, const std::vector & cutoff_hyps, - const std::vector & descriptor_settings) + const std::vector & descriptor_settings, + int descriptor_index) : DescriptorCalculator(radial_basis, cutoff_function, radial_hyps, - cutoff_hyps, descriptor_settings){} + cutoff_hyps, descriptor_settings, descriptor_index){} void B1_Calculator :: compute(const LocalEnvironment & env){ // Initialize single bond vectors. @@ -122,7 +125,8 @@ void B1_Calculator :: compute(const LocalEnvironment & env){ // Compute single bond vector. single_bond_sum_env(single_bond_vals, single_bond_force_dervs, single_bond_stress_dervs, radial_pointer, - cutoff_pointer, env, env.cutoff, N, + cutoff_pointer, env, + env.many_body_cutoffs[descriptor_index], N, lmax, radial_hyps, cutoff_hyps); // Set B1 values. @@ -137,9 +141,10 @@ B2_Calculator :: B2_Calculator(const std::string & radial_basis, const std::string & cutoff_function, const std::vector & radial_hyps, const std::vector & cutoff_hyps, - const std::vector & descriptor_settings) + const std::vector & descriptor_settings, + int descriptor_index) : DescriptorCalculator(radial_basis, cutoff_function, radial_hyps, - cutoff_hyps, descriptor_settings){} + cutoff_hyps, descriptor_settings, descriptor_index){} void B2_Calculator :: compute(const LocalEnvironment & env){ // Initialize single bond vectors. @@ -161,7 +166,8 @@ void B2_Calculator :: compute(const LocalEnvironment & env){ // Compute single bond vector. single_bond_sum_env(single_bond_vals, single_bond_force_dervs, single_bond_stress_dervs, radial_pointer, - cutoff_pointer, env, env.cutoff, N, + cutoff_pointer, env, + env.many_body_cutoffs[descriptor_index], N, lmax, radial_hyps, cutoff_hyps); // Initialize B2 vectors. diff --git a/src/descriptor.h b/src/descriptor.h index 85b4fb7d3..58ea0b76c 100644 --- a/src/descriptor.h +++ b/src/descriptor.h @@ -4,9 +4,9 @@ #include #include -// Descriptor calculator. class LocalEnvironment; +// Descriptor calculator. class DescriptorCalculator{ protected: void (*radial_pointer)(double *, double *, double, int, @@ -19,6 +19,7 @@ class DescriptorCalculator{ std::string radial_basis, cutoff_function; std::vector radial_hyps, cutoff_hyps; std::vector descriptor_settings; + int descriptor_index; DescriptorCalculator(); @@ -26,7 +27,7 @@ class DescriptorCalculator{ const std::string & radial_basis, const std::string & cutoff_function, const std::vector & radial_hyps, const std::vector & cutoff_hyps, - const std::vector & descriptor_settings); + const std::vector & descriptor_settings, int descriptor_index); virtual void compute(const LocalEnvironment & env) = 0; @@ -51,7 +52,8 @@ class B1_Calculator : public DescriptorCalculator{ const std::string & cutoff_function, const std::vector & radial_hyps, const std::vector & cutoff_hyps, - const std::vector & descriptor_settings); + const std::vector & descriptor_settings, + int descriptor_index); void compute(const LocalEnvironment & env); }; @@ -65,7 +67,8 @@ class B2_Calculator : public DescriptorCalculator{ const std::string & cutoff_function, const std::vector & radial_hyps, const std::vector & cutoff_hyps, - const std::vector & descriptor_settings); + const std::vector & descriptor_settings, + int descriptor_index); void compute(const LocalEnvironment & env); }; diff --git a/src/kernels.cpp b/src/kernels.cpp index b3364df25..699cebf41 100644 --- a/src/kernels.cpp +++ b/src/kernels.cpp @@ -41,26 +41,30 @@ double TwoBodyKernel :: env_env(const LocalEnvironment & env1, const LocalEnvironment & env2){ double kern = 0; double ri, rj, fi, fj, rdiff; - int e1, e2; + int e1, e2, ind1, ind2; - double cut1 = env1.cutoff; - double cut2 = env2.cutoff; + double cut1 = env1.n_body_cutoffs[0]; + double cut2 = env2.n_body_cutoffs[0]; double rcut_vals_1[2]; double rcut_vals_2[2]; int c1 = env1.central_species; int c2 = env2.central_species; + std::vector inds1 = env1.n_body_indices[0]; + std::vector inds2 = env2.n_body_indices[0]; - for (int m = 0; m < env1.rs.size(); m ++){ - ri = env1.rs[m]; - e1 = env1.environment_species[m]; + for (int m = 0; m < inds1.size(); m ++){ + ind1 = inds1[m]; + ri = env1.rs[ind1]; + e1 = env1.environment_species[ind1]; (*cutoff_pointer)(rcut_vals_1, ri, cut1, cutoff_hyps); fi = rcut_vals_1[0]; - for (int n = 0; n < env2.rs.size(); n ++){ - e2 = env2.environment_species[n]; + for (int n = 0; n < inds2.size(); n ++){ + ind2 = inds2[n]; + e2 = env2.environment_species[ind2]; // Proceed only if the pairs match. if ((c1 == c2 && e1 == e2) || (c1 == e2 && c2 == e1)){ - rj = env2.rs[n]; + rj = env2.rs[ind2]; (*cutoff_pointer)(rcut_vals_2, rj, cut2, cutoff_hyps); fj = rcut_vals_2[0]; rdiff = ri - rj; @@ -80,42 +84,48 @@ Eigen::VectorXd TwoBodyKernel :: env_struc(const LocalEnvironment & env1, double ri, rj, fi, fj, fdj, rdiff, c1, c2, c3, c4, fx, fy, fz, xval, yval, zval, xrel, yrel, zrel, en_kern; - int e1, e2, cent2; + int e1, e2, cent2, ind1, ind2; int cent1 = env1.central_species; - double cut1 = env1.cutoff; - double cut2 = struc1.cutoff; + double cut1 = env1.n_body_cutoffs[0]; + double cut2 = struc1.n_body_cutoffs[0]; double rcut_vals_1[2]; double rcut_vals_2[2]; + std::vector inds1 = env1.n_body_indices[0]; + std::vector inds2; + double vol_inv = 1 / struc1.volume; LocalEnvironment env_curr; for (int i = 0; i < struc1.noa; i ++){ env_curr = struc1.local_environments[i]; + inds2 = env_curr.n_body_indices[0]; cent2 = env_curr.central_species; - for (int m = 0; m < env1.rs.size(); m ++){ - ri = env1.rs[m]; + for (int m = 0; m < inds1.size(); m ++){ + ind1 = inds1[m]; + ri = env1.rs[ind1]; (*cutoff_pointer)(rcut_vals_1, ri, cut1, cutoff_hyps); fi = rcut_vals_1[0]; - e1 = env1.environment_species[m]; + e1 = env1.environment_species[ind1]; - for (int n = 0; n < env_curr.rs.size(); n ++){ - e2 = env_curr.environment_species[n]; + for (int n = 0; n < inds2.size(); n ++){ + ind2 = inds2[n]; + e2 = env_curr.environment_species[ind2]; // Proceed only if the pairs match. if ((cent1 == cent2 && e1 == e2) || (cent1 == e2 && cent2 == e1)){ - rj = env_curr.rs[n]; + rj = env_curr.rs[ind2]; rdiff = ri - rj; - xval = env_curr.xs[n]; - yval = env_curr.ys[n]; - zval = env_curr.zs[n]; - xrel = env_curr.xrel[n]; - yrel = env_curr.yrel[n]; - zrel = env_curr.zrel[n]; + xval = env_curr.xs[ind2]; + yval = env_curr.ys[ind2]; + zval = env_curr.zs[ind2]; + xrel = env_curr.xrel[ind2]; + yrel = env_curr.yrel[ind2]; + zrel = env_curr.zrel[ind2]; (*cutoff_pointer)(rcut_vals_2, rj, cut2, cutoff_hyps); fj = rcut_vals_2[0]; @@ -498,11 +508,13 @@ void ThreeBodyKernel :: update_kernel_vector(Eigen::VectorXd & kernel_vector, DotProductKernel :: DotProductKernel() {}; -DotProductKernel :: DotProductKernel(double sigma, double power){ +DotProductKernel :: DotProductKernel(double sigma, double power, + int descriptor_index){ this->sigma = sigma; sig2 = sigma * sigma; this->power = power; + this -> descriptor_index = descriptor_index; kernel_hyperparameters = std::vector {sigma, power}; } @@ -512,9 +524,10 @@ double DotProductKernel :: env_env(const LocalEnvironment & env1, // Central species must be the same to give a nonzero kernel. if (env1.central_species != env2.central_species) return 0; - double dot = env1.descriptor_vals.dot(env2.descriptor_vals); - double d1 = env1.descriptor_norm; - double d2 = env2.descriptor_norm; + double dot = env1.descriptor_vals[descriptor_index] + .dot(env2.descriptor_vals[descriptor_index]); + double d1 = env1.descriptor_norm[descriptor_index]; + double d2 = env2.descriptor_norm[descriptor_index]; return sig2 * pow(dot / (d1 * d2), power); } @@ -527,7 +540,7 @@ Eigen::VectorXd DotProductKernel // Account for edge case where d1 is zero. double empty_thresh = 1e-8; - double d1 = env1.descriptor_norm; + double d1 = env1.descriptor_norm[descriptor_index]; if (d1 < empty_thresh) return kern_vec; double en_kern = 0; @@ -546,26 +559,29 @@ Eigen::VectorXd DotProductKernel if (env1.central_species != env_curr.central_species) continue; // Check that d2 is nonzero. - d2 = env_curr.descriptor_norm; + d2 = env_curr.descriptor_norm[descriptor_index]; if (d2 < empty_thresh) continue; d2_cubed = d2 * d2 * d2; // Energy kernel - dot_val = env1.descriptor_vals.dot(env_curr.descriptor_vals); + dot_val = env1.descriptor_vals[descriptor_index] + .dot(env_curr.descriptor_vals[descriptor_index]); norm_dot = dot_val / (d1 * d2); en_kern += pow(norm_dot, power); // Force kernel - force_dot = env_curr.descriptor_force_dervs * env1.descriptor_vals; + force_dot = env_curr.descriptor_force_dervs[descriptor_index] * + env1.descriptor_vals[descriptor_index]; f1 = (force_dot / (d1 * d2)) - - (dot_val * env_curr.force_dot / (d2_cubed * d1)); + (dot_val * env_curr.force_dot[descriptor_index] / (d2_cubed * d1)); dval = power * pow(norm_dot, power - 1); force_kern += dval * f1; // Stress kernel - stress_dot = env_curr.descriptor_stress_dervs * env1.descriptor_vals; + stress_dot = env_curr.descriptor_stress_dervs[descriptor_index] * + env1.descriptor_vals[descriptor_index]; s1 = (stress_dot / (d1 * d2)) - - (dot_val * env_curr.stress_dot /(d2_cubed * d1)); + (dot_val * env_curr.stress_dot[descriptor_index] /(d2_cubed * d1)); stress_kern += dval * s1; } diff --git a/src/kernels.h b/src/kernels.h index caf0bb506..f52ac62e5 100644 --- a/src/kernels.h +++ b/src/kernels.h @@ -27,10 +27,11 @@ class Kernel{ class DotProductKernel : public Kernel{ public: double sigma, sig2, power; + int descriptor_index; DotProductKernel(); - DotProductKernel(double sigma, double power); + DotProductKernel(double sigma, double power, int descriptor_index); double env_env(const LocalEnvironment & env1, const LocalEnvironment & env2); diff --git a/src/local_environment.cpp b/src/local_environment.cpp index 9fddd49b0..1630a9a7e 100644 --- a/src/local_environment.cpp +++ b/src/local_environment.cpp @@ -2,6 +2,8 @@ #include #include +// TODO: carefully review and edit changes + LocalEnvironment :: LocalEnvironment(){} LocalEnvironment :: LocalEnvironment(const Structure & structure, int atom, @@ -36,32 +38,38 @@ LocalEnvironment :: LocalEnvironment(const Structure & structure, int atom, this->zrel = zrel; } +// n-body LocalEnvironment :: LocalEnvironment(const Structure & structure, int atom, - double cutoff, std::vector nested_cutoffs) + double cutoff, std::vector n_body_cutoffs) : LocalEnvironment(structure, atom, cutoff){ - this->nested_cutoffs = nested_cutoffs; - compute_nested_environment(); + this->n_body_cutoffs = n_body_cutoffs; + compute_indices(); } -// TODO: base "compute descriptor" on many body cutoff. +// many-body LocalEnvironment :: LocalEnvironment(const Structure & structure, int atom, - double cutoff, DescriptorCalculator * descriptor_calculator) + double cutoff, std::vector many_body_cutoffs, + std::vector descriptor_calculators) : LocalEnvironment(structure, atom, cutoff){ - this->descriptor_calculator = descriptor_calculator; - compute_descriptor(); + this->many_body_cutoffs = many_body_cutoffs; + this->descriptor_calculators = descriptor_calculators; + compute_descriptors(); } +// n-body + many-body LocalEnvironment :: LocalEnvironment(const Structure & structure, int atom, - double cutoff, std::vector nested_cutoffs, - DescriptorCalculator * descriptor_calculator) + double cutoff, std::vector n_body_cutoffs, + std::vector many_body_cutoffs, + std::vector descriptor_calculators) : LocalEnvironment(structure, atom, cutoff){ - this->nested_cutoffs = nested_cutoffs; - this->descriptor_calculator = descriptor_calculator; - compute_nested_environment(); - compute_descriptor(); + this->n_body_cutoffs = n_body_cutoffs; + this->many_body_cutoffs = many_body_cutoffs; + this->descriptor_calculators = descriptor_calculators; + compute_indices(); + compute_descriptors(); } void LocalEnvironment :: compute_environment( @@ -174,25 +182,48 @@ void LocalEnvironment :: compute_environment( delete [] dists; delete [] xvals; delete [] yvals; delete [] zvals; } -void LocalEnvironment :: compute_nested_environment(){ - - double two_body_cutoff = nested_cutoffs[0]; - double three_body_cutoff = nested_cutoffs[1]; - double many_body_cutoff = nested_cutoffs[2]; +void LocalEnvironment :: compute_indices(){ int no_atoms = rs.size(); std::vector three_body_inds; + // Initialize a list of lists storing atom indices. + int n_cutoffs = n_body_cutoffs.size(); + std::vector empty; + for (int i = 0; i < n_cutoffs; i ++){ + n_body_indices.push_back(empty); + } + + int n_mb_cutoffs = many_body_cutoffs.size(); + for (int i = 0; i < n_mb_cutoffs; i ++){ + many_body_indices.push_back(empty); + } + // Store indices of atoms inside the 2-, 3-, and many-body cutoff spheres. + double current_cutoff; for (int i = 0; i < no_atoms; i ++){ double r_curr = rs[i]; - if (r_curr <= two_body_cutoff) two_body_indices.push_back(i); - if (r_curr <= three_body_cutoff) three_body_inds.push_back(i); - if (r_curr <= many_body_cutoff) many_body_indices.push_back(i); + + // Store n-body index + for (int j = 0; j < n_cutoffs; j ++){ + current_cutoff = n_body_cutoffs[j]; + if (r_curr < current_cutoff){ + n_body_indices[j].push_back(i); + } + } + + // Store body index + for (int j = 0; j < n_mb_cutoffs; j ++){ + current_cutoff = many_body_cutoffs[j]; + if (r_curr < current_cutoff){ + many_body_indices[j].push_back(i); + } + } } // Store triplets if the 3-body cutoff is nonzero. - if (three_body_cutoff > 0){ + if (n_cutoffs > 1){ + double three_body_cutoff = n_body_cutoffs[1]; double cross_bond_dist, x1, y1, z1, x2, y2, z2, x_diff, y_diff, z_diff; int ind1, ind2; std::vector triplet = std::vector {0, 0}; @@ -219,13 +250,19 @@ void LocalEnvironment :: compute_nested_environment(){ } } -void LocalEnvironment :: compute_descriptor(){ - descriptor_calculator->compute(*this); - descriptor_vals = descriptor_calculator->descriptor_vals; - descriptor_force_dervs = descriptor_calculator->descriptor_force_dervs; - descriptor_stress_dervs = descriptor_calculator->descriptor_stress_dervs; - - descriptor_norm = sqrt(descriptor_vals.dot(descriptor_vals)); - force_dot = descriptor_force_dervs * descriptor_vals; - stress_dot = descriptor_stress_dervs * descriptor_vals; +void LocalEnvironment :: compute_descriptors(){ + int n_calculators = descriptor_calculators.size(); + for (int i = 0; i < n_calculators; i ++){ + descriptor_calculators[i]->compute(*this); + descriptor_vals.push_back(descriptor_calculators[i]->descriptor_vals); + descriptor_force_dervs.push_back( + descriptor_calculators[i]->descriptor_force_dervs); + descriptor_stress_dervs.push_back( + descriptor_calculators[i]->descriptor_stress_dervs); + + descriptor_norm.push_back(sqrt( + descriptor_vals[i].dot(descriptor_vals[i]))); + force_dot.push_back(descriptor_force_dervs[i] * descriptor_vals[i]); + stress_dot.push_back(descriptor_stress_dervs[i] * descriptor_vals[i]); + } } diff --git a/src/local_environment.h b/src/local_environment.h index acb48517f..62643a73d 100644 --- a/src/local_environment.h +++ b/src/local_environment.h @@ -8,35 +8,47 @@ class LocalEnvironment{ public: std::vector environment_indices, environment_species, - neighbor_list, two_body_indices, many_body_indices; - std::vector> three_body_indices; + neighbor_list; int central_index, central_species, noa, sweep; - std::vector rs, xs, ys, zs, xrel, yrel, zrel, - nested_cutoffs, cross_bond_dists; - double cutoff, structure_volume, descriptor_norm; + std::vector rs, xs, ys, zs, xrel, yrel, zrel; + double cutoff, structure_volume; + + // Store cutoffs for each kernel and indices of atoms inside each cutoff sphere. + std::vector n_body_cutoffs, many_body_cutoffs; + std::vector> n_body_indices, many_body_indices; + + // Triplet indices and cross bond distances are stored only if the 3-body kernel is used + std::vector> three_body_indices; + std::vector cross_bond_dists; - DescriptorCalculator * descriptor_calculator; - Eigen::VectorXd descriptor_vals; - Eigen::MatrixXd descriptor_force_dervs, descriptor_stress_dervs, - force_dot, stress_dot; + // Store descriptor calculators for each many body cutoff. + std::vector descriptor_calculators; + std::vector descriptor_vals; + std::vector descriptor_force_dervs, + descriptor_stress_dervs, force_dot, stress_dot; + std::vector descriptor_norm; LocalEnvironment(); LocalEnvironment(const Structure & structure, int atom, double cutoff); - // If nested cutoffs are given, compute 2-, 3-, and many-body indices. + // n-body LocalEnvironment(const Structure & structure, int atom, - double cutoff, std::vector nested_cutoffs); + double cutoff, std::vector n_body_cutoffs); - // If a descriptor calculator is given as well, compute and store the - // many-body descriptor. + // many-body LocalEnvironment(const Structure & structure, int atom, double cutoff, - DescriptorCalculator * descriptor_calculator); + std::vector many_body_cutoffs, + std::vector + descriptor_calculators); + // n-body + many-body LocalEnvironment(const Structure & structure, int atom, double cutoff, - std::vector nested_cutoffs, - DescriptorCalculator * descriptor_calculator); + std::vector n_body_cutoffs, + std::vector many_body_cutoffs, + std::vector + descriptor_calculator); static void compute_environment(const Structure & structure, int noa, int atom, @@ -52,8 +64,8 @@ class LocalEnvironment{ std::vector & yrel, std::vector & zrel); - void compute_nested_environment(); - void compute_descriptor(); + void compute_indices(); + void compute_descriptors(); }; #endif diff --git a/src/structure.cpp b/src/structure.cpp index 7fd118c47..5c154e83a 100644 --- a/src/structure.cpp +++ b/src/structure.cpp @@ -74,6 +74,10 @@ double Structure :: get_max_cutoff(){ return max_cutoff; } +// ---------------------------------------------------------------------------- +// structure descriptor +// ---------------------------------------------------------------------------- + StructureDescriptor :: StructureDescriptor(){} StructureDescriptor :: StructureDescriptor(const Eigen::MatrixXd & cell, @@ -85,37 +89,44 @@ StructureDescriptor :: StructureDescriptor(const Eigen::MatrixXd & cell, this->compute_environments(); } +// n-body StructureDescriptor :: StructureDescriptor(const Eigen::MatrixXd & cell, const std::vector & species, const Eigen::MatrixXd & positions, - double cutoff, std::vector nested_cutoffs) + double cutoff, std::vector n_body_cutoffs) : Structure(cell, species, positions){ this->cutoff = cutoff; - this->nested_cutoffs = nested_cutoffs; + this->n_body_cutoffs = n_body_cutoffs; this->compute_nested_environments(); } +// many-body StructureDescriptor :: StructureDescriptor(const Eigen::MatrixXd & cell, const std::vector & species, const Eigen::MatrixXd & positions, double cutoff, - DescriptorCalculator * descriptor_calculator) + std::vector many_body_cutoffs, + std::vector + descriptor_calculators) : Structure(cell, species, positions){ - this->descriptor_calculator = descriptor_calculator; + this->descriptor_calculators = descriptor_calculators; this->cutoff = cutoff; + this->many_body_cutoffs = many_body_cutoffs; this->compute_descriptors(); } StructureDescriptor :: StructureDescriptor(const Eigen::MatrixXd & cell, const std::vector & species, const Eigen::MatrixXd & positions, - double cutoff, std::vector nested_cutoffs, - DescriptorCalculator * descriptor_calculator) + double cutoff, std::vector n_body_cutoffs, + std::vector many_body_cutoffs, + std::vector descriptor_calculators) : Structure(cell, species, positions){ - this->descriptor_calculator = descriptor_calculator; + this->descriptor_calculators = descriptor_calculators; this->cutoff = cutoff; - this->nested_cutoffs = nested_cutoffs; + this->n_body_cutoffs = n_body_cutoffs; + this->many_body_cutoffs = many_body_cutoffs; this->nested_descriptors(); } @@ -134,7 +145,7 @@ void StructureDescriptor :: compute_nested_environments(){ LocalEnvironment env; for (int i = 0; i < noa; i ++){ - env = LocalEnvironment(*this, i, cutoff, nested_cutoffs); + env = LocalEnvironment(*this, i, cutoff, n_body_cutoffs); local_environments.push_back(env); } } @@ -144,7 +155,8 @@ void StructureDescriptor :: compute_descriptors(){ LocalEnvironment env; for (int i = 0; i < noa; i ++){ - env = LocalEnvironment(*this, i, cutoff, descriptor_calculator); + env = LocalEnvironment(*this, i, cutoff, many_body_cutoffs, + descriptor_calculators); local_environments.push_back(env); } } @@ -154,8 +166,8 @@ void StructureDescriptor :: nested_descriptors(){ LocalEnvironment env; for (int i = 0; i < noa; i ++){ - env = LocalEnvironment(*this, i, cutoff, nested_cutoffs, - descriptor_calculator); + env = LocalEnvironment(*this, i, cutoff, n_body_cutoffs, + many_body_cutoffs, descriptor_calculators); local_environments.push_back(env); } } diff --git a/src/structure.h b/src/structure.h index 0ea448ab9..2037111b5 100644 --- a/src/structure.h +++ b/src/structure.h @@ -30,10 +30,11 @@ class Structure{ // Structure descriptor. Stores the atomic environments in a structure. class StructureDescriptor : public Structure{ public: - DescriptorCalculator * descriptor_calculator; + std::vector descriptor_calculators; std::vector local_environments; double cutoff; - std::vector nested_cutoffs; + std::vector n_body_cutoffs; + std::vector many_body_cutoffs; // Make structure labels empty by default. Eigen::VectorXd energy; @@ -47,22 +48,28 @@ class StructureDescriptor : public Structure{ const Eigen::MatrixXd & positions, double cutoff); + // n-body StructureDescriptor(const Eigen::MatrixXd & cell, const std::vector & species, const Eigen::MatrixXd & positions, - double cutoff, std::vector nested_cutoffs); + double cutoff, std::vector n_body_cutoffs); + // many-body StructureDescriptor(const Eigen::MatrixXd & cell, const std::vector & species, - const Eigen::MatrixXd & positions, - double cutoff, - DescriptorCalculator * descriptor_calculator); + const Eigen::MatrixXd & positions, double cutoff, + std::vector many_body_cutoffs, + std::vector + descriptor_calculators); + // n-body + many-body StructureDescriptor(const Eigen::MatrixXd & cell, const std::vector & species, const Eigen::MatrixXd & positions, - double cutoff, std::vector nested_cutoffs, - DescriptorCalculator * descriptor_calculator); + double cutoff, std::vector n_body_cutoffs, + std::vector many_body_cutoffs, + std::vector + descriptor_calculators); void compute_environments(); void compute_nested_environments(); diff --git a/tests/test_descriptor.cpp b/tests/test_descriptor.cpp index 2dc1af55f..ec1f8b6cf 100644 --- a/tests/test_descriptor.cpp +++ b/tests/test_descriptor.cpp @@ -30,6 +30,7 @@ class DescriptorTest : public ::testing::Test{ std::vector descriptor_settings {5, 10, 10}; double rcut = 3; + std::vector many_body_cutoffs {rcut}; // Choose arbitrary rotation angles. double xrot = 1.28; @@ -64,6 +65,7 @@ class DescriptorTest : public ::testing::Test{ // Set up descriptor calculator. std::string radial_string = "chebyshev"; std::string cutoff_string = "quadratic"; + int descriptor_index = 0; DescriptorTest(){ // Define rotation matrices. @@ -96,6 +98,7 @@ class DescriptorTest : public ::testing::Test{ env1 = LocalEnvironment(struc1, 0, rcut); env2 = LocalEnvironment(struc2, 0, rcut); + env1.many_body_cutoffs = env2.many_body_cutoffs = many_body_cutoffs; single_bond_vals = Eigen::VectorXd::Zero(no_descriptors); force_dervs = Eigen::MatrixXd::Zero(noa * 3, no_descriptors); @@ -103,11 +106,13 @@ class DescriptorTest : public ::testing::Test{ // Create B1 and B2 descriptor calculators. desc1 = B1_Calculator(radial_string, cutoff_string, - radial_hyps, cutoff_hyps, descriptor_settings); + radial_hyps, cutoff_hyps, descriptor_settings, + descriptor_index); desc2 = desc3 = desc1; desc4 = B2_Calculator(radial_string, cutoff_string, - radial_hyps, cutoff_hyps, descriptor_settings); + radial_hyps, cutoff_hyps, descriptor_settings, + descriptor_index); desc5 = desc6 = desc4; } @@ -171,6 +176,7 @@ TEST_F(DescriptorTest, CentTest){ positions_3(0, m) += delta; struc3 = Structure(cell, species, positions_3); env3 = LocalEnvironment(struc3, 0, rcut); + env3.many_body_cutoffs = many_body_cutoffs; desc3.compute(env3); // Check derivatives. @@ -194,6 +200,7 @@ TEST_F(DescriptorTest, CentTest){ positions_3(0, m) += delta; struc3 = Structure(cell, species, positions_3); env3 = LocalEnvironment(struc3, 0, rcut); + env3.many_body_cutoffs = many_body_cutoffs; desc6.compute(env3); // Check derivatives. @@ -221,6 +228,7 @@ TEST_F(DescriptorTest, EnvTest){ positions_3(p, m) += delta; struc3 = Structure(cell, species, positions_3); env3 = LocalEnvironment(struc3, 0, rcut); + env3.many_body_cutoffs = many_body_cutoffs; desc3.compute(env3); // Check derivatives. @@ -247,6 +255,7 @@ TEST_F(DescriptorTest, EnvTest){ positions_3(p, m) += delta; struc3 = Structure(cell, species, positions_3); env3 = LocalEnvironment(struc3, 0, rcut); + env3.many_body_cutoffs = many_body_cutoffs; desc6.compute(env3); // Check derivatives. @@ -286,6 +295,7 @@ TEST_F(DescriptorTest, StressTest){ struc2 = Structure(cell_2, species, positions_2); env2 = LocalEnvironment(struc2, 0, rcut); + env2.many_body_cutoffs = many_body_cutoffs; desc2.compute(env2); // Check stress derivatives. @@ -324,6 +334,7 @@ TEST_F(DescriptorTest, StressTest){ struc2 = Structure(cell_2, species, positions_2); env2 = LocalEnvironment(struc2, 0, rcut); + env2.many_body_cutoffs = many_body_cutoffs; desc5.compute(env2); // Check stress derivatives. diff --git a/tests/test_kernels.cpp b/tests/test_kernels.cpp index 0c931fc0a..8e74984f4 100644 --- a/tests/test_kernels.cpp +++ b/tests/test_kernels.cpp @@ -27,9 +27,12 @@ class KernelTest : public ::testing::Test{ std::vector radial_hyps {0, 5}; std::vector cutoff_hyps; std::vector descriptor_settings {2, 5, 5}; - double cutoff = 5; - std::vector nested_cutoffs {5, 5}; + double cutoff = 6; + std::vector nested_cutoffs {5, 4}; + std::vector many_body_cutoffs {3}; B2_Calculator desc1; + std::vector descriptor_calculators; + int descriptor_index = 0; // kernel double signal_variance = 2; @@ -58,15 +61,17 @@ class KernelTest : public ::testing::Test{ 0.35585787, -0.87190223, 0.06770428; desc1 = B2_Calculator(radial_string, cutoff_string, - radial_hyps, cutoff_hyps, descriptor_settings); + radial_hyps, cutoff_hyps, descriptor_settings, 0); + descriptor_calculators.push_back(&desc1); + test_struc = StructureDescriptor(cell, species, positions, cutoff, - nested_cutoffs, &desc1); + nested_cutoffs, many_body_cutoffs, descriptor_calculators); test_struc_2 = StructureDescriptor(cell, species, positions_2, cutoff, - nested_cutoffs, &desc1); + nested_cutoffs, many_body_cutoffs, descriptor_calculators); test_env = test_struc_2.local_environments[0]; - kernel = DotProductKernel(signal_variance, power); + kernel = DotProductKernel(signal_variance, power, descriptor_index); two_body_kernel = TwoBodyKernel(signal_variance, length_scale, cutoff_string, cutoff_hyps); three_body_kernel = ThreeBodyKernel(signal_variance, length_scale, @@ -96,8 +101,9 @@ TEST_F(KernelTest, ForceTest){ for (int m = 0; m < 3; m ++){ positions_3 = positions; positions_3(p, m) += delta; - test_struc_3 = StructureDescriptor(cell, species, positions_3, - cutoff, nested_cutoffs, &desc1); + test_struc_3 = StructureDescriptor(cell, species, positions_3, + cutoff, nested_cutoffs, many_body_cutoffs, + descriptor_calculators); kern_pert = kernel.env_struc(test_env, test_struc_3); fin_val = -(kern_pert(0) - kern_vec(0)) / delta; exact_val = kern_vec(1 + 3 * p + m); @@ -120,7 +126,8 @@ TEST_F(KernelTest, TwoBodyForceTest){ positions_3 = positions; positions_3(p, m) += delta; test_struc_3 = StructureDescriptor(cell, species, positions_3, - cutoff, nested_cutoffs, &desc1); + cutoff, nested_cutoffs, many_body_cutoffs, + descriptor_calculators); kern_pert = two_body_kernel.env_struc(test_env, test_struc_3); fin_val = -(kern_pert(0) - kern_vec_2(0)) / delta; exact_val = kern_vec_2(1 + 3 * p + m); @@ -142,7 +149,8 @@ TEST_F(KernelTest, ThreeBodyForceTest){ positions_3 = positions; positions_3(p, m) += delta; test_struc_3 = StructureDescriptor(cell, species, positions_3, - cutoff, nested_cutoffs, &desc1); + cutoff, nested_cutoffs, many_body_cutoffs, + descriptor_calculators); kern_pert = three_body_kernel.env_struc(test_env, test_struc_3); fin_val = -(kern_pert(0) - kern_vec_3(0)) / delta; exact_val = kern_vec_3(1 + 3 * p + m); @@ -174,7 +182,8 @@ TEST_F(KernelTest, StressTest){ } test_struc_3 = StructureDescriptor(cell_2, species, positions_3, - cutoff, nested_cutoffs, &desc1); + cutoff, nested_cutoffs, many_body_cutoffs, + descriptor_calculators); kern_pert = kernel.env_struc(test_env, test_struc_3); fin_val = -(kern_pert(0) - kern_vec(0)) / delta; @@ -211,7 +220,8 @@ TEST_F(KernelTest, TwoBodyStressTest){ } test_struc_3 = StructureDescriptor(cell_2, species, positions_3, - cutoff, nested_cutoffs, &desc1); + cutoff, nested_cutoffs, many_body_cutoffs, + descriptor_calculators); kern_pert = two_body_kernel.env_struc(test_env, test_struc_3); fin_val = -(kern_pert(0) - kern_vec_2(0)) / delta; @@ -248,7 +258,8 @@ TEST_F(KernelTest, ThreeBodyStressTest){ } test_struc_3 = StructureDescriptor(cell_2, species, positions_3, - cutoff, nested_cutoffs, &desc1); + cutoff, nested_cutoffs, many_body_cutoffs, + descriptor_calculators); kern_pert = three_body_kernel.env_struc(test_env, test_struc_3); fin_val = -(kern_pert(0) - kern_vec_3(0)) / delta; diff --git a/tests/test_local_environment.cpp b/tests/test_local_environment.cpp index b7245459d..dbcdc1a91 100644 --- a/tests/test_local_environment.cpp +++ b/tests/test_local_environment.cpp @@ -13,6 +13,7 @@ class EnvironmentTest : public :: testing :: Test{ std::vector species {0, 1, 2, 3, 4}; Eigen::MatrixXd positions{5, 3}; B2_Calculator desc1; + std::vector descriptor_calculators; StructureDescriptor test_struc; int atom; LocalEnvironment test_env; @@ -22,8 +23,10 @@ class EnvironmentTest : public :: testing :: Test{ std::vector radial_hyps {0, 5}; std::vector cutoff_hyps; std::vector descriptor_settings {5, 5, 5}; + int descriptor_index = 0; std::vector nested_cutoffs {3, 3, 3}; double cutoff = 3; + std::vector many_body_cutoffs {cutoff}; EnvironmentTest(){ cell << 1.3, 0.5, 0.8, @@ -37,13 +40,14 @@ class EnvironmentTest : public :: testing :: Test{ 3.2, 1.1, 3.3; desc1 = B2_Calculator(radial_string, cutoff_string, - radial_hyps, cutoff_hyps, descriptor_settings); + radial_hyps, cutoff_hyps, descriptor_settings, descriptor_index); + descriptor_calculators.push_back(&desc1); test_struc = StructureDescriptor(cell, species, positions, cutoff, - &desc1); + many_body_cutoffs, descriptor_calculators); atom = 0; test_env = LocalEnvironment(test_struc, atom, cutoff, - nested_cutoffs, &desc1); + nested_cutoffs, descriptor_calculators); } }; @@ -76,14 +80,14 @@ TEST_F(EnvironmentTest, DotTest){ // Calculate the descriptor norm the old fashioned way. double norm_val = 0; double val_curr; - int no_desc = test_env.descriptor_vals.rows(); + int no_desc = test_env.descriptor_vals[0].rows(); for (int i = 0; i < no_desc; i++){ - val_curr = test_env.descriptor_vals(i); + val_curr = test_env.descriptor_vals[0](i); norm_val += val_curr * val_curr; } norm_val = sqrt(norm_val); - EXPECT_NEAR(norm_val, test_env.descriptor_norm, THRESHOLD); + EXPECT_NEAR(norm_val, test_env.descriptor_norm[0], THRESHOLD); } // TEST_F(EnvironmentTest, NestedTest){ diff --git a/tests/test_sparse_gp.cpp b/tests/test_sparse_gp.cpp index ddfd36188..cc45e4ec9 100644 --- a/tests/test_sparse_gp.cpp +++ b/tests/test_sparse_gp.cpp @@ -29,6 +29,7 @@ class SparseTest : public ::testing::Test{ double signal_variance = 1; double length_scale = 1; double power = 1; + int descriptor_index = 0; DotProductKernel kernel; TwoBodyKernel two_body_kernel; ThreeBodyKernel three_body_kernel; @@ -48,7 +49,7 @@ class SparseTest : public ::testing::Test{ stresses = Eigen::VectorXd::Random(6); desc1 = B2_Calculator(radial_string, cutoff_string, - radial_hyps, cutoff_hyps, descriptor_settings); + radial_hyps, cutoff_hyps, descriptor_settings, 0); test_struc = StructureDescriptor(cell, species, positions, cutoff, nested_cutoffs); @@ -60,7 +61,8 @@ class SparseTest : public ::testing::Test{ cutoff_string, cutoff_hyps); three_body_kernel = ThreeBodyKernel(signal_variance, length_scale, cutoff_string, cutoff_hyps); - many_body_kernel = DotProductKernel(signal_variance, power); + many_body_kernel = DotProductKernel(signal_variance, power, + descriptor_index); // kernels = // std::vector {&two_body_kernel, &three_body_kernel, @@ -102,50 +104,50 @@ TEST_F(SparseTest, UpdateK){ EXPECT_EQ(kern_val, sparse_gp.Kuu(0,0)); } -// TEST_F(SparseTest, Predict){ -// // Eigen::initParallel(); - -// double sigma_e = 0.1; -// double sigma_f = 0.01; -// double sigma_s = 1; - -// SparseGP sparse_gp = SparseGP(kernels, sigma_e, sigma_f, sigma_s); -// LocalEnvironment env1 = test_struc.local_environments[0]; -// LocalEnvironment env2 = test_struc.local_environments[1]; -// sparse_gp.add_sparse_environment(env1); -// sparse_gp.add_sparse_environment(env2); - -// test_struc.stresses = Eigen::VectorXd {}; - -// sparse_gp.add_training_structure(test_struc); -// sparse_gp.update_alpha(); -// Eigen::VectorXd pred_vals; - -// // predict in parallel -// auto t1 = std::chrono::high_resolution_clock::now(); -// for (int i = 0; i < 1000; i ++){ -// pred_vals = sparse_gp.predict(test_struc); -// } -// auto t2 = std::chrono::high_resolution_clock::now(); -// auto tot_time = -// std::chrono::duration_cast(t2-t1).count(); -// // std::cout << tot_time << std::endl; - -// // // predict in serial -// // t1 = std::chrono::high_resolution_clock::now(); -// // for (int i = 0; i < 1000; i ++){ -// // Eigen::VectorXd pred_vals = sparse_gp.predict_serial(test_struc); -// // } -// // t2 = std::chrono::high_resolution_clock::now(); -// // tot_time = -// // std::chrono::duration_cast(t2-t1).count(); -// // std::cout << tot_time << std::endl; - -// // std::cout << "predicted values:" << std::endl; -// // std::cout << pred_vals << std::endl; - -// // std::cout << sparse_gp.Kuu << std::endl; -// // std::cout << sparse_gp.Kuf << std::endl; -// // std::cout << sparse_gp.Sigma << std::endl; -// // std::cout << sparse_gp.alpha << std::endl; -// } \ No newline at end of file +TEST_F(SparseTest, Predict){ + // Eigen::initParallel(); + + double sigma_e = 0.1; + double sigma_f = 0.01; + double sigma_s = 1; + + SparseGP sparse_gp = SparseGP(kernels, sigma_e, sigma_f, sigma_s); + LocalEnvironment env1 = test_struc.local_environments[0]; + LocalEnvironment env2 = test_struc.local_environments[1]; + sparse_gp.add_sparse_environment(env1); + sparse_gp.add_sparse_environment(env2); + + test_struc.stresses = Eigen::VectorXd {}; + + sparse_gp.add_training_structure(test_struc); + sparse_gp.update_alpha(); + Eigen::VectorXd pred_vals; + + // predict in parallel + auto t1 = std::chrono::high_resolution_clock::now(); + for (int i = 0; i < 1000; i ++){ + pred_vals = sparse_gp.predict(test_struc); + } + auto t2 = std::chrono::high_resolution_clock::now(); + auto tot_time = + std::chrono::duration_cast(t2-t1).count(); + std::cout << tot_time << std::endl; + + // // predict in serial + // t1 = std::chrono::high_resolution_clock::now(); + // for (int i = 0; i < 1000; i ++){ + // Eigen::VectorXd pred_vals = sparse_gp.predict_serial(test_struc); + // } + // t2 = std::chrono::high_resolution_clock::now(); + // tot_time = + // std::chrono::duration_cast(t2-t1).count(); + // std::cout << tot_time << std::endl; + + // std::cout << "predicted values:" << std::endl; + // std::cout << pred_vals << std::endl; + + // std::cout << sparse_gp.Kuu << std::endl; + // std::cout << sparse_gp.Kuf << std::endl; + // std::cout << sparse_gp.Sigma << std::endl; + // std::cout << sparse_gp.alpha << std::endl; +} \ No newline at end of file diff --git a/tests/test_structure.cpp b/tests/test_structure.cpp index 6f6ed6c3a..f4973761c 100644 --- a/tests/test_structure.cpp +++ b/tests/test_structure.cpp @@ -12,6 +12,7 @@ class StructureTest : public ::testing::Test{ std::vector species {0, 1, 2, 3, 4}; Eigen::MatrixXd positions{5, 3}; B2_Calculator desc1; + std::vector descriptor_calculators; StructureDescriptor test_struc; std::string radial_string = "chebyshev"; @@ -19,7 +20,9 @@ class StructureTest : public ::testing::Test{ std::vector radial_hyps {0, 5}; std::vector cutoff_hyps; std::vector descriptor_settings {5, 5, 5}; + int descriptor_index = 0; double cutoff = 3; + std::vector many_body_cutoffs {cutoff}; StructureTest(){ cell << 4.0, 0.5, 0.8, @@ -33,9 +36,10 @@ class StructureTest : public ::testing::Test{ 3.2, 1.1, 3.3; desc1 = B2_Calculator(radial_string, cutoff_string, - radial_hyps, cutoff_hyps, descriptor_settings); + radial_hyps, cutoff_hyps, descriptor_settings, descriptor_index); + descriptor_calculators.push_back(&desc1); test_struc = StructureDescriptor(cell, species, positions, cutoff, - &desc1); + many_body_cutoffs, descriptor_calculators); } }; @@ -66,16 +70,17 @@ TEST_F(StructureTest, StructureDescriptor){ LocalEnvironment env; for (int i = 0; i < test_struc.species.size(); i ++){ env = LocalEnvironment(test_struc, i, cutoff); + env.many_body_cutoffs = many_body_cutoffs; desc1.compute(env); for (int j = 0; j < desc1.descriptor_vals.size(); j ++){ EXPECT_EQ(desc1.descriptor_vals(j), test_struc.local_environments[i] - .descriptor_vals(j)); + .descriptor_vals[0](j)); for (int k = 0; k < test_struc.species.size(); k ++){ EXPECT_EQ(desc1.descriptor_force_dervs(k, j), test_struc.local_environments[i] - .descriptor_force_dervs(k, j)); + .descriptor_force_dervs[0](k, j)); } } } From ce85df0345c79ddd849633fd81b40fd3f2d730dd Mon Sep 17 00:00:00 2001 From: Jonathan Vandermause Date: Sat, 7 Mar 2020 16:46:19 -0500 Subject: [PATCH 18/24] use n body cutoff in 3-body kernel --- src/kernels.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/kernels.cpp b/src/kernels.cpp index 699cebf41..b31398218 100644 --- a/src/kernels.cpp +++ b/src/kernels.cpp @@ -205,8 +205,8 @@ double ThreeBodyKernel :: env_env(const LocalEnvironment & env1, r22, r23, r31, r32, r33, p1, p2, p3, p4, p5, p6; int i1, i2, j1, j2, ei1, ei2, ej1, ej2; - double cut1 = env1.cutoff; - double cut2 = env2.cutoff; + double cut1 = env1.n_body_cutoffs[1]; + double cut2 = env2.n_body_cutoffs[1]; double rcut_vals_i1[2], rcut_vals_i2[2], rcut_vals_i3[2], rcut_vals_j1[2], rcut_vals_j2[2], rcut_vals_j3[2]; int c1 = env1.central_species; @@ -324,8 +324,8 @@ Eigen::VectorXd ThreeBodyKernel :: env_struc(const LocalEnvironment & env1, LocalEnvironment env2; double vol_inv = 1 / struc1.volume; - double cut1 = env1.cutoff; - double cut2 = struc1.cutoff; + double cut1 = env1.n_body_cutoffs[1]; + double cut2 = struc1.n_body_cutoffs[1]; double rcut_vals_i1[2], rcut_vals_i2[2], rcut_vals_i3[2], rcut_vals_j1[2], rcut_vals_j2[2], rcut_vals_j3[2]; int c1 = env1.central_species; From f7bfcc33abc445ea10c326c3f7e8fe28d6152fc0 Mon Sep 17 00:00:00 2001 From: Jonathan Vandermause Date: Sat, 7 Mar 2020 17:27:22 -0500 Subject: [PATCH 19/24] loop over mb indices only when performing single bond sums --- src/descriptor.cpp | 6 ++---- src/local_environment.cpp | 1 + src/single_bond.cpp | 10 ++++++---- src/single_bond.h | 2 +- tests/test_bond_env.cpp | 22 ++++++++++++++++------ tests/test_descriptor.cpp | 10 +++++++++- tests/test_structure.cpp | 4 ++-- 7 files changed, 37 insertions(+), 18 deletions(-) diff --git a/src/descriptor.cpp b/src/descriptor.cpp index b7b4c5366..2702b9c7f 100644 --- a/src/descriptor.cpp +++ b/src/descriptor.cpp @@ -125,8 +125,7 @@ void B1_Calculator :: compute(const LocalEnvironment & env){ // Compute single bond vector. single_bond_sum_env(single_bond_vals, single_bond_force_dervs, single_bond_stress_dervs, radial_pointer, - cutoff_pointer, env, - env.many_body_cutoffs[descriptor_index], N, + cutoff_pointer, env, descriptor_index, N, lmax, radial_hyps, cutoff_hyps); // Set B1 values. @@ -166,8 +165,7 @@ void B2_Calculator :: compute(const LocalEnvironment & env){ // Compute single bond vector. single_bond_sum_env(single_bond_vals, single_bond_force_dervs, single_bond_stress_dervs, radial_pointer, - cutoff_pointer, env, - env.many_body_cutoffs[descriptor_index], N, + cutoff_pointer, env, descriptor_index, N, lmax, radial_hyps, cutoff_hyps); // Initialize B2 vectors. diff --git a/src/local_environment.cpp b/src/local_environment.cpp index 1630a9a7e..fd6809456 100644 --- a/src/local_environment.cpp +++ b/src/local_environment.cpp @@ -55,6 +55,7 @@ LocalEnvironment :: LocalEnvironment(const Structure & structure, int atom, this->many_body_cutoffs = many_body_cutoffs; this->descriptor_calculators = descriptor_calculators; + compute_indices(); compute_descriptors(); } diff --git a/src/single_bond.cpp b/src/single_bond.cpp index 8e03b8071..3f9c45847 100644 --- a/src/single_bond.cpp +++ b/src/single_bond.cpp @@ -95,15 +95,17 @@ void single_bond_sum_env( Eigen::MatrixXd & force_dervs, Eigen::MatrixXd & stress_dervs, void (*basis_function)(double *, double *, double, int, vector), void (*cutoff_function)(double *, double, double, vector), - const LocalEnvironment & env, double rcut, int N, int lmax, + const LocalEnvironment & env, int descriptor_index, int N, int lmax, const vector & radial_hyps, const vector & cutoff_hyps){ - int noa = env.rs.size(); + int noa = env.many_body_indices[descriptor_index].size(); int cent_ind = env.central_index; double x, y, z, r; - int s, env_ind; + int s, env_ind, atom_index; + double rcut = env.many_body_cutoffs[descriptor_index]; - for (int atom_index = 0; atom_index < noa; atom_index ++){ + for (int n = 0; n < noa; n ++){ + atom_index = env.many_body_indices[descriptor_index][n]; x = env.xs[atom_index]; y = env.ys[atom_index]; z = env.zs[atom_index]; diff --git a/src/single_bond.h b/src/single_bond.h index c6243922e..6da343ee8 100644 --- a/src/single_bond.h +++ b/src/single_bond.h @@ -25,7 +25,7 @@ void single_bond_sum_env( void (*basis_function)(double *, double *, double, int, std::vector), void (*cutoff_function)(double *, double, double, std::vector), - const LocalEnvironment & env, double rcut, int N, int lmax, + const LocalEnvironment & env, int descriptor_index, int N, int lmax, const std::vector & radial_hyps, const std::vector & cutoff_hyps); diff --git a/tests/test_bond_env.cpp b/tests/test_bond_env.cpp index 0fda4cdcc..89255fdd8 100644 --- a/tests/test_bond_env.cpp +++ b/tests/test_bond_env.cpp @@ -4,6 +4,7 @@ #include "structure.h" #include "cutoffs.h" #include "radial.h" +#include "descriptor.h" #include #include #include @@ -23,6 +24,7 @@ class BondEnv : public ::testing::Test{ LocalEnvironment env1, env2; double rcut = 3; + std::vector many_body_cutoffs {rcut}; // Prepare cutoff. std::vector cutoff_hyps; @@ -62,6 +64,8 @@ class BondEnv : public ::testing::Test{ struc1 = Structure(cell, species, positions_1); env1 = LocalEnvironment(struc1, 0, rcut); + env1.many_body_cutoffs = many_body_cutoffs; + env1.compute_indices(); single_bond_vals =Eigen::VectorXd::Zero(no_descriptors); force_dervs = Eigen::MatrixXd::Zero(noa * 3, no_descriptors); @@ -72,7 +76,7 @@ class BondEnv : public ::testing::Test{ TEST_F(BondEnv, CentTest){ single_bond_sum_env(single_bond_vals, force_dervs, stress_dervs, - basis_function, cutoff_function, env1, rcut, N, lmax, + basis_function, cutoff_function, env1, 0, N, lmax, radial_hyps, cutoff_hyps); // Perturb the coordinates of the central atom. @@ -81,6 +85,8 @@ TEST_F(BondEnv, CentTest){ positions_2(0, m) += delta; struc2 = Structure(cell, species, positions_2); env2 = LocalEnvironment(struc2, 0, rcut); + env2.many_body_cutoffs = many_body_cutoffs; + env2.compute_indices(); // Initialize matrices. single_bond_vals_2 = Eigen::VectorXd::Zero(no_descriptors); @@ -88,7 +94,7 @@ TEST_F(BondEnv, CentTest){ stress_dervs_2 = Eigen::MatrixXd::Zero(6, no_descriptors); single_bond_sum_env(single_bond_vals_2, force_dervs_2, stress_dervs_2, - basis_function, cutoff_function, env2, rcut, N, lmax, + basis_function, cutoff_function, env2, 0, N, lmax, radial_hyps, cutoff_hyps); double finite_diff, exact, diff; @@ -107,7 +113,7 @@ TEST_F(BondEnv, CentTest){ TEST_F(BondEnv, EnvTest){ single_bond_sum_env(single_bond_vals, force_dervs, stress_dervs, - basis_function, cutoff_function, env1, rcut, N, lmax, + basis_function, cutoff_function, env1, 0, N, lmax, radial_hyps, cutoff_hyps); double finite_diff, exact, diff; @@ -120,6 +126,8 @@ TEST_F(BondEnv, EnvTest){ positions_2(p, m) += delta; struc2 = Structure(cell, species, positions_2); env2 = LocalEnvironment(struc2, 0, rcut); + env2.many_body_cutoffs = many_body_cutoffs; + env2.compute_indices(); // Initialize matrices. single_bond_vals_2 = Eigen::VectorXd::Zero(no_descriptors); @@ -127,7 +135,7 @@ TEST_F(BondEnv, EnvTest){ stress_dervs_2 = Eigen::MatrixXd::Zero(6, no_descriptors); single_bond_sum_env(single_bond_vals_2, force_dervs_2, - stress_dervs_2, basis_function, cutoff_function, env2, rcut, + stress_dervs_2, basis_function, cutoff_function, env2, 0, N, lmax, radial_hyps, cutoff_hyps); // Check derivatives. @@ -148,7 +156,7 @@ TEST_F(BondEnv, StressTest){ double tolerance = 1e-5; single_bond_sum_env(single_bond_vals, force_dervs, stress_dervs, - basis_function, cutoff_function, env1, rcut, N, lmax, + basis_function, cutoff_function, env1, 0, N, lmax, radial_hyps, cutoff_hyps); // Test all 6 independent strains (xx, xy, xz, yy, yz, zz). @@ -167,6 +175,8 @@ TEST_F(BondEnv, StressTest){ struc2 = Structure(cell_2, species, positions_2); env2 = LocalEnvironment(struc2, 0, rcut); + env2.many_body_cutoffs = many_body_cutoffs; + env2.compute_indices(); single_bond_vals_2 = Eigen::VectorXd::Zero(no_descriptors); force_dervs_2 = Eigen::MatrixXd::Zero(noa * 3, no_descriptors); @@ -174,7 +184,7 @@ TEST_F(BondEnv, StressTest){ // Calculate descriptors. single_bond_sum_env(single_bond_vals_2, force_dervs_2, - stress_dervs_2, basis_function, cutoff_function, env2, rcut, + stress_dervs_2, basis_function, cutoff_function, env2, 0, N, lmax, radial_hyps, cutoff_hyps); // Check stress derivatives. diff --git a/tests/test_descriptor.cpp b/tests/test_descriptor.cpp index ec1f8b6cf..0bfd80c0d 100644 --- a/tests/test_descriptor.cpp +++ b/tests/test_descriptor.cpp @@ -99,6 +99,8 @@ class DescriptorTest : public ::testing::Test{ env1 = LocalEnvironment(struc1, 0, rcut); env2 = LocalEnvironment(struc2, 0, rcut); env1.many_body_cutoffs = env2.many_body_cutoffs = many_body_cutoffs; + env1.compute_indices(); + env2.compute_indices(); single_bond_vals = Eigen::VectorXd::Zero(no_descriptors); force_dervs = Eigen::MatrixXd::Zero(noa * 3, no_descriptors); @@ -177,6 +179,7 @@ TEST_F(DescriptorTest, CentTest){ struc3 = Structure(cell, species, positions_3); env3 = LocalEnvironment(struc3, 0, rcut); env3.many_body_cutoffs = many_body_cutoffs; + env3.compute_indices(); desc3.compute(env3); // Check derivatives. @@ -201,6 +204,7 @@ TEST_F(DescriptorTest, CentTest){ struc3 = Structure(cell, species, positions_3); env3 = LocalEnvironment(struc3, 0, rcut); env3.many_body_cutoffs = many_body_cutoffs; + env3.compute_indices(); desc6.compute(env3); // Check derivatives. @@ -228,7 +232,8 @@ TEST_F(DescriptorTest, EnvTest){ positions_3(p, m) += delta; struc3 = Structure(cell, species, positions_3); env3 = LocalEnvironment(struc3, 0, rcut); - env3.many_body_cutoffs = many_body_cutoffs; + env3.many_body_cutoffs = many_body_cutoffs; + env3.compute_indices(); desc3.compute(env3); // Check derivatives. @@ -256,6 +261,7 @@ TEST_F(DescriptorTest, EnvTest){ struc3 = Structure(cell, species, positions_3); env3 = LocalEnvironment(struc3, 0, rcut); env3.many_body_cutoffs = many_body_cutoffs; + env3.compute_indices(); desc6.compute(env3); // Check derivatives. @@ -296,6 +302,7 @@ TEST_F(DescriptorTest, StressTest){ struc2 = Structure(cell_2, species, positions_2); env2 = LocalEnvironment(struc2, 0, rcut); env2.many_body_cutoffs = many_body_cutoffs; + env2.compute_indices(); desc2.compute(env2); // Check stress derivatives. @@ -335,6 +342,7 @@ TEST_F(DescriptorTest, StressTest){ struc2 = Structure(cell_2, species, positions_2); env2 = LocalEnvironment(struc2, 0, rcut); env2.many_body_cutoffs = many_body_cutoffs; + env2.compute_indices(); desc5.compute(env2); // Check stress derivatives. diff --git a/tests/test_structure.cpp b/tests/test_structure.cpp index f4973761c..b36f881fa 100644 --- a/tests/test_structure.cpp +++ b/tests/test_structure.cpp @@ -69,8 +69,8 @@ TEST_F(StructureTest, StructureDescriptor){ // Check that structure descriptors match environment descriptors. LocalEnvironment env; for (int i = 0; i < test_struc.species.size(); i ++){ - env = LocalEnvironment(test_struc, i, cutoff); - env.many_body_cutoffs = many_body_cutoffs; + env = LocalEnvironment(test_struc, i, cutoff, many_body_cutoffs, + descriptor_calculators); desc1.compute(env); for (int j = 0; j < desc1.descriptor_vals.size(); j ++){ From 3c05ecf4665b5a9fd99fa97d39ec07ce8cb12bbe Mon Sep 17 00:00:00 2001 From: Jonathan Vandermause Date: Sat, 7 Mar 2020 17:29:41 -0500 Subject: [PATCH 20/24] add descriptor todo --- src/descriptor.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/descriptor.cpp b/src/descriptor.cpp index 2702b9c7f..5b9ae5ae9 100644 --- a/src/descriptor.cpp +++ b/src/descriptor.cpp @@ -70,6 +70,7 @@ for (int n1 = 0; n1 < no_radial; n1 ++){ single_bond_vals(n2_l); // Store force derivatives. + // TODO: loop over many body indices, not entire neighbor list for (int atom_index = 0; atom_index < neigh_size; atom_index ++){ env_ind = env.neighbor_list[atom_index]; From e48ca45ae74b6c63bdf5557a645f58c1c5aa52ac Mon Sep 17 00:00:00 2001 From: Jonathan Vandermause Date: Sun, 8 Mar 2020 15:08:27 -0400 Subject: [PATCH 21/24] fix b2 binding --- src/ace_binding.cpp | 2 +- src/sparse_gp.h | 4 ++++ 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/src/ace_binding.cpp b/src/ace_binding.cpp index c93fde854..3a353f169 100644 --- a/src/ace_binding.cpp +++ b/src/ace_binding.cpp @@ -101,7 +101,7 @@ PYBIND11_MODULE(ace, m){ py::class_(m, "B2_Calculator") .def(py::init &, const std::vector &, - const std::vector &>()); + const std::vector &, int>()); // Kernel functions py::class_(m, "Kernel") diff --git a/src/sparse_gp.h b/src/sparse_gp.h index 79a100965..1fceb9c94 100644 --- a/src/sparse_gp.h +++ b/src/sparse_gp.h @@ -22,6 +22,8 @@ class SparseGP{ double sigma_e, sigma_f, sigma_s; + double log_marginal_likelihood, data_fit, complexity_penalty; + SparseGP(); SparseGP(std::vector, double sigma_e, double sigma_f, @@ -36,6 +38,8 @@ class SparseGP{ void update_alpha(); + void compute_likelihood(); + Eigen::VectorXd predict(StructureDescriptor test_structure); Eigen::VectorXd predict_serial(StructureDescriptor test_structure); }; From 93d3c61b3e88d5a628dfe7990182baccf3d11b70 Mon Sep 17 00:00:00 2001 From: Jonathan Vandermause Date: Sun, 8 Mar 2020 15:19:02 -0400 Subject: [PATCH 22/24] add kuu inverse to sparse gp --- src/sparse_gp.cpp | 2 ++ src/sparse_gp.h | 5 ++++- 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/src/sparse_gp.cpp b/src/sparse_gp.cpp index 31c5a0a19..47fb2e451 100644 --- a/src/sparse_gp.cpp +++ b/src/sparse_gp.cpp @@ -330,7 +330,9 @@ void SparseGP :: add_training_structure(StructureDescriptor training_structure){ void SparseGP::update_alpha(){ Eigen::MatrixXd sigma_inv = Kuu + Kuf * noise_matrix * Kuf.transpose(); + // TODO: Use Woodbury identity to perform inversion once. Sigma = sigma_inv.inverse(); + Kuu_inverse = Kuu.inverse(); alpha = Sigma * Kuf * noise_matrix * y; } diff --git a/src/sparse_gp.h b/src/sparse_gp.h index 1fceb9c94..e89e1453a 100644 --- a/src/sparse_gp.h +++ b/src/sparse_gp.h @@ -10,7 +10,7 @@ class SparseGP{ public: - Eigen::MatrixXd Kuu, Kuf, Sigma, noise_matrix; + Eigen::MatrixXd Kuu, Kuf, Sigma, Kuu_inverse, noise_matrix; Eigen::VectorXd y, alpha, hyperparameters, noise; std::vector kernels; @@ -42,6 +42,9 @@ class SparseGP{ Eigen::VectorXd predict(StructureDescriptor test_structure); Eigen::VectorXd predict_serial(StructureDescriptor test_structure); + + void predict_distribution(StructureDescriptor test_structure, + Eigen::VectorXd & mean_vector, Eigen::VectorXd & std_vector); }; #endif \ No newline at end of file From 462284c5eddd1e1b7c48414c7381011854e66590 Mon Sep 17 00:00:00 2001 From: Jonathan Vandermause Date: Sun, 8 Mar 2020 15:47:07 -0400 Subject: [PATCH 23/24] add env struc partial method to 3-body kernel --- src/kernels.cpp | 265 +++++++++++++++++++++++++----------------------- src/kernels.h | 9 ++ src/sparse_gp.h | 5 +- 3 files changed, 151 insertions(+), 128 deletions(-) diff --git a/src/kernels.cpp b/src/kernels.cpp index b31398218..36751c35a 100644 --- a/src/kernels.cpp +++ b/src/kernels.cpp @@ -307,8 +307,9 @@ double ThreeBodyKernel :: env_env(const LocalEnvironment & env1, return sig2 * kern; } -Eigen::VectorXd ThreeBodyKernel :: env_struc(const LocalEnvironment & env1, - const StructureDescriptor & struc1){ +Eigen::VectorXd ThreeBodyKernel :: env_struc_partial( + const LocalEnvironment & env1, const StructureDescriptor & struc1, + int atom){ int no_elements = 1 + 3 * struc1.noa + 6; Eigen::VectorXd kernel_vector = @@ -319,9 +320,9 @@ Eigen::VectorXd ThreeBodyKernel :: env_struc(const LocalEnvironment & env1, p2, p3, p4, xval1, xval2, xrel1, xrel2, yval1, yval2, yrel1, yrel2, zval1, zval2, zrel1, zrel2, fdjx1, fdjx2, fdjy1, fdjy2, fdjz1, fdjz2, fx1, fx2, fy1, fy2, fz1, fz2; - int i1, i2, j1, j2, ei1, ei2, ej1, ej2, c2; + int i1, i2, j1, j2, ei1, ei2, ej1, ej2; - LocalEnvironment env2; + LocalEnvironment env2 = struc1.local_environments[atom]; double vol_inv = 1 / struc1.volume; double cut1 = env1.n_body_cutoffs[1]; @@ -329,134 +330,130 @@ Eigen::VectorXd ThreeBodyKernel :: env_struc(const LocalEnvironment & env1, double rcut_vals_i1[2], rcut_vals_i2[2], rcut_vals_i3[2], rcut_vals_j1[2], rcut_vals_j2[2], rcut_vals_j3[2]; int c1 = env1.central_species; + int c2 = env2.central_species; + + for (int m = 0; m < env1.three_body_indices.size(); m ++){ + i1 = env1.three_body_indices[m][0]; + i2 = env1.three_body_indices[m][1]; - for (int i = 0; i < struc1.noa; i ++){ - env2 = struc1.local_environments[i]; - c2 = env2.central_species; - - for (int m = 0; m < env1.three_body_indices.size(); m ++){ - i1 = env1.three_body_indices[m][0]; - i2 = env1.three_body_indices[m][1]; - - ri1 = env1.rs[i1]; - ri2 = env1.rs[i2]; - ri3 = env1.cross_bond_dists[m]; - - ei1 = env1.environment_species[i1]; - ei2 = env1.environment_species[i2]; - - (*cutoff_pointer)(rcut_vals_i1, ri1, cut1, cutoff_hyps); - (*cutoff_pointer)(rcut_vals_i2, ri2, cut1, cutoff_hyps); - (*cutoff_pointer)(rcut_vals_i3, ri3, cut1, cutoff_hyps); - - fi1 = rcut_vals_i1[0]; - fi2 = rcut_vals_i2[0]; - fi3 = rcut_vals_i3[0]; - fi = fi1 * fi2 * fi3; - - for (int n = 0; n < env2.cross_bond_dists.size(); n ++){ - j1 = env2.three_body_indices[n][0]; - j2 = env2.three_body_indices[n][1]; - - rj1 = env2.rs[j1]; - rj2 = env2.rs[j2]; - rj3 = env2.cross_bond_dists[n]; - - ej1 = env2.environment_species[j1]; - ej2 = env2.environment_species[j2]; - - xval1 = env2.xs[j1]; - yval1 = env2.ys[j1]; - zval1 = env2.zs[j1]; - xrel1 = env2.xrel[j1]; - yrel1 = env2.yrel[j1]; - zrel1 = env2.zrel[j1]; - - xval2 = env2.xs[j2]; - yval2 = env2.ys[j2]; - zval2 = env2.zs[j2]; - xrel2 = env2.xrel[j2]; - yrel2 = env2.yrel[j2]; - zrel2 = env2.zrel[j2]; - - (*cutoff_pointer)(rcut_vals_j1, rj1, cut2, cutoff_hyps); - (*cutoff_pointer)(rcut_vals_j2, rj2, cut2, cutoff_hyps); - (*cutoff_pointer)(rcut_vals_j3, rj3, cut2, cutoff_hyps); - - fj1 = rcut_vals_j1[0]; - fdj1 = rcut_vals_j1[1]; - fj2 = rcut_vals_j2[0]; - fdj2 = rcut_vals_j2[1]; - fj3 = rcut_vals_j3[0]; - fj = fj1 * fj2 * fj3; - - fdjx1 = xrel1 * fdj1 * fj2 * fj3; - fdjy1 = yrel1 * fdj1 * fj2 * fj3; - fdjz1 = zrel1 * fdj1 * fj2 * fj3; - fdjx2 = xrel2 * fj1 * fdj2 * fj3; - fdjy2 = yrel2 * fj1 * fdj2 * fj3; - fdjz2 = zrel2 * fj1 * fdj2 * fj3; - - r11 = ri1 - rj1; - r12 = ri1 - rj2; - r13 = ri1 - rj3; - r21 = ri2 - rj1; - r22 = ri2 - rj2; - r23 = ri2 - rj3; - r31 = ri3 - rj1; - r32 = ri3 - rj2; - r33 = ri3 - rj3; - - // Sum over six permutations. - if (c1 == c2){ - if (ei1 == ej1 && ei2 == ej2){ - update_kernel_vector(kernel_vector, no_elements, i, - vol_inv, r11, r22, r33, fi, fj, - fdjx1, fdjx2, fdjy1, fdjy2, fdjz1, fdjz2, - xrel1, xval1, xrel2, xval2, yrel1, yval1, yrel2, - yval2, zrel1, zval1, zrel2, zval2); - } - if (ei1 == ej2 && ei2 == ej1){ - update_kernel_vector(kernel_vector, no_elements, i, - vol_inv, r21, r12, r33, fi, fj, - fdjx1, fdjx2, fdjy1, fdjy2, fdjz1, fdjz2, - xrel1, xval1, xrel2, xval2, yrel1, yval1, yrel2, - yval2, zrel1, zval1, zrel2, zval2); - } + ri1 = env1.rs[i1]; + ri2 = env1.rs[i2]; + ri3 = env1.cross_bond_dists[m]; + + ei1 = env1.environment_species[i1]; + ei2 = env1.environment_species[i2]; + + (*cutoff_pointer)(rcut_vals_i1, ri1, cut1, cutoff_hyps); + (*cutoff_pointer)(rcut_vals_i2, ri2, cut1, cutoff_hyps); + (*cutoff_pointer)(rcut_vals_i3, ri3, cut1, cutoff_hyps); + + fi1 = rcut_vals_i1[0]; + fi2 = rcut_vals_i2[0]; + fi3 = rcut_vals_i3[0]; + fi = fi1 * fi2 * fi3; + + for (int n = 0; n < env2.cross_bond_dists.size(); n ++){ + j1 = env2.three_body_indices[n][0]; + j2 = env2.three_body_indices[n][1]; + + rj1 = env2.rs[j1]; + rj2 = env2.rs[j2]; + rj3 = env2.cross_bond_dists[n]; + + ej1 = env2.environment_species[j1]; + ej2 = env2.environment_species[j2]; + + xval1 = env2.xs[j1]; + yval1 = env2.ys[j1]; + zval1 = env2.zs[j1]; + xrel1 = env2.xrel[j1]; + yrel1 = env2.yrel[j1]; + zrel1 = env2.zrel[j1]; + + xval2 = env2.xs[j2]; + yval2 = env2.ys[j2]; + zval2 = env2.zs[j2]; + xrel2 = env2.xrel[j2]; + yrel2 = env2.yrel[j2]; + zrel2 = env2.zrel[j2]; + + (*cutoff_pointer)(rcut_vals_j1, rj1, cut2, cutoff_hyps); + (*cutoff_pointer)(rcut_vals_j2, rj2, cut2, cutoff_hyps); + (*cutoff_pointer)(rcut_vals_j3, rj3, cut2, cutoff_hyps); + + fj1 = rcut_vals_j1[0]; + fdj1 = rcut_vals_j1[1]; + fj2 = rcut_vals_j2[0]; + fdj2 = rcut_vals_j2[1]; + fj3 = rcut_vals_j3[0]; + fj = fj1 * fj2 * fj3; + + fdjx1 = xrel1 * fdj1 * fj2 * fj3; + fdjy1 = yrel1 * fdj1 * fj2 * fj3; + fdjz1 = zrel1 * fdj1 * fj2 * fj3; + fdjx2 = xrel2 * fj1 * fdj2 * fj3; + fdjy2 = yrel2 * fj1 * fdj2 * fj3; + fdjz2 = zrel2 * fj1 * fdj2 * fj3; + + r11 = ri1 - rj1; + r12 = ri1 - rj2; + r13 = ri1 - rj3; + r21 = ri2 - rj1; + r22 = ri2 - rj2; + r23 = ri2 - rj3; + r31 = ri3 - rj1; + r32 = ri3 - rj2; + r33 = ri3 - rj3; + + // Sum over six permutations. + if (c1 == c2){ + if (ei1 == ej1 && ei2 == ej2){ + update_kernel_vector(kernel_vector, no_elements, atom, + vol_inv, r11, r22, r33, fi, fj, + fdjx1, fdjx2, fdjy1, fdjy2, fdjz1, fdjz2, + xrel1, xval1, xrel2, xval2, yrel1, yval1, yrel2, + yval2, zrel1, zval1, zrel2, zval2); + } + if (ei1 == ej2 && ei2 == ej1){ + update_kernel_vector(kernel_vector, no_elements, atom, + vol_inv, r21, r12, r33, fi, fj, + fdjx1, fdjx2, fdjy1, fdjy2, fdjz1, fdjz2, + xrel1, xval1, xrel2, xval2, yrel1, yval1, yrel2, + yval2, zrel1, zval1, zrel2, zval2); } + } - if (c1 == ej1){ - if (ei1 == ej2 && ei2 == c2){ - update_kernel_vector(kernel_vector, no_elements, i, - vol_inv, r21, r32, r13, fi, fj, - fdjx1, fdjx2, fdjy1, fdjy2, fdjz1, fdjz2, - xrel1, xval1, xrel2, xval2, yrel1, yval1, yrel2, - yval2, zrel1, zval1, zrel2, zval2); - } - if (ei1 == c2 && ei2 == ej2){ - update_kernel_vector(kernel_vector, no_elements, i, - vol_inv, r11, r32, r23, fi, fj, - fdjx1, fdjx2, fdjy1, fdjy2, fdjz1, fdjz2, - xrel1, xval1, xrel2, xval2, yrel1, yval1, yrel2, - yval2, zrel1, zval1, zrel2, zval2); - } + if (c1 == ej1){ + if (ei1 == ej2 && ei2 == c2){ + update_kernel_vector(kernel_vector, no_elements, atom, + vol_inv, r21, r32, r13, fi, fj, + fdjx1, fdjx2, fdjy1, fdjy2, fdjz1, fdjz2, + xrel1, xval1, xrel2, xval2, yrel1, yval1, yrel2, + yval2, zrel1, zval1, zrel2, zval2); + } + if (ei1 == c2 && ei2 == ej2){ + update_kernel_vector(kernel_vector, no_elements, atom, + vol_inv, r11, r32, r23, fi, fj, + fdjx1, fdjx2, fdjy1, fdjy2, fdjz1, fdjz2, + xrel1, xval1, xrel2, xval2, yrel1, yval1, yrel2, + yval2, zrel1, zval1, zrel2, zval2); } + } - if (c1 == ej2){ - if (ei1 == ej1 && ei2 == c2){ - update_kernel_vector(kernel_vector, no_elements, i, - vol_inv, r31, r22, r13, fi, fj, - fdjx1, fdjx2, fdjy1, fdjy2, fdjz1, fdjz2, - xrel1, xval1, xrel2, xval2, yrel1, yval1, yrel2, - yval2, zrel1, zval1, zrel2, zval2); - } - if (ei1 == c2 && ei2 == ej1){ - update_kernel_vector(kernel_vector, no_elements, i, - vol_inv, r31, r12, r23, fi, fj, - fdjx1, fdjx2, fdjy1, fdjy2, fdjz1, fdjz2, - xrel1, xval1, xrel2, xval2, yrel1, yval1, yrel2, - yval2, zrel1, zval1, zrel2, zval2); - } + if (c1 == ej2){ + if (ei1 == ej1 && ei2 == c2){ + update_kernel_vector(kernel_vector, no_elements, atom, + vol_inv, r31, r22, r13, fi, fj, + fdjx1, fdjx2, fdjy1, fdjy2, fdjz1, fdjz2, + xrel1, xval1, xrel2, xval2, yrel1, yval1, yrel2, + yval2, zrel1, zval1, zrel2, zval2); + } + if (ei1 == c2 && ei2 == ej1){ + update_kernel_vector(kernel_vector, no_elements, atom, + vol_inv, r31, r12, r23, fi, fj, + fdjx1, fdjx2, fdjy1, fdjy2, fdjz1, fdjz2, + xrel1, xval1, xrel2, xval2, yrel1, yval1, yrel2, + yval2, zrel1, zval1, zrel2, zval2); } } } @@ -465,6 +462,20 @@ Eigen::VectorXd ThreeBodyKernel :: env_struc(const LocalEnvironment & env1, return kernel_vector; } +Eigen::VectorXd ThreeBodyKernel :: env_struc(const LocalEnvironment & env1, + const StructureDescriptor & struc1){ + + int no_elements = 1 + 3 * struc1.noa + 6; + Eigen::VectorXd kernel_vector = + Eigen::VectorXd::Zero(no_elements); + + for (int i = 0; i < struc1.noa; i ++){ + kernel_vector += env_struc_partial(env1, struc1, i); + } + + return kernel_vector; +} + void ThreeBodyKernel :: update_kernel_vector(Eigen::VectorXd & kernel_vector, int no_elements, int i, double vol_inv, double r11, double r22, double r33, diff --git a/src/kernels.h b/src/kernels.h index f52ac62e5..feaef7630 100644 --- a/src/kernels.h +++ b/src/kernels.h @@ -38,6 +38,9 @@ class DotProductKernel : public Kernel{ Eigen::VectorXd env_struc(const LocalEnvironment & env1, const StructureDescriptor & struc1); + + Eigen::VectorXd env_struc_partial(const LocalEnvironment & env1, + const StructureDescriptor & struc1, int atom); }; @@ -58,6 +61,9 @@ class TwoBodyKernel : public Kernel{ Eigen::VectorXd env_struc(const LocalEnvironment & env1, const StructureDescriptor & struc1); + + Eigen::VectorXd env_struc_partial(const LocalEnvironment & env1, + const StructureDescriptor & struc1, int atom); }; class ThreeBodyKernel : public Kernel{ @@ -78,6 +84,9 @@ class ThreeBodyKernel : public Kernel{ Eigen::VectorXd env_struc(const LocalEnvironment & env1, const StructureDescriptor & struc1); + Eigen::VectorXd env_struc_partial(const LocalEnvironment & env1, + const StructureDescriptor & struc1, int atom); + void update_kernel_vector(Eigen::VectorXd & kernel_vector, int no_elements, int i, double vol_inv, double r11, double r22, double r33, diff --git a/src/sparse_gp.h b/src/sparse_gp.h index e89e1453a..2623c77d2 100644 --- a/src/sparse_gp.h +++ b/src/sparse_gp.h @@ -43,7 +43,10 @@ class SparseGP{ Eigen::VectorXd predict(StructureDescriptor test_structure); Eigen::VectorXd predict_serial(StructureDescriptor test_structure); - void predict_distribution(StructureDescriptor test_structure, + void predict_DTC(StructureDescriptor test_structure, + Eigen::VectorXd & mean_vector, Eigen::VectorXd & std_vector); + + void predict_SOR(StructureDescriptor test_structure, Eigen::VectorXd & mean_vector, Eigen::VectorXd & std_vector); }; From 8613c5cd3c9e4e03b401326b9c785c05ac41c696 Mon Sep 17 00:00:00 2001 From: Jonathan Vandermause Date: Sun, 8 Mar 2020 16:16:11 -0400 Subject: [PATCH 24/24] add env struc partial methods to 2- and many-body kernels --- src/kernels.cpp | 243 ++++++++++++++++++++++++++---------------------- src/kernels.h | 25 ++--- 2 files changed, 147 insertions(+), 121 deletions(-) diff --git a/src/kernels.cpp b/src/kernels.cpp index 36751c35a..902f0248e 100644 --- a/src/kernels.cpp +++ b/src/kernels.cpp @@ -75,9 +75,10 @@ double TwoBodyKernel :: env_env(const LocalEnvironment & env1, return sig2 * kern; } -Eigen::VectorXd TwoBodyKernel :: env_struc(const LocalEnvironment & env1, - const StructureDescriptor & struc1){ - +Eigen::VectorXd TwoBodyKernel :: env_struc_partial( + const LocalEnvironment & env1, const StructureDescriptor & struc1, + int atom){ + int no_elements = 1 + 3 * struc1.noa + 6; Eigen::VectorXd kernel_vector = Eigen::VectorXd::Zero(no_elements); @@ -92,80 +93,89 @@ Eigen::VectorXd TwoBodyKernel :: env_struc(const LocalEnvironment & env1, double rcut_vals_1[2]; double rcut_vals_2[2]; + double vol_inv = 1 / struc1.volume; + + LocalEnvironment env_curr = struc1.local_environments[atom]; std::vector inds1 = env1.n_body_indices[0]; - std::vector inds2; + std::vector inds2 = env_curr.n_body_indices[0]; + cent2 = env_curr.central_species; - double vol_inv = 1 / struc1.volume; + for (int m = 0; m < inds1.size(); m ++){ + ind1 = inds1[m]; + ri = env1.rs[ind1]; + (*cutoff_pointer)(rcut_vals_1, ri, cut1, cutoff_hyps); + fi = rcut_vals_1[0]; + e1 = env1.environment_species[ind1]; - LocalEnvironment env_curr; - for (int i = 0; i < struc1.noa; i ++){ - env_curr = struc1.local_environments[i]; - inds2 = env_curr.n_body_indices[0]; - cent2 = env_curr.central_species; - - for (int m = 0; m < inds1.size(); m ++){ - ind1 = inds1[m]; - ri = env1.rs[ind1]; - (*cutoff_pointer)(rcut_vals_1, ri, cut1, cutoff_hyps); - fi = rcut_vals_1[0]; - e1 = env1.environment_species[ind1]; - - for (int n = 0; n < inds2.size(); n ++){ - ind2 = inds2[n]; - e2 = env_curr.environment_species[ind2]; - - // Proceed only if the pairs match. - if ((cent1 == cent2 && e1 == e2) || (cent1 == e2 && - cent2 == e1)){ - rj = env_curr.rs[ind2]; - rdiff = ri - rj; - - xval = env_curr.xs[ind2]; - yval = env_curr.ys[ind2]; - zval = env_curr.zs[ind2]; - xrel = env_curr.xrel[ind2]; - yrel = env_curr.yrel[ind2]; - zrel = env_curr.zrel[ind2]; - - (*cutoff_pointer)(rcut_vals_2, rj, cut2, cutoff_hyps); - fj = rcut_vals_2[0]; - fdj = rcut_vals_2[1]; - - // energy kernel - c1 = rdiff * rdiff; - c2 = exp(-c1 * ls1); - kernel_vector(0) += sig2 * fi * fj * c2 / 2; - - // helper constants - c3 = c2 * ls2 * fi * fj * rdiff; - c4 = c2 * fi * fdj; - - // fx + exx, exy, exz stress components - fx = xrel * c3 + c4 * xrel; - kernel_vector(1 + 3 * i) += sig2 * fx; - kernel_vector(no_elements - 6) -= - sig2 * fx * xval * vol_inv / 2; - kernel_vector(no_elements - 5) -= - sig2 * fx * yval * vol_inv / 2; - kernel_vector(no_elements - 4) -= - sig2 * fx * zval * vol_inv / 2; - - // fy + eyy, eyz stress components - fy = yrel * c3 + c4 * yrel; - kernel_vector(2 + 3 * i) += sig2 * fy; - kernel_vector(no_elements - 3) -= - sig2 * fy * yval * vol_inv / 2; - kernel_vector(no_elements - 2) -= - sig2 * fy * zval * vol_inv / 2; - - // fz + ezz stress component - fz = zrel * c3 + c4 * zrel; - kernel_vector(3 + 3 * i) += sig2 * fz; - kernel_vector(no_elements - 1) -= - sig2 * fz * zval * vol_inv / 2; - } + for (int n = 0; n < inds2.size(); n ++){ + ind2 = inds2[n]; + e2 = env_curr.environment_species[ind2]; + + // Proceed only if the pairs match. + if ((cent1 == cent2 && e1 == e2) || (cent1 == e2 && + cent2 == e1)){ + rj = env_curr.rs[ind2]; + rdiff = ri - rj; + + xval = env_curr.xs[ind2]; + yval = env_curr.ys[ind2]; + zval = env_curr.zs[ind2]; + xrel = env_curr.xrel[ind2]; + yrel = env_curr.yrel[ind2]; + zrel = env_curr.zrel[ind2]; + + (*cutoff_pointer)(rcut_vals_2, rj, cut2, cutoff_hyps); + fj = rcut_vals_2[0]; + fdj = rcut_vals_2[1]; + + // energy kernel + c1 = rdiff * rdiff; + c2 = exp(-c1 * ls1); + kernel_vector(0) += sig2 * fi * fj * c2 / 2; + + // helper constants + c3 = c2 * ls2 * fi * fj * rdiff; + c4 = c2 * fi * fdj; + + // fx + exx, exy, exz stress components + fx = xrel * c3 + c4 * xrel; + kernel_vector(1 + 3 * atom) += sig2 * fx; + kernel_vector(no_elements - 6) -= + sig2 * fx * xval * vol_inv / 2; + kernel_vector(no_elements - 5) -= + sig2 * fx * yval * vol_inv / 2; + kernel_vector(no_elements - 4) -= + sig2 * fx * zval * vol_inv / 2; + + // fy + eyy, eyz stress components + fy = yrel * c3 + c4 * yrel; + kernel_vector(2 + 3 * atom) += sig2 * fy; + kernel_vector(no_elements - 3) -= + sig2 * fy * yval * vol_inv / 2; + kernel_vector(no_elements - 2) -= + sig2 * fy * zval * vol_inv / 2; + + // fz + ezz stress component + fz = zrel * c3 + c4 * zrel; + kernel_vector(3 + 3 * atom) += sig2 * fz; + kernel_vector(no_elements - 1) -= + sig2 * fz * zval * vol_inv / 2; } - } + } + } + + return kernel_vector; +} + +Eigen::VectorXd TwoBodyKernel :: env_struc(const LocalEnvironment & env1, + const StructureDescriptor & struc1){ + + int no_elements = 1 + 3 * struc1.noa + 6; + Eigen::VectorXd kernel_vector = + Eigen::VectorXd::Zero(no_elements); + + for (int i = 0; i < struc1.noa; i ++){ + kernel_vector += env_struc_partial(env1, struc1, i); } return kernel_vector; @@ -543,9 +553,9 @@ double DotProductKernel :: env_env(const LocalEnvironment & env1, return sig2 * pow(dot / (d1 * d2), power); } -Eigen::VectorXd DotProductKernel - :: env_struc(const LocalEnvironment & env1, - const StructureDescriptor & struc1){ +Eigen::VectorXd DotProductKernel :: env_struc_partial( + const LocalEnvironment & env1, const StructureDescriptor & struc1, + int atom){ Eigen::VectorXd kern_vec = Eigen::VectorXd::Zero(1 + struc1.noa * 3 + 6); @@ -561,43 +571,56 @@ Eigen::VectorXd DotProductKernel Eigen::VectorXd force_dot, stress_dot, f1, s1; const double vol_inv = 1 / struc1.volume; double dot_val, d2, norm_dot, dval, d2_cubed; - LocalEnvironment env_curr; - - for (int i = 0; i < struc1.noa; i ++){ - env_curr = struc1.local_environments[i]; - - // Check that the environments have the same central species. - if (env1.central_species != env_curr.central_species) continue; - - // Check that d2 is nonzero. - d2 = env_curr.descriptor_norm[descriptor_index]; - if (d2 < empty_thresh) continue; - d2_cubed = d2 * d2 * d2; - - // Energy kernel - dot_val = env1.descriptor_vals[descriptor_index] - .dot(env_curr.descriptor_vals[descriptor_index]); - norm_dot = dot_val / (d1 * d2); - en_kern += pow(norm_dot, power); - - // Force kernel - force_dot = env_curr.descriptor_force_dervs[descriptor_index] * - env1.descriptor_vals[descriptor_index]; - f1 = (force_dot / (d1 * d2)) - - (dot_val * env_curr.force_dot[descriptor_index] / (d2_cubed * d1)); - dval = power * pow(norm_dot, power - 1); - force_kern += dval * f1; - - // Stress kernel - stress_dot = env_curr.descriptor_stress_dervs[descriptor_index] * - env1.descriptor_vals[descriptor_index]; - s1 = (stress_dot / (d1 * d2)) - - (dot_val * env_curr.stress_dot[descriptor_index] /(d2_cubed * d1)); - stress_kern += dval * s1; - } + LocalEnvironment env_curr = struc1.local_environments[atom]; + + // Check that the environments have the same central species. + if (env1.central_species != env_curr.central_species){ + return kern_vec; + }; + + // Check that d2 is nonzero. + d2 = env_curr.descriptor_norm[descriptor_index]; + if (d2 < empty_thresh){ + return kern_vec; + }; + d2_cubed = d2 * d2 * d2; + + // Energy kernel + dot_val = env1.descriptor_vals[descriptor_index] + .dot(env_curr.descriptor_vals[descriptor_index]); + norm_dot = dot_val / (d1 * d2); + en_kern += pow(norm_dot, power); + + // Force kernel + force_dot = env_curr.descriptor_force_dervs[descriptor_index] * + env1.descriptor_vals[descriptor_index]; + f1 = (force_dot / (d1 * d2)) - + (dot_val * env_curr.force_dot[descriptor_index] / (d2_cubed * d1)); + dval = power * pow(norm_dot, power - 1); + force_kern += dval * f1; + + // Stress kernel + stress_dot = env_curr.descriptor_stress_dervs[descriptor_index] * + env1.descriptor_vals[descriptor_index]; + s1 = (stress_dot / (d1 * d2)) - + (dot_val * env_curr.stress_dot[descriptor_index] /(d2_cubed * d1)); + stress_kern += dval * s1; kern_vec(0) = sig2 * en_kern; kern_vec.segment(1, struc1.noa * 3) = -sig2 * force_kern; kern_vec.tail(6) = -sig2 * stress_kern * vol_inv; return kern_vec; } + +Eigen::VectorXd DotProductKernel + :: env_struc(const LocalEnvironment & env1, + const StructureDescriptor & struc1){ + + Eigen::VectorXd kern_vec = Eigen::VectorXd::Zero(1 + struc1.noa * 3 + 6); + + for (int i = 0; i < struc1.noa; i ++){ + kern_vec += env_struc_partial(env1, struc1, i); + } + + return kern_vec; +} diff --git a/src/kernels.h b/src/kernels.h index feaef7630..6e7792600 100644 --- a/src/kernels.h +++ b/src/kernels.h @@ -18,6 +18,10 @@ class Kernel{ virtual double env_env(const LocalEnvironment & env1, const LocalEnvironment & env2) = 0; + virtual Eigen::VectorXd env_struc_partial( + const LocalEnvironment & env1, const StructureDescriptor & struc1, + int atom) = 0; + virtual Eigen::VectorXd env_struc(const LocalEnvironment & env1, const StructureDescriptor & struc1) = 0; @@ -35,13 +39,12 @@ class DotProductKernel : public Kernel{ double env_env(const LocalEnvironment & env1, const LocalEnvironment & env2); - - Eigen::VectorXd env_struc(const LocalEnvironment & env1, - const StructureDescriptor & struc1); Eigen::VectorXd env_struc_partial(const LocalEnvironment & env1, const StructureDescriptor & struc1, int atom); + Eigen::VectorXd env_struc(const LocalEnvironment & env1, + const StructureDescriptor & struc1); }; class TwoBodyKernel : public Kernel{ @@ -58,12 +61,12 @@ class TwoBodyKernel : public Kernel{ double env_env(const LocalEnvironment & env1, const LocalEnvironment & env2); - - Eigen::VectorXd env_struc(const LocalEnvironment & env1, - const StructureDescriptor & struc1); - + Eigen::VectorXd env_struc_partial(const LocalEnvironment & env1, const StructureDescriptor & struc1, int atom); + + Eigen::VectorXd env_struc(const LocalEnvironment & env1, + const StructureDescriptor & struc1); }; class ThreeBodyKernel : public Kernel{ @@ -81,12 +84,12 @@ class ThreeBodyKernel : public Kernel{ double env_env(const LocalEnvironment & env1, const LocalEnvironment & env2); - Eigen::VectorXd env_struc(const LocalEnvironment & env1, - const StructureDescriptor & struc1); - Eigen::VectorXd env_struc_partial(const LocalEnvironment & env1, const StructureDescriptor & struc1, int atom); + Eigen::VectorXd env_struc(const LocalEnvironment & env1, + const StructureDescriptor & struc1); + void update_kernel_vector(Eigen::VectorXd & kernel_vector, int no_elements, int i, double vol_inv, double r11, double r22, double r33, @@ -97,4 +100,4 @@ class ThreeBodyKernel : public Kernel{ double zrel1, double zval1, double zrel2, double zval2); }; -#endif \ No newline at end of file +#endif