From 66e4e01870b5816f056fc578bc75f9c78367a78c Mon Sep 17 00:00:00 2001
From: kra568 <kra568@login1.head.smc.salomon.it4i.cz>
Date: Sat, 12 Oct 2019 20:34:23 +0200
Subject: [PATCH] [ENH] improved the acsf example to take into account a wide
 array of cmd line parameters

---
 run_test.sh           |  22 ++-
 src/examples/acsf.cpp | 358 ++++++++++++++++++++----------------------
 2 files changed, 189 insertions(+), 191 deletions(-)

diff --git a/run_test.sh b/run_test.sh
index 0f3de9b9..f3d5a54e 100755
--- a/run_test.sh
+++ b/run_test.sh
@@ -3,11 +3,27 @@ export OMP_NUM_THREADS=1
 #./build_scripts/compile_salomon.sh
 source ./build_scripts/load_salomon_modules.inc
 cd ./build/examples
-test_type=0
-niteration=100
+
+niteration=10
 batch_size=0
 tolerance=1e-1
 nlayers=2
 layers_cardinality="5 1"
-mpirun -np 24 ./acsf $test_type $niteration $batch_size $tolerance $nlayers $layers_cardinality > dimer_0001.txt
+
+source_data_file="/home/kra568/lib4neuro/data/lj_all.xyz"
+target_net_file="/home/kra568/lib4neuro/nets/net_dimer_5_1.l4n"
+
+nsym_g2=1
+g2_cutoffs="13.23"
+g2_extensions="0.2"
+g2_shifts="0"
+
+nsym_g5=0
+g5_cutoffs=""
+g5_extensions=""
+g5_shifts=""
+g5_angles=""
+
+
+mpirun -np 24 ./acsf $niteration $batch_size $tolerance $source_data_file $target_net_file $nlayers $layers_cardinality $nsym_g2 $g2_cutoffs $g2_extensions $g2_shifts $nsym_g5 $g5_cutoffs $g5_extensions $g5_shifts $g5_angles > dimer_0001.txt
 cd ../.. 
diff --git a/src/examples/acsf.cpp b/src/examples/acsf.cpp
index e28c9a89..a212e05a 100644
--- a/src/examples/acsf.cpp
+++ b/src/examples/acsf.cpp
@@ -129,41 +129,58 @@ double optimize_via_NelderMead(l4n::NeuralNetwork& net, l4n::ErrorFunction& ef)
 
 }
 
-void dimer_test( unsigned int niters, unsigned int batch_size, double tol, std::vector<unsigned int> net_complx ){
+void dimer_test(
+    unsigned int niters,
+    unsigned int batch_size,
+    double tol,
+    std::vector<unsigned int> net_complx,
+    const char* src_file,
+    const char* tgt_net_file,
+    const std::vector<double> &g2_cutoff_coefficients,
+    const std::vector<double> &g2_extensions,
+    const std::vector<double> &g2_shifts,
+    const std::vector<double> &g5_cutoff_coefficients,
+    const std::vector<double> &g5_extensions,
+    const std::vector<bool> &g5_shiftsmax,
+    const std::vector<double> &g5_angles
+){
 
     /* Specify cutoff functions */
-    l4n::CutoffFunction2 cutoff2(13.23);
-
-    /* Specify symmetry functions */
-    l4n::G2 sym_f1(&cutoff2, 0.2, 0);
-    l4n::G2 sym_f2(&cutoff2, 0.2, 3);
-    l4n::G2 sym_f3(&cutoff2, 0.2, 6);
-    l4n::G2 sym_f4(&cutoff2, 0.2, 9);
-    l4n::G2 sym_f5(&cutoff2, 0.2, 12);
-    l4n::G2 sym_f6(&cutoff2, 0.2, 0.04);
-    l4n::G2 sym_f7(&cutoff2, 0.2, 0.04);
-
-    l4n::G5 sym_f8(&cutoff2, 0.05, 1, 0.15);
-    l4n::G5 sym_f9(&cutoff2, 0.05, 1, 0.55);
-    l4n::G5 sym_f10(&cutoff2, 0.05, 1, 0.85);
-
-    l4n::G5 sym_f11(&cutoff2, 1, -1, 0.9);
-    l4n::G5 sym_f12(&cutoff2, 1.1, -1, 0.9);
-    l4n::G5 sym_f13(&cutoff2, 1.2, -1, 0.9);
-    l4n::G5 sym_f14(&cutoff2, 1.3, -1, 0.9);
-    l4n::G5 sym_f15(&cutoff2, 1.4, -1, 0.9);
-    l4n::G5 sym_f16(&cutoff2, 1.5, -1, 0.9);
-    l4n::G5 sym_f17(&cutoff2, 1.6, -1, 0.9);
-    l4n::G5 sym_f18(&cutoff2, 1.7, -1, 0.9);
-
-    std::vector<l4n::SymmetryFunction*> helium_sym_funcs = {&sym_f1};
+    std::vector<l4n::CutoffFunction2> g2_cutofffunctions;
+    for( auto el: g2_cutoff_coefficients ){
+        g2_cutofffunctions.emplace_back( l4n::CutoffFunction2(el) );
+    }
+    std::vector<l4n::CutoffFunction2> g5_cutofffunctions;
+    for( auto el: g5_cutoff_coefficients ){
+        g5_cutofffunctions.emplace_back( l4n::CutoffFunction2(el) );
+    }
+
+    /* the used symmetry functions */
+    std::vector<l4n::SymmetryFunction*> helium_sym_funcs(g2_extensions.size() + g5_extensions.size() );
+    /* Specify G2 symmetry functions */
+    std::vector<l4n::G2> G2functions;
+    for( size_t i = 0; i < g2_extensions.size(); ++i ){
+        G2functions.emplace_back( l4n::G2(&g2_cutofffunctions[ i ], g2_extensions[ i ], g2_shifts[ i ] ) );
+        helium_sym_funcs[ i ] = &G2functions[ i ];
+    }
+
+    /* Specify G2 symmetry functions */
+    std::vector<l4n::G5> G5functions;
+    auto shift = G2functions.size();
+    for( size_t i = 0; i < g5_extensions.size(); ++i ){
+        G5functions.emplace_back( l4n::G5(&g5_cutofffunctions[ i ], g5_extensions[ i ], g5_shiftsmax[ i ], g5_angles[ i ] ) );
+        helium_sym_funcs[ i + shift ] = &G5functions[ i ];
+    }
+
+
+
 
     l4n::Element helium = l4n::Element("He",helium_sym_funcs);
     std::unordered_map<l4n::ELEMENT_SYMBOL, l4n::Element*> elements;
     elements[l4n::ELEMENT_SYMBOL::He] = &helium;
 
     /* Read data */
-    l4n::XYZReader reader("/home/kra568/lib4neuro/data/lj_all.xyz", true);
+    l4n::XYZReader reader(src_file, true);
     reader.read();
 
     if( l4n::mpi_rank == 0 ){
@@ -220,8 +237,10 @@ void dimer_test( unsigned int niters, unsigned int batch_size, double tol, std::
     MPI_Allreduce( MPI_IN_PLACE, &max_error_rel, 1, MPI_DOUBLE, MPI_MAX, l4n::mpi_active_comm );
 
     if( l4n::mpi_rank == 0 ){
+        std::cout << "*****************************************************" << std::endl;
         std::cout << "maximal absolute error: " << max_error_abs << std::endl;
         std::cout << "maximal relative error: " << 100 * max_error_rel << " %" << std::endl;
+        net.save_text(tgt_net_file);
     }
 
 }
@@ -230,196 +249,159 @@ void dimer_test( unsigned int niters, unsigned int batch_size, double tol, std::
 int main( int argc, char** argv ) {
     MPI_INIT
 
-    unsigned int testt = 0;
-
     unsigned int niters = 1000;
     unsigned int batch_size = 0;
     double tol = 1e-1;
-    unsigned int nlayers = 2;
+    char *src_file = "/mem.txt";
+    char *tgt_net_file = "/mem.txt";
+    unsigned int nlayers = 0;
     std::vector<unsigned int> net_complx = {2, 1};
 
+    size_t id = 1;
+    if( argc > id ){
+        niters = atoi( argv[id] );
+    }
+    id++;
 
-    if( argc > 1 ){
-        testt = atoi( argv[1] );
+    if( argc > id ){
+        batch_size = atoi( argv[id] );
     }
-    if( argc > 2 ){
-        niters = atoi( argv[2] );
+    id++;
+
+    if( argc > id ){
+        tol = atof( argv[id] );
     }
-    if( argc > 3 ){
-        batch_size = atoi( argv[3] );
+    id++;
+
+    if( argc > id ){
+        src_file = argv[id];
     }
-    if( argc > 4 ){
-        tol = atof( argv[4] );
+    id++;
+
+    if( argc > id ){
+        tgt_net_file = argv[id];
     }
-    if( argc > 5 ){
-        nlayers = atoi( argv[5] );
+    id++;
+
+    if( argc > id ){
+        nlayers = atoi( argv[id] );
         net_complx.resize(nlayers);
     }
-    for( unsigned int lidx = 6; lidx < 6 + nlayers; ++lidx ){
-        net_complx[lidx - 6] = atoi( argv[lidx] );
+    id++;
+
+    for( unsigned int lidx = id; lidx < id + nlayers; ++lidx ){
+        net_complx[lidx - id] = atoi( argv[lidx] );
+    }
+    id += nlayers;
+
+    std::vector<double> g2_cutoff_coefficients = {};
+    std::vector<double> g2_extensions = {};
+    std::vector<double> g2_shifts = {};
+    unsigned int ng2 = 0;
+    if( argc > id ){
+        ng2 = atoi( argv[id] );
+        g2_extensions.resize(ng2);
+        g2_shifts.resize(ng2);
+        g2_cutoff_coefficients.resize(ng2);
     }
+    id++;
 
+    for( unsigned int sidx = id; sidx < id + ng2; ++sidx ){
+        g2_cutoff_coefficients[sidx - id] = atof( argv[sidx] );
+    }
+    id += ng2;
+    for( unsigned int sidx = id; sidx < id + ng2; ++sidx ){
+        g2_extensions[sidx - id] = atof( argv[sidx] );
+    }
+    id += ng2;
+    for( unsigned int sidx = id; sidx < id + ng2; ++sidx ){
+        g2_shifts[sidx - id] = atof( argv[sidx] );
+    }
+    id += ng2;
+
+    std::vector<double> g5_cutoff_coefficients = {};
+    std::vector<double> g5_extensions = {};
+    std::vector<bool> g5_shifts = {};
+    std::vector<double> g5_angles = {};
+    unsigned int ng5 = 0;
+    if( argc > id ){
+        ng5 = atoi( argv[id] );
+        g5_extensions.resize(ng5);
+        g5_shifts.resize(ng5);
+        g5_angles.resize(ng5);
+        g5_cutoff_coefficients.resize(ng5);
+    }
+    id++;
+
+    for( unsigned int sidx = id; sidx < id + ng5; ++sidx ){
+        g5_cutoff_coefficients[sidx - id] = atof( argv[sidx] );
+    }
+    id += ng5;
+    for( unsigned int sidx = id; sidx < id + ng5; ++sidx ){
+        g5_extensions[sidx - id] = atof( argv[sidx] );
+    }
+    id += ng5;
+    for( unsigned int sidx = id; sidx < id + ng5; ++sidx ){
+        g5_shifts[sidx - id] = (bool)atoi( argv[sidx] );
+    }
+    id += ng5;
+    for( unsigned int sidx = id; sidx < id + ng5; ++sidx ){
+        g5_angles[sidx - id] = atof( argv[sidx] );
+    }
+    id += ng5;
+
+    nlayers = net_complx.size();
+    ng2 = g2_extensions.size();
+    ng5 = g5_extensions.size();
 
     if( l4n::mpi_rank == 0 ){
-        std::cout << "Test type: " << testt << std::endl;
+
         std::cout << "Maximal number of iterations: " << niters << std::endl;
         std::cout << "Batch size: " << batch_size << std::endl;
         std::cout << "Learning tolerance: " << tol << std::endl;
+
         std::cout << "# of network layers: " << nlayers << std::endl;
         std::cout << "# of neurons in the layers: ";
         for( auto el: net_complx ){
             std::cout << el << " ";
         }
         std::cout << std::endl;
+        if( ng2 > 0 ){
+            std::cout << "***********************************************************" << std::endl;
+            std::cout << "number of G2 symmetry functions: " << ng2 << std::endl;
+            for( size_t i = 0; i < ng2; ++i ){
+                std::cout << " f" << i << ": cutoff " << g2_cutoff_coefficients[ i ] << ", extension " << g2_extensions[ i ] << ", shift " << g2_shifts[ i ] << std::endl;
+            }
+        }
+        if( ng5 > 0 ){
+            std::cout << "***********************************************************" << std::endl;
+            std::cout << "number of G5 symmetry functions: " << ng5 << std::endl;
+            for( size_t i = 0; i < ng5; ++i ){
+                std::cout << " f" << i << ": cutoff " << g5_cutoff_coefficients[ i ] << ", extension " << g5_extensions[ i ] << ", shift " << g5_shifts[ i ] << ", angle " << g5_angles[ i ] << std::endl;
+            }
+        }
+        std::cout << "***********************************************************" << std::endl;
+        std::cout << "Data file: " << src_file << std::endl;
+        std::cout << "Network will be saved in: " << tgt_net_file << std::endl;
         std::cout << "***********************************************************" << std::endl << std::endl;
     }
     try{
 
-        if( testt == 0 ){
-            dimer_test( niters, batch_size, tol, net_complx );
-        }
-
-//        /* Specify cutoff functions */
-//        l4n::CutoffFunction2 cutoff2(13.23);
-//
-//        /* Specify symmetry functions */
-//        l4n::G2 sym_f1(&cutoff2, 0.2, 0);
-//        l4n::G2 sym_f2(&cutoff2, 0.2, 3);
-//        l4n::G2 sym_f3(&cutoff2, 0.2, 6);
-//        l4n::G2 sym_f4(&cutoff2, 0.2, 9);
-//        l4n::G2 sym_f5(&cutoff2, 0.2, 12);
-//        l4n::G2 sym_f6(&cutoff2, 0.2, 0.04);
-//        l4n::G2 sym_f7(&cutoff2, 0.2, 0.04);
-//
-//        l4n::G5 sym_f8(&cutoff2, 0.05, 1, 0.15);
-//        l4n::G5 sym_f9(&cutoff2, 0.05, 1, 0.55);
-//        l4n::G5 sym_f10(&cutoff2, 0.05, 1, 0.85);
-//
-//        l4n::G5 sym_f11(&cutoff2, 1, -1, 0.9);
-//        l4n::G5 sym_f12(&cutoff2, 1.1, -1, 0.9);
-//        l4n::G5 sym_f13(&cutoff2, 1.2, -1, 0.9);
-//        l4n::G5 sym_f14(&cutoff2, 1.3, -1, 0.9);
-//        l4n::G5 sym_f15(&cutoff2, 1.4, -1, 0.9);
-//        l4n::G5 sym_f16(&cutoff2, 1.5, -1, 0.9);
-//        l4n::G5 sym_f17(&cutoff2, 1.6, -1, 0.9);
-//        l4n::G5 sym_f18(&cutoff2, 1.7, -1, 0.9);
-//
-//        std::vector<l4n::SymmetryFunction*> helium_sym_funcs = {&sym_f1
-//                                                                // &sym_f2
-//                                                                // &sym_f3,
-//                                                                // &sym_f4,
-//                                                                // &sym_f5,
-////                                                                &sym_f6,
-////                                                                &sym_f7,
-//                                                                // &sym_f8,
-//                                                                // &sym_f9,
-//                                                                // &sym_f10
-////                                                                &sym_f11,
-////                                                                &sym_f12};
-////                                                                &sym_f13,
-////                                                                &sym_f14,
-////                                                                &sym_f15,
-////                                                                &sym_f16,
-////                                                                &sym_f17,
-////                                                                &sym_f18
-//    };
-//
-//        l4n::Element helium = l4n::Element("He",
-//                                           helium_sym_funcs);
-//        std::unordered_map<l4n::ELEMENT_SYMBOL, l4n::Element*> elements;
-//        elements[l4n::ELEMENT_SYMBOL::He] = &helium;
-//
-//        /* Read data */
-//        // l4n::XYZReader reader("/home/kra568/lib4neuro/data/HE4+T0.xyz", true);
-////        l4n::XYZReader reader("/home/kra568/lib4neuro/data/HE21+T0.xyz", true);
-//       l4n::XYZReader reader("/home/kra568/lib4neuro/data/lj_all.xyz", true);
-//        reader.read();
-//
-//		if( l4n::mpi_rank == 0 ){
-//			std::cout << "Finished reading data" << std::endl;
-//		}
-//
-//        std::shared_ptr<l4n::DataSet> ds = reader.get_acsf_data_set(elements);
-////        std::cout << "ahoj" << std::endl;
-//
-//        /* Create a neural network */
-//        std::unordered_map<l4n::ELEMENT_SYMBOL, std::vector<unsigned int>> n_hidden_neurons;
-//        n_hidden_neurons[l4n::ELEMENT_SYMBOL::He] = {20, 20, 1};
-//
-//        std::unordered_map<l4n::ELEMENT_SYMBOL, std::vector<l4n::NEURON_TYPE>> type_hidden_neurons;
-//        type_hidden_neurons[l4n::ELEMENT_SYMBOL::He] = {l4n::NEURON_TYPE::LOGISTIC, l4n::NEURON_TYPE::LOGISTIC, l4n::NEURON_TYPE::LINEAR};
-//
-//        l4n::ACSFNeuralNetwork net(elements, *reader.get_element_list(), reader.contains_charge(), n_hidden_neurons, type_hidden_neurons);
-//
-//        l4n::MSE mse(&net, ds.get(), false);
-//
-//        net.randomize_parameters();
-//
-////        for(size_t i = 0; i < ds->get_data()->at(0).first.size(); i++) {
-////            std::cout << ds->get_data()->at(0).first.at(i) << " ";
-////            if(i % 2 == 1) {
-////                std::cout << std::endl;
-////            }
-////        }
-////        std::cout << "----" << std::endl;
-//
-//       // l4n::ACSFParametersOptimizer param_optim(&mse, &reader);
-//       // std::vector<l4n::SYMMETRY_FUNCTION_PARAMETER> fitted_params = {l4n::SYMMETRY_FUNCTION_PARAMETER::EXTENSION,
-//                                                                      // l4n::SYMMETRY_FUNCTION_PARAMETER::SHIFT_MAX,
-//                                                                      // l4n::SYMMETRY_FUNCTION_PARAMETER::SHIFT,
-//                                                                      // l4n::SYMMETRY_FUNCTION_PARAMETER::ANGULAR_RESOLUTION};
-//
-//		// for( auto el:helium_sym_funcs ){
-//			// el->print_parameters();
-//		// }
-//       // param_optim.fit_ACSF_parameters(fitted_params, false);
-//		// for( auto el:helium_sym_funcs ){
-//			// el->print_parameters();
-//		// }
-////
-////        for(size_t i = 0; i < mse.get_dataset()->get_data()->at(0).first.size(); i++) {
-////            std::cout << mse.get_dataset()->get_data()->at(0).first.at(i) << " ";
-////            if(i % 2 == 1) {
-////                std::cout << std::endl;
-////            }
-////        }
-////        std::cout << "----" << std::endl;
-//
-//
-////        optimize_via_particle_swarm(net, mse);
-////        optimize_via_NelderMead(net, mse);
-//
-//        double err1 = optimize_via_LBMQ(net, mse);
-//        double err2 = optimize_via_gradient_descent(net, mse);
-//
-//        /* Print fit comparison with real data */
-//        std::vector<double> output;
-//        output.resize(1);
-//
-//		double max_error_abs = 0.0, max_error_rel = 0.0;
-//
-//        for(auto e : *mse.get_dataset()->get_data()) {
-//			net.eval_single(e.first, output);
-//			double error_= std::abs(e.second.at(0) - output.at(0));
-//			double error_rel = error_ / std::abs(output.at(0));
-//
-//			if( max_error_abs < error_ ){
-//				max_error_abs = error_;
-//			}
-//
-//			if( error_rel > max_error_rel ){
-//				max_error_rel = error_rel;
-//			}
-//        }
-//
-//		MPI_Allreduce( MPI_IN_PLACE, &max_error_abs, 1, MPI_DOUBLE, MPI_MAX, l4n::mpi_active_comm );
-//		MPI_Allreduce( MPI_IN_PLACE, &max_error_rel, 1, MPI_DOUBLE, MPI_MAX, l4n::mpi_active_comm );
-//
-//		if( l4n::mpi_rank == 0 ){
-//			std::cout << "maximal absolute error: " << max_error_abs << std::endl;
-//			std::cout << "maximal relative error: " << 100 * max_error_rel << " %" << std::endl;
-//		}
-
+        dimer_test(
+            niters,
+            batch_size,
+            tol,
+            net_complx,
+            src_file,
+            tgt_net_file,
+            g2_cutoff_coefficients,
+            g2_extensions,
+            g2_shifts,
+            g2_cutoff_coefficients,
+            g5_extensions,
+            g5_shifts,
+            g5_angles
+        );
 
     } catch (const std::exception& e) {
         std::cerr << e.what() << std::endl;
-- 
GitLab