From 212e4d44b3ad0fd5afa796333de3bea7dcc820ce Mon Sep 17 00:00:00 2001 From: Danilo Horta Date: Thu, 18 Aug 2022 12:17:15 +0100 Subject: [PATCH] LU instead of L --- .gitignore | 1 + .pre-commit-config.yaml | 10 +--- glimix_core/__init__.py | 2 +- glimix_core/_ep/ep.py | 24 ++++---- glimix_core/_ep/linear_kernel.py | 6 +- glimix_core/_ep/posterior.py | 25 +++++---- glimix_core/_ep/posterior_linear_kernel.py | 20 +++---- glimix_core/glmm/test/test_glmm_expfam.py | 4 -- glimix_core/lmm/_kron2sum.py | 32 +++++------ glimix_core/lmm/_kron2sum_scan.py | 4 +- pyproject.toml | 44 +++++++++++++++ setup.cfg | 64 ---------------------- setup.py | 4 -- 13 files changed, 103 insertions(+), 137 deletions(-) create mode 100644 pyproject.toml delete mode 100644 setup.py diff --git a/.gitignore b/.gitignore index 2bb22482..74021076 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,4 @@ +poetry.lock .pytest_cache/ **/__pycache__/ .vscode diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index c7ccf70e..36115bf1 100755 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -1,19 +1,15 @@ repos: - repo: https://github.com/ambv/black - rev: 21.5b1 + rev: 22.6.0 hooks: - id: black language_version: python3 - repo: https://github.com/timothycrosley/isort - rev: 5.8.0 + rev: 5.10.1 hooks: - id: isort - - repo: https://gitlab.com/pycqa/flake8 - rev: 3.9.2 - hooks: - - id: flake8 - repo: https://github.com/pre-commit/pre-commit-hooks - rev: v4.0.1 + rev: v4.3.0 hooks: - id: check-added-large-files - id: trailing-whitespace diff --git a/glimix_core/__init__.py b/glimix_core/__init__.py index 255a85ad..c6d15e29 100644 --- a/glimix_core/__init__.py +++ b/glimix_core/__init__.py @@ -6,7 +6,7 @@ from . import cov, example, ggp, glmm, gp, lik, link, lmm, mean, random from ._testit import test -__version__ = "3.1.12" +__version__ = "3.1.13" __all__ = [ "__version__", diff --git a/glimix_core/_ep/ep.py b/glimix_core/_ep/ep.py index aefdafa2..4e03008c 100644 --- a/glimix_core/_ep/ep.py +++ b/glimix_core/_ep/ep.py @@ -122,14 +122,14 @@ def cav(self): return self._cav def lml(self): - from numpy_sugar.linalg import cho_solve + from numpy_sugar.linalg import lu_slogdet, lu_solve if self._cache["lml"] is not None: return self._cache["lml"] self._update() - L = self._posterior.L() + LU = self._posterior.LU() Q, S = self._posterior.cov["QS"] Q = Q[0] ttau = self._site.tau @@ -141,13 +141,13 @@ def lml(self): TS = ttau + ctau lml = [ - -log(L.diagonal()).sum(), + -lu_slogdet(LU)[1], -0.5 * sum(log(S)), # lml += 0.5 * sum(log(ttau)), - +0.5 * dot(teta, dot(Q, cho_solve(L, dot(Q.T, teta)))), + +0.5 * dot(teta, dot(Q, lu_solve(LU, dot(Q.T, teta)))), -0.5 * dot(teta, teta / TS), +dot(m, teta) - 0.5 * dot(m, ttau * m), - -0.5 * dot(m * ttau, dot(Q, cho_solve(L, dot(Q.T, 2 * teta - ttau * m)))), + -0.5 * dot(m * ttau, dot(Q, lu_solve(LU, dot(Q.T, 2 * teta - ttau * m)))), +sum(self._moments["log_zeroth"]), +0.5 * sum(log(TS)), # lml -= 0.5 * sum(log(ttau)), @@ -163,11 +163,11 @@ def lml(self): return lml def lml_derivative_over_cov(self, dQS): - from numpy_sugar.linalg import cho_solve, ddot, dotd + from numpy_sugar.linalg import ddot, dotd, lu_solve self._update() - L = self._posterior.L() + LU = self._posterior.LU() Q = self._posterior.cov["QS"][0][0] ttau = self._site.tau teta = self._site.eta @@ -175,7 +175,7 @@ def lml_derivative_over_cov(self, dQS): diff = teta - ttau * self._posterior.mean v0 = dot(dQS[0][0], dQS[1] * dot(dQS[0][0].T, diff)) - v1 = ttau * dot(Q, cho_solve(L, dot(Q.T, diff))) + v1 = ttau * dot(Q, lu_solve(LU, dot(Q.T, diff))) dlml = 0.5 * dot(diff, v0) dlml -= dot(v0, v1) dlml += 0.5 * dot(v1, dot(dQS[0][0], dQS[1] * dot(dQS[0][0].T, v1))) @@ -183,17 +183,17 @@ def lml_derivative_over_cov(self, dQS): diag = dotd(dQS[0][0], dqs) dlml -= 0.5 * sum(ttau * diag) - tmp = cho_solve(L, dot(ddot(Q.T, ttau, left=False), dQS[0][0])) + tmp = lu_solve(LU, dot(ddot(Q.T, ttau, left=False), dQS[0][0])) dlml += 0.5 * sum(ttau * dotd(Q, dot(tmp, dqs))) return dlml def lml_derivative_over_mean(self, dm): - from numpy_sugar.linalg import cho_solve + from numpy_sugar.linalg import lu_solve self._update() - L = self._posterior.L() + LU = self._posterior.LU() Q = self._posterior.cov["QS"][0][0] ttau = self._site.tau teta = self._site.eta @@ -201,7 +201,7 @@ def lml_derivative_over_mean(self, dm): diff = teta - ttau * self._posterior.mean dlml = dot(diff, dm) - dlml -= dot(diff, dot(Q, cho_solve(L, dot(Q.T, (ttau * dm.T).T)))) + dlml -= dot(diff, dot(Q, lu_solve(LU, dot(Q.T, (ttau * dm.T).T)))) return dlml diff --git a/glimix_core/_ep/linear_kernel.py b/glimix_core/_ep/linear_kernel.py index 0f13547b..55b294fa 100644 --- a/glimix_core/_ep/linear_kernel.py +++ b/glimix_core/_ep/linear_kernel.py @@ -25,12 +25,14 @@ def __init__(self, nsites, rtol=1.49e-05, atol=1.49e-08): ) def lml(self): + from numpy_sugar.linalg import lu_slogdet + if self._cache["lml"] is not None: return self._cache["lml"] self._update() - L = self._posterior.L() + LU = self._posterior.LU() LQt = self._posterior.LQt() cov = self._posterior.cov Q = cov["QS"][0][0] @@ -49,7 +51,7 @@ def lml(self): dif = 2 * A * teta - A * ttau * m lml = [ - -log(L.diagonal()).sum(), + -lu_slogdet(LU)[1], -0.5 * sum(log(s * S)), +0.5 * sum(log(A)), # lml += 0.5 * sum(log(ttau)), diff --git a/glimix_core/_ep/posterior.py b/glimix_core/_ep/posterior.py index c7092de2..3358f5a0 100644 --- a/glimix_core/_ep/posterior.py +++ b/glimix_core/_ep/posterior.py @@ -1,4 +1,5 @@ -from numpy import dot, empty, empty_like, sqrt, sum as npsum +from numpy import dot, empty, empty_like, sqrt +from numpy import sum as npsum class Posterior(object): @@ -23,7 +24,7 @@ def __init__(self, site): self._NxR_data = None self._RxN_data = None self._RxR_data = None - self._L_cache = None + self._LU_cache = None self._LQt_cache = None self._QS_cache = None self._AQ_cache = None @@ -53,7 +54,7 @@ def _RxR(self): return self._RxR_data def _flush_cache(self): - self._L_cache = None + self._LU_cache = None self._LQt_cache = None self._QS_cache = None self._AQ_cache = None @@ -104,8 +105,8 @@ def cov(self, v): self._initialize() self._cov = v - def L(self): - r"""Cholesky decomposition of :math:`\mathrm B`. + def LU(self): + r"""LU factor of :math:`\mathrm B`. .. math:: @@ -113,20 +114,20 @@ def L(self): + \mathrm{S}^{-1} """ from numpy_sugar.linalg import ddot, sum2diag - from scipy.linalg import cho_factor + from scipy.linalg import lu_factor - if self._L_cache is not None: - return self._L_cache + if self._LU_cache is not None: + return self._LU_cache Q = self._cov["QS"][0][0] S = self._cov["QS"][1] B = dot(Q.T, ddot(self._site.tau, Q, left=True)) sum2diag(B, 1.0 / S, out=B) - self._L_cache = cho_factor(B, lower=True)[0] - return self._L_cache + self._LU_cache = lu_factor(B, overwrite_a=True, check_finite=False) + return self._LU_cache def LQt(self): - from numpy_sugar.linalg import cho_solve + from numpy_sugar.linalg import lu_solve if self._LQt_cache is not None: return self._LQt_cache @@ -134,7 +135,7 @@ def LQt(self): L = self.L() Q = self._cov["QS"][0][0] - self._LQt_cache = cho_solve(L, Q.T) + self._LQt_cache = lu_solve(L, Q.T) return self._LQt_cache def QS(self): diff --git a/glimix_core/_ep/posterior_linear_kernel.py b/glimix_core/_ep/posterior_linear_kernel.py index 033f05fe..ff9f319c 100644 --- a/glimix_core/_ep/posterior_linear_kernel.py +++ b/glimix_core/_ep/posterior_linear_kernel.py @@ -3,13 +3,6 @@ from .posterior import Posterior -def _cho_factor(B): - from scipy.linalg import cho_factor - - B = cho_factor(B, overwrite_a=True, lower=True, check_finite=False)[0] - return B - - class PosteriorLinearKernel(Posterior): """ EP posterior. @@ -62,8 +55,8 @@ def QSQtATQLQtA(self): self._QSQtATQLQtA_cache = dotd(QS, dot(dot(ddot(AQ.T, T), Q), LQtA)) return self._QSQtATQLQtA_cache - def L(self): - r"""Cholesky decomposition of :math:`\mathrm B`. + def LU(self): + r"""LU factor of :math:`\mathrm B`. .. math:: @@ -71,9 +64,10 @@ def L(self): + \mathrm{S}^{-1} """ from numpy_sugar.linalg import ddot, sum2diag + from scipy.linalg import lu_factor - if self._L_cache is not None: - return self._L_cache + if self._LU_cache is not None: + return self._LU_cache s = self._cov["scale"] d = self._cov["delta"] @@ -84,8 +78,8 @@ def L(self): B = dot(Q.T, self._NxR, out=self._RxR) B *= 1 - d sum2diag(B, 1.0 / S / s, out=B) - self._L_cache = _cho_factor(B) - return self._L_cache + self._LU_cache = lu_factor(B, overwrite_a=True, check_finite=False) + return self._LU_cache def update(self): from numpy_sugar.linalg import ddot, dotd diff --git a/glimix_core/glmm/test/test_glmm_expfam.py b/glimix_core/glmm/test/test_glmm_expfam.py index be85fab7..919e7028 100644 --- a/glimix_core/glmm/test/test_glmm_expfam.py +++ b/glimix_core/glmm/test/test_glmm_expfam.py @@ -16,7 +16,6 @@ from numpy.random import RandomState from numpy.testing import assert_, assert_allclose from numpy_sugar.linalg import economic_qs, economic_qs_linear -from pandas import DataFrame from glimix_core.example import linear_eye_cov from glimix_core.glmm import GLMMExpFam, GLMMNormal @@ -606,9 +605,6 @@ def test_glmmexpfam_poisson(): offset = ones(n) * random.randn() age = random.randint(16, 75, n) M = stack((offset, age), axis=1) - M = DataFrame(stack([offset, age], axis=1), columns=["offset", "age"]) - M["sample"] = [f"sample{i}" for i in range(n)] - M = M.set_index("sample") # genetic variants G = random.randn(n, 4) diff --git a/glimix_core/lmm/_kron2sum.py b/glimix_core/lmm/_kron2sum.py index e16e767b..f35e23b5 100644 --- a/glimix_core/lmm/_kron2sum.py +++ b/glimix_core/lmm/_kron2sum.py @@ -196,7 +196,7 @@ def C0(self): C0 : ndarray C₀. """ - return self._cov.C0.value() / (self._G_norm ** 2) + return self._cov.C0.value() / (self._G_norm**2) @property def C1(self): @@ -419,8 +419,8 @@ def _XY(self): @property def _terms(self): - from numpy_sugar.linalg import ddot, sum2diag - from scipy.linalg import cho_factor, cho_solve + from numpy_sugar.linalg import ddot, lu_solve, sum2diag + from scipy.linalg import lu_factor if self._cache["terms"] is not None: return self._cache["terms"] @@ -439,7 +439,7 @@ def _terms(self): Z = kron(L0.T @ WL0, self._GG) Z = sum2diag(Z, 1) - Lz = cho_factor(Z, lower=True) + Lz = lu_factor(Z, check_finite=False) # 𝐲ᵀR⁻¹𝐲 = vec(YW)ᵀ𝐲 yRiy = (YW * self._Y).sum() @@ -452,8 +452,8 @@ def _terms(self): # MᵀR⁻¹𝐲 = vec(XᵀYWA) MRiy = vec(self._XY @ WA) - ZiXRiM = cho_solve(Lz, XRiM) - ZiXRiy = cho_solve(Lz, XRiy) + ZiXRiM = lu_solve(Lz, XRiM) + ZiXRiy = lu_solve(Lz, XRiy) MRiXZiXRiy = ZiXRiM.T @ XRiy MRiXZiXRiM = XRiM.T @ ZiXRiM @@ -461,8 +461,8 @@ def _terms(self): yKiy = yRiy - XRiy @ ZiXRiy MKiy = MRiy - MRiXZiXRiy H = MRiM - MRiXZiXRiM - Lh = cho_factor(H) - b = cho_solve(Lh, MKiy) + Lh = lu_factor(H, check_finite=False) + b = lu_solve(Lh, MKiy) B = unvec(b, (self.ncovariates, -1)) self._mean.B = B XRim = XRiM @ b @@ -571,7 +571,7 @@ def _lml_gradient(self): C1.Lu : ndarray Gradient of the log of the marginal likelihood over C₁ parameters. """ - from scipy.linalg import cho_solve + from numpy_sugar.linalg import lu_solve terms = self._terms dC0 = self._cov.C0.gradient()["Lu"] @@ -621,10 +621,10 @@ def _lml_gradient(self): yR0y = vec(_mdot(self._GY, WdC0)).T @ vec(self._GY @ W) yR1y = (YW.T * _mdot(self._Y, WdC1).T).T.sum(axis=(0, 1)) - ZiXR0X = cho_solve(Lz, XR0X) - ZiXR1X = cho_solve(Lz, XR1X) - ZiXR0y = cho_solve(Lz, XR0y) - ZiXR1y = cho_solve(Lz, XR1y) + ZiXR0X = lu_solve(Lz, XR0X) + ZiXR1X = lu_solve(Lz, XR1X) + ZiXR0y = lu_solve(Lz, XR0y) + ZiXR1y = lu_solve(Lz, XR1y) # Mᵀ𝕂y = Mᵀ𝓡𝐲 - (MᵀR⁻¹X)Z⁻¹(Xᵀ𝓡𝐲) - (Mᵀ𝓡X)Z⁻¹(XᵀR⁻¹𝐲) # + (MᵀR⁻¹X)Z⁻¹(Xᵀ𝓡X)Z⁻¹(XᵀR⁻¹𝐲) @@ -655,7 +655,7 @@ def _lml_gradient(self): XRim = XRiM @ b MRim = MRiM @ b - db = {"C0.Lu": cho_solve(Lh, MK0m - MK0y), "C1.Lu": cho_solve(Lh, MK1m - MK1y)} + db = {"C0.Lu": lu_solve(Lh, MK0m - MK0y), "C1.Lu": lu_solve(Lh, MK1m - MK1y)} grad = { "C0.Lu": -trace(WdC0) * self._trGG + trace(ZiXR0X), @@ -663,8 +663,8 @@ def _lml_gradient(self): } if self._restricted: - grad["C0.Lu"] += cho_solve(Lh, MK0M).diagonal().sum(1) - grad["C1.Lu"] += cho_solve(Lh, MK1M).diagonal().sum(1) + grad["C0.Lu"] += lu_solve(Lh, MK0M).diagonal().sum(1) + grad["C1.Lu"] += lu_solve(Lh, MK1M).diagonal().sum(1) mKiM = MRim.T - XRim.T @ ZiXRiM yKiM = MRiy.T - XRiy.T @ ZiXRiM diff --git a/glimix_core/lmm/_kron2sum_scan.py b/glimix_core/lmm/_kron2sum_scan.py index bc9bc778..17221e5a 100644 --- a/glimix_core/lmm/_kron2sum_scan.py +++ b/glimix_core/lmm/_kron2sum_scan.py @@ -173,7 +173,7 @@ def scan(self, A1, X1): from numpy import empty from numpy.linalg import multi_dot from numpy_sugar import epsilon, is_all_finite - from scipy.linalg import cho_solve + from numpy_sugar.linalg import lu_solve A1 = asarray(A1, float) X1 = asarray(X1, float) @@ -206,7 +206,7 @@ def scan(self, A1, X1): M1Riy = vec(multi_dot([X1.T, self._Y, A1W.T])) XRiM1 = kron(self._WL0.T @ A1, GX1) - ZiXRiM1 = cho_solve(self._Lz, XRiM1) + ZiXRiM1 = lu_solve(self._Lz, XRiM1) MRiXZiXRiM1 = self._XRiM.T @ ZiXRiM1 M1RiXZiXRiM1 = XRiM1.T @ ZiXRiM1 diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 00000000..1e682f06 --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,44 @@ +[tool.poetry] +name = "glimix-core" +version = "3.1.13" +description = "Fast inference over mean and covariance parameters for Generalised Linear Mixed Models" + +license = "MIT" + +authors = ["Danilo Horta "] + +readme = "README.md" + +repository = "https://github.com/limix/glimix-core" +homepage = "https://github.com/limix/glimix-core" + +keywords = ["glmm", "lmm"] + +classifiers = ["License :: OSI Approved :: MIT License"] + +[tool.poetry.dependencies] +brent-search = "*" +liknorm = ">=1.2.7" +ndarray-listener = "*" +numpy = "*" +numpy-sugar = "*" +optimix = "*" +pytest = "*" +pytest-doctestplus = "*" +python = ">=3.7,<3.11" +scipy = "*" +tqdm = "*" + +[tool.poetry.dev-dependencies] +black = "*" +isort = "*" +pyright = "*" +pytest = "*" +pytest-doctestplus = "*" + +[build-system] +requires = ["poetry_core>=1.0.0"] +build-backend = "poetry.core.masonry.api" + +[tool.isort] +profile = "black" diff --git a/setup.cfg b/setup.cfg index a70fe5db..1890bcd6 100644 --- a/setup.cfg +++ b/setup.cfg @@ -1,52 +1,6 @@ -[metadata] -author = Danilo Horta -author_email = horta@ebi.ac.uk -classifiers = - Development Status :: 5 - Production/Stable - License :: OSI Approved :: MIT License - Operating System :: OS Independent - Programming Language :: Python -description = Fast inference over mean and covariance parameters for Generalised Linear Mixed Models -download_url = https://github.com/limix/glimix-core -keywords = function, optimisation -license = MIT -long_description = file: README.md -long_description_content_type = text/markdown -maintainer = Danilo Horta -platforms = Windows, MacOS, Linux -maintainer_email = horta@ebi.ac.uk -name = glimix-core -url = https://github.com/limix/glimix-core -version = attr: version.get - -[options] -zip_safe = True -include_package_data = True -packages = find: -python_requires = >= 3.7 -setup_requires = - pytest-runner>=5.2 -install_requires = - brent-search>=2.0.1 - cachetools>=2.1.0 - liknorm>=1.2.4 - ndarray-listener>=2.0.1 - numpy-sugar>=1.5.1 - numpy>=1.14.3 - optimix>=3.0.4 - pandas>=1.0.4 - pytest-doctestplus>=0.8.0 - pytest>=5.4.3 - scipy>=1.0.1 - tqdm>=4.46.1 - -[aliases] -test = pytest - [tool:pytest] addopts = --doctest-plus - --doctest-modules --doctest-glob='*.rst' --ignore="setup.py" --ignore="doc/conf.py" @@ -56,21 +10,3 @@ doctest_plus_atol = 1e-03 doctest_plus_rtol = 1e-03 doctest_rst = enabled norecursedirs = .eggs .git *.egg-info build .ropeproject .undodir - -[tool:isort] -multi_line_output=3 -include_trailing_comma=True -force_grid_wrap=0 -combine_as_imports=True -line_length=88 - -[pylint] -disable = redefined-builtin,R0915 - -[flake8] -ignore = E501 E741 E203 W503 W0212 W0622 R0915 - -[rstcheck] -ignore_substitutions = today, version -ignore_directives = plot, autofunction, command-output, autmodule, automodule, autoclass, autoattribute, automethod, doctest -ignore_messages = Error in "math" directive diff --git a/setup.py b/setup.py deleted file mode 100644 index 7f1a1763..00000000 --- a/setup.py +++ /dev/null @@ -1,4 +0,0 @@ -from setuptools import setup - -if __name__ == "__main__": - setup()