<commit_before>#!/usr/bin/env python3
# pylint: disable=redefined-outer-name

"""" Gaussian Process Regression model

Implementation is based on Algorithm 2.1 (pg. 19) of
"Gaussian Processes for Machine Learning" by <NAME> Williams

<NAME>, <NAME>
"""

import math

import numpy as np
from scipy.linalg import solve_triangular
from scipy.optimize import minimize
from typing import List

from env import ChemicalEnvironment, two_body, two_body_py, three_body,\
    three_body_py
from struc import Structure


class GaussianProcess:
    """ Gaussian Process Regression Model """

    def __init__(self, kernel: str):
        """Initialize GP parameters and training data.

        :param kernel: covariance / kernel function to be used
        :type kernel: str
        """

        # gp kernel and hyperparameters
        self.kernel_type = kernel
        self.length_scale = 1
        self.sigma_n = 1
        self.sigma_f = 0.01

        # quantities used in GPR algorithm
        self.k_mat = None
        self.l_mat = None
        self.alpha = None

        # training set
        self.training_data = []
        self.training_labels = []
        self.training_labels_np = np.empty(0,)

    def kernel(self, env1: ChemicalEnvironment, env2: ChemicalEnvironment,
               d1: int, d2: int, sig: float, length_scale: float):
        """Evaluate specified kernel.

        :param env1: first chemical environment to compare
        :type evn1: ChemicalEnvironment
        :param env2: second chemical environment to compare
        :type env2: ChemicalEnvironment
        :param d1:
        :type d1: int
        :param d2:
        :type d2: int
        :param sig: signal variance
        :type sig: float
        :param length_scale: length scale of kernel
        :type: length_scale: float
        """

        if self.kernel_type == 'two_body':
            return two_body(env1, env2, d1, d2, sig, length_scale)

        elif self.kernel_type == 'two_body_py':
            return two_body_py(env1, env2, d1, d2, sig, length_scale)

        elif self.kernel_type == 'three_body':
            return three_body(env1, env2, d1, d2, sig, length_scale)

        elif self.kernel_type == 'three_body_py':
            return three_body_py(env1, env2, d1, d2, sig, length_scale)

        else:
            raise ValueError('{} is not a valid kernel'.format(self.kernel))

    def update_db(self, struc: Structure, forces: list):
        """Given structure and forces, add to training set.

        :param struc: structure to add to db
        :type struc: Structure
        :param forces: list of corresponding forces to add to db
        :type forces: list<float>
        """

        noa = len(struc.positions)

        for atom in range(noa):
            env_curr = ChemicalEnvironment(struc, atom)
            forces_curr = forces[atom]

            self.training_data.append(env_curr)
            self.training_labels.append(forces_curr)

        # create numpy array of training labels
        self.training_labels_np = self.force_list_to_np(self.training_labels)

    @staticmethod
    def force_list_to_np(forces: list) -> np.ndarray:
        """ Convert list of forces to numpy array of forces.

        :param forces: list of forces to convert
        :type forces: list<float>
        :return: numpy array forces
        :rtype: np.ndarray
        """
        forces_np = []

        for m in range(len(forces)):
            for n in range(3):
                forces_np.append(forces[m][n])

        forces_np = np.array(forces_np)

        return forces_np

    def train(self):
        """ Train Gaussian Process model on training data. """

        # optimize hyperparameters
        self.opt_hyper()

        # following: Algorithm 2.1 (pg. 19) of
        # "Gaussian Processes for Machine Learning" by <NAME> Williams
        self.set_kernel(sigma_f=self.sigma_f,
                        length_scale=self.length_scale,
                        sigma_n=self.sigma_n)

        # get alpha and likelihood
        self.set_alpha()

    def opt_hyper(self):
        """
        Optimize hyperparameters of GP by minimizing minus log likelihood. """
        # initial guess
        x_0 = np.array([self.sigma_f, self.length_scale, self.sigma_n])

        # nelder-mead optimization
        args = (self,)
        res = minimize(minus_like_hyp, x_0, args,
                       method='nelder-mead',
                       options={'xtol': 1e-8, 'disp': False})

        self.sigma_f = res.x[0]
        self.length_scale = res.x[1]
        self.sigma_n = res.x[2]
<commit_msg>add changes to make H2 example work<commit_after>#!/usr/bin/env python3
# pylint: disable=redefined-outer-name

"""" Gaussian Process Regression model

Implementation is based on Algorithm 2.1 (pg. 19) of
"Gaussian Processes for Machine Learning" by <NAME> Williams

<NAME>, <NAME>
"""

import math

import numpy as np
from scipy.linalg import solve_triangular
from scipy.optimize import minimize
from typing import List

from env import ChemicalEnvironment, two_body, two_body_py, three_body,\
    three_body_py
from struc import Structure


class GaussianProcess:
    """ Gaussian Process Regression Model """

    def __init__(self, kernel: str):
        """Initialize GP parameters and training data.

        :param kernel: covariance / kernel function to be used
        :type kernel: str
        """

        # gp kernel and hyperparameters
        self.kernel_type = kernel
        self.length_scale = 1
        self.sigma_n = 1
        self.sigma_f = 0.01

        # quantities used in GPR algorithm
        self.k_mat = None
        self.l_mat = None
        self.alpha = None

        # training set
        self.training_data = []
        self.training_labels = []
        self.training_labels_np = np.empty(0,)

    def kernel(self, env1: ChemicalEnvironment, env2: ChemicalEnvironment,
               d1: int, d2: int, sig: float, length_scale: float):
        """Evaluate specified kernel.

        :param env1: first chemical environment to compare
        :type evn1: ChemicalEnvironment
        :param env2: second chemical environment to compare
        :type env2: ChemicalEnvironment
        :param d1:
        :type d1: int
        :param d2:
        :type d2: int
        :param sig: signal variance
        :type sig: float
        :param length_scale: length scale of kernel
        :type: length_scale: float
        """

        if self.kernel_type == 'two_body':
            return two_body(env1, env2, d1, d2, sig, length_scale)

        elif self.kernel_type == 'two_body_py':
            return two_body_py(env1, env2, d1, d2, sig, length_scale)

        elif self.kernel_type == 'three_body':
            return three_body(env1, env2, d1, d2, sig, length_scale)

        elif self.kernel_type == 'three_body_py':
            return three_body_py(env1, env2, d1, d2, sig, length_scale)

        else:
            raise ValueError('{} is not a valid kernel'.format(self.kernel))

    def update_db(self, struc: Structure, forces: list):
        """Given structure and forces, add to training set.

        :param struc: structure to add to db
        :type struc: Structure
        :param forces: list of corresponding forces to add to db
        :type forces: list<float>
        """

        noa = len(struc.positions)

        for atom in range(noa):
            env_curr = ChemicalEnvironment(struc, atom)
            forces_curr = forces[atom]

            self.training_data.append(env_curr)
            self.training_labels.append(forces_curr)

        # create numpy array of training labels
        self.training_labels_np = self.force_list_to_np(self.training_labels)

    @staticmethod
    def force_list_to_np(forces: list) -> np.ndarray:
        """ Convert list of forces to numpy array of forces.

        :param forces: list of forces to convert
        :type forces: list<float>
        :return: numpy array forces
        :rtype: np.ndarray
        """
        forces_np = []

        for m in range(len(forces)):
            for n in range(3):
                forces_np.append(forces[m][n])

        forces_np = np.array(forces_np)

        return forces_np

    def train(self):
        """ Train Gaussian Process model on training data. """

        # optimize hyperparameters
        self.opt_hyper()

        # following: Algorithm 2.1 (pg. 19) of
        # "Gaussian Processes for Machine Learning" by <NAME> Williams
        self.set_kernel(sigma_f=self.sigma_f,
                        length_scale=self.length_scale,
                        sigma_n=self.sigma_n)

        # get alpha and likelihood
        self.set_alpha()

    def opt_hyper(self):
        """
        Optimize hyperparameters of GP by minimizing minus log likelihood. """
        # initial guess
        self.sigma_f = 1
        self.length_scale = 1
        self.sigma_n = 0.01

        x_0 = np.array([self.sigma_f, self.length_scale, self.sigma_n])

        # nelder-mead optimization
        args = (self,)
        res = minimize(minus_like_hyp, x_0, args,
                       method='nelder-mead',
                       options={'xtol': 1e-8, 'disp': False})

        self.sigma_f = res.x[0]
        self.length_scale = res.x[1]
        self.sigma_n = res.x[2]

{'fraction_non_alphanumeric': 0.08970493215593367, 'fraction_numerical': 0.011630411371957786, 'mean_word_length': 2.194702442380461, 'pattern_counts': {'":': 0, '<': 15, '<?xml version=': 0, '>': 17, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3256954', 'n_tokens_mistral': 2790, 'n_tokens_neox': 2630, 'n_words': 804}