SymmetryFunction.cpp

//
// Created by martin on 02.08.19.
//
#define _USE_MATH_DEFINES

#include <cmath>

#include "../exceptions.h"
#include "../message.h"

#include "SymmetryFunction.h"
#include "SymmetryFunctionSerialization.h"

BOOST_CLASS_EXPORT_IMPLEMENT(lib4neuro::CutoffFunction)
BOOST_CLASS_EXPORT_IMPLEMENT(lib4neuro::CutoffFunction1)
BOOST_CLASS_EXPORT_IMPLEMENT(lib4neuro::CutoffFunction2)
BOOST_CLASS_EXPORT_IMPLEMENT(lib4neuro::SymmetryFunction)
BOOST_CLASS_EXPORT_IMPLEMENT(lib4neuro::G1)
BOOST_CLASS_EXPORT_IMPLEMENT(lib4neuro::G2)


lib4neuro::CutoffFunction::CutoffFunction(double radius) {
    this->radius = radius;
}

double lib4neuro::CutoffFunction::get_radius() {
    return this->radius;
}

lib4neuro::CutoffFunction::CutoffFunction() {}

double lib4neuro::CutoffFunction1::eval(double ri) {
    return (ri > this->radius) ? 0 : 0.5 * ( cos( (M_PI * ri) / this->radius) + 1 );
}

lib4neuro::CutoffFunction1::CutoffFunction1(double radius) : CutoffFunction(radius) {}

lib4neuro::CutoffFunction1::CutoffFunction1() {}

//unsigned int lib4neuro::SymmetryFunction::get_n_particles() {
//    return this->n_particles;
//}

lib4neuro::SymmetryFunction::SymmetryFunction(lib4neuro::CutoffFunction* cutoff_func) {
    this->cutoff_func = cutoff_func;
}

double lib4neuro::SymmetryFunction::euclid_distance(std::vector<double> v1,
                                                    std::vector<double> v2) {
    double distance = 0;
    for(size_t j = 0; j < v1.size(); j++) {
        distance += (v1.at(j)-v2.at(j))*(v1.at(j)-v2.at(j));
    }
    return std::sqrt(distance);
}

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;

    for(size_t i = 0; i < cartesian_coord.size(); i++) {
        if(i == particle_ind) {
            continue;
        }

        /* Compute Euclid distance of i-th particle with the others */
        distance = this->euclid_distance(cartesian_coord.at(particle_ind), cartesian_coord.at(i));

        /* Call cutoff function */
        sum += this->cutoff_func->eval(distance);
    }

    return sum;
}

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) {
    double sum = 0;
    double distance;

    for(size_t i = 0; i < cartesian_coord.size(); i++) {
        if(i == particle_ind) {
            continue;
        }

        /* Compute Euclid distance of i-th particle with the others */
        distance = this->euclid_distance(cartesian_coord.at(particle_ind).second, cartesian_coord.at(i).second);

        /* Call cutoff function */
        sum += this->cutoff_func->eval(distance);
    }

    return sum;
}

lib4neuro::G1::G1() {}

lib4neuro::CutoffFunction2::CutoffFunction2(double radius) : CutoffFunction(radius) {}

double lib4neuro::CutoffFunction2::eval(double ri) {
    return (ri > this->radius) ? 0 : std::pow(tanh(1 - ri/this->radius),3);
}

lib4neuro::CutoffFunction2::CutoffFunction2() {}

lib4neuro::G2::G2(lib4neuro::CutoffFunction* cutoff_func,
                  double extension, double shift) : SymmetryFunction(cutoff_func) {
//    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,
                           std::vector<std::pair<ELEMENT_SYMBOL, std::vector<double>>>& cartesian_coord) {
    double sum = 0;
    double distance;

    double extension = this->parameters.at(SYMMETRY_FUNCTION_PARAMETER::EXTENSION);
    double shift = this->parameters.at(SYMMETRY_FUNCTION_PARAMETER::SHIFT);

//    std::cout << "particle " << particle_ind << "(";
//    for( auto el: cartesian_coord.at(particle_ind).second ){
//        std::cout << el << " ";
//    }
//    std::cout << ")";
    for(size_t i = 0; i < cartesian_coord.size(); i++) {
        if(i == particle_ind) {
            continue;
        }
//        std::cout << " comparing with particle " << i << "(";
//        for( auto el: cartesian_coord.at(i).second ){
//            std::cout << el << " ";
//        }

        /* Compute Euclid distance of i-th particle with the others */
        distance = this->euclid_distance(cartesian_coord.at(particle_ind).second, cartesian_coord.at(i).second);
//        std::cout << ") distance = " << distance << " partial contribution " << std::exp(-extension * (distance - shift) * (distance - shift)) * this->cutoff_func->eval(distance) << std::endl;

        /* Call cutoff function */
        sum +=  std::exp(-extension * (distance - shift) * (distance - shift)) * this->cutoff_func->eval(distance);
    }
//    COUT_INFO( "result " << sum);

    return sum;
}

double lib4neuro::G2::eval(unsigned int particle_ind,
                           std::vector<std::vector<double>> cartesian_coord) {
    THROW_NOT_IMPLEMENTED_ERROR();
    return 0;
}

lib4neuro::G2::G2() {}

lib4neuro::G3::G3(lib4neuro::CutoffFunction* cutoff_func,
                  double period_length) : SymmetryFunction(cutoff_func) {
//    this->period_length = period_length;
    this->parameters[SYMMETRY_FUNCTION_PARAMETER::PERIOD_LENGTH] = period_length;
//    this->type = SYMMETRY_FUNCTION_TYPE::G3;
}

lib4neuro::G3::G3() {}

double lib4neuro::G3::eval(unsigned int particle_ind,
                           std::vector<std::pair<ELEMENT_SYMBOL, std::vector<double>>>& cartesian_coord) {
    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;
        }

        /* Compute Euclid distance of i-th particle with the others */
        distance = this->euclid_distance(cartesian_coord.at(particle_ind).second, cartesian_coord.at(i).second);

        /* Call cutoff function */
        sum +=  cos(period_length*distance) * this->cutoff_func->eval(distance);
    }

    return sum;
}

double lib4neuro::G3::eval(unsigned int particle_ind,
                           std::vector<std::vector<double>> cartesian_coord) {
    THROW_NOT_IMPLEMENTED_ERROR();
    return 0;
}

lib4neuro::G4::G4(lib4neuro::CutoffFunction* cutoff_func,
                  double extension,
                  int 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, 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() {}

double lib4neuro::G4::eval(unsigned int particle_ind,
                           std::vector<std::vector<double>> cartesian_coord) {
    THROW_NOT_IMPLEMENTED_ERROR();
    return 0;
}

double lib4neuro::G4::eval(unsigned int particle_ind,
                           std::vector<std::pair<ELEMENT_SYMBOL, std::vector<double>>>& cartesian_coord) {
    double sum = 0;
    double distance1;
    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;
        }

        for(size_t j = 0; j < cartesian_coord.size(); j++) {
            if(j == particle_ind || i == j) {
                continue;
            }

            /* Compute Euclid distances*/
            distance1 = this->euclid_distance(cartesian_coord.at(particle_ind).second, cartesian_coord.at(i).second);
            distance2 = this->euclid_distance(cartesian_coord.at(particle_ind).second, cartesian_coord.at(j).second);
            distance3 = this->euclid_distance(cartesian_coord.at(i).second, cartesian_coord.at(j).second);

            /* Call cutoff function */
            int lambda = shift_max;
            std::vector<double> v1;
            std::vector<double> v2;
            double scalar_product = 0;
            for(size_t k = 0; k < cartesian_coord.at(particle_ind).second.size(); k++) {
                v1.emplace_back(cartesian_coord.at(particle_ind).second.at(k) - cartesian_coord.at(i).second.at(k));
                v2.emplace_back(cartesian_coord.at(particle_ind).second.at(k) - cartesian_coord.at(j).second.at(k));
                scalar_product += v1.at(k)*v2.at(k);
            }

            double theta = acos(scalar_product / (distance1 * distance2));
            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-angular_resolution)*sum;
}

lib4neuro::G5::G5(lib4neuro::CutoffFunction* cutoff_func,
                  double extension,
                  bool shift_max,
                  double angular_resolution) : SymmetryFunction(cutoff_func) {
    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() {}

double lib4neuro::G5::eval(unsigned int particle_ind,
                           std::vector<std::vector<double>> cartesian_coord) {
    THROW_NOT_IMPLEMENTED_ERROR();
    return 0;
}

double lib4neuro::G5::eval(unsigned int particle_ind,
                           std::vector<std::pair<ELEMENT_SYMBOL, std::vector<double>>>& cartesian_coord) {
    double sum = 0;
    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;
        }

        for(size_t j = 0; j < cartesian_coord.size(); j++) {
            if(j == particle_ind || i == j) {
                continue;
            }

            /* Compute Euclid distances*/
            distance1 = this->euclid_distance(cartesian_coord.at(particle_ind).second, cartesian_coord.at(i).second);
            distance2 = this->euclid_distance(cartesian_coord.at(particle_ind).second, cartesian_coord.at(j).second);

            /* Call cutoff function */
            int lambda = (shift_max) ? 1 : -1;
            std::vector<double> v1;
            std::vector<double> v2;
            double scalar_product = 0;
            for(size_t k = 0; k < cartesian_coord.at(particle_ind).second.size(); k++) {
                v1.emplace_back(cartesian_coord.at(particle_ind).second.at(k) - cartesian_coord.at(i).second.at(k));
                v2.emplace_back(cartesian_coord.at(particle_ind).second.at(k) - cartesian_coord.at(j).second.at(k));
                scalar_product += v1.at(k)*v2.at(k);
            }

            double theta = acos(scalar_product / (distance1 * 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-angular_resolution)*sum;
}