Remove unused MvNormalLogp Op

Author: ricardoV94

pymc/distributions/dist_math.py

@@ -28,21 +28,16 @@
 import scipy.linalg
 import scipy.stats
 
-from pytensor.compile.builders import OpFromGraph
 from pytensor.graph.basic import Apply, Variable
 from pytensor.graph.op import Op
 from pytensor.scalar import UnaryScalarOp, upgrade_to_float_no_complex
 from pytensor.tensor import gammaln
 from pytensor.tensor.elemwise import Elemwise
-from pytensor.tensor.slinalg import Cholesky, SolveTriangular
 
 from pymc.distributions.shape_utils import to_tuple
 from pymc.logprob.utils import CheckParameterValue
 from pymc.pytensorf import floatX
 
-solve_lower = SolveTriangular(lower=True)
-solve_upper = SolveTriangular(lower=False)
-
 f = floatX
 c = -0.5 * np.log(2.0 * np.pi)
 _beta_clip_values = {
@@ -242,69 +237,6 @@ def log_normal(x, mean, **kwargs):
     return f(c) - pt.log(pt.abs(std)) - (x - mean) ** 2 / (2.0 * std**2)
 
 
-def MvNormalLogp():
-    """Compute the log pdf of a multivariate normal distribution.
-
-    This should be used in MvNormal.logp once Theano#5908 is released.
-
-    Parameters
-    ----------
-    cov: pt.matrix
-        The covariance matrix.
-    delta: pt.matrix
-        Array of deviations from the mean.
-    """
-    cov = pt.matrix("cov")
-    cov.tag.test_value = floatX(np.eye(3))
-    delta = pt.matrix("delta")
-    delta.tag.test_value = floatX(np.zeros((2, 3)))
-
-    cholesky = Cholesky(lower=True, on_error="nan")
-
-    n, k = delta.shape
-    n, k = f(n), f(k)
-    chol_cov = cholesky(cov)
-    diag = pt.diag(chol_cov)
-    ok = pt.all(diag > 0)
-
-    chol_cov = pt.switch(ok, chol_cov, pt.fill(chol_cov, 1))
-    delta_trans = solve_lower(chol_cov, delta.T).T
-
-    result = n * k * pt.log(f(2) * np.pi)
-    result += f(2) * n * pt.sum(pt.log(diag))
-    result += (delta_trans ** f(2)).sum()
-    result = f(-0.5) * result
-    logp = pt.switch(ok, result, -np.inf)
-
-    def dlogp(inputs, gradients):
-        (g_logp,) = gradients
-        cov, delta = inputs
-
-        g_logp.tag.test_value = floatX(1.0)
-        n, k = delta.shape
-
-        chol_cov = cholesky(cov)
-        diag = pt.diag(chol_cov)
-        ok = pt.all(diag > 0)
-
-        chol_cov = pt.switch(ok, chol_cov, pt.fill(chol_cov, 1))
-        delta_trans = solve_lower(chol_cov, delta.T).T
-
-        inner = n * pt.eye(k) - pt.dot(delta_trans.T, delta_trans)
-        g_cov = solve_upper(chol_cov.T, inner)
-        g_cov = solve_upper(chol_cov.T, g_cov.T)
-
-        tau_delta = solve_upper(chol_cov.T, delta_trans.T)
-        g_delta = tau_delta.T
-
-        g_cov = pt.switch(ok, g_cov, -np.nan)
-        g_delta = pt.switch(ok, g_delta, -np.nan)
-
-        return [-0.5 * g_cov * g_logp, -g_delta * g_logp]
-
-    return OpFromGraph([cov, delta], [logp], grad_overrides=dlogp, inline=True)
-
-
 class SplineWrapper(Op):
     """
     Creates an PyTensor operation from scipy.interpolate.UnivariateSpline

pymc/distributions/multivariate.py

@@ -155,8 +155,6 @@ def quaddist_parse(value, mu, cov, mat_type="cov"):
         onedim = False
 
     delta = value - mu
-    # Use this when Theano#5908 is released.
-    # return MvNormalLogp()(self.cov, delta)
     chol_cov = cholesky(cov)
     if mat_type != "tau":
         dist, logdet, ok = quaddist_chol(delta, chol_cov)

tests/distributions/test_dist_math.py

@@ -26,7 +26,6 @@
 
 from pymc.distributions import Discrete
 from pymc.distributions.dist_math import (
-    MvNormalLogp,
     SplineWrapper,
     check_parameters,
     clipped_beta_rvs,
@@ -123,65 +122,6 @@ def test_multinomial_check_parameters():
     )
 
 
-class TestMvNormalLogp:
-    def test_logp(self):
-        np.random.seed(42)
-
-        chol_val = floatX(np.array([[1, 0.9], [0, 2]]))
-        cov_val = floatX(np.dot(chol_val, chol_val.T))
-        cov = pt.matrix("cov")
-        cov.tag.test_value = cov_val
-        delta_val = floatX(np.random.randn(5, 2))
-        delta = pt.matrix("delta")
-        delta.tag.test_value = delta_val
-        expect = stats.multivariate_normal(mean=np.zeros(2), cov=cov_val)
-        expect = expect.logpdf(delta_val).sum()
-        logp = MvNormalLogp()(cov, delta)
-        logp_f = pytensor.function([cov, delta], logp)
-        logp = logp_f(cov_val, delta_val)
-        npt.assert_allclose(logp, expect)
-
-    @pytensor.config.change_flags(compute_test_value="ignore")
-    def test_grad(self):
-        np.random.seed(42)
-
-        def func(chol_vec, delta):
-            chol = pt.stack(
-                [
-                    pt.stack([pt.exp(0.1 * chol_vec[0]), 0]),
-                    pt.stack([chol_vec[1], 2 * pt.exp(chol_vec[2])]),
-                ]
-            )
-            cov = pt.dot(chol, chol.T)
-            return MvNormalLogp()(cov, delta)
-
-        chol_vec_val = floatX(np.array([0.5, 1.0, -0.1]))
-
-        delta_val = floatX(np.random.randn(1, 2))
-        verify_grad(func, [chol_vec_val, delta_val])
-
-        delta_val = floatX(np.random.randn(5, 2))
-        verify_grad(func, [chol_vec_val, delta_val])
-
-    @pytensor.config.change_flags(compute_test_value="ignore")
-    def test_hessian(self):
-        chol_vec = pt.vector("chol_vec")
-        chol_vec.tag.test_value = floatX(np.array([0.1, 2, 3]))
-        chol = pt.stack(
-            [
-                pt.stack([pt.exp(0.1 * chol_vec[0]), 0]),
-                pt.stack([chol_vec[1], 2 * pt.exp(chol_vec[2])]),
-            ]
-        )
-        cov = pt.dot(chol, chol.T)
-        delta = pt.matrix("delta")
-        delta.tag.test_value = floatX(np.ones((5, 2)))
-        logp = MvNormalLogp()(cov, delta)
-        g_cov, g_delta = pt.grad(logp, [cov, delta])
-        # TODO: What's the test?  Something needs to be asserted.
-        pt.grad(g_delta.sum() + g_cov.sum(), [delta, cov])
-
-
 class TestSplineWrapper:
     @pytensor.config.change_flags(compute_test_value="ignore")
     def test_grad(self):