<commit_before># -*- coding: UTF-8 -*-
"""
    A collection of methods for sampling from continuous models in pgmpy
"""
from __future__ import division
from math import sqrt

import numpy as np

from pgmpy.utils import _check_1d_array_object, _check_length_equal
from pgmpy.sampling import LeapFrog, BaseSimulateHamiltonianDynamics, BaseGradLogPDF, _return_samples


class HamiltonianMC(object):
    """
    Class for performing sampling using simple
    Hamiltonian Monte Carlo

    Parameters:
    -----------
    model: An instance pgmpy.models
        Model from which sampling has to be done

    grad_log_pdf: A subclass of pgmpy.inference.continuous.BaseGradLogPDF, defaults to None
        A class to find log and gradient of log distribution for a given assignment
        If None, then will use model.get_gradient_log_pdf

    simulate_dynamics: A subclass of pgmpy.inference.continuous.BaseSimulateHamiltonianDynamics
        A class to propose future values of momentum and position in time by simulating
        Hamiltonian Dynamics

    Public Methods:
    ---------------
    sample()
    generate_sample()

    Example:
    --------
    >>> from pgmpy.sampling import HamiltonianMC as HMC, LeapFrog, GradLogPDFGaussian
    >>> from pgmpy.factors.continuous import JointGaussianDistribution as JGD
    >>> import numpy as np
    >>> mean = np.array([-3, 4])
    >>> covariance = np.array([[3, 0.7], [0.7, 5]])
    >>> model = JGD(['x', 'y'], mean, covariance)
    >>> sampler = HMC(model=model, grad_log_pdf=GradLogPDFGaussian, simulate_dynamics=LeapFrog)
    >>> samples = sampler.sample(initial_pos=np.array([1, 1]), num_samples = 10000,
    ...                          trajectory_length=2, stepsize=0.4, return_type='recarray')
    >>> samples
    rec.array([(1.0, 1.0), (-3.1861687131079086, 3.7940994520145654),
     (-1.6920542547310844, 6.347410703806017), ...,
     (-1.8093621120575312, 5.940860883943261),
     (0.3933248026088032, 6.3853098838119235),
     (-0.8654072934719572, 6.023803629334816)],
              dtype=[('x', '<f8'), ('y', '<f8')])

    >>> samples = np.array([samples[var_name] for var_name in model.variables])
    >>> np.cov(samples)
    array([[ 3.0352818 ,  0.71379304],
           [ 0.71379304,  4.91776713]])
    >>> sampler.accepted_proposals
    9932.0
    >>> sampler.acceptance_rate
    0.9932

    References
    ----------
    R.Neal. Handbook of Markov Chain Monte Carlo,
    chapter 5: MCMC Using Hamiltonian Dynamics.
    CRC Press, 2011.
    """

    def __init__(self, model, grad_log_pdf, simulate_dynamics=LeapFrog):

        if not issubclass(grad_log_pdf, BaseGradLogPDF):
            raise TypeError("grad_log_pdf must be an instance of " +
                            "pgmpy.inference.base_continuous.BaseGradLogPDF")

        if not issubclass(simulate_dynamics, BaseSimulateHamiltonianDynamics):
            raise TypeError("split_time must be an instance of " +
                            "pgmpy.inference.base_continuous.BaseSimulateHamiltonianDynamics")

        self.model = model
        self.grad_log_pdf = grad_log_pdf
        self.simulate_dynamics = simulate_dynamics
        self.accepted_proposals = 0.0
        self.acceptance_rate = 0

    def _acceptance_prob(self, position, position_bar, momentum, momentum_bar):
        """
        Returns the acceptance probability for given new position(position) and momentum
        """

        # Parameters to help in evaluating Joint distribution P(position, momentum)
        _, logp = self.grad_log_pdf(position, self.model).get_gradient_log_pdf()
        _, logp_bar = self.grad_log_pdf(position_bar, self.model).get_gradient_log_pdf()

        # acceptance_prob = P(position_bar, momentum_bar)/ P(position, momentum)
        potential_change = logp_bar - logp  # Negative change
        kinetic_change = 0.5 * np.float(np.dot(momentum_bar.T, momentum_bar) - np.dot(momentum.T, momentum))

        # acceptance probability
        return np.exp(potential_change - kinetic_change)

    def _find_reasonable_stepsize(self, position, stepsize_app=1):
        """
        Method for choosing initial value of stepsize

        References
        -----------
        <NAME>, <NAME>, The No-U-Turn Sampler: Adaptively
        Setting Path Lengths in Hamiltonian Monte Carlo. Journal of
        Machine Learning Research 15 (2014) 1351-1381
        Algorithm 4 : Heuristic for choosing an initial value of epsilon
        """
        # momentum = N(0, I)
        momentum = np.reshape(np.random.normal(0, 1, len(position)), position.shape)

        # Take a single step in time
        position_bar, momentum_bar, _ =\
            self.simulate_dynamics(self.model, position, momentum,
                                   stepsize_app, self.grad_log_pdf).get_proposed_values()

        acceptance_prob = self._acceptance_prob(position, position_bar, momentum, momentum_bar)

        # a = 2I[acceptance_prob] -1
        a = 2 * (acceptance_prob > 0.5) - 1

        condition = (acceptance_prob ** a) > (2 ** (-a))

        while condition:
            stepsize_app = (2 ** a) * stepsize_app

            position_bar, momentum_bar, _ =\
                self.simulate_dynamics(self.model, position, momentum,
                                       stepsize_app, self.grad_log_pdf).get_proposed_values()

            acceptance_prob = self._acceptance_prob(position, position_bar, momentum, momentum_bar)

            condition = (acceptance_prob ** a) > (2 ** (-a))

        return stepsize_app

    def _sample(self, position, trajectory_length, stepsize, lsteps=None):
        """
        Runs a single sampling iteration to return a sample
        """
        # Resampling momentum
        momentum = np.reshape(np.random.normal(0, 1, len(position)), position.shape)

        # position_m here will be the previous sampled value of position
        position_bar, momentum_bar = position.copy(), momentum

        # Number of steps L to simulate dynamics
        if lsteps is None:
            lsteps = int(max(1, round(trajectory_length / stepsize, 0)))<commit_msg>Tempp orary fix for numpy 0.12:<commit_after># -*- coding: UTF-8 -*-
"""
    A collection of methods for sampling from continuous models in pgmpy
"""
from __future__ import division
from math import sqrt

import numpy as np

from pgmpy.utils import _check_1d_array_object, _check_length_equal
from pgmpy.sampling import LeapFrog, BaseSimulateHamiltonianDynamics, BaseGradLogPDF, _return_samples


class HamiltonianMC(object):
    """
    Class for performing sampling using simple
    Hamiltonian Monte Carlo

    Parameters:
    -----------
    model: An instance pgmpy.models
        Model from which sampling has to be done

    grad_log_pdf: A subclass of pgmpy.inference.continuous.BaseGradLogPDF, defaults to None
        A class to find log and gradient of log distribution for a given assignment
        If None, then will use model.get_gradient_log_pdf

    simulate_dynamics: A subclass of pgmpy.inference.continuous.BaseSimulateHamiltonianDynamics
        A class to propose future values of momentum and position in time by simulating
        Hamiltonian Dynamics

    Public Methods:
    ---------------
    sample()
    generate_sample()

    Example:
    --------
    >>> from pgmpy.sampling import HamiltonianMC as HMC, LeapFrog, GradLogPDFGaussian
    >>> from pgmpy.factors.continuous import JointGaussianDistribution as JGD
    >>> import numpy as np
    >>> mean = np.array([-3, 4])
    >>> covariance = np.array([[3, 0.7], [0.7, 5]])
    >>> model = JGD(['x', 'y'], mean, covariance)
    >>> sampler = HMC(model=model, grad_log_pdf=GradLogPDFGaussian, simulate_dynamics=LeapFrog)
    >>> samples = sampler.sample(initial_pos=np.array([1, 1]), num_samples = 10000,
    ...                          trajectory_length=2, stepsize=0.4, return_type='recarray')
    >>> samples
    rec.array([(1.0, 1.0), (-3.1861687131079086, 3.7940994520145654),
     (-1.6920542547310844, 6.347410703806017), ...,
     (-1.8093621120575312, 5.940860883943261),
     (0.3933248026088032, 6.3853098838119235),
     (-0.8654072934719572, 6.023803629334816)],
              dtype=[('x', '<f8'), ('y', '<f8')])

    >>> samples = np.array([samples[var_name] for var_name in model.variables])
    >>> np.cov(samples)
    array([[ 3.0352818 ,  0.71379304],
           [ 0.71379304,  4.91776713]])
    >>> sampler.accepted_proposals
    9932.0
    >>> sampler.acceptance_rate
    0.9932

    References
    ----------
    R.Neal. Handbook of Markov Chain Monte Carlo,
    chapter 5: MCMC Using Hamiltonian Dynamics.
    CRC Press, 2011.
    """

    def __init__(self, model, grad_log_pdf, simulate_dynamics=LeapFrog):

        if not issubclass(grad_log_pdf, BaseGradLogPDF):
            raise TypeError("grad_log_pdf must be an instance of " +
                            "pgmpy.inference.base_continuous.BaseGradLogPDF")

        if not issubclass(simulate_dynamics, BaseSimulateHamiltonianDynamics):
            raise TypeError("split_time must be an instance of " +
                            "pgmpy.inference.base_continuous.BaseSimulateHamiltonianDynamics")

        self.model = model
        self.grad_log_pdf = grad_log_pdf
        self.simulate_dynamics = simulate_dynamics
        self.accepted_proposals = 0.0
        self.acceptance_rate = 0

    def _acceptance_prob(self, position, position_bar, momentum, momentum_bar):
        """
        Returns the acceptance probability for given new position(position) and momentum
        """

        # Parameters to help in evaluating Joint distribution P(position, momentum)
        _, logp = self.grad_log_pdf(position, self.model).get_gradient_log_pdf()
        _, logp_bar = self.grad_log_pdf(position_bar, self.model).get_gradient_log_pdf()

        # acceptance_prob = P(position_bar, momentum_bar)/ P(position, momentum)
        potential_change = logp_bar - logp  # Negative change
        kinetic_change = 0.5 * np.float(np.dot(momentum_bar.T, momentum_bar) - np.dot(momentum.T, momentum))

        # acceptance probability
        return np.exp(potential_change - kinetic_change)

    def _get_condition(self, acceptance_prob, a):
        """
        Temporary method to fix issue in numpy 0.12 #852
        """
        if a == 1:
            return (acceptance_prob ** a) > (1/(2**a))
        else:
            return (1/(acceptance_prob ** a)) > (2**(-a))

    def _find_reasonable_stepsize(self, position, stepsize_app=1):
        """
        Method for choosing initial value of stepsize

        References
        -----------
        <NAME>, <NAME>, The No-U-Turn Sampler: Adaptively
        Setting Path Lengths in Hamiltonian Monte Carlo. Journal of
        Machine Learning Research 15 (2014) 1351-1381
        Algorithm 4 : Heuristic for choosing an initial value of epsilon
        """
        # momentum = N(0, I)
        momentum = np.reshape(np.random.normal(0, 1, len(position)), position.shape)

        # Take a single step in time
        position_bar, momentum_bar, _ =\
            self.simulate_dynamics(self.model, position, momentum,
                                   stepsize_app, self.grad_log_pdf).get_proposed_values()

        acceptance_prob = self._acceptance_prob(position, position_bar, momentum, momentum_bar)

        # a = 2I[acceptance_prob] -1
        a = 2 * (acceptance_prob > 0.5) - 1

        condition = self._get_condition(acceptance_prob, a)

        while condition:
            stepsize_app = (2 ** a) * stepsize_app

            position_bar, momentum_bar, _ =\
                self.simulate_dynamics(self.model, position, momentum,
                                       stepsize_app, self.grad_log_pdf).get_proposed_values()

            acceptance_prob = self._acceptance_prob(position, position_bar, momentum, momentum_bar)

            condition = self._get_condition(acceptance_prob, a)

        return stepsize_app

    def _sample(self, position, trajectory_length, stepsize, lsteps=None):
        """
        Runs a single sampling iteration to return a sample
        """
        # Resampling momentum
        momentum = np.reshape(np.random.normal(0, 1, len(position)), position.shape)

        # position_m here will be the previous sampled value of position
        position_bar, momentum_bar = position.copy(), momentum

        # Number of steps L to simulate dynamics
        if lsteps is None:
            lsteps = int(max(1, round(trajectory_length / stepsize, 0)))
{'fraction_non_alphanumeric': 0.10347589300016018, 'fraction_numerical': 0.04581130866570559, 'mean_word_length': 2.6384032634032635, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 91, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1164863', 'n_tokens_mistral': 4111, 'n_tokens_neox': 3606, 'n_words': 943}