Repository: mir-group/flare_pp (language: python)

Issue:
Pulling master into my dev branch 



Discussion / hints:

Fix patch:
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 0e480101..4e194174 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -6,7 +6,7 @@ if(${CMAKE_VERSION} VERSION_LESS 3.12)
 endif()
 
 # set the project name
-project(ACE_Tests VERSION 1.0
+project(flare_pp VERSION 1.0
 	    DESCRIPTION "High performance implementation of the atomic cluster expansion."
 	    LANGUAGES CXX)
 
@@ -29,14 +29,11 @@ set(CMAKE_CXX_STANDARD_REQUIRED True)
 add_subdirectory(src)
 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
-		    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)
+add_subdirectory(pybind11)
+pybind11_add_module(ace SHARED src/ace_binding.cpp)
+target_link_libraries(ace PRIVATE flare_pp)
+target_include_directories(ace PRIVATE ${ACE_INCLUDE_DIR})
diff --git a/timing/learning_curve.py b/scripts/learning_curve.py
similarity index 100%
rename from timing/learning_curve.py
rename to scripts/learning_curve.py
diff --git a/scripts/n_body.py b/scripts/n_body.py
new file mode 100644
index 00000000..757692b2
--- /dev/null
+++ b/scripts/n_body.py
@@ -0,0 +1,115 @@
+import argparse
+import os
+import pickle
+
+import numpy as np
+
+import ace
+
+# Gets data dir from input
+parser = argparse.ArgumentParser()
+parser.add_argument('--data', type=str, required=True)
+args = parser.parse_args()
+
+# create two and three body kernels
+kernel_2b = ace.TwoBodyKernel(0.02, 0.60, "quadratic", [])
+kernel_3b = ace.ThreeBodyKernel(0.013, 0.2, "quadratic", [])
+
+# Defines the cutoff radii
+cutoff = 10.0
+nested_cutoffs = [10.0, 4.]
+
+# create sparse gp
+sigma_e = 0.1
+sigma_f = 0.1
+sigma_s = 0.1
+sparse_gp = ace.SparseGP([kernel_2b, kernel_3b],
+                         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
+data_path = args.data
+file_base = 'extract_run'
+file_tail = '.pkl'
+file_nos = [0]
+snapshots = range(0, 100, 50)
+test_snapshot = 300
+
+# create test structure
+test_no = 0
+data_file = '{}{}{}'.format(file_base, test_no, file_tail)
+open_file = open(os.path.join(data_path, 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)
+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(os.path.join(data_path, 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)
+        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 i in sparse_inds:
+            sparse_gp.add_sparse_environment(train_struc.local_environments[i])
+
+        # 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)
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 7aaadba5..9dca1db2 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -1,9 +1,3 @@
-# make library
-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)
+add_subdirectory(flare_pp)
 
-find_package(OpenMP REQUIRED)
-target_link_libraries(Ace PRIVATE OpenMP::OpenMP_CXX)
-
-set(ACE_INCLUDE_DIR "${ACE_ROOT_DIR}/src" PARENT_SCOPE)
+set(ACE_INCLUDE_DIR ${ACE_INCLUDE_DIR} PARENT_SCOPE)
diff --git a/src/flare_pp/CMakeLists.txt b/src/flare_pp/CMakeLists.txt
new file mode 100644
index 00000000..f501b993
--- /dev/null
+++ b/src/flare_pp/CMakeLists.txt
@@ -0,0 +1,11 @@
+# make library
+add_library(flare_pp SHARED 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)
+if(OpenMP_CXX_FOUND)
+    target_link_libraries(flare_pp PUBLIC OpenMP::OpenMP_CXX)
+endif()
+
+set(ACE_INCLUDE_DIR "${ACE_ROOT_DIR}/src/flare_pp" PARENT_SCOPE)
diff --git a/src/cutoffs.cpp b/src/flare_pp/cutoffs.cpp
similarity index 100%
rename from src/cutoffs.cpp
rename to src/flare_pp/cutoffs.cpp
diff --git a/src/cutoffs.h b/src/flare_pp/cutoffs.h
similarity index 100%
rename from src/cutoffs.h
rename to src/flare_pp/cutoffs.h
diff --git a/src/descriptor.cpp b/src/flare_pp/descriptor.cpp
similarity index 92%
rename from src/descriptor.cpp
rename to src/flare_pp/descriptor.cpp
index 5b9ae5ae..d437cd03 100644
--- a/src/descriptor.cpp
+++ b/src/flare_pp/descriptor.cpp
@@ -6,6 +6,7 @@
 #include <cmath>
 #include <iostream>
 
+
 DescriptorCalculator::DescriptorCalculator(){}
 
 DescriptorCalculator::DescriptorCalculator(
@@ -38,6 +39,7 @@ DescriptorCalculator::DescriptorCalculator(
 
 }
 
+
 void B2_descriptor(
 Eigen::VectorXd & B2_vals,
 Eigen::MatrixXd & B2_force_dervs,
@@ -51,18 +53,25 @@ int neigh_size = env.neighbor_list.size();
 int cent_ind = env.central_index;
 int no_radial = nos * N;
 int no_harmonics = (lmax + 1) * (lmax + 1);
-int n1_ind, n2_ind, l_ind, n1_l, n2_l, env_ind;
-int counter = 0;
 
-for (int n1 = 0; n1 < no_radial; n1 ++){
-    n1_ind = n1 * no_harmonics;
+
+int n1_l, n2_l, env_ind;
+int counter;
+int n1_count;
+int n2_count;
+
+for (int n1 = no_radial-1; n1 >= 0; n1 --){
+
     for (int n2 = n1; n2 < no_radial; n2 ++){
-        n2_ind = n2 * no_harmonics;
-        l_ind = 0;
+
         for (int l = 0; l < (lmax + 1); l ++){
             for (int m = 0; m < (2 * l + 1); m ++){
-                n1_l = n1_ind + l_ind;
-                n2_l = n2_ind + l_ind;
+                n1_l = n1 * no_harmonics + (l * l + m);
+                n2_l = n2 * no_harmonics + (l * l + m);
+
+                n1_count = (n1 * (2*no_radial - n1 + 1)) / 2;
+                n2_count = n2 - n1;
+                counter = l + (n1_count + n2_count) * (lmax+1);
 
                 // Store B2 value.
                 B2_vals(counter) +=
@@ -91,12 +100,11 @@ for (int n1 = 0; n1 < no_radial; n1 ++){
                         single_bond_stress_dervs(p, n1_l) *
                         single_bond_vals(n2_l);
                 }
-                l_ind ++;
             }
-            counter ++;
         }
     }
 }
+
 };
 
 B1_Calculator :: B1_Calculator(){}
@@ -154,8 +162,9 @@ void B2_Calculator :: compute(const LocalEnvironment & env){
     int no_radial = nos * N;
 
     int no_harmonics = (lmax + 1) * (lmax + 1);
-    int no_bond = nos * N * no_harmonics; 
-    int no_descriptors = no_radial * (no_radial + 1) * (lmax + 1) / 2;
+
+    int no_bond = no_radial * no_harmonics;
+    int no_descriptors = (no_radial * (no_radial + 1) / 2) * (lmax + 1);
 
     single_bond_vals = Eigen::VectorXd::Zero(no_bond);
     single_bond_force_dervs =
@@ -175,7 +184,7 @@ void B2_Calculator :: compute(const LocalEnvironment & env){
         Eigen::MatrixXd::Zero(env.noa * 3, no_descriptors);
     descriptor_stress_dervs =
         Eigen::MatrixXd::Zero(6, no_descriptors);
-    
+
     B2_descriptor(descriptor_vals, descriptor_force_dervs,
                   descriptor_stress_dervs,
                   single_bond_vals, single_bond_force_dervs,
diff --git a/src/descriptor.h b/src/flare_pp/descriptor.h
similarity index 100%
rename from src/descriptor.h
rename to src/flare_pp/descriptor.h
diff --git a/src/kernels.cpp b/src/flare_pp/kernels.cpp
similarity index 100%
rename from src/kernels.cpp
rename to src/flare_pp/kernels.cpp
diff --git a/src/kernels.h b/src/flare_pp/kernels.h
similarity index 100%
rename from src/kernels.h
rename to src/flare_pp/kernels.h
diff --git a/src/local_environment.cpp b/src/flare_pp/local_environment.cpp
similarity index 100%
rename from src/local_environment.cpp
rename to src/flare_pp/local_environment.cpp
diff --git a/src/local_environment.h b/src/flare_pp/local_environment.h
similarity index 100%
rename from src/local_environment.h
rename to src/flare_pp/local_environment.h
diff --git a/src/radial.cpp b/src/flare_pp/radial.cpp
similarity index 100%
rename from src/radial.cpp
rename to src/flare_pp/radial.cpp
diff --git a/src/radial.h b/src/flare_pp/radial.h
similarity index 100%
rename from src/radial.h
rename to src/flare_pp/radial.h
diff --git a/src/single_bond.cpp b/src/flare_pp/single_bond.cpp
similarity index 100%
rename from src/single_bond.cpp
rename to src/flare_pp/single_bond.cpp
diff --git a/src/single_bond.h b/src/flare_pp/single_bond.h
similarity index 100%
rename from src/single_bond.h
rename to src/flare_pp/single_bond.h
diff --git a/src/sparse_gp.cpp b/src/flare_pp/sparse_gp.cpp
similarity index 100%
rename from src/sparse_gp.cpp
rename to src/flare_pp/sparse_gp.cpp
diff --git a/src/sparse_gp.h b/src/flare_pp/sparse_gp.h
similarity index 100%
rename from src/sparse_gp.h
rename to src/flare_pp/sparse_gp.h
diff --git a/src/structure.cpp b/src/flare_pp/structure.cpp
similarity index 100%
rename from src/structure.cpp
rename to src/flare_pp/structure.cpp
diff --git a/src/structure.h b/src/flare_pp/structure.h
similarity index 100%
rename from src/structure.h
rename to src/flare_pp/structure.h
diff --git a/src/y_grad.cpp b/src/flare_pp/y_grad.cpp
similarity index 100%
rename from src/y_grad.cpp
rename to src/flare_pp/y_grad.cpp
diff --git a/src/y_grad.h b/src/flare_pp/y_grad.h
similarity index 100%
rename from src/y_grad.h
rename to src/flare_pp/y_grad.h
diff --git a/timing/CMakeLists.txt b/timing/CMakeLists.txt
index 7c8f4821..22f2c0f2 100644
--- a/timing/CMakeLists.txt
+++ b/timing/CMakeLists.txt
@@ -1,4 +1,9 @@
 # make executable
 add_executable(time time_single_bond.cpp)
 target_include_directories(time PUBLIC ${ACE_INCLUDE_DIR})
-target_link_libraries(time PUBLIC Ace)
+target_link_libraries(time PUBLIC flare_pp)
+
+
+add_executable(benchmark benchmark_B2.cpp)
+target_include_directories(benchmark PUBLIC ${ACE_INCLUDE_DIR})
+target_link_libraries(benchmark PUBLIC flare_pp)
diff --git a/timing/benchmark_B2.cpp b/timing/benchmark_B2.cpp
new file mode 100644
index 00000000..470195bd
--- /dev/null
+++ b/timing/benchmark_B2.cpp
@@ -0,0 +1,94 @@
+#include <chrono>
+#include <iostream>
+#include <cmath>
+
+#include <Eigen/Dense>
+
+#include "structure.h"
+#include "local_environment.h"
+
+
+void construct_structure(Structure** struc) {
+    int noa = 50;
+    int nos = 5;
+
+    std::vector<int> species = std::vector<int>(noa);
+    for (int i = 0; i < noa; i++) {
+        species[i] = rand() % nos;
+    }
+
+    Eigen::MatrixXd cell {3, 3};
+    Eigen::MatrixXd positions_1 = 10 * Eigen::MatrixXd::Random(noa, 3);
+
+
+    cell << 10, 0.52, 0.12,
+           -0.93, 10, 0.32,
+            0.1, -0.2, 10;
+
+    (*struc) = new Structure(cell, species, positions_1);
+}
+
+void deconstruct_structure(Structure* struc) {
+    delete struc;
+}
+
+void construct_environment(LocalEnvironment** env, Structure* struc) {
+    double rcut = 10;
+    (*env) = new LocalEnvironment(*struc, 0, rcut);
+}
+
+void deconstruct_environment(LocalEnvironment* env) {
+    delete env;
+}
+
+void construct_calculator(B2_Calculator** calc) {
+    std::string radial_string = "chebyshev";
+    std::string cutoff_string = "quadratic";
+
+    std::vector<double> cutoff_hyps;
+
+    double first_gauss = 0;
+    double final_gauss = 10;
+    std::vector<double> radial_hyps = {first_gauss, final_gauss};
+
+    std::vector<int> descriptor_settings {5, 10, 10};
+
+    (*calc) = new B2_Calculator(radial_string, cutoff_string, radial_hyps,
+                            cutoff_hyps, descriptor_settings, 0);
+}
+
+void deconstruct_calculator(B2_Calculator* calc) {
+    delete calc;
+}
+
+int main() {
+    Structure* struc;
+    LocalEnvironment* env;
+    B2_Calculator* calc;
+    double duration = 0;
+    int num_reps = 20;
+
+    construct_calculator(&calc);
+
+
+    for (int i = 0; i < num_reps; i++) {
+        construct_structure(&struc);
+        construct_environment(&env, struc);
+
+        std::chrono::high_resolution_clock::time_point t1 = std::chrono::high_resolution_clock::now();
+
+        calc->compute(*env);
+
+        std::chrono::high_resolution_clock::time_point t2 = std::chrono::high_resolution_clock::now();
+        duration += (double) std::chrono::duration_cast<std::chrono::milliseconds>( t2 - t1 ).count();
+
+        deconstruct_structure(struc);
+        deconstruct_environment(env);
+    }
+
+
+    std::cout << "time: " << duration << " ms" << std::endl;
+
+    deconstruct_calculator(calc);
+    return 0;
+}
diff --git a/timing/time_single_bond.cpp b/timing/time_single_bond.cpp
index 171eca06..6c3c6833 100644
--- a/timing/time_single_bond.cpp
+++ b/timing/time_single_bond.cpp
@@ -1,9 +1,11 @@
-#include "y_grad.h"
-#include "radial.h"
-#include "cutoffs.h"
 #include <chrono>
 #include <iostream>
 #include <cmath>
+
+#include "y_grad.h"
+#include "radial.h"
+#include "cutoffs.h"
+
 using namespace std;
 
 int main(){
@@ -25,8 +27,8 @@ int main(){
 
     // Prepare cutoff.
     double rcut = 7;
-    std::vector<double> cutoff_hyps;
-    void (*cutoff_function)(double *, double, double, std::vector<double>) =
+    vector<double> cutoff_hyps;
+    void (*cutoff_function)(double *, double, double, vector<double>) =
         cos_cutoff;
 
     // Prepare spherical harmonics.
@@ -38,9 +40,9 @@ int main(){
     double first_gauss = 1;
     double final_gauss = 6;
     int N = 10;
-    std::vector<double> radial_hyps = {sigma, first_gauss, final_gauss};
+    vector<double> radial_hyps = {sigma, first_gauss, final_gauss};
     void (*basis_function)(double *, double *, double, int,
-        std::vector<double>) = equispaced_gaussians;
+        vector<double>) = equispaced_gaussians;
 
     // Time radial basis.
     double * g = new double[N];
@@ -48,20 +50,20 @@ int main(){
     double * gy = new double[N];
     double * gz = new double[N];
 
-    auto t1 = std::chrono::high_resolution_clock::now();
+    auto t1 = chrono::high_resolution_clock::now();
 
-    int reps = 10000;
+    int reps = 1000000;
     for (int n = 0; n < reps; n++){
-    calculate_radial(g, gx, gy, gz, basis_function, cutoff_function, 
+    calculate_radial(g, gx, gy, gz, basis_function, cutoff_function,
                      x, y, z, r, rcut, N, radial_hyps, cutoff_hyps);
     }
 
-    auto t2 = std::chrono::high_resolution_clock::now();
-    auto tot_time = 
-        std::chrono::duration_cast<std::chrono::microseconds>(t2-t1).count();
+    auto t2 = chrono::high_resolution_clock::now();
+    auto tot_time =
+        (double) chrono::duration_cast<chrono::microseconds>(t2-t1).count();
     auto mean_time = tot_time / reps;
-    std::cout << "calculating the radial basis took "
-              << mean_time << " microseconds\n";
+    cout << "calculating the radial basis took "
+              << mean_time << " microseconds" << endl;
 
     delete [] g; delete [] gx; delete [] gy; delete [] gz;
 
@@ -71,15 +73,15 @@ int main(){
     vector<double> hy = vector<double>(number_of_harmonics, 0);
     vector<double> hz = vector<double>(number_of_harmonics, 0);
 
-    t1 = std::chrono::high_resolution_clock::now();
+    t1 = chrono::high_resolution_clock::now();
     for (int n = 0; n < reps; n++){
     get_Y(h, hx, hy, hz, x, y, z, lmax);
     }
-    t2 = std::chrono::high_resolution_clock::now();
+    t2 = chrono::high_resolution_clock::now();
 
-    tot_time = 
-        std::chrono::duration_cast<std::chrono::microseconds>(t2-t1).count();
+    tot_time =
+        (double) chrono::duration_cast<chrono::microseconds>(t2-t1).count();
     mean_time = tot_time / reps;
-    std::cout << "calculating the spherical harmonics took "
-              << mean_time << " microseconds\n";
+    cout << "calculating the spherical harmonics took "
+              << mean_time << " microseconds" << endl;
 }


Test patch:
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 4411bc4c..f2e766a1 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -3,4 +3,4 @@ add_executable(tests test_main.cpp test_y_grad.cpp test_radial.cpp
                test_structure.cpp test_local_environment.cpp test_bond_env.cpp
 	       test_descriptor.cpp test_kernels.cpp test_sparse_gp.cpp)
 target_include_directories(tests PUBLIC ${ACE_INCLUDE_DIR})
-target_link_libraries(tests PUBLIC gtest gtest_main Ace)
+target_link_libraries(tests PUBLIC gtest gtest_main flare_pp)