<commit_before>from __future__ import absolute_import, division, unicode_literals

import logging

from liknorm import LikNormMachine
from numpy import (asarray, clip, concatenate, dot, exp, inf, log, ndarray,
                   zeros)
from numpy.linalg import LinAlgError

from numpy_sugar import epsilon
from optimix import Function, Scalar, Vector
from optimix.optimize import BadSolutionError

from ..ep import EP
from ..ep.posterior import PosteriorLinearKernel


class GLMM(EP, Function):
    r"""Expectation Propagation for Generalised Gaussian Processes.

    Let

    .. math::

        \mathrm Q \mathrm E \mathrm Q^{\intercal}
        = \mathrm G\mathrm G^{\intercal}

    be the eigen decomposition of the random effect's covariance.
    It turns out that the prior covariance of GLMM can be described as

    .. math::

        \mathrm Q s((1-\delta)\mathrm E
        + \delta\mathrm I) \mathrm Q^{\intercal}.

    This means that :math:`\mathrm Q` does not change during inference, and
    this fact is used here to speed-up EP inference for GLMM.

    Args:
        y (array_like): outcome variable.
        lik_name (str): likelihood name.
        mean (:class:`optimix.Function`): mean function.
                                          (Refer to :doc:`mean`.)
        cov (:class:`optimix.Function`): covariance function.
                                         (Refer to :doc:`cov`.)

    Example
    -------

    .. doctest::

        >>> from numpy import dot, sqrt, zeros
        >>> from numpy.random import RandomState
        >>>
        >>> from numpy_sugar.linalg import economic_qs
        >>>
        >>> from glimix_core.glmm import GLMM
        >>>
        >>> random = RandomState(0)
        >>> nsamples = 50
        >>>
        >>> X = random.randn(50, 2)
        >>> G = random.randn(50, 100)
        >>> K = dot(G, G.T)
        >>> ntrials = random.randint(1, 100, nsamples)
        >>> z = dot(G, random.randn(100)) / sqrt(100)
        >>>
        >>> successes = zeros(len(ntrials), int)
        >>> for i in range(len(ntrials)):
        ...     successes[i] = sum(z[i] + 0.2 * random.randn(ntrials[i]) > 0)
        >>>
        >>> y = (successes, ntrials)
        >>>
        >>> QS = economic_qs(K)
        >>> glmm = GLMM(y, 'binomial', X, QS)
        >>> print('Before: %.4f' % glmm.feed().value())
        Before: -95.1854
        >>> glmm.feed().maximize(progress=False)
        >>> print('After: %.2f' % glmm.feed().value())
        After: -92.24
    """

    def __init__(self, y, lik_name, X, QS):
        super(GLMM, self).__init__(PosteriorLinearKernel)
        Function.__init__(
            self,
            beta=Vector(zeros(X.shape[1])),
            logscale=Scalar(0.0),
            logitdelta=Scalar(0.0))

        self._logger = logging.getLogger(__name__)

        logscale = self.variables()['logscale']
        logscale.bounds = (-5, 4.0)
        logscale.listen(self._clear_cache)

        logitdelta = self.variables()['logitdelta']
        logitdelta.bounds = (-inf, +inf)
        logitdelta.listen(self._clear_cache)

        self.variables()['beta'].listen(self._clear_cache)

        if isinstance(y, list):
            y = tuple(y)
        elif not isinstance(y, tuple):
            y = (y, )

        self._y = tuple([asarray(i, float) for i in y])
        self._X = X

        if not isinstance(QS, tuple):
            raise ValueError("QS must be a tuple.")

        if not isinstance(QS[0], tuple):
            raise ValueError("QS[0] must be a tuple.")

        self._QS = QS

        self._machine = LikNormMachine(lik_name, 500)
        self._initialized = False
        self.set_nodata()
<commit_msg>extreme scale allowed<commit_after>from __future__ import absolute_import, division, unicode_literals

import logging

from liknorm import LikNormMachine
from numpy import (asarray, clip, concatenate, dot, exp, inf, log, ndarray,
                   zeros)
from numpy.linalg import LinAlgError

from numpy_sugar import epsilon
from optimix import Function, Scalar, Vector
from optimix.optimize import BadSolutionError

from ..ep import EP
from ..ep.posterior import PosteriorLinearKernel


class GLMM(EP, Function):
    r"""Expectation Propagation for Generalised Gaussian Processes.

    Let

    .. math::

        \mathrm Q \mathrm E \mathrm Q^{\intercal}
        = \mathrm G\mathrm G^{\intercal}

    be the eigen decomposition of the random effect's covariance.
    It turns out that the prior covariance of GLMM can be described as

    .. math::

        \mathrm Q s((1-\delta)\mathrm E
        + \delta\mathrm I) \mathrm Q^{\intercal}.

    This means that :math:`\mathrm Q` does not change during inference, and
    this fact is used here to speed-up EP inference for GLMM.

    Args:
        y (array_like): outcome variable.
        lik_name (str): likelihood name.
        mean (:class:`optimix.Function`): mean function.
                                          (Refer to :doc:`mean`.)
        cov (:class:`optimix.Function`): covariance function.
                                         (Refer to :doc:`cov`.)

    Example
    -------

    .. doctest::

        >>> from numpy import dot, sqrt, zeros
        >>> from numpy.random import RandomState
        >>>
        >>> from numpy_sugar.linalg import economic_qs
        >>>
        >>> from glimix_core.glmm import GLMM
        >>>
        >>> random = RandomState(0)
        >>> nsamples = 50
        >>>
        >>> X = random.randn(50, 2)
        >>> G = random.randn(50, 100)
        >>> K = dot(G, G.T)
        >>> ntrials = random.randint(1, 100, nsamples)
        >>> z = dot(G, random.randn(100)) / sqrt(100)
        >>>
        >>> successes = zeros(len(ntrials), int)
        >>> for i in range(len(ntrials)):
        ...     successes[i] = sum(z[i] + 0.2 * random.randn(ntrials[i]) > 0)
        >>>
        >>> y = (successes, ntrials)
        >>>
        >>> QS = economic_qs(K)
        >>> glmm = GLMM(y, 'binomial', X, QS)
        >>> print('Before: %.4f' % glmm.feed().value())
        Before: -95.1854
        >>> glmm.feed().maximize(progress=False)
        >>> print('After: %.2f' % glmm.feed().value())
        After: -92.24
    """

    def __init__(self, y, lik_name, X, QS):
        super(GLMM, self).__init__(PosteriorLinearKernel)
        Function.__init__(
            self,
            beta=Vector(zeros(X.shape[1])),
            logscale=Scalar(0.0),
            logitdelta=Scalar(0.0))

        self._logger = logging.getLogger(__name__)

        logscale = self.variables()['logscale']
        logscale.bounds = (-400., 400.)
        # logscale.bounds = (-5., 4.)
        logscale.listen(self._clear_cache)

        logitdelta = self.variables()['logitdelta']
        logitdelta.bounds = (-inf, +inf)
        logitdelta.listen(self._clear_cache)

        self.variables()['beta'].listen(self._clear_cache)

        if isinstance(y, list):
            y = tuple(y)
        elif not isinstance(y, tuple):
            y = (y, )

        self._y = tuple([asarray(i, float) for i in y])
        self._X = X

        if not isinstance(QS, tuple):
            raise ValueError("QS must be a tuple.")

        if not isinstance(QS[0], tuple):
            raise ValueError("QS[0] must be a tuple.")

        self._QS = QS

        self._machine = LikNormMachine(lik_name, 500)
        self._initialized = False
        self.set_nodata()

{'fraction_non_alphanumeric': 0.13344138303619665, 'fraction_numerical': 0.014181523500810372, 'mean_word_length': 2.0410677618069815, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 161, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4591919', 'n_tokens_mistral': 2342, 'n_tokens_neox': 2184, 'n_words': 597}