diff --git a/include/4neuro.h b/include/4neuro.h index 2ccc19415ff0a5eb8910422951ce84c14994c881..c963a773017eb617fbbc8a1074e8ae38ae973615 100644 --- a/include/4neuro.h +++ b/include/4neuro.h @@ -33,6 +33,7 @@ #include "../src/constants.h" #include "../src/Coordinates/coordinates.h" #include "../src/LearningMethods/NelderMead.h" +#include "../src/SymmetryFunction/ACSFParametersOptimizer.h" // Abbreviate lib4neuro namespace to l4n namespace l4n = lib4neuro; diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 575da7ef454ffc396f5d4bc2eb1fd3c0a40d30e1..865ac0bd50b8c578574db06064b1e8c78847050b 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -120,6 +120,7 @@ IF("${BUILD_LIB}" STREQUAL "yes") Reader/Reader.cpp Coordinates/coordinates.cpp LearningMethods/NelderMead.cpp + SymmetryFunction/ACSFParametersOptimizer.cpp ) # Detect Threading library diff --git a/src/ErrorFunction/ErrorFunctions.cpp b/src/ErrorFunction/ErrorFunctions.cpp index 45139360e8460a2ccd5e6b35b8a5885fb5494589..7e67dc32d87bbf58bc7dbbcc1ccac1f5a0fab293 100644 --- a/src/ErrorFunction/ErrorFunctions.cpp +++ b/src/ErrorFunction/ErrorFunctions.cpp @@ -15,6 +15,14 @@ namespace lib4neuro { return this->dimension; } + std::vector<NeuralNetwork*>& ErrorFunction::get_nets() { + return nets; + } + + DataSet* ErrorFunction::get_dataset() const { + return ds; + } + void MSE::divide_data_train_test(double percent_test) { size_t ds_size = this->ds->get_n_elements(); diff --git a/src/ErrorFunction/ErrorFunctions.h b/src/ErrorFunction/ErrorFunctions.h index f15e10d2cfbfbfe0b282526a7ee5854b10a271bd..d3aa5126d749e40f89017b0c19503ea8d5a42ca7 100644 --- a/src/ErrorFunction/ErrorFunctions.h +++ b/src/ErrorFunction/ErrorFunctions.h @@ -232,6 +232,10 @@ namespace lib4neuro { */ virtual void randomize_parameters(double scaling) = 0; + [[nodiscard]] std::vector<NeuralNetwork*>& get_nets(); + + [[nodiscard]] DataSet* get_dataset() const; + protected: /** diff --git a/src/LearningMethods/NelderMead.cpp b/src/LearningMethods/NelderMead.cpp index 15db89df126ce00a76133dc41ae9393860e7e3c9..8c43c2a85d21a4d46e352653c1a0573d243dcec0 100644 --- a/src/LearningMethods/NelderMead.cpp +++ b/src/LearningMethods/NelderMead.cpp @@ -141,7 +141,7 @@ void lib4neuro::NelderMead::optimize(lib4neuro::ErrorFunction& ef, ef.set_parameters(this->optimal_parameters); /* Termination condition */ - + //TODO /* Compute centroid */ arma::Col<double> centroid(problem_dim); diff --git a/src/Network/ACSFNeuralNetwork.cpp b/src/Network/ACSFNeuralNetwork.cpp index c2231666ceeb97bd0dd41ddbc0760d5b9ef8235e..9ab44255af60f1801a68aa2db5ad165c5ebc9f2b 100644 --- a/src/Network/ACSFNeuralNetwork.cpp +++ b/src/Network/ACSFNeuralNetwork.cpp @@ -7,6 +7,7 @@ #include "../exceptions.h" #include "../settings.h" #include "ACSFNeuralNetwork.h" +#include "../ErrorFunction/ErrorFunctions.h" lib4neuro::ACSFNeuralNetwork::ACSFNeuralNetwork(std::unordered_map<ELEMENT_SYMBOL, Element*>& elements, std::vector<ELEMENT_SYMBOL>& elements_list, @@ -24,6 +25,10 @@ lib4neuro::ACSFNeuralNetwork::ACSFNeuralNetwork(std::unordered_map<ELEMENT_SYMBO } } + /* Save info about elements */ + this->elements = &elements; + this->elements_list = &elements_list; + /* Construct the neural network */ std::vector<size_t> inputs; @@ -46,7 +51,7 @@ lib4neuro::ACSFNeuralNetwork::ACSFNeuralNetwork(std::unordered_map<ELEMENT_SYMBO /* Create input neurons for sub-net */ std::shared_ptr<NeuronLinear> inp_n; - for(size_t j = 0; j < elements[elements_list.at(i)]->getSymmetryFunctions().size(); j++) { + for(size_t j = 0; j < elements[elements_list.at(i)]->getSymmetryFunctions()->size(); j++) { inp_n = std::make_shared<NeuronLinear>(); last_neuron_idx = this->add_neuron(inp_n, BIAS_TYPE::NO_BIAS); previous_layer.emplace_back(last_neuron_idx); @@ -136,52 +141,21 @@ lib4neuro::ACSFNeuralNetwork::ACSFNeuralNetwork(std::unordered_map<ELEMENT_SYMBO new_layer.clear(); } - /* Create hidden layers in sub-net */ -// std::vector<unsigned int> n_neurons = n_hidden_neurons[elements_list.at(i)]; -// std::vector<NEURON_TYPE> types = type_hidden_neurons[elements_list.at(i)]; -// for(size_t j = 0; j < n_neurons.size(); j++) { /* Iterate over hidden layers */ -// /* Create hidden neurons */ -// for(size_t k = 0; k < n_neurons.at(j); k++) { -// std::shared_ptr<Neuron> hid_n; -// switch(types.at(j)) { -// case NEURON_TYPE::LOGISTIC: { -// hid_n = std::make_shared<NeuronLogistic>(); -// break; -// } -// case NEURON_TYPE::BINARY: { -// hid_n = std::make_shared<NeuronBinary>(); -// break; -// } -// case NEURON_TYPE::CONSTANT: { -// hid_n = std::make_shared<NeuronConstant>(); -// break; -// } -// case NEURON_TYPE::LINEAR: { -// hid_n = std::make_shared<NeuronLinear>(); -// break; -// } -// } -// -// neuron_idx = this->add_neuron(hid_n, BIAS_TYPE::NEXT_BIAS); -// new_layer.emplace_back(neuron_idx); -// -// /* Connect hidden neuron to the previous layer */ -// for(auto prev_n : previous_layer) { -// this->add_connection_simple(prev_n, neuron_idx, SIMPLE_CONNECTION_TYPE::NEXT_WEIGHT); -// } -// previous_layer = new_layer; -// new_layer.clear(); -// } -// } - /* Create output neurons for sub-net */ for(auto prev_n : previous_layer) { this->add_connection_constant(prev_n, outputs[ 0 ], 1.0); } - } /* Specify network inputs and outputs */ this->specify_input_neurons(inputs); this->specify_output_neurons(outputs); } + +std::unordered_map<lib4neuro::ELEMENT_SYMBOL, lib4neuro::Element*>* lib4neuro::ACSFNeuralNetwork::get_elements() { + return this->elements; +} + +std::vector<lib4neuro::ELEMENT_SYMBOL>* lib4neuro::ACSFNeuralNetwork::get_elements_list() { + return this->elements_list; +} diff --git a/src/Network/ACSFNeuralNetwork.h b/src/Network/ACSFNeuralNetwork.h index 2d6fa571d5924a089530b0dd35a1b2beab05cbf8..071c58fb2b867ef6ac0942ac14b9daaac43814d3 100644 --- a/src/Network/ACSFNeuralNetwork.h +++ b/src/Network/ACSFNeuralNetwork.h @@ -9,15 +9,25 @@ #include "NeuralNetwork.h" #include "../DataSet/DataSet.h" +#include "../SymmetryFunction/SymmetryFunction.h" namespace lib4neuro { class ACSFNeuralNetwork : public NeuralNetwork { + private: + std::unordered_map<ELEMENT_SYMBOL, Element*>* elements; + + std::vector<ELEMENT_SYMBOL>* elements_list; + public: LIB4NEURO_API explicit ACSFNeuralNetwork(std::unordered_map<ELEMENT_SYMBOL, Element*>& elements, std::vector<ELEMENT_SYMBOL>& elements_list, bool with_charge, std::unordered_map<ELEMENT_SYMBOL, std::vector<unsigned int >> n_hidden_neurons, std::unordered_map<ELEMENT_SYMBOL, std::vector<NEURON_TYPE>> type_hidden_neurons); + + LIB4NEURO_API std::unordered_map<ELEMENT_SYMBOL, Element*>* get_elements(); + + LIB4NEURO_API std::vector<ELEMENT_SYMBOL>* get_elements_list(); }; } diff --git a/src/Network/NeuralNetwork.cpp b/src/Network/NeuralNetwork.cpp index a492ddeb05c70f0eb209e2fbf2bba39eabe7f151..338b1595724f4a2fd2566c90a40b982ada589ffb 100644 --- a/src/Network/NeuralNetwork.cpp +++ b/src/Network/NeuralNetwork.cpp @@ -271,7 +271,7 @@ namespace lib4neuro { double potential, bias; int bias_idx; - this->copy_parameter_space(custom_weights_and_biases); + this->copy_parameter_space(custom_weights_and_biases); // TODO rewrite, so the original parameters are not edited! this->analyze_layer_structure(); @@ -292,6 +292,7 @@ namespace lib4neuro { for (auto layer: this->neuron_layers_feedforward) { /* we iterate through all neurons in this layer and propagate the signal to the neighboring neurons */ +//#pragma omp parallel for collapse(3) for (auto si: *layer) { bias = 0.0; bias_idx = this->neuron_bias_indices.at(si); @@ -1057,5 +1058,13 @@ namespace lib4neuro { return std::make_pair(*std::min_element(this->connection_weights.begin(), this->connection_weights.end()), *std::max_element(this->connection_weights.begin(), this->connection_weights.end())); } + + size_t NeuralNetwork::get_input_neurons_number() { + return this->input_neuron_indices.size(); + } + + size_t NeuralNetwork::get_output_neurons_number() { + return this->output_neuron_indices.size(); + } } diff --git a/src/Network/NeuralNetwork.h b/src/Network/NeuralNetwork.h index f6c1b157a1c8605b17679e185df1a5b30c883194..303e5f67affdc7f57b298e6016e83602eebed172 100644 --- a/src/Network/NeuralNetwork.h +++ b/src/Network/NeuralNetwork.h @@ -494,6 +494,10 @@ namespace lib4neuro { //TODO WHY IS THIS HERE? LIB4NEURO_API void set_normalization_strategy_instance(NormalizationStrategy* ns); + LIB4NEURO_API size_t get_input_neurons_number(); + + LIB4NEURO_API size_t get_output_neurons_number(); + }; // class NeuralNetwork class FullyConnectedFFN : public NeuralNetwork { diff --git a/src/Reader/XYZReader.cpp b/src/Reader/XYZReader.cpp index aa68f95bba8f2d9454ad3dc438d285ce7f75b65d..9214f415b3ed794931090ffb6693edded928442c 100644 --- a/src/Reader/XYZReader.cpp +++ b/src/Reader/XYZReader.cpp @@ -150,7 +150,7 @@ void lib4neuro::XYZReader::transform_input_to_acsf(std::unordered_map<ELEMENT_SY double coord_val; for(size_t i = 0; i < particles.size(); i++) { /* Iterate over all the particles */ single_coords_check.clear(); - for (auto sym_func : element_description[particles.at(i).first]->getSymmetryFunctions()) { + for (auto sym_func : *element_description[particles.at(i).first]->getSymmetryFunctions()) { coord_val = sym_func->eval(i, particles); acsf_coords.emplace_back(coord_val); single_coords_check.emplace_back(coord_val); @@ -180,9 +180,9 @@ void lib4neuro::XYZReader::transform_input_to_acsf(std::unordered_map<ELEMENT_SY THROW_RUNTIME_ERROR("Not all descriptors are unique with currently specified symmetry functions!"); } - if(this->acsf_data_set == nullptr) { +// if(this->acsf_data_set == nullptr) { this->acsf_data_set = std::make_shared<DataSet>(); - } +// } this->acsf_data_set->add_data_pair(acsf_coords, configuration.second); acsf_coords.clear(); @@ -208,9 +208,9 @@ lib4neuro::ELEMENT_SYMBOL lib4neuro::XYZReader::get_element_symbol(std::string s std::shared_ptr<lib4neuro::DataSet> lib4neuro::XYZReader::get_acsf_data_set(std::unordered_map<ELEMENT_SYMBOL, Element*>& element_description) { - if(acsf_data_set == nullptr) { +// if(acsf_data_set == nullptr) { this->transform_input_to_acsf(element_description); - } +// } return this->acsf_data_set; } diff --git a/src/SymmetryFunction/ACSFParametersOptimizer.cpp b/src/SymmetryFunction/ACSFParametersOptimizer.cpp new file mode 100644 index 0000000000000000000000000000000000000000..a83cbca9b7b87753b2e9cc97ee5a37f48e214154 --- /dev/null +++ b/src/SymmetryFunction/ACSFParametersOptimizer.cpp @@ -0,0 +1,148 @@ +// +// Created by martin on 10.09.19. +// + +#include <algorithm> +#include <random> +#include <unordered_map> + +#include "Reader/XYZReader.h" +#include "../exceptions.h" +#include "ACSFParametersOptimizer.h" +#include "../Network/ACSFNeuralNetwork.h" +#include "../LearningMethods/GradientDescentBB.h" + +lib4neuro::ACSFParametersOptimizer::ACSFParametersOptimizer(lib4neuro::ErrorFunction* ef, XYZReader* reader) { + this->error_func = ef; + + if(this->error_func->get_nets().size() > 1) { + THROW_LOGIC_ERROR("Only single network can be passed as an argument!"); + } + + this->net = dynamic_cast<ACSFNeuralNetwork*>(this->error_func->get_nets().at(0)); + this->data_set = this->error_func->get_dataset(); + this->reader = reader; +} + +void lib4neuro::ACSFParametersOptimizer::fit_ACSF_parameters(std::vector<SYMMETRY_FUNCTION_PARAMETER>& fitted_params, + bool random_init_params) { + + /* Seed the random numbers' generator */ + std::random_device rd{}; + std::mt19937 gen{rd()}; + + DataSet* ds; + + GradientDescentBB gd(1e-6, + 1000, + 10000); + + /* Construct and train a new linear neural network to evaluate the chosen symmetry functions' parameters */ + std::vector<unsigned int> neuron_numbers = {static_cast<unsigned int>(this->net->get_input_neurons_number()), 1}; + std::vector<NEURON_TYPE> type_hidden_neurons = {}; + FullyConnectedFFN nn(&neuron_numbers, &type_hidden_neurons); + + /* Initial parameters' values */ + std::unordered_map<SYMMETRY_FUNCTION_PARAMETER, double> init_param_values; + std::normal_distribution<> d{1, 10}; // TODO make probability distribution choice more flexible by adding function parameters + std::uniform_real_distribution<> d_uni(0, 2); + std::vector<short int> shift_max_vals = {-1, 1}; + if(random_init_params) { + /* Iterate over element's description */ + for(auto el : *this->net->get_elements()) { + for(auto func : *el.second->getSymmetryFunctions()) { + for (auto param : *func->get_parameters()) { + /* Check, if the parameter is supposed to be optimized */ + if (std::find(fitted_params.begin(), + fitted_params.end(), + param.first) != fitted_params.end()) { + if (param.first == SYMMETRY_FUNCTION_PARAMETER::SHIFT_MAX) { + func->set_parameter(param.first, + shift_max_vals.at((int) (d_uni(gen)))); + } else { + func->set_parameter(param.first, + d(gen)); + } + } + } + } + } + } + + double previous_error; + double current_error; + double delta_error; + + double temperature = 3e-4; //TODO make flexible by adding function parameters + + ds = this->reader->get_acsf_data_set(*this->net->get_elements()).get(); + MSE mse_init(&nn, ds); + previous_error = mse_init.eval(); + + // TODO what is the termination condition? + /* Iteratively improve the parameters */ + for(unsigned int i = 0; i < 10; i++) { + /* Temperature iterations */ + for(unsigned int j = 0; j < 10; j++) { + /* Change the selected parameters in ACSFs */ + for (auto el : *this->net->get_elements()) { + for (auto func : *el.second->getSymmetryFunctions()) { + for (auto param : *func->get_parameters()) { + /* Check, if the parameter is supposed to be optimized */ + if (std::find(fitted_params.begin(), + fitted_params.end(), + param.first) != fitted_params.end()) { + if (param.first == SYMMETRY_FUNCTION_PARAMETER::SHIFT_MAX) { + func->set_parameter(param.first, + shift_max_vals.at((int) (d_uni(gen)))); + } else { + std::normal_distribution<> dd{param.second, 10}; + func->set_parameter(param.first, + dd(gen)); + } + } + } + } + } + + /* Re-compute coordinates for the new set of ACSFs */ + try { + ds = this->reader->get_acsf_data_set(*this->net->get_elements()).get(); + } catch(const std::runtime_error& err) { + continue; + } + + /* Train the network */ + MSE mse(&nn, + ds); + + gd.optimize(mse); + + current_error = mse.eval(); + delta_error = current_error - previous_error; + + std::bernoulli_distribution bernoulli_d(std::exp(-delta_error / temperature)); + if (!(delta_error < 0 || (delta_error > 0 && (bool)bernoulli_d(gen)))) { + /* If conditions are not met, revert the parameter changes */ + for (auto el : *this->net->get_elements()) { + for (auto func : *el.second->getSymmetryFunctions()) { + for (auto param : *func->get_parameters()) { + /* Check, if the parameter is supposed to be optimized */ + if (std::find(fitted_params.begin(), + fitted_params.end(), + param.first) != fitted_params.end()) { + func->revert_parameter_change(param.first); + } + } + } + } + } + + previous_error = current_error; + + std::cout << "it: " << i*j+j << " Delta error: " << delta_error << std::endl; + } + + temperature *= 0.98; + } +} diff --git a/src/SymmetryFunction/ACSFParametersOptimizer.h b/src/SymmetryFunction/ACSFParametersOptimizer.h new file mode 100644 index 0000000000000000000000000000000000000000..c0e8bd4213dd575b3bf275f0d2c1b1f07316d4be --- /dev/null +++ b/src/SymmetryFunction/ACSFParametersOptimizer.h @@ -0,0 +1,32 @@ +// +// Created by martin on 10.09.19. +// + +#ifndef LIB4NEURO_ACSFPARAMETERSOPTIMIZER_H +#define LIB4NEURO_ACSFPARAMETERSOPTIMIZER_H + +#include "../settings.h" +#include "../ErrorFunction/ErrorFunctions.h" +#include "../SymmetryFunction/SymmetryFunction.h" +#include "../Network/ACSFNeuralNetwork.h" + +namespace lib4neuro { + class ACSFParametersOptimizer { + private: + ErrorFunction* error_func; + + ACSFNeuralNetwork* net; + + DataSet* data_set; + + XYZReader* reader; // TODO implement for Reader class in general + + public: + LIB4NEURO_API explicit ACSFParametersOptimizer(ErrorFunction* ef, XYZReader* reader); + + LIB4NEURO_API void fit_ACSF_parameters(std::vector<SYMMETRY_FUNCTION_PARAMETER>& fitted_params, + bool random_init_params); + }; +} + +#endif //LIB4NEURO_ACSFPARAMETERSOPTIMIZER_H diff --git a/src/SymmetryFunction/Element.cpp b/src/SymmetryFunction/Element.cpp index a9986ee89accc7ce4f48ddb3734275832bc52ede..1876ee679337f2ede845047090d7a0fb5625e8b5 100644 --- a/src/SymmetryFunction/Element.cpp +++ b/src/SymmetryFunction/Element.cpp @@ -15,8 +15,8 @@ const std::string& lib4neuro::Element::getElementSymbol() const { // return n_particles; //} -std::vector<lib4neuro::SymmetryFunction*>& lib4neuro::Element::getSymmetryFunctions() { - return symmetry_functions; +std::vector<lib4neuro::SymmetryFunction*>* lib4neuro::Element::getSymmetryFunctions() { + return &symmetry_functions; } //const std::vector<unsigned int>& lib4neuro::Element::getNHiddenNeurons() const { diff --git a/src/SymmetryFunction/Element.h b/src/SymmetryFunction/Element.h index da726901b631defee3e9f7f90e9e4e6fb84c92e1..3bb1d849767d9415f22f5b41fa133709e60caf96 100644 --- a/src/SymmetryFunction/Element.h +++ b/src/SymmetryFunction/Element.h @@ -29,7 +29,7 @@ namespace lib4neuro { // LIB4NEURO_API [[nodiscard]] unsigned int getNParticles() const; - LIB4NEURO_API [[nodiscard]] std::vector<SymmetryFunction*>& getSymmetryFunctions(); + LIB4NEURO_API [[nodiscard]] std::vector<SymmetryFunction*>* getSymmetryFunctions(); // LIB4NEURO_API [[nodiscard]] const std::vector<unsigned int>& getNHiddenNeurons() const; // diff --git a/src/SymmetryFunction/SymmetryFunction.cpp b/src/SymmetryFunction/SymmetryFunction.cpp index 0d7d74d80f65c5678048062ed8d3fcfc4f9a1b21..33e60cf29809737c77eb78d5ab4a603209c2dcd5 100644 --- a/src/SymmetryFunction/SymmetryFunction.cpp +++ b/src/SymmetryFunction/SymmetryFunction.cpp @@ -54,6 +54,25 @@ double lib4neuro::SymmetryFunction::euclid_distance(std::vector<double> v1, lib4neuro::SymmetryFunction::SymmetryFunction() {} +std::unordered_map<lib4neuro::SYMMETRY_FUNCTION_PARAMETER, double>* lib4neuro::SymmetryFunction::get_parameters() { + return &this->parameters; +} + +void lib4neuro::SymmetryFunction::set_parameter(lib4neuro::SYMMETRY_FUNCTION_PARAMETER parameter, double value) { + if(this->parameters.find(parameter) != this->parameters.end()) { + this->previous_parameters[parameter] = this->parameters.at(parameter); + } + this->parameters[parameter] = value; +} + +void lib4neuro::SymmetryFunction::revert_parameter_change(SYMMETRY_FUNCTION_PARAMETER parameter) { + this->parameters.at(parameter) = this->previous_parameters.at(parameter); +} + +//lib4neuro::SYMMETRY_FUNCTION_TYPE lib4neuro::SymmetryFunction::get_type() { +// return this->type; +//} + double lib4neuro::G1::eval(unsigned int particle_ind, std::vector<std::vector<double>> cartesian_coord) { double sum = 0; double distance; @@ -73,7 +92,10 @@ double lib4neuro::G1::eval(unsigned int particle_ind, std::vector<std::vector<do return sum; } -lib4neuro::G1::G1(lib4neuro::CutoffFunction* cutoff_func) : SymmetryFunction(cutoff_func) {} +lib4neuro::G1::G1(lib4neuro::CutoffFunction* cutoff_func) : SymmetryFunction(cutoff_func) { + this->parameters = {}; +// this->type = SYMMETRY_FUNCTION_TYPE::G1; +} double lib4neuro::G1::eval(unsigned int particle_ind, std::vector<std::pair<ELEMENT_SYMBOL, std::vector<double>>>& cartesian_coord) { @@ -107,8 +129,12 @@ lib4neuro::CutoffFunction2::CutoffFunction2() {} lib4neuro::G2::G2(lib4neuro::CutoffFunction* cutoff_func, double extension, double shift) : SymmetryFunction(cutoff_func) { - this->extension = extension; - this->shift = shift; +// this->extension = extension; +// this->shift = shift; +// this->parameters = {SYMMETRY_FUNCTION_PARAMETER::EXTENSION, SYMMETRY_FUNCTION_PARAMETER::SHIFT}; + this->parameters[SYMMETRY_FUNCTION_PARAMETER::EXTENSION]= extension; + this->parameters[SYMMETRY_FUNCTION_PARAMETER::SHIFT] = shift; +// this->type = SYMMETRY_FUNCTION_TYPE::G2; } double lib4neuro::G2::eval(unsigned int particle_ind, @@ -116,6 +142,9 @@ double lib4neuro::G2::eval(unsigned int particle_ind, double sum = 0; double distance; + double extension = this->parameters.at(SYMMETRY_FUNCTION_PARAMETER::EXTENSION); + double shift = this->parameters.at(SYMMETRY_FUNCTION_PARAMETER::SHIFT); + for(size_t i = 0; i < cartesian_coord.size(); i++) { if(i == particle_ind) { continue; @@ -125,7 +154,7 @@ double lib4neuro::G2::eval(unsigned int particle_ind, distance = this->euclid_distance(cartesian_coord.at(particle_ind).second, cartesian_coord.at(i).second); /* Call cutoff function */ - sum += std::exp(-this->extension * (distance - this->shift) * (distance - this->shift)) * this->cutoff_func->eval(distance); + sum += std::exp(-extension * (distance - shift) * (distance - shift)) * this->cutoff_func->eval(distance); } return sum; @@ -141,7 +170,9 @@ lib4neuro::G2::G2() {} lib4neuro::G3::G3(lib4neuro::CutoffFunction* cutoff_func, double period_length) : SymmetryFunction(cutoff_func) { - this->period_length = period_length; +// this->period_length = period_length; + this->parameters[SYMMETRY_FUNCTION_PARAMETER::PERIOD_LENGTH] = period_length; +// this->type = SYMMETRY_FUNCTION_TYPE::G3; } lib4neuro::G3::G3() {} @@ -151,6 +182,8 @@ double lib4neuro::G3::eval(unsigned int particle_ind, double sum = 0; double distance; + double period_length = this->parameters.at(SYMMETRY_FUNCTION_PARAMETER::PERIOD_LENGTH); + for(size_t i = 0; i < cartesian_coord.size(); i++) { if(i == particle_ind) { continue; @@ -160,7 +193,7 @@ double lib4neuro::G3::eval(unsigned int particle_ind, distance = this->euclid_distance(cartesian_coord.at(particle_ind).second, cartesian_coord.at(i).second); /* Call cutoff function */ - sum += cos(this->period_length*distance) * this->cutoff_func->eval(distance); + sum += cos(period_length*distance) * this->cutoff_func->eval(distance); } return sum; @@ -174,11 +207,21 @@ double lib4neuro::G3::eval(unsigned int particle_ind, lib4neuro::G4::G4(lib4neuro::CutoffFunction* cutoff_func, double extension, - bool shift_max, + int shift_max, double angular_resolution) : SymmetryFunction(cutoff_func) { - this->extension = extension; - this->shift_max = shift_max; - this->angular_resolution = angular_resolution; +// this->extension = extension; +// this->shift_max = shift_max; +// this->angular_resolution = angular_resolution; +// this->parameters = {SYMMETRY_FUNCTION_PARAMETER::EXTENSION, SYMMETRY_FUNCTION_PARAMETER::ANGULAR_RESOLUTION}; +// this->type = SYMMETRY_FUNCTION_TYPE::G4; + + if(shift_max != 1 && shift_max != -1) { + THROW_LOGIC_ERROR("shift_max parameter can be only 1 or -1!"); + } + + this->parameters[SYMMETRY_FUNCTION_PARAMETER::EXTENSION] = extension; + this->parameters[SYMMETRY_FUNCTION_PARAMETER::SHIFT_MAX] = shift_max; + this->parameters[SYMMETRY_FUNCTION_PARAMETER::ANGULAR_RESOLUTION] = angular_resolution; } lib4neuro::G4::G4() {} @@ -196,6 +239,10 @@ double lib4neuro::G4::eval(unsigned int particle_ind, double distance2; double distance3; + double angular_resolution = this->parameters.at(SYMMETRY_FUNCTION_PARAMETER::ANGULAR_RESOLUTION); + double shift_max = this->parameters.at(SYMMETRY_FUNCTION_PARAMETER::SHIFT_MAX); + double extension = this->parameters.at(SYMMETRY_FUNCTION_PARAMETER::EXTENSION); + for(size_t i = 0; i < cartesian_coord.size(); i++) { if(i == particle_ind) { continue; @@ -212,7 +259,7 @@ double lib4neuro::G4::eval(unsigned int particle_ind, distance3 = this->euclid_distance(cartesian_coord.at(i).second, cartesian_coord.at(j).second); /* Call cutoff function */ - int lambda = (this->shift_max) ? 1 : -1; + int lambda = shift_max; std::vector<double> v1; std::vector<double> v2; double scalar_product = 0; @@ -223,24 +270,30 @@ double lib4neuro::G4::eval(unsigned int particle_ind, } double theta = acos(scalar_product / (distance1 * distance2)); - sum += pow(1 + lambda * cos(theta), this->angular_resolution) - * exp(-this->extension * (distance1*distance1 + distance2*distance2 + distance3*distance3)) + sum += pow(1 + lambda * cos(theta), angular_resolution) + * exp(-extension * (distance1*distance1 + distance2*distance2 + distance3*distance3)) * this->cutoff_func->eval(distance1) * this->cutoff_func->eval(distance2) * this->cutoff_func->eval(distance3); } } - return pow(2, 1-this->angular_resolution)*sum; + return pow(2, 1-angular_resolution)*sum; } lib4neuro::G5::G5(lib4neuro::CutoffFunction* cutoff_func, double extension, bool shift_max, double angular_resolution) : SymmetryFunction(cutoff_func) { - this->extension = extension; - this->shift_max = shift_max; - this->angular_resolution = angular_resolution; + this->parameters[SYMMETRY_FUNCTION_PARAMETER::EXTENSION] = extension; + this->parameters[SYMMETRY_FUNCTION_PARAMETER::SHIFT_MAX] = shift_max; + this->parameters[SYMMETRY_FUNCTION_PARAMETER::ANGULAR_RESOLUTION] = angular_resolution; + +// this->extension = extension; +// this->shift_max = shift_max; +// this->angular_resolution = angular_resolution; +// this->parameters = {SYMMETRY_FUNCTION_PARAMETER::EXTENSION, SYMMETRY_FUNCTION_PARAMETER::ANGULAR_RESOLUTION}; +// this->type = SYMMETRY_FUNCTION_TYPE::G5; } lib4neuro::G5::G5() {} @@ -257,6 +310,10 @@ double lib4neuro::G5::eval(unsigned int particle_ind, double distance1; double distance2; + double angular_resolution = this->parameters.at(SYMMETRY_FUNCTION_PARAMETER::ANGULAR_RESOLUTION); + double shift_max = this->parameters.at(SYMMETRY_FUNCTION_PARAMETER::SHIFT_MAX); + double extension = this->parameters.at(SYMMETRY_FUNCTION_PARAMETER::EXTENSION); + for(size_t i = 0; i < cartesian_coord.size(); i++) { if(i == particle_ind) { continue; @@ -272,7 +329,7 @@ double lib4neuro::G5::eval(unsigned int particle_ind, distance2 = this->euclid_distance(cartesian_coord.at(particle_ind).second, cartesian_coord.at(j).second); /* Call cutoff function */ - int lambda = (this->shift_max) ? 1 : -1; + int lambda = (shift_max) ? 1 : -1; std::vector<double> v1; std::vector<double> v2; double scalar_product = 0; @@ -283,18 +340,12 @@ double lib4neuro::G5::eval(unsigned int particle_ind, } double theta = acos(scalar_product / (distance1 * distance2)); - sum += pow(1 + lambda * cos(theta), this->angular_resolution) - * exp(-this->extension * (distance1*distance1 + distance2*distance2)) + sum += pow(1 + lambda * cos(theta), angular_resolution) + * exp(-extension * (distance1*distance1 + distance2*distance2)) * this->cutoff_func->eval(distance1) * this->cutoff_func->eval(distance2); } } - return pow(2, 1-this->angular_resolution)*sum; + return pow(2, 1-angular_resolution)*sum; } - - - - - - diff --git a/src/SymmetryFunction/SymmetryFunction.h b/src/SymmetryFunction/SymmetryFunction.h index 8e267a53d860e447ebdb1912faecb33f44517c58..d756b32f4c86fca3a1c3f96e9fd4d7f3d9a7b8ce 100644 --- a/src/SymmetryFunction/SymmetryFunction.h +++ b/src/SymmetryFunction/SymmetryFunction.h @@ -5,6 +5,7 @@ #ifndef LIB4NEURO_SYMMETRYFUNCTION_H #define LIB4NEURO_SYMMETRYFUNCTION_H +#include <unordered_map> #include <vector> #include "../settings.h" @@ -17,6 +18,19 @@ namespace lib4neuro { Fl,Mc,Lv,Ts,Og }; + enum class SYMMETRY_FUNCTION_PARAMETER { + //radius, // TODO make available also fitting of cutoff functions + EXTENSION, + SHIFT, + ANGULAR_RESOLUTION, + SHIFT_MAX, + PERIOD_LENGTH + }; + +// enum class SYMMETRY_FUNCTION_TYPE { +// G1, G2, G3, G4, G5 +// }; + class CutoffFunction { public: LIB4NEURO_API explicit CutoffFunction(double radius); @@ -81,15 +95,27 @@ namespace lib4neuro { std::vector<std::pair<ELEMENT_SYMBOL, std::vector<double>>>& cartesian_coord) = 0; - struct access; + LIB4NEURO_API std::unordered_map<SYMMETRY_FUNCTION_PARAMETER, double>* get_parameters(); + + LIB4NEURO_API void set_parameter(SYMMETRY_FUNCTION_PARAMETER parameter, double value); -// LIB4NEURO_API unsigned int get_n_particles(); + LIB4NEURO_API void revert_parameter_change(SYMMETRY_FUNCTION_PARAMETER parameter); + +// LIB4NEURO_API SYMMETRY_FUNCTION_TYPE get_type(); + + struct access; protected: -// unsigned int n_particles; +// std::vector<SYMMETRY_FUNCTION_PARAMETER> parameters; CutoffFunction* cutoff_func; + std::unordered_map<SYMMETRY_FUNCTION_PARAMETER, double> parameters; + + std::unordered_map<SYMMETRY_FUNCTION_PARAMETER, double> previous_parameters; + +// SYMMETRY_FUNCTION_TYPE type; + double euclid_distance(std::vector<double> v1, std::vector<double> v2); }; @@ -128,9 +154,9 @@ namespace lib4neuro { struct access; - protected: - double extension; - double shift; +// protected: +// double extension; +// double shift; }; class G3 : public SymmetryFunction { @@ -150,15 +176,13 @@ namespace lib4neuro { struct access; - protected: - double period_length; }; class G4 : public SymmetryFunction { public: LIB4NEURO_API explicit G4(CutoffFunction* cutoff_func, double extension, - bool shift_max, + int shift_max, double angular_resolution); /** @@ -173,11 +197,6 @@ namespace lib4neuro { std::vector<double>>>& cartesian_coord) override; struct access; - - protected: - double extension; - bool shift_max; - double angular_resolution; }; class G5 : public SymmetryFunction { @@ -200,10 +219,6 @@ namespace lib4neuro { struct access; - protected: - double extension; - bool shift_max; - double angular_resolution; }; } diff --git a/src/SymmetryFunction/SymmetryFunctionSerialization.h b/src/SymmetryFunction/SymmetryFunctionSerialization.h index 89e4a36ca0cc926fe901858810831e007f0784fb..7879969d946a3270de342b83e749e0d14c16ad4f 100644 --- a/src/SymmetryFunction/SymmetryFunctionSerialization.h +++ b/src/SymmetryFunction/SymmetryFunctionSerialization.h @@ -46,8 +46,6 @@ struct lib4neuro::G2::access { lib4neuro::G2& sf, const unsigned int version) { ar & boost::serialization::base_object<lib4neuro::SymmetryFunction>(sf); - ar & sf.shift; - ar & sf.extension; } }; @@ -56,7 +54,6 @@ struct lib4neuro::CutoffFunction::access { static void serialize(Archive& ar, lib4neuro::CutoffFunction& f, const unsigned int version) { - ar & f.radius; } }; diff --git a/src/examples/dev_sandbox.cpp b/src/examples/dev_sandbox.cpp index 53fcb8433743483597b3fbd464bd72424b45d308..135f5b7c82b30af1e08cfd79e46ebd088d57279e 100644 --- a/src/examples/dev_sandbox.cpp +++ b/src/examples/dev_sandbox.cpp @@ -111,7 +111,7 @@ double optimize_via_LBMQ(l4n::NeuralNetwork& net, } double optimize_via_NelderMead(l4n::NeuralNetwork& net, l4n::ErrorFunction& ef) { - l4n::NelderMead nm(10e-4, 10e-7, 5, 50); + l4n::NelderMead nm(10e-4, 10e-7, 500, 150); nm.optimize(ef); net.copy_parameter_space(nm.get_parameters()); @@ -132,7 +132,7 @@ int main() { /* Specify cutoff functions */ // l4n::CutoffFunction1 cutoff1(10.1); - l4n::CutoffFunction2 cutoff1(8); + l4n::CutoffFunction2 cutoff2(8); //l4n::CutoffFunction2 cutoff2(25); // l4n::CutoffFunction2 cutoff2(15.2); // l4n::CutoffFunction2 cutoff4(10.3); @@ -140,12 +140,13 @@ int main() { // l4n::CutoffFunction2 cutoff6(11); /* Specify symmetry functions */ -// l4n::G1 sym_f1(&cutoff1); - l4n::G2 sym_f2(&cutoff1, 2.09, 0.8); - l4n::G2 sym_f3(&cutoff1, 0.01, 0.04); -// l4n::G2 sym_f4(&cutoff2, 0.02, 0.04); -// l4n::G2 sym_f5(&cutoff2, 2.09, 0.04); - + l4n::G2 sym_f1(&cutoff2, 1.9, 0.7); + l4n::G2 sym_f2(&cutoff2, 2.09, 0.8); + l4n::G2 sym_f3(&cutoff2, 0.01, 0.04); + l4n::G2 sym_f4(&cutoff2, 0.02, 0.04); + l4n::G2 sym_f5(&cutoff2, 2.09, 0.04); + l4n::G2 sym_f6(&cutoff2, 2.09, 0.04); + l4n::G2 sym_f7(&cutoff2, 2.09, 0.04); // l4n::G3 sym_f4(&cutoff4, 0.3); // l4n::G4 sym_f5(&cutoff5, 0.05, true, 0.05); @@ -153,6 +154,13 @@ int main() { // l4n::G4 sym_f7(&cutoff6, 0.5, true, 0.05); // l4n::G4 sym_f8(&cutoff6, 0.5, false, 0.05); + l4n::G5 sym_f8(&cutoff2, 2.1, -1, 0.9); + l4n::G5 sym_f9(&cutoff2, 2.1, -1, 0.9); + l4n::G5 sym_f10(&cutoff2, 2.1, -1, 0.9); + l4n::G5 sym_f11(&cutoff2, 2.1, -1, 0.9); + l4n::G5 sym_f12(&cutoff2, 2.1, -1, 0.9); + l4n::G5 sym_f13(&cutoff2, 2.1, -1, 0.9); + std::vector<l4n::SymmetryFunction*> helium_sym_funcs = {&sym_f2, &sym_f3}; //, &sym_f4, &sym_f5}; //, &sym_f6, &sym_f7, &sym_f8}; l4n::Element helium = l4n::Element("He", @@ -170,23 +178,28 @@ int main() { /* 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, 1}; + 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::LINEAR}; + 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()); net.randomize_parameters(); + + l4n::ACSFParametersOptimizer param_optim(&mse, &reader); + std::vector<l4n::SYMMETRY_FUNCTION_PARAMETER> fitted_params = {l4n::SYMMETRY_FUNCTION_PARAMETER::EXTENSION}; + param_optim.fit_ACSF_parameters(fitted_params, true); + // optimize_via_particle_swarm(net, mse); - optimize_via_NelderMead(net, mse); +// optimize_via_NelderMead(net, mse); // double err1 = optimize_via_LBMQ(net, mse); - double err2 = optimize_via_gradient_descent(net, mse); +// double err2 = optimize_via_gradient_descent(net, mse); std::cout << "Weights: " << net.get_min_max_weight().first << " " << net.get_min_max_weight().second << std::endl; @@ -194,17 +207,17 @@ int main() { std::vector<double> output; output.resize(1); - for(auto e : *ds->get_data()) { - for(unsigned int i = 0; i < e.first.size(); i++) { - std::cout << e.first.at(i) << " "; - if(i % 2 == 1) { - std::cout << std::endl; - } - } - std::cout << "OUTS (DS, predict): " << e.second.at(0) << " "; - net.eval_single(e.first, output); - std::cout << output.at(0) << std::endl; - } +// for(auto e : *ds->get_data()) { +// for(unsigned int i = 0; i < e.first.size(); i++) { +// std::cout << e.first.at(i) << " "; +// if(i % 2 == 1) { +// std::cout << std::endl; +// } +// } +// std::cout << "OUTS (DS, predict): " << e.second.at(0) << " "; +// net.eval_single(e.first, output); +// std::cout << output.at(0) << std::endl; +// } } catch (const std::exception& e) { std::cerr << e.what() << std::endl;