# Source code for glimix_core.ggp._expfam

import warnings

from liknorm import LikNormMachine
from numpy import ascontiguousarray, sign
from numpy.linalg import LinAlgError
from optimix import Function

from .._ep import EP
from .._util import check_outcome

[docs]class ExpFamGP(Function):
r"""Expectation Propagation for Generalised Gaussian Processes.

Parameters
----------
y : array_like
Outcome variable.
lik_name : str
Likelihood name.
mean : function
Mean function. (Refer to :doc:mean.)
cov : function
Covariance function. (Refer to :doc:cov.)

Example
-------

.. doctest::

>>> from numpy.random import RandomState
>>>
>>> from glimix_core.example import offset_mean
>>> from glimix_core.example import linear_eye_cov
>>> from glimix_core.ggp import ExpFamGP
>>> from glimix_core.lik import BernoulliProdLik
>>> from glimix_core.random import GGPSampler
>>>
>>> random = RandomState(1)
>>>
>>> mean = offset_mean()
>>> cov = linear_eye_cov()
>>>
>>> y = GGPSampler(lik, mean, cov).sample(random)
>>>
>>> ggp = ExpFamGP(y, 'bernoulli', mean, cov)
>>> print('Before: %.4f' % ggp.lml())
Before: -6.9802
>>> ggp.fit(verbose=False)
>>> print('After: %.2f' % ggp.lml())
After: -6.55
"""

def __init__(self, y, lik, mean, cov):
if isinstance(y, tuple):
n = len(y[0])
else:
n = len(y)
Function.__init__(self, "ExpFamGP", composite=[mean, cov])

if not isinstance(lik, (tuple, list)):
lik = (lik,)

self._lik = (lik[0].lower(),) + tuple(ascontiguousarray(i) for i in lik[1:])
self._y = check_outcome(y, self._lik)

self._mean = mean
self._cov = cov
self._ep = EP(n)
self._ep.set_compute_moments(self.compute_moments)

self._mean_value = None
self._cov_value = None

self._machine = LikNormMachine(lik[0], 500)

def __del__(self):
if hasattr(self, "_machine"):
self._machine.finish()

[docs]    def fit(self, verbose=True, factr=1e5, pgtol=1e-7):
r"""Maximise the marginal likelihood.

Parameters
----------
verbose : bool
True for progress output; False otherwise.
Defaults to True.
factr : float, optional
The iteration stops when
(f^k - f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= factr * eps, where eps is
the machine precision.
pgtol : float, optional
The iteration will stop when max{|proj g_i | i = 1, ..., n} <= pgtol
where pg_i is the i-th component of the projected gradient.

Notes
-----
Please, refer to :func:scipy.optimize.fmin_l_bfgs_b for further information
about factr and pgtol.
"""
self._maximize(verbose=verbose, factr=factr, pgtol=pgtol)

[docs]    def lml(self):
r"""Log of the marginal likelihood.

Returns
-------
float
:math:\log p(\mathbf y)
"""
return self.value()

def compute_moments(self, eta, tau, moments):
y = (self._y,) + self._lik[1:]
self._machine.moments(y, eta, tau, moments)

def value(self):
from numpy_sugar import epsilon
from numpy_sugar.linalg import economic_qs

mean = self._mean.value()
cov = self._cov.value()
try:
self._ep.set_prior(mean, dict(QS=economic_qs(cov)))
lml = self._ep.lml()
except (ValueError, LinAlgError) as e:
warnings.warn(str(e), RuntimeWarning)
lml = -1 / epsilon.small
return lml

from numpy_sugar import epsilon
from numpy_sugar.linalg import economic_qs

mean = self._mean.value()
cov = self._cov.value()

try:
self._ep.set_prior(mean, dict(QS=economic_qs(cov)))

for n, g in iter(gmean.items()):