From 20f0e20b313c2a20efa4abd34ed7879635d02bd8 Mon Sep 17 00:00:00 2001 From: Rutger van Haasteren Date: Fri, 21 Apr 2023 17:55:05 +0200 Subject: [PATCH 01/20] First addition of FastShermanMorrison for testing --- enterprise/fastshermanmorrison/__init__.py | 1 + .../cython_fastshermanmorrison.pyx | 945 ++++++++++++++++++ .../fastshermanmorrison/fastshermanmorrison.c | 233 +++++ .../fastshermanmorrison.py | 163 +++ setup.py | 12 +- 5 files changed, 1353 insertions(+), 1 deletion(-) create mode 100644 enterprise/fastshermanmorrison/__init__.py create mode 100644 enterprise/fastshermanmorrison/cython_fastshermanmorrison.pyx create mode 100644 enterprise/fastshermanmorrison/fastshermanmorrison.c create mode 100644 enterprise/fastshermanmorrison/fastshermanmorrison.py diff --git a/enterprise/fastshermanmorrison/__init__.py b/enterprise/fastshermanmorrison/__init__.py new file mode 100644 index 00000000..3037395f --- /dev/null +++ b/enterprise/fastshermanmorrison/__init__.py @@ -0,0 +1 @@ +from . import fastshermanmorrison diff --git a/enterprise/fastshermanmorrison/cython_fastshermanmorrison.pyx b/enterprise/fastshermanmorrison/cython_fastshermanmorrison.pyx new file mode 100644 index 00000000..c8a3b8c2 --- /dev/null +++ b/enterprise/fastshermanmorrison/cython_fastshermanmorrison.pyx @@ -0,0 +1,945 @@ +cimport numpy as np +import numpy as np +np.import_array() + +from libc.math cimport log, sqrt +import cython +from scipy.linalg.cython_blas cimport dgemm, dger, dgemv + +cdef public void dgemm_(char *transa, char *transb, int *m, int *n, int *k, + double *alpha, double *a, int *lda, double *b, + int *ldb, double *beta, double *c, int *ldc): + """Public dgemm that can be used in the external C code""" + dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc) + +cdef public void dgemv_(char *trans, int *m, int *n, double *alpha, double *a, + int *lda, double *x, int *incx, double *beta, + double *y, int *incy): + """Public dgemv that can be used in the external C code""" + dgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy) + +cdef public void dger_(int *m, int *n, double *alpha, double *x, int *incx, + double *y, int *incy, double *a, int *lda): + """Public dger that can be used in the external C code""" + dger(m, n, alpha, x, incx, y, incy, a, lda) + + +cdef extern from "fastshermanmorrison.c": + void blas_block_shermor_2D_asym( + int n_Z1_rows, + int n_Z1_cols, + int n_Z1_row_major, + double *pd_Z1, + int n_Z2_cols, + int n_Z2_row_major, + double *pd_Z2, + double *pd_Nvec, + int n_J_rows, + double *pd_Jvec, + int *pn_Uinds, + double *pd_ZNZ, + double *pd_Jldet + ) + +@cython.boundscheck(False) +@cython.wraparound(False) +def cython_block_shermor_0D( \ + np.ndarray[np.double_t,ndim=1] r, \ + np.ndarray[np.double_t,ndim=1] Nvec, \ + np.ndarray[np.double_t,ndim=1] Jvec, \ + np.ndarray[np.int_t,ndim=2] Uinds): + """ + Sherman-Morrison block-inversion for Jitter (Cythonized) + + @param r: The timing residuals, array (n) + @param Nvec: The white noise amplitude, array (n) + @param Jvec: The jitter amplitude, array (k) + @param Uinds: The start/finish indices for the jitter blocks (k x 2) + + For this version, the residuals need to be sorted properly so that all the + blocks are continuous in memory. Here, there are n residuals, and k jitter + parameters. + """ + cdef unsigned int cc, ii, cols = len(Jvec) + cdef double Jldet=0.0, ji, beta, nir, nisum + cdef np.ndarray[np.double_t,ndim=1] ni = np.empty(len(r), 'd') + cdef np.ndarray[np.double_t,ndim=1] Nx = r / Nvec + + ni = 1.0 / Nvec + + for cc in range(cols): + if Jvec[cc] > 0.0: + ji = 1.0 / Jvec[cc] + + nir = 0.0 + nisum = 0.0 + for ii in range(Uinds[cc,0],Uinds[cc,1]): + nisum += ni[ii] + nir += r[ii]*ni[ii] + + beta = 1.0 / (nisum + ji) + + for ii in range(Uinds[cc,0],Uinds[cc,1]): + Nx[ii] -= beta * nir * ni[ii] + + return Nx + +@cython.boundscheck(False) +@cython.wraparound(False) +def cython_block_shermor_0D_ld( \ + np.ndarray[np.double_t,ndim=1] r, \ + np.ndarray[np.double_t,ndim=1] Nvec, \ + np.ndarray[np.double_t,ndim=1] Jvec, \ + np.ndarray[np.int_t,ndim=2] Uinds): + """ + Sherman-Morrison block-inversion for Jitter (Cythonized) + + @param r: The timing residuals, array (n) + @param Nvec: The white noise amplitude, array (n) + @param Jvec: The jitter amplitude, array (k) + @param Uinds: The start/finish indices for the jitter blocks (k x 2) + + For this version, the residuals need to be sorted properly so that all the + blocks are continuous in memory. Here, there are n residuals, and k jitter + parameters. + """ + cdef unsigned int cc, ii, rows = len(r), cols = len(Jvec) + cdef double Jldet=0.0, ji, beta, nir, nisum + cdef np.ndarray[np.double_t,ndim=1] ni = np.empty(len(r), 'd') + cdef np.ndarray[np.double_t,ndim=1] Nx = r / Nvec + + ni = 1.0 / Nvec + + for cc in range(rows): + Jldet += log(Nvec[cc]) + + for cc in range(cols): + if Jvec[cc] > 0.0: + ji = 1.0 / Jvec[cc] + + nir = 0.0 + nisum = 0.0 + for ii in range(Uinds[cc,0],Uinds[cc,1]): + nisum += ni[ii] + nir += r[ii]*ni[ii] + + beta = 1.0 / (nisum + ji) + Jldet += log(Jvec[cc]) - log(beta) + + for ii in range(Uinds[cc,0],Uinds[cc,1]): + Nx[ii] -= beta * nir * ni[ii] + + return Jldet, Nx + + +def python_block_shermor_1D(r, Nvec, Jvec, Uinds): + """ + Sherman-Morrison block-inversion for Jitter + + @param r: The timing residuals, array (n) + @param Nvec: The white noise amplitude, array (n) + @param Jvec: The jitter amplitude, array (k) + @param Uinds: The start/finish indices for the jitter blocks (k x 2) + + For this version, the residuals need to be sorted properly so that all the + blocks are continuous in memory. Here, there are n residuals, and k jitter + parameters. + """ + ni = 1.0 / Nvec + Jldet = np.einsum('i->', np.log(Nvec)) + xNx = np.dot(r, r * ni) + + for cc, jv in enumerate(Jvec): + if jv > 0.0: + rblock = r[Uinds[cc,0]:Uinds[cc,1]] + niblock = ni[Uinds[cc,0]:Uinds[cc,1]] + + beta = 1.0 / (np.einsum('i->', niblock)+1.0/jv) + xNx -= beta * np.dot(rblock, niblock)**2 + Jldet += np.log(jv) - np.log(beta) + + return Jldet, xNx + +@cython.boundscheck(False) +@cython.wraparound(False) +def cython_block_shermor_1D( \ + np.ndarray[np.double_t,ndim=1] r, \ + np.ndarray[np.double_t,ndim=1] Nvec, \ + np.ndarray[np.double_t,ndim=1] Jvec, \ + np.ndarray[np.int_t,ndim=2] Uinds): + """ + Sherman-Morrison block-inversion for Jitter (Cythonized) + + @param r: The timing residuals, array (n) + @param Nvec: The white noise amplitude, array (n) + @param Jvec: The jitter amplitude, array (k) + @param Uinds: The start/finish indices for the jitter blocks (k x 2) + + For this version, the residuals need to be sorted properly so that all the + blocks are continuous in memory. Here, there are n residuals, and k jitter + parameters. + """ + cdef unsigned int cc, ii, rows = len(r), cols = len(Jvec) + cdef double Jldet=0.0, ji, beta, xNx=0.0, nir, nisum + cdef np.ndarray[np.double_t,ndim=1] ni = np.empty(rows, 'd') + + ni = 1.0 / Nvec + + for cc in range(rows): + Jldet += log(Nvec[cc]) + xNx += r[cc]*r[cc]*ni[cc] + + for cc in range(cols): + if Jvec[cc] > 0.0: + ji = 1.0 / Jvec[cc] + + nir = 0.0 + nisum = 0.0 + for ii in range(Uinds[cc,0],Uinds[cc,1]): + nisum += ni[ii] + nir += r[ii]*ni[ii] + + beta = 1.0 / (nisum + ji) + Jldet += log(Jvec[cc]) - log(beta) + xNx -= beta * nir * nir + + return Jldet, xNx + +@cython.boundscheck(False) +@cython.wraparound(False) +def cython_block_shermor_1D1( \ + np.ndarray[np.double_t,ndim=1] x, \ + np.ndarray[np.double_t,ndim=1] y, \ + np.ndarray[np.double_t,ndim=1] Nvec, \ + np.ndarray[np.double_t,ndim=1] Jvec, \ + np.ndarray[np.int_t,ndim=2] Uinds): + """ + Sherman-Morrison block-inversion for Jitter (Cythonized) + + @param r: The timing residuals, array (n) + @param Nvec: The white noise amplitude, array (n) + @param Jvec: The jitter amplitude, array (k) + @param Uinds: The start/finish indices for the jitter blocks (k x 2) + + For this version, the residuals need to be sorted properly so that all the + blocks are continuous in memory. Here, there are n residuals, and k jitter + parameters. + """ + cdef unsigned int cc, ii, rows = len(x), cols = len(Jvec) + cdef double Jldet=0.0, ji, beta, yNx=0.0, nix, niy, nisum + cdef np.ndarray[np.double_t,ndim=1] ni = np.empty(rows, 'd') + + ni = 1.0 / Nvec + + for cc in range(rows): + Jldet += log(Nvec[cc]) + yNx += y[cc]*x[cc]*ni[cc] + + for cc in range(cols): + if Jvec[cc] > 0.0: + ji = 1.0 / Jvec[cc] + + nix = 0.0 + niy = 0.0 + nisum = 0.0 + for ii in range(Uinds[cc,0],Uinds[cc,1]): + nisum += ni[ii] + nix += x[ii]*ni[ii] + niy += y[ii]*ni[ii] + + beta = 1.0 / (nisum + ji) + Jldet += log(Jvec[cc]) - log(beta) + yNx -= beta * nix * niy + + return Jldet, yNx + + +def python_block_shermor_2D(Z, Nvec, Jvec, Uinds): + """ + Sherman-Morrison block-inversion for Jitter, ZNiZ + + @param Z: The design matrix, array (n x m) + @param Nvec: The white noise amplitude, array (n) + @param Jvec: The jitter amplitude, array (k) + @param Uinds: The start/finish indices for the jitter blocks (k x 2) + + For this version, the residuals need to be sorted properly so that all the + blocks are continuous in memory. Here, there are n residuals, and k jitter + parameters. + + N = D + U*J*U.T + calculate: log(det(N)), Z.T * N^-1 * Z + """ + ni = 1.0 / Nvec + Jldet = np.einsum('i->', np.log(Nvec)) + zNz = np.dot(Z.T*ni, Z) + + for cc, jv in enumerate(Jvec): + if jv > 0.0: + Zblock = Z[Uinds[cc,0]:Uinds[cc,1], :] + niblock = ni[Uinds[cc,0]:Uinds[cc,1]] + + beta = 1.0 / (np.einsum('i->', niblock)+1.0/jv) + zn = np.dot(niblock, Zblock) + zNz -= beta * np.outer(zn.T, zn) + Jldet += np.log(jv) - np.log(beta) + + return Jldet, zNz + +@cython.boundscheck(False) +@cython.wraparound(False) +def cython_block_shermor_2D( \ + np.ndarray[np.double_t,ndim=2] Z, \ + np.ndarray[np.double_t,ndim=1] Nvec, \ + np.ndarray[np.double_t,ndim=1] Jvec, \ + np.ndarray[np.int_t,ndim=2] Uinds): + """ + Sherman-Morrison block-inversion for Jitter (Cythonized) + + @param Z: The design matrix, array (n x m) + @param Nvec: The white noise amplitude, array (n) + @param Jvec: The jitter amplitude, array (k) + @param Uinds: The start/finish indices for the jitter blocks (k x 2) + + For this version, the residuals need to be sorted properly so that all the + blocks are continuous in memory. Here, there are n residuals, and k jitter + parameters. + + N = D + U*J*U.T + calculate: log(det(N)), Z.T * N^-1 * Z + """ + cdef unsigned int cc, ii, rows = len(Nvec), cols = len(Jvec) + cdef double Jldet=0.0, ji, beta, nir, nisum + cdef np.ndarray[np.double_t,ndim=1] ni = np.empty(len(Nvec), 'd') + cdef np.ndarray[np.double_t,ndim=2] zNz + + ni = 1.0 / Nvec + zNz = np.dot(Z.T*ni, Z) + + for cc in range(rows): + Jldet += log(Nvec[cc]) + + for cc in range(cols): + if Jvec[cc] > 0.0: + Zblock = Z[Uinds[cc,0]:Uinds[cc,1], :] + niblock = ni[Uinds[cc,0]:Uinds[cc,1]] + + nisum = 0.0 + for ii in range(len(niblock)): + nisum += niblock[ii] + + beta = 1.0 / (nisum+1.0/Jvec[cc]) + Jldet += log(Jvec[cc]) - log(beta) + zn = np.dot(niblock, Zblock) + zNz -= beta * np.outer(zn.T, zn) + + return Jldet, zNz + +def python_block_shermor_2D_asymm(Z1, Z2, Nvec, Jvec, Uinds): + """ + Sherman-Morrison block-inversion for Jitter, ZNiZ + + @param Z: The design matrix, array (n x m) + @param Z2: The second design matrix, array (n x m2) + @param Nvec: The white noise amplitude, array (n) + @param Jvec: The jitter amplitude, array (k) + @param Uinds: The start/finish indices for the jitter blocks (k x 2) + + For this version, the residuals need to be sorted properly so that all the + blocks are continuous in memory. Here, there are n residuals, and k jitter + parameters. + + N = D + U*J*U.T + calculate: log(det(N)), Z.T * N^-1 * Z + """ + ni = 1.0 / Nvec + Jldet = np.einsum('i->', np.log(Nvec)) + zNz = np.dot(Z1.T*ni, Z2) + + for cc, jv in enumerate(Jvec): + if jv > 0.0: + Zblock1 = Z1[Uinds[cc,0]:Uinds[cc,1], :] + Zblock2 = Z2[Uinds[cc,0]:Uinds[cc,1], :] + niblock = ni[Uinds[cc,0]:Uinds[cc,1]] + + beta = 1.0 / (np.einsum('i->', niblock)+1.0/jv) + zn1 = np.dot(niblock, Zblock1) + zn2 = np.dot(niblock, Zblock2) + zNz -= beta * np.outer(zn1.T, zn2) + Jldet += np.log(jv) - np.log(beta) + + return Jldet, zNz + + +@cython.boundscheck(False) +@cython.wraparound(False) +def cython_block_shermor_2D_asymm( + np.ndarray[np.double_t,ndim=2] Z1, + np.ndarray[np.double_t,ndim=2] Z2, + np.ndarray[np.double_t,ndim=1] Nvec, + np.ndarray[np.double_t,ndim=1] Jvec, + np.ndarray[np.int_t,ndim=2] Uinds): + """ + Sherman-Morrison block-inversion for Jitter, ZNiZ + + @param Z: The design matrix, array (n x m) + @param Z2: The second design matrix, array (n x m2) + @param Nvec: The white noise amplitude, array (n) + @param Jvec: The jitter amplitude, array (k) + @param Uinds: The start/finish indices for the jitter blocks (k x 2) + + For this version, the residuals need to be sorted properly so that all the + blocks are continuous in memory. Here, there are n residuals, and k jitter + parameters. + + N = D + U*J*U.T + calculate: log(det(N)), Z.T * N^-1 * Z + """ + cdef unsigned int cc, ii, rows = len(Nvec), cols = len(Jvec) + cdef double Jldet=0.0, ji, beta, nir, nisum + cdef np.ndarray[np.double_t,ndim=1] ni = np.empty(len(Nvec), 'd') + cdef np.ndarray[np.double_t,ndim=2] zNz + + print("WARNING: cython_block_shermor_2D_asymm is deprecated.") + print(" use cython_blas_block_shermor_2D_asymm") + + ni = 1.0 / Nvec + for cc in range(rows): + Jldet += log(Nvec[cc]) + zNz = np.dot(Z1.T*ni, Z2) + + for cc in range(cols): + if Jvec[cc] > 0.0: + Zblock1 = Z1[Uinds[cc,0]:Uinds[cc,1], :] + Zblock2 = Z2[Uinds[cc,0]:Uinds[cc,1], :] + niblock = ni[Uinds[cc,0]:Uinds[cc,1]] + + nisum = 0.0 + for ii in range(len(niblock)): + nisum += niblock[ii] + + beta = 1.0 / (nisum+1.0/Jvec[cc]) + zn1 = np.dot(niblock, Zblock1) + zn2 = np.dot(niblock, Zblock2) + + zNz -= beta * np.outer(zn1.T, zn2) + Jldet += log(Jvec[cc]) - log(beta) + + return Jldet, zNz + + +def python_draw_ecor(r, Nvec, Jvec, Uinds): + """ + Given Jvec, draw new epoch-averaged residuals + + @param r: The timing residuals, array (n) + @param Nvec: The white noise amplitude, array (n) + @param Jvec: The jitter amplitude, array (k) + @param Uinds: The start/finish indices for the jitter blocks (k x 2) + + For this version, the residuals need to be sorted properly so that all the + blocks are continuous in memory. Here, there are n residuals, and k jitter + parameters. + + N = D + U*J*U.T + calculate: Norm(0, sqrt(J)) + (U^T * D^{-1} * U)^{-1}U.T D^{-1} r + """ + + rv = np.random.randn(len(Jvec)) * np.sqrt(Jvec) + ni = 1.0 / Nvec + + for cc in range(len(Jvec)): + rblock = r[Uinds[cc,0]:Uinds[cc,1]] + niblock = ni[Uinds[cc,0]:Uinds[cc,1]] + beta = 1.0 / np.einsum('i->', niblock) + + rv[cc] += beta * np.dot(rblock, niblock) + + return rv + +@cython.boundscheck(False) +@cython.wraparound(False) +def cython_draw_ecor( \ + np.ndarray[np.double_t,ndim=1] r, \ + np.ndarray[np.double_t,ndim=1] Nvec, \ + np.ndarray[np.double_t,ndim=1] Jvec, \ + np.ndarray[np.int_t,ndim=2] Uinds): + """ + Given Jvec, draw new epoch-averaged residuals + + @param r: The timing residuals, array (n) + @param Nvec: The white noise amplitude, array (n) + @param Jvec: The jitter amplitude, array (k) + @param Uinds: The start/finish indices for the jitter blocks (k x 2) + + For this version, the residuals need to be sorted properly so that all the + blocks are continuous in memory. Here, there are n residuals, and k jitter + parameters. + + N = D + U*J*U.T + calculate: Norm(0, sqrt(J)) + (U^T * D^{-1} * U)^{-1}U.T D^{-1} r + """ + cdef unsigned int cc, ii, rows = len(r), cols = len(Jvec) + cdef double ji, nir, nisum + cdef np.ndarray[np.double_t,ndim=1] ni = np.empty(rows, 'd') + cdef np.ndarray[np.double_t,ndim=1] rv = np.random.randn(cols) + + for cc in range(cols): + rv[cc] *= sqrt(Jvec[cc]) + + ni = 1.0 / Nvec + + for cc in range(cols): + ji = 1.0 / Jvec[cc] + + nir = 0.0 + nisum = 0.0 + for ii in range(Uinds[cc,0],Uinds[cc,1]): + nisum += ni[ii] + nir += r[ii]*ni[ii] + + rv[cc] += nir / nisum + + return rv + + +@cython.boundscheck(False) +@cython.wraparound(False) +def cython_shermor_draw_ecor( \ + np.ndarray[np.double_t,ndim=1] r, \ + np.ndarray[np.double_t,ndim=1] Nvec, \ + np.ndarray[np.double_t,ndim=1] Jvec, \ + np.ndarray[np.int_t,ndim=2] Uinds): + """ + Do both the Sherman-Morrison block-inversion for Jitter, + and the draw of the ecor parameters together (Cythonized) + + @param r: The timing residuals, array (n) + @param Nvec: The white noise amplitude, array (n) + @param Jvec: The jitter amplitude, array (k) + @param Uinds: The start/finish indices for the jitter blocks (k x 2) + + For this version, the residuals need to be sorted properly so that all the + blocks are continuous in memory. Here, there are n residuals, and k jitter + parameters. + + N = D + U*J*U.T + calculate: r.T * N^-1 * r, log(det(N)), Norm(0, sqrt(J)) + (U^T * D^{-1} * U)^{-1}U.T D^{-1} r + """ + cdef unsigned int cc, ii, rows = len(r), cols = len(Jvec) + cdef double Jldet=0.0, ji, beta, xNx=0.0, nir, nisum + cdef np.ndarray[np.double_t,ndim=1] ni = np.empty(rows, 'd') + cdef np.ndarray[np.double_t,ndim=1] rv = np.random.randn(cols) + + ni = 1.0 / Nvec + + for cc in range(cols): + rv[cc] *= sqrt(Jvec[cc]) + + for cc in range(rows): + Jldet += log(Nvec[cc]) + xNx += r[cc]*r[cc]*ni[cc] + + for cc in range(cols): + nir = 0.0 + nisum = 0.0 + for ii in range(Uinds[cc,0],Uinds[cc,1]): + nisum += ni[ii] + nir += r[ii]*ni[ii] + + rv[cc] += nir / nisum + + if Jvec[cc] > 0.0: + ji = 1.0 / Jvec[cc] + + beta = 1.0 / (nisum + ji) + Jldet += log(Jvec[cc]) - log(beta) + xNx -= beta * nir * nir + + return Jldet, xNx, rv + + +@cython.boundscheck(False) +@cython.wraparound(False) +def cython_update_ea_residuals( \ + np.ndarray[np.double_t,ndim=1] gibbsresiduals, \ + np.ndarray[np.double_t,ndim=1] gibbssubresiduals, \ + np.ndarray[np.double_t,ndim=1] eat, \ + np.ndarray[np.int_t,ndim=2] Uinds): + """ + Given epoch-averaged residuals, update the residuals, and the subtracted + residuals, so that these can be further processed by the other conditional + probability density functions. + + @param gibbsresiduals: The timing residuals, array (n) + @param gibbssubresiduals: The white noise amplitude, array (n) + @param eat: epoch averaged residuals (k) + @param Uinds: The start/finish indices for the jitter blocks + (k x 2) + + """ + cdef unsigned int k = Uinds.shape[0], ii, cc + + for cc in range(Uinds.shape[0]): + for ii in range(Uinds[cc,0],Uinds[cc,1]): + gibbssubresiduals[ii] += eat[cc] + gibbsresiduals[ii] -= eat[cc] + + return gibbsresiduals, gibbssubresiduals + + +@cython.boundscheck(False) +@cython.wraparound(False) +def cython_Uj(np.ndarray[np.double_t,ndim=1] j, \ + np.ndarray[np.int_t,ndim=2] Uinds, nobs): + """ + Given epoch-averaged residuals (j), get the residuals. + Used in 'updateDetSources' + + @param j: epoch averaged residuals (k) + @param Uinds: The start/finish indices for the jitter blocks + (k x 2) + @param nobs: Number of observations (length return vector) + + """ + cdef unsigned int k = Uinds.shape[0], ii, cc + cdef np.ndarray[np.double_t,ndim=1] Uj = np.zeros(nobs, 'd') + + for cc in range(k): + for ii in range(Uinds[cc,0],Uinds[cc,1]): + Uj[ii] += j[cc] + + return Uj + +@cython.boundscheck(False) +@cython.wraparound(False) +def cython_UTx(np.ndarray[np.double_t,ndim=1] x, \ + np.ndarray[np.int_t,ndim=2] Uinds): + """ + Given residuals (x), get np.dot(U.T, x) + Used in 'updateDetSources' + + @param j: epoch averaged residuals (k) + @param Uinds: The start/finish indices for the jitter blocks + (k x 2) + + """ + cdef unsigned int k = Uinds.shape[0], ii, cc + cdef np.ndarray[np.double_t,ndim=1] UTx = np.zeros(k, 'd') + + for cc in range(k): + for ii in range(Uinds[cc,0],Uinds[cc,1]): + UTx[cc] += x[ii] + + return UTx + +@cython.boundscheck(False) +@cython.wraparound(False) +def cython_logdet_dN( \ + np.ndarray[np.double_t,ndim=1] Nvec, \ + np.ndarray[np.double_t,ndim=1] Jvec, \ + np.ndarray[np.double_t,ndim=1] dNvec, \ + np.ndarray[np.int_t,ndim=2] Uinds): + """ + Sherman-Morrison block-inversion for Jitter (Cythonized) + + Calculates Trace(N^{-1} dN/dNp), where: + - N^{-1} is the ecorr-include N inverse + - dN/dNp is the diagonal derivate of N wrt Np + + @param Nvec: The white noise amplitude, array (n) + @param Jvec: The jitter amplitude, array (k) + @param dNvec: The white noise derivative, array (n) + @param Uinds: The start/finish indices for the jitter blocks (k x 2) + + For this version, the residuals need to be sorted properly so that all the + blocks are continuous in memory. Here, there are n residuals, and k jitter + parameters. + """ + cdef unsigned int cc, ii, rows = len(Nvec), cols = len(Jvec) + cdef double tr=0.0, ji, nisum, Nnisum + cdef np.ndarray[np.double_t,ndim=1] ni = np.empty(rows, 'd') + cdef np.ndarray[np.double_t,ndim=1] Nni = np.empty(rows, 'd') + + ni = 1.0 / Nvec + Nni = dNvec / Nvec**2 + + for cc in range(rows): + tr += dNvec[cc] * ni[cc] + + for cc in range(cols): + if Jvec[cc] > 0.0: + ji = 1.0 / Jvec[cc] + + nisum = 0.0 + Nnisum = 0.0 + for ii in range(Uinds[cc,0],Uinds[cc,1]): + nisum += ni[ii] + Nnisum += Nni[ii] + + tr -= Nnisum / (nisum + ji) + + return tr + +@cython.boundscheck(False) +@cython.wraparound(False) +def cython_logdet_dJ( \ + np.ndarray[np.double_t,ndim=1] Nvec, \ + np.ndarray[np.double_t,ndim=1] Jvec, \ + np.ndarray[np.double_t,ndim=1] dJvec, \ + np.ndarray[np.int_t,ndim=2] Uinds): + """ + Sherman-Morrison block-inversion for Jitter (Cythonized) + + Calculates Trace(N^{-1} dN/dJp), where: + - N^{-1} is the ecorr-include N inverse + - dN/dJp = U dJ/dJp U^{T}, with dJ/dJp the diagnal derivative of J wrt + Jp + + @param Nvec: The white noise amplitude, array (n) + @param Jvec: The jitter amplitude, array (k) + @param dJvec: The jitter derivative, array (k) + @param Uinds: The start/finish indices for the jitter blocks (k x 2) + + For this version, the residuals need to be sorted properly so that all the + blocks are continuous in memory. Here, there are n residuals, and k jitter + parameters. + """ + cdef unsigned int cc, ii, rows = len(Nvec), cols = len(Jvec) + cdef double dJldet=0.0, ji, beta, nisum + cdef np.ndarray[np.double_t,ndim=1] ni = np.empty(rows, 'd') + + ni = 1.0 / Nvec + + for cc in range(cols): + if Jvec[cc] > 0.0: + ji = 1.0 / Jvec[cc] + + nisum = 0.0 + for ii in range(Uinds[cc,0],Uinds[cc,1]): + nisum += ni[ii] + + beta = 1.0 / (nisum + ji) + + dJldet += dJvec[cc]*(nisum - beta*nisum**2) + + return dJldet + +@cython.boundscheck(False) +@cython.wraparound(False) +def cython_logdet_dN_dN( \ + np.ndarray[np.double_t,ndim=1] Nvec, \ + np.ndarray[np.double_t,ndim=1] Jvec, \ + np.ndarray[np.double_t,ndim=1] dNvec1, \ + np.ndarray[np.double_t,ndim=1] dNvec2, \ + np.ndarray[np.int_t,ndim=2] Uinds): + """ + Sherman-Morrison block-inversion for Jitter (Cythonized) + + Calculates Trace(N^{-1} dN/dNp1 N^{-1} dN/dNp2), where: + - N^{-1} is the ecorr-include N inverse + - dN/dNpx is the diagonal derivate of N wrt Npx + + @param Nvec: The white noise amplitude, array (n) + @param Jvec: The jitter amplitude, array (k) + @param dNvec1: The white noise derivative, array (n) + @param dNvec2: The white noise derivative, array (n) + @param Uinds: The start/finish indices for the jitter blocks (k x 2) + + For this version, the residuals need to be sorted properly so that all the + blocks are continuous in memory. Here, there are n residuals, and k jitter + parameters. + """ + cdef unsigned int cc, ii, rows = len(Nvec), cols = len(Jvec) + cdef double tr=0.0, ji, nisum, Nnisum1, Nnisum2, NniNnisum, beta + cdef np.ndarray[np.double_t,ndim=1] ni = np.empty(rows, 'd') + cdef np.ndarray[np.double_t,ndim=1] Nni1 = np.empty(rows, 'd') + cdef np.ndarray[np.double_t,ndim=1] Nni2 = np.empty(rows, 'd') + + ni = 1.0 / Nvec + Nni1 = dNvec1 / Nvec**2 + Nni2 = dNvec2 / Nvec**2 + + for cc in range(rows): + tr += dNvec1[cc] * dNvec2[cc] * ni[cc]**2 + + for cc in range(cols): + if Jvec[cc] > 0.0: + ji = 1.0 / Jvec[cc] + + nisum = 0.0 + Nnisum1 = 0.0 + Nnisum2 = 0.0 + NniNnisum = 0.0 + for ii in range(Uinds[cc,0],Uinds[cc,1]): + nisum += ni[ii] + Nnisum1 += Nni1[ii] + Nnisum2 += Nni2[ii] + NniNnisum += Nni1[ii]*Nni2[ii]*Nvec[ii] + + beta = 1.0 / (nisum + ji) + + tr += Nnisum1 * Nnisum2 * beta**2 + tr -= 2 * NniNnisum * beta + + return tr + +@cython.boundscheck(False) +@cython.wraparound(False) +def cython_logdet_dN_dJ( \ + np.ndarray[np.double_t,ndim=1] Nvec, \ + np.ndarray[np.double_t,ndim=1] Jvec, \ + np.ndarray[np.double_t,ndim=1] dNvec, \ + np.ndarray[np.double_t,ndim=1] dJvec, \ + np.ndarray[np.int_t,ndim=2] Uinds): + """ + Sherman-Morrison block-inversion for Jitter (Cythonized) + + Calculates Trace(N^{-1} dN/dNp N^{-1} dN/dJp), where: + - N^{-1} is the ecorr-include N inverse + - dN/dNp is the diagonal derivate of N wrt Np + - dN/dJp = U dJ/dJp U^{T}, with dJ/dJp the diagnal derivative of J wrt + Jp + + @param Nvec: The white noise amplitude, array (n) + @param Jvec: The jitter amplitude, array (k) + @param dNvec: The white noise derivative, array (n) + @param dJvec: The white noise ecor derivative, array (k) + @param Uinds: The start/finish indices for the jitter blocks (k x 2) + + For this version, the residuals need to be sorted properly so that all the + blocks are continuous in memory. Here, there are n residuals, and k jitter + parameters. + """ + cdef unsigned int cc, ii, rows = len(Nvec), cols = len(Jvec) + cdef double tr=0.0, ji, nisum, Nnisum, beta + cdef np.ndarray[np.double_t,ndim=1] ni = np.empty(rows, 'd') + cdef np.ndarray[np.double_t,ndim=1] Nni = np.empty(rows, 'd') + + ni = 1.0 / Nvec + Nni = dNvec / Nvec**2 + + for cc in range(cols): + if Jvec[cc] > 0.0: + ji = 1.0 / Jvec[cc] + + nisum = 0.0 + Nnisum = 0.0 + for ii in range(Uinds[cc,0],Uinds[cc,1]): + nisum += ni[ii] + Nnisum += Nni[ii] + + beta = 1.0 / (nisum + ji) + + tr += Nnisum * dJvec[cc] + tr -= 2 * nisum * dJvec[cc] * Nnisum * beta + tr += Nnisum * nisum**2 * dJvec[cc] *beta**2 + + return tr + +@cython.boundscheck(False) +@cython.wraparound(False) +def cython_logdet_dJ_dJ( \ + np.ndarray[np.double_t,ndim=1] Nvec, \ + np.ndarray[np.double_t,ndim=1] Jvec, \ + np.ndarray[np.double_t,ndim=1] dJvec1, \ + np.ndarray[np.double_t,ndim=1] dJvec2, \ + np.ndarray[np.int_t,ndim=2] Uinds): + """ + Sherman-Morrison block-inversion for Jitter (Cythonized) + + Calculates Trace(N^{-1} dN/dJp1 N^{-1} dN/dJp2), where: + - N^{-1} is the ecorr-include N inverse + - dN/dJpx = U dJ/dJpx U^{T}, with dJ/dJpx the diagnal derivative of J wrt + Jpx + + @param Nvec: The white noise amplitude, array (n) + @param Jvec: The jitter amplitude, array (k) + @param dJvec1: The white noise derivative, array (k) + @param dJvec2: The white noise derivative, array (k) + @param Uinds: The start/finish indices for the jitter blocks (k x 2) + + For this version, the residuals need to be sorted properly so that all the + blocks are continuous in memory. Here, there are n residuals, and k jitter + parameters. + """ + cdef unsigned int cc, ii, rows = len(Nvec), cols = len(Jvec) + cdef double tr=0.0, ji, nisum, beta + cdef np.ndarray[np.double_t,ndim=1] ni = np.empty(rows, 'd') + + ni = 1.0 / Nvec + + for cc in range(cols): + if Jvec[cc] > 0.0: + ji = 1.0 / Jvec[cc] + + nisum = 0.0 + for ii in range(Uinds[cc,0],Uinds[cc,1]): + nisum += ni[ii] + + beta = 1.0 / (nisum + ji) + + tr += dJvec1[cc] * dJvec2[cc] * nisum**2 + tr -= 2 * dJvec1[cc] * dJvec2[cc] * beta * nisum**3 + tr += dJvec1[cc] * dJvec2[cc] * beta**2 * nisum**4 + + return tr + + +@cython.boundscheck(False) +@cython.wraparound(False) +cpdef double c_blas_block_shermor_2D_asymm( + np.ndarray[np.double_t,ndim=2] Z1, + np.ndarray[np.double_t,ndim=2] Z2, + np.ndarray[np.double_t,ndim=1] Nvec, + np.ndarray[np.double_t,ndim=1] Jvec, + Uinds, #Need to copy, because Cython can't always use numpy integer arrays. Bah + np.ndarray[np.double_t,ndim=2] ZNZ, + ): + cdef double d_Jldet + + cdef int n_Z1_rows = Z1.shape[0] + cdef int n_Z1_cols = Z1.shape[1] + cdef int n_Z2_cols = Z2.shape[1] + cdef int n_J_rows = len(Jvec) + cdef int n_Z1_row_major, n_Z2_row_major + + n_Z1_row_major = 0 if Z1.flags['F_CONTIGUOUS'] else 1 + n_Z2_row_major = 0 if Z2.flags['F_CONTIGUOUS'] else 1 + + # Hack, because somehow we can't pass integer arrays? + # This is a tiny bit of overhead right here + cdef int [:] Uinds_new + + # This makes it into a C-ordered (Row-Major) array we can pass + Uinds_new = np.ascontiguousarray(Uinds.flatten(), dtype=np.dtype("i")) + + assert Z1.shape[0] == Z2.shape[0] + assert Z1.shape[0] == len(Nvec) + assert Uinds.shape[0] == len(Jvec) + + blas_block_shermor_2D_asym( + n_Z1_rows, + n_Z1_cols, + n_Z1_row_major, + &Z1[0,0], + n_Z2_cols, + n_Z2_row_major, + &Z2[0,0], + &Nvec[0], + n_J_rows, + &Jvec[0], + &Uinds_new[0], + &ZNZ[0,0], + &d_Jldet) + + return d_Jldet + +def cython_blas_block_shermor_2D_asymm(Z1, Z2, Nvec, Jvec, Uinds): + """Wrapper for the C/Cython code""" + + ZNZ = np.zeros((Z1.shape[1], Z2.shape[1]), order='F') + Jldet = c_blas_block_shermor_2D_asymm(Z1, Z2, Nvec, Jvec, np.array(Uinds, order="C"), ZNZ) + + return Jldet, ZNZ + diff --git a/enterprise/fastshermanmorrison/fastshermanmorrison.c b/enterprise/fastshermanmorrison/fastshermanmorrison.c new file mode 100644 index 00000000..f6f13c9d --- /dev/null +++ b/enterprise/fastshermanmorrison/fastshermanmorrison.c @@ -0,0 +1,233 @@ +/* cython_fastshermor.c + * + * Rutger van Haasteren, April 19 2023, Hannover + * + */ + +#include +#include +#include + + +extern void dgemm_(char *transa, char *transb, int *m, int *n, int *k, + double *alpha, double *a, int *lda, double *b, + int *ldb, double *beta, double *c, int *ldc); +extern void dgemv_(char *trans, int *m, int *n, double *alpha, double *a, + int *lda, double *x, int *incx, double *beta, + double *y, int *incy); +extern void dger_(int *m, int *n, double *alpha, double *x, int *incx, + double *y, int *incy, double *a, int *lda); + +static void blas_block_shermor_2D_asym( + int n_Z1_rows, + int n_Z1_cols, + int n_Z1_row_major, + double *pd_Z1, + int n_Z2_cols, + int n_Z2_row_major, + double *pd_Z2, + double *pd_Nvec, + int n_J_rows, + double *pd_Jvec, + int *pn_Uinds, + double *pd_ZNZ, + double *pd_Jldet + ) { + /* C implementation of python_block_shermor_2D_asym, because the python + * overhead is large + * + * parameters + * ---------- + * + * :param n_Z1_rows: Number of rows of Z1 + * :param n_Z1_cols: Number of columns of Z1 + * :param n_Z1_row_major: 1 if Z1 is Row-Major, 0 if Column-Major + * :param pd_Z1: The Z1 matrix + * :param n_Z2_cols: Number of columns of Z2 + * :param n_Z2_row_major: 1 if Z2 is Row-Major, 0 if Column-Major + * :param pd_Z2: The Z2 matrix + * :param pd_Nvec: The Nvec vector + * :param n_J_rows: The number of Jvec elements + * :param pd_Jvec: The Jvec vector + * :param pn_Uinds: The matrix of quantization indices (Row-Major) + * :param pd_ZNZ: The return value of ZNZ (Column-Major) + * :param pd_Jldet: The return value of log(det(J)) + */ + + double d_galpha=1.0, d_gbeta=0.0, d_nisum=0.0, d_beta; + double *pd_Z1ni, *pd_ZNZ_add, *pd_ni, *pd_zn1, *pd_zn2; + int cc, i, j, m, n, k, lda, ldb, ldc, n_jblock, n_jblock_i, n_index; + char *transa, *transb; + + pd_Z1ni = malloc(n_Z1_rows*n_Z1_cols * sizeof(double)); + pd_ZNZ_add = calloc(n_Z1_rows*n_Z1_cols, sizeof(double)); + pd_ni = malloc(n_Z1_rows * sizeof(double)); + pd_zn1 = calloc(n_Z1_cols, sizeof(double)); + pd_zn2 = calloc(n_Z2_cols, sizeof(double)); + + /* openmp this? */ + for(i=0; i 0.0) { + + /* Note: pn_Uinds is row-major */ + d_nisum = 0.0; + n_jblock_i = pn_Uinds[2*cc]; + n_jblock = pn_Uinds[2*cc+1] - pn_Uinds[2*cc]; + for(i=pn_Uinds[2*cc]; i 1: + rblock = x[slc] + niblock = 1 / self._nvec[slc] + beta = 1.0 / (np.einsum("i->", niblock) + 1.0 / jv) + Nx[slc] -= beta * np.dot(niblock, rblock) * niblock + return Nx + + def _solve_1D1(self, x, y): + """Solves :math:`y^T N^{-1}x`, where :math:`x` and + :math:`y` are vectors. + """ + + Nx = x / self._nvec + yNx = np.dot(y, Nx) + for slc, jv in zip(self._slices, self._jvec): + if slc.stop - slc.start > 1: + xblock = x[slc] + yblock = y[slc] + niblock = 1 / self._nvec[slc] + beta = 1.0 / (np.einsum("i->", niblock) + 1.0 / jv) + yNx -= beta * np.dot(niblock, xblock) * np.dot(niblock, yblock) + return yNx + + def _solve_2D2(self, X, Z): + """Solves :math:`Z^T N^{-1}X`, where :math:`X` + and :math:`Z` are 2-d arrays. + """ + + ZNX = np.dot(Z.T / self._nvec, X) + for slc, jv in zip(self._slices, self._jvec): + if slc.stop - slc.start > 1: + Zblock = Z[slc, :] + Xblock = X[slc, :] + niblock = 1 / self._nvec[slc] + beta = 1.0 / (np.einsum("i->", niblock) + 1.0 / jv) + zn = np.dot(niblock, Zblock) + xn = np.dot(niblock, Xblock) + ZNX -= beta * np.outer(zn.T, xn) + return ZNX + + def _get_logdet(self): + """Returns log determinant of :math:`N+UJU^{T}` where :math:`U` + is a quantization matrix. + """ + logdet = np.einsum("i->", np.log(self._nvec)) + for slc, jv in zip(self._slices, self._jvec): + if slc.stop - slc.start > 1: + niblock = 1 / self._nvec[slc] + beta = 1.0 / (np.einsum("i->", niblock) + 1.0 / jv) + logdet += np.log(jv) - np.log(beta) + return logdet + + def solve(self, other, left_array=None, logdet=False): + if other.ndim == 1: + if left_array is None: + ret = self._solve_D1(other) + elif left_array is not None and left_array.ndim == 1: + ret = self._solve_1D1(other, left_array) + elif left_array is not None and left_array.ndim == 2: + ret = np.dot(left_array.T, self._solve_D1(other)) + else: + raise TypeError + elif other.ndim == 2: + if left_array is None: + raise NotImplementedError("ShermanMorrison does not implement _solve_D2") + elif left_array is not None and left_array.ndim == 2: + ret = self._solve_2D2(other, left_array) + elif left_array is not None and left_array.ndim == 1: + ret = np.dot(other.T, self._solve_D1(left_array)) + else: + raise TypeError + else: + raise TypeError + + return (ret, self._get_logdet()) if logdet else ret + +class FastShermanMorrison(ShermanMorrison): + """Custom container class for Sherman-morrison array inversion.""" + + def __init__(self, jvec, slices, nvec=0.0): + self._uinds = np.vstack([[slc.start, slc.stop] for slc in slices]) + super().__init__(jvec, slices, nvec=nvec) + + def __add__(self, other): + nvec = self._nvec + other + return FastShermanMorrison(self._jvec, self._slices, nvec) + + def _solve_D1(self, x): + """Solves :math:`N^{-1}x` where :math:`x` is a vector.""" + return cfsm.cython_block_shermor_0D(x, self._nvec, self._jvec, self._uinds) + + def _solve_1D1(self, x, y): + """Solves :math:`y^T N^{-1}x`, where :math:`x` and + :math:`y` are vectors. + """ + logJdet, yNx = cfsm.cython_block_shermor_1D1(x, y, self._nvec, self._jvec, self._uinds) + return yNx + + def _solve_2D2(self, X, Z): + """Solves :math:`Z^T N^{-1}X`, where :math:`X` + and :math:`Z` are 2-d arrays. + """ + logJdet, ZNX = cfsm.cython_blas_block_shermor_2D_asymm(Z, X, self._nvec, self._jvec, self._uinds) + return ZNX + + def _get_logdet(self): + """Returns log determinant of :math:`N+UJU^{T}` where :math:`U` + is a quantization matrix. + """ + logJdet, xNx = cfsm.cython_block_shermor_1D(np.zeros_like(self._nvec), self._nvec, self._jvec, self._uinds) + return logJdet + + def solve(self, other, left_array=None, logdet=False): + if other.ndim == 1: + if left_array is None: + ret = self._solve_D1(other) + elif left_array is not None and left_array.ndim == 1: + ret = self._solve_1D1(other, left_array) + elif left_array is not None and left_array.ndim == 2: + ret = np.dot(left_array.T, self._solve_D1(other)) + else: + raise TypeError + elif other.ndim == 2: + if left_array is None: + raise NotImplementedError("ShermanMorrison does not implement _solve_D2") + elif left_array is not None and left_array.ndim == 2: + ret = self._solve_2D2(other, left_array) + elif left_array is not None and left_array.ndim == 1: + ret = np.dot(other.T, self._solve_D1(left_array)) + else: + raise TypeError + else: + raise TypeError + + return (ret, self._get_logdet()) if logdet else ret + diff --git a/setup.py b/setup.py index 2e5993a3..a04085da 100644 --- a/setup.py +++ b/setup.py @@ -2,6 +2,8 @@ # -*- coding: utf-8 -*- from setuptools import setup +from setuptools import Extension +from Cython.Build import cythonize with open("README.md", encoding="utf-8") as readme_file: readme = readme_file.read() @@ -14,6 +16,14 @@ "scikit-sparse>=0.4.5", "pint-pulsar>=0.8.3", "libstempo>=2.4.4", + "cython>=0.29.34", +] + +ext_modules=[ + Extension('enterprise.fastshermanmorrison.cython_fastshermanmorrison', + ['enterprise/fastshermanmorrison/cython_fastshermanmorrison.pyx'], + include_dirs = [numpy.get_include(), 'fastshermanmorrison/'], + extra_compile_args=["-O2", "-fno-wrapv"]) # 50% more efficient! ] test_requirements = [] @@ -26,7 +36,7 @@ author="Justin A. Ellis", author_email="justin.ellis18@gmail.com", url="https://github.com/nanograv/enterprise", - packages=["enterprise", "enterprise.signals"], + packages=["enterprise", "enterprise.signals", "enterprise.fastshermanmorrison"], package_dir={"enterprise": "enterprise"}, include_package_data=True, package_data={"enterprise": ["datafiles/*", "datafiles/ephemeris/*", "datafiles/ng9/*", "datafiles/mdc_open1/*"]}, From 312093b7128c9591b4ca0378603e6a3bb8853ec5 Mon Sep 17 00:00:00 2001 From: Rutger van Haasteren Date: Fri, 21 Apr 2023 18:11:19 +0200 Subject: [PATCH 02/20] Added numpy and cython extensions to setup.py --- setup.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/setup.py b/setup.py index a04085da..fb2efe85 100644 --- a/setup.py +++ b/setup.py @@ -1,6 +1,7 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- +import numpy from setuptools import setup from setuptools import Extension from Cython.Build import cythonize @@ -61,4 +62,5 @@ ], test_suite="tests", tests_require=test_requirements, + ext_modules = cythonize(ext_modules) ) From 1b536bc39f3068fed0fa83a38f7f3429a9984669 Mon Sep 17 00:00:00 2001 From: Rutger van Haasteren Date: Mon, 24 Apr 2023 10:40:44 +0200 Subject: [PATCH 03/20] Updated docstrings and method selector for fast-sherman-morrison --- .../fastshermanmorrison.py | 107 +----------------- enterprise/signals/white_signals.py | 20 +++- 2 files changed, 21 insertions(+), 106 deletions(-) diff --git a/enterprise/fastshermanmorrison/fastshermanmorrison.py b/enterprise/fastshermanmorrison/fastshermanmorrison.py index db558403..5a374e24 100644 --- a/enterprise/fastshermanmorrison/fastshermanmorrison.py +++ b/enterprise/fastshermanmorrison/fastshermanmorrison.py @@ -1,108 +1,9 @@ import numpy as np from . import cython_fastshermanmorrison as cfsm +from enterprise.signals import signal_base -class ShermanMorrison(object): - """Custom container class for Sherman-morrison array inversion.""" - - def __init__(self, jvec, slices, nvec=0.0): - self._jvec = jvec - self._slices = slices - self._nvec = nvec - - def __add__(self, other): - nvec = self._nvec + other - return ShermanMorrison(self._jvec, self._slices, nvec) - - # hacky way to fix adding 0 - def __radd__(self, other): - if other == 0: - return self.__add__(other) - else: - raise TypeError - - def _solve_D1(self, x): - """Solves :math:`N^{-1}x` where :math:`x` is a vector.""" - - Nx = x / self._nvec - for slc, jv in zip(self._slices, self._jvec): - if slc.stop - slc.start > 1: - rblock = x[slc] - niblock = 1 / self._nvec[slc] - beta = 1.0 / (np.einsum("i->", niblock) + 1.0 / jv) - Nx[slc] -= beta * np.dot(niblock, rblock) * niblock - return Nx - - def _solve_1D1(self, x, y): - """Solves :math:`y^T N^{-1}x`, where :math:`x` and - :math:`y` are vectors. - """ - - Nx = x / self._nvec - yNx = np.dot(y, Nx) - for slc, jv in zip(self._slices, self._jvec): - if slc.stop - slc.start > 1: - xblock = x[slc] - yblock = y[slc] - niblock = 1 / self._nvec[slc] - beta = 1.0 / (np.einsum("i->", niblock) + 1.0 / jv) - yNx -= beta * np.dot(niblock, xblock) * np.dot(niblock, yblock) - return yNx - - def _solve_2D2(self, X, Z): - """Solves :math:`Z^T N^{-1}X`, where :math:`X` - and :math:`Z` are 2-d arrays. - """ - - ZNX = np.dot(Z.T / self._nvec, X) - for slc, jv in zip(self._slices, self._jvec): - if slc.stop - slc.start > 1: - Zblock = Z[slc, :] - Xblock = X[slc, :] - niblock = 1 / self._nvec[slc] - beta = 1.0 / (np.einsum("i->", niblock) + 1.0 / jv) - zn = np.dot(niblock, Zblock) - xn = np.dot(niblock, Xblock) - ZNX -= beta * np.outer(zn.T, xn) - return ZNX - - def _get_logdet(self): - """Returns log determinant of :math:`N+UJU^{T}` where :math:`U` - is a quantization matrix. - """ - logdet = np.einsum("i->", np.log(self._nvec)) - for slc, jv in zip(self._slices, self._jvec): - if slc.stop - slc.start > 1: - niblock = 1 / self._nvec[slc] - beta = 1.0 / (np.einsum("i->", niblock) + 1.0 / jv) - logdet += np.log(jv) - np.log(beta) - return logdet - - def solve(self, other, left_array=None, logdet=False): - if other.ndim == 1: - if left_array is None: - ret = self._solve_D1(other) - elif left_array is not None and left_array.ndim == 1: - ret = self._solve_1D1(other, left_array) - elif left_array is not None and left_array.ndim == 2: - ret = np.dot(left_array.T, self._solve_D1(other)) - else: - raise TypeError - elif other.ndim == 2: - if left_array is None: - raise NotImplementedError("ShermanMorrison does not implement _solve_D2") - elif left_array is not None and left_array.ndim == 2: - ret = self._solve_2D2(other, left_array) - elif left_array is not None and left_array.ndim == 1: - ret = np.dot(other.T, self._solve_D1(left_array)) - else: - raise TypeError - else: - raise TypeError - - return (ret, self._get_logdet()) if logdet else ret - -class FastShermanMorrison(ShermanMorrison): - """Custom container class for Sherman-morrison array inversion.""" +class FastShermanMorrison(signal_base.ShermanMorrison): + """Custom container class for Fast-Sherman-morrison array inversion.""" def __init__(self, jvec, slices, nvec=0.0): self._uinds = np.vstack([[slc.start, slc.stop] for slc in slices]) @@ -149,7 +50,7 @@ def solve(self, other, left_array=None, logdet=False): raise TypeError elif other.ndim == 2: if left_array is None: - raise NotImplementedError("ShermanMorrison does not implement _solve_D2") + raise NotImplementedError("FastShermanMorrison does not implement _solve_D2") elif left_array is not None and left_array.ndim == 2: ret = self._solve_2D2(other, left_array) elif left_array is not None and left_array.ndim == 1: diff --git a/enterprise/signals/white_signals.py b/enterprise/signals/white_signals.py index a7d794b0..77e077fb 100644 --- a/enterprise/signals/white_signals.py +++ b/enterprise/signals/white_signals.py @@ -8,6 +8,7 @@ import scipy.sparse from enterprise.signals import parameter, selections, signal_base, utils +from enterprise.fastshermanmorrison import fastshermanmorrison from enterprise.signals.parameter import function from enterprise.signals.selections import Selection @@ -114,7 +115,7 @@ def EquadNoise(*args, **kwargs): def EcorrKernelNoise( log10_ecorr=parameter.Uniform(-10, -5), selection=Selection(selections.no_selection), - method="sherman-morrison", + method="fast-sherman-morrison", name="", ): r"""Class factory for ECORR type noise. @@ -123,7 +124,8 @@ def EcorrKernelNoise( :param selection: ``Selection`` object specifying masks for backends, time segments, etc. :param method: Method for computing noise covariance matrix. - Options include `sherman-morrison`, `sparse`, and `block` + Options include `fast-sherman-morrison`, `sherman-morrison`, `sparse`, + and `block` :return: ``EcorrKernelNoise`` class. @@ -140,6 +142,12 @@ def EcorrKernelNoise( In this signal implementation we offer three methods of performing these matrix operations: + fast-sherman-morrison + Uses the `Sherman-Morrison`_ forumla to compute the matrix + inverse and other matrix operations. **Note:** This method can only + be used for covariances that make up ECorrKernelNoise, :math:`uv^T`. + This version is Cython optimized. + sherman-morrison Uses the `Sherman-Morrison`_ forumla to compute the matrix inverse and other matrix operations. **Note:** This method can only @@ -166,7 +174,7 @@ def EcorrKernelNoise( """ - if method not in ["sherman-morrison", "block", "sparse"]: + if method not in ["fast-sherman-morrison", "sherman-morrison", "block", "sparse"]: msg = "EcorrKernelNoise does not support method: {}".format(method) raise TypeError(msg) @@ -210,6 +218,8 @@ def ndiag_params(self): def get_ndiag(self, params): if method == "sherman-morrison": return self._get_ndiag_sherman_morrison(params) + elif method == "fast-sherman-morrison": + return self._get_ndiag_sparse(params) elif method == "sparse": return self._get_ndiag_sparse(params) elif method == "block": @@ -238,6 +248,10 @@ def _get_ndiag_sherman_morrison(self, params): slices, jvec = self._get_jvecs(params) return signal_base.ShermanMorrison(jvec, slices) + def _get_ndiag_fast_sherman_morrison(self, params): + slices, jvec = self._get_jvecs(params) + return fastshermanmorrison.FastShermanMorrison(jvec, slices) + def _get_ndiag_block(self, params): slices, jvec = self._get_jvecs(params) blocks = [] From db81f64e77eba74f633d35315831ed803da3e91a Mon Sep 17 00:00:00 2001 From: Rutger van Haasteren Date: Mon, 24 Apr 2023 16:28:20 +0200 Subject: [PATCH 04/20] Added numpy and cython to requirements_dev.txt --- requirements_dev.txt | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/requirements_dev.txt b/requirements_dev.txt index 73ae670a..86d0cbb2 100644 --- a/requirements_dev.txt +++ b/requirements_dev.txt @@ -14,4 +14,6 @@ sphinx-rtd-theme>=0.4.0 pytest-cov>=2.7.0 coverage-conditional-plugin>=0.4.0 jupyter>=1.0.0 -build==0.3.1.post1 \ No newline at end of file +build==0.3.1.post1 +numpy>=1.16.3 +cython>=0.29.34 From d8e1a3559ccaedc47aa140a1dd6243d713e717ee Mon Sep 17 00:00:00 2001 From: Rutger van Haasteren Date: Mon, 24 Apr 2023 17:44:02 +0200 Subject: [PATCH 05/20] Moved FastShermanMorrison to external package --- enterprise/fastshermanmorrison/__init__.py | 1 - .../cython_fastshermanmorrison.pyx | 945 ------------------ .../fastshermanmorrison/fastshermanmorrison.c | 233 ----- .../fastshermanmorrison.py | 64 -- enterprise/signals/white_signals.py | 18 +- requirements_dev.txt | 2 - setup.py | 12 - 7 files changed, 16 insertions(+), 1259 deletions(-) delete mode 100644 enterprise/fastshermanmorrison/__init__.py delete mode 100644 enterprise/fastshermanmorrison/cython_fastshermanmorrison.pyx delete mode 100644 enterprise/fastshermanmorrison/fastshermanmorrison.c delete mode 100644 enterprise/fastshermanmorrison/fastshermanmorrison.py diff --git a/enterprise/fastshermanmorrison/__init__.py b/enterprise/fastshermanmorrison/__init__.py deleted file mode 100644 index 3037395f..00000000 --- a/enterprise/fastshermanmorrison/__init__.py +++ /dev/null @@ -1 +0,0 @@ -from . import fastshermanmorrison diff --git a/enterprise/fastshermanmorrison/cython_fastshermanmorrison.pyx b/enterprise/fastshermanmorrison/cython_fastshermanmorrison.pyx deleted file mode 100644 index c8a3b8c2..00000000 --- a/enterprise/fastshermanmorrison/cython_fastshermanmorrison.pyx +++ /dev/null @@ -1,945 +0,0 @@ -cimport numpy as np -import numpy as np -np.import_array() - -from libc.math cimport log, sqrt -import cython -from scipy.linalg.cython_blas cimport dgemm, dger, dgemv - -cdef public void dgemm_(char *transa, char *transb, int *m, int *n, int *k, - double *alpha, double *a, int *lda, double *b, - int *ldb, double *beta, double *c, int *ldc): - """Public dgemm that can be used in the external C code""" - dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc) - -cdef public void dgemv_(char *trans, int *m, int *n, double *alpha, double *a, - int *lda, double *x, int *incx, double *beta, - double *y, int *incy): - """Public dgemv that can be used in the external C code""" - dgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy) - -cdef public void dger_(int *m, int *n, double *alpha, double *x, int *incx, - double *y, int *incy, double *a, int *lda): - """Public dger that can be used in the external C code""" - dger(m, n, alpha, x, incx, y, incy, a, lda) - - -cdef extern from "fastshermanmorrison.c": - void blas_block_shermor_2D_asym( - int n_Z1_rows, - int n_Z1_cols, - int n_Z1_row_major, - double *pd_Z1, - int n_Z2_cols, - int n_Z2_row_major, - double *pd_Z2, - double *pd_Nvec, - int n_J_rows, - double *pd_Jvec, - int *pn_Uinds, - double *pd_ZNZ, - double *pd_Jldet - ) - -@cython.boundscheck(False) -@cython.wraparound(False) -def cython_block_shermor_0D( \ - np.ndarray[np.double_t,ndim=1] r, \ - np.ndarray[np.double_t,ndim=1] Nvec, \ - np.ndarray[np.double_t,ndim=1] Jvec, \ - np.ndarray[np.int_t,ndim=2] Uinds): - """ - Sherman-Morrison block-inversion for Jitter (Cythonized) - - @param r: The timing residuals, array (n) - @param Nvec: The white noise amplitude, array (n) - @param Jvec: The jitter amplitude, array (k) - @param Uinds: The start/finish indices for the jitter blocks (k x 2) - - For this version, the residuals need to be sorted properly so that all the - blocks are continuous in memory. Here, there are n residuals, and k jitter - parameters. - """ - cdef unsigned int cc, ii, cols = len(Jvec) - cdef double Jldet=0.0, ji, beta, nir, nisum - cdef np.ndarray[np.double_t,ndim=1] ni = np.empty(len(r), 'd') - cdef np.ndarray[np.double_t,ndim=1] Nx = r / Nvec - - ni = 1.0 / Nvec - - for cc in range(cols): - if Jvec[cc] > 0.0: - ji = 1.0 / Jvec[cc] - - nir = 0.0 - nisum = 0.0 - for ii in range(Uinds[cc,0],Uinds[cc,1]): - nisum += ni[ii] - nir += r[ii]*ni[ii] - - beta = 1.0 / (nisum + ji) - - for ii in range(Uinds[cc,0],Uinds[cc,1]): - Nx[ii] -= beta * nir * ni[ii] - - return Nx - -@cython.boundscheck(False) -@cython.wraparound(False) -def cython_block_shermor_0D_ld( \ - np.ndarray[np.double_t,ndim=1] r, \ - np.ndarray[np.double_t,ndim=1] Nvec, \ - np.ndarray[np.double_t,ndim=1] Jvec, \ - np.ndarray[np.int_t,ndim=2] Uinds): - """ - Sherman-Morrison block-inversion for Jitter (Cythonized) - - @param r: The timing residuals, array (n) - @param Nvec: The white noise amplitude, array (n) - @param Jvec: The jitter amplitude, array (k) - @param Uinds: The start/finish indices for the jitter blocks (k x 2) - - For this version, the residuals need to be sorted properly so that all the - blocks are continuous in memory. Here, there are n residuals, and k jitter - parameters. - """ - cdef unsigned int cc, ii, rows = len(r), cols = len(Jvec) - cdef double Jldet=0.0, ji, beta, nir, nisum - cdef np.ndarray[np.double_t,ndim=1] ni = np.empty(len(r), 'd') - cdef np.ndarray[np.double_t,ndim=1] Nx = r / Nvec - - ni = 1.0 / Nvec - - for cc in range(rows): - Jldet += log(Nvec[cc]) - - for cc in range(cols): - if Jvec[cc] > 0.0: - ji = 1.0 / Jvec[cc] - - nir = 0.0 - nisum = 0.0 - for ii in range(Uinds[cc,0],Uinds[cc,1]): - nisum += ni[ii] - nir += r[ii]*ni[ii] - - beta = 1.0 / (nisum + ji) - Jldet += log(Jvec[cc]) - log(beta) - - for ii in range(Uinds[cc,0],Uinds[cc,1]): - Nx[ii] -= beta * nir * ni[ii] - - return Jldet, Nx - - -def python_block_shermor_1D(r, Nvec, Jvec, Uinds): - """ - Sherman-Morrison block-inversion for Jitter - - @param r: The timing residuals, array (n) - @param Nvec: The white noise amplitude, array (n) - @param Jvec: The jitter amplitude, array (k) - @param Uinds: The start/finish indices for the jitter blocks (k x 2) - - For this version, the residuals need to be sorted properly so that all the - blocks are continuous in memory. Here, there are n residuals, and k jitter - parameters. - """ - ni = 1.0 / Nvec - Jldet = np.einsum('i->', np.log(Nvec)) - xNx = np.dot(r, r * ni) - - for cc, jv in enumerate(Jvec): - if jv > 0.0: - rblock = r[Uinds[cc,0]:Uinds[cc,1]] - niblock = ni[Uinds[cc,0]:Uinds[cc,1]] - - beta = 1.0 / (np.einsum('i->', niblock)+1.0/jv) - xNx -= beta * np.dot(rblock, niblock)**2 - Jldet += np.log(jv) - np.log(beta) - - return Jldet, xNx - -@cython.boundscheck(False) -@cython.wraparound(False) -def cython_block_shermor_1D( \ - np.ndarray[np.double_t,ndim=1] r, \ - np.ndarray[np.double_t,ndim=1] Nvec, \ - np.ndarray[np.double_t,ndim=1] Jvec, \ - np.ndarray[np.int_t,ndim=2] Uinds): - """ - Sherman-Morrison block-inversion for Jitter (Cythonized) - - @param r: The timing residuals, array (n) - @param Nvec: The white noise amplitude, array (n) - @param Jvec: The jitter amplitude, array (k) - @param Uinds: The start/finish indices for the jitter blocks (k x 2) - - For this version, the residuals need to be sorted properly so that all the - blocks are continuous in memory. Here, there are n residuals, and k jitter - parameters. - """ - cdef unsigned int cc, ii, rows = len(r), cols = len(Jvec) - cdef double Jldet=0.0, ji, beta, xNx=0.0, nir, nisum - cdef np.ndarray[np.double_t,ndim=1] ni = np.empty(rows, 'd') - - ni = 1.0 / Nvec - - for cc in range(rows): - Jldet += log(Nvec[cc]) - xNx += r[cc]*r[cc]*ni[cc] - - for cc in range(cols): - if Jvec[cc] > 0.0: - ji = 1.0 / Jvec[cc] - - nir = 0.0 - nisum = 0.0 - for ii in range(Uinds[cc,0],Uinds[cc,1]): - nisum += ni[ii] - nir += r[ii]*ni[ii] - - beta = 1.0 / (nisum + ji) - Jldet += log(Jvec[cc]) - log(beta) - xNx -= beta * nir * nir - - return Jldet, xNx - -@cython.boundscheck(False) -@cython.wraparound(False) -def cython_block_shermor_1D1( \ - np.ndarray[np.double_t,ndim=1] x, \ - np.ndarray[np.double_t,ndim=1] y, \ - np.ndarray[np.double_t,ndim=1] Nvec, \ - np.ndarray[np.double_t,ndim=1] Jvec, \ - np.ndarray[np.int_t,ndim=2] Uinds): - """ - Sherman-Morrison block-inversion for Jitter (Cythonized) - - @param r: The timing residuals, array (n) - @param Nvec: The white noise amplitude, array (n) - @param Jvec: The jitter amplitude, array (k) - @param Uinds: The start/finish indices for the jitter blocks (k x 2) - - For this version, the residuals need to be sorted properly so that all the - blocks are continuous in memory. Here, there are n residuals, and k jitter - parameters. - """ - cdef unsigned int cc, ii, rows = len(x), cols = len(Jvec) - cdef double Jldet=0.0, ji, beta, yNx=0.0, nix, niy, nisum - cdef np.ndarray[np.double_t,ndim=1] ni = np.empty(rows, 'd') - - ni = 1.0 / Nvec - - for cc in range(rows): - Jldet += log(Nvec[cc]) - yNx += y[cc]*x[cc]*ni[cc] - - for cc in range(cols): - if Jvec[cc] > 0.0: - ji = 1.0 / Jvec[cc] - - nix = 0.0 - niy = 0.0 - nisum = 0.0 - for ii in range(Uinds[cc,0],Uinds[cc,1]): - nisum += ni[ii] - nix += x[ii]*ni[ii] - niy += y[ii]*ni[ii] - - beta = 1.0 / (nisum + ji) - Jldet += log(Jvec[cc]) - log(beta) - yNx -= beta * nix * niy - - return Jldet, yNx - - -def python_block_shermor_2D(Z, Nvec, Jvec, Uinds): - """ - Sherman-Morrison block-inversion for Jitter, ZNiZ - - @param Z: The design matrix, array (n x m) - @param Nvec: The white noise amplitude, array (n) - @param Jvec: The jitter amplitude, array (k) - @param Uinds: The start/finish indices for the jitter blocks (k x 2) - - For this version, the residuals need to be sorted properly so that all the - blocks are continuous in memory. Here, there are n residuals, and k jitter - parameters. - - N = D + U*J*U.T - calculate: log(det(N)), Z.T * N^-1 * Z - """ - ni = 1.0 / Nvec - Jldet = np.einsum('i->', np.log(Nvec)) - zNz = np.dot(Z.T*ni, Z) - - for cc, jv in enumerate(Jvec): - if jv > 0.0: - Zblock = Z[Uinds[cc,0]:Uinds[cc,1], :] - niblock = ni[Uinds[cc,0]:Uinds[cc,1]] - - beta = 1.0 / (np.einsum('i->', niblock)+1.0/jv) - zn = np.dot(niblock, Zblock) - zNz -= beta * np.outer(zn.T, zn) - Jldet += np.log(jv) - np.log(beta) - - return Jldet, zNz - -@cython.boundscheck(False) -@cython.wraparound(False) -def cython_block_shermor_2D( \ - np.ndarray[np.double_t,ndim=2] Z, \ - np.ndarray[np.double_t,ndim=1] Nvec, \ - np.ndarray[np.double_t,ndim=1] Jvec, \ - np.ndarray[np.int_t,ndim=2] Uinds): - """ - Sherman-Morrison block-inversion for Jitter (Cythonized) - - @param Z: The design matrix, array (n x m) - @param Nvec: The white noise amplitude, array (n) - @param Jvec: The jitter amplitude, array (k) - @param Uinds: The start/finish indices for the jitter blocks (k x 2) - - For this version, the residuals need to be sorted properly so that all the - blocks are continuous in memory. Here, there are n residuals, and k jitter - parameters. - - N = D + U*J*U.T - calculate: log(det(N)), Z.T * N^-1 * Z - """ - cdef unsigned int cc, ii, rows = len(Nvec), cols = len(Jvec) - cdef double Jldet=0.0, ji, beta, nir, nisum - cdef np.ndarray[np.double_t,ndim=1] ni = np.empty(len(Nvec), 'd') - cdef np.ndarray[np.double_t,ndim=2] zNz - - ni = 1.0 / Nvec - zNz = np.dot(Z.T*ni, Z) - - for cc in range(rows): - Jldet += log(Nvec[cc]) - - for cc in range(cols): - if Jvec[cc] > 0.0: - Zblock = Z[Uinds[cc,0]:Uinds[cc,1], :] - niblock = ni[Uinds[cc,0]:Uinds[cc,1]] - - nisum = 0.0 - for ii in range(len(niblock)): - nisum += niblock[ii] - - beta = 1.0 / (nisum+1.0/Jvec[cc]) - Jldet += log(Jvec[cc]) - log(beta) - zn = np.dot(niblock, Zblock) - zNz -= beta * np.outer(zn.T, zn) - - return Jldet, zNz - -def python_block_shermor_2D_asymm(Z1, Z2, Nvec, Jvec, Uinds): - """ - Sherman-Morrison block-inversion for Jitter, ZNiZ - - @param Z: The design matrix, array (n x m) - @param Z2: The second design matrix, array (n x m2) - @param Nvec: The white noise amplitude, array (n) - @param Jvec: The jitter amplitude, array (k) - @param Uinds: The start/finish indices for the jitter blocks (k x 2) - - For this version, the residuals need to be sorted properly so that all the - blocks are continuous in memory. Here, there are n residuals, and k jitter - parameters. - - N = D + U*J*U.T - calculate: log(det(N)), Z.T * N^-1 * Z - """ - ni = 1.0 / Nvec - Jldet = np.einsum('i->', np.log(Nvec)) - zNz = np.dot(Z1.T*ni, Z2) - - for cc, jv in enumerate(Jvec): - if jv > 0.0: - Zblock1 = Z1[Uinds[cc,0]:Uinds[cc,1], :] - Zblock2 = Z2[Uinds[cc,0]:Uinds[cc,1], :] - niblock = ni[Uinds[cc,0]:Uinds[cc,1]] - - beta = 1.0 / (np.einsum('i->', niblock)+1.0/jv) - zn1 = np.dot(niblock, Zblock1) - zn2 = np.dot(niblock, Zblock2) - zNz -= beta * np.outer(zn1.T, zn2) - Jldet += np.log(jv) - np.log(beta) - - return Jldet, zNz - - -@cython.boundscheck(False) -@cython.wraparound(False) -def cython_block_shermor_2D_asymm( - np.ndarray[np.double_t,ndim=2] Z1, - np.ndarray[np.double_t,ndim=2] Z2, - np.ndarray[np.double_t,ndim=1] Nvec, - np.ndarray[np.double_t,ndim=1] Jvec, - np.ndarray[np.int_t,ndim=2] Uinds): - """ - Sherman-Morrison block-inversion for Jitter, ZNiZ - - @param Z: The design matrix, array (n x m) - @param Z2: The second design matrix, array (n x m2) - @param Nvec: The white noise amplitude, array (n) - @param Jvec: The jitter amplitude, array (k) - @param Uinds: The start/finish indices for the jitter blocks (k x 2) - - For this version, the residuals need to be sorted properly so that all the - blocks are continuous in memory. Here, there are n residuals, and k jitter - parameters. - - N = D + U*J*U.T - calculate: log(det(N)), Z.T * N^-1 * Z - """ - cdef unsigned int cc, ii, rows = len(Nvec), cols = len(Jvec) - cdef double Jldet=0.0, ji, beta, nir, nisum - cdef np.ndarray[np.double_t,ndim=1] ni = np.empty(len(Nvec), 'd') - cdef np.ndarray[np.double_t,ndim=2] zNz - - print("WARNING: cython_block_shermor_2D_asymm is deprecated.") - print(" use cython_blas_block_shermor_2D_asymm") - - ni = 1.0 / Nvec - for cc in range(rows): - Jldet += log(Nvec[cc]) - zNz = np.dot(Z1.T*ni, Z2) - - for cc in range(cols): - if Jvec[cc] > 0.0: - Zblock1 = Z1[Uinds[cc,0]:Uinds[cc,1], :] - Zblock2 = Z2[Uinds[cc,0]:Uinds[cc,1], :] - niblock = ni[Uinds[cc,0]:Uinds[cc,1]] - - nisum = 0.0 - for ii in range(len(niblock)): - nisum += niblock[ii] - - beta = 1.0 / (nisum+1.0/Jvec[cc]) - zn1 = np.dot(niblock, Zblock1) - zn2 = np.dot(niblock, Zblock2) - - zNz -= beta * np.outer(zn1.T, zn2) - Jldet += log(Jvec[cc]) - log(beta) - - return Jldet, zNz - - -def python_draw_ecor(r, Nvec, Jvec, Uinds): - """ - Given Jvec, draw new epoch-averaged residuals - - @param r: The timing residuals, array (n) - @param Nvec: The white noise amplitude, array (n) - @param Jvec: The jitter amplitude, array (k) - @param Uinds: The start/finish indices for the jitter blocks (k x 2) - - For this version, the residuals need to be sorted properly so that all the - blocks are continuous in memory. Here, there are n residuals, and k jitter - parameters. - - N = D + U*J*U.T - calculate: Norm(0, sqrt(J)) + (U^T * D^{-1} * U)^{-1}U.T D^{-1} r - """ - - rv = np.random.randn(len(Jvec)) * np.sqrt(Jvec) - ni = 1.0 / Nvec - - for cc in range(len(Jvec)): - rblock = r[Uinds[cc,0]:Uinds[cc,1]] - niblock = ni[Uinds[cc,0]:Uinds[cc,1]] - beta = 1.0 / np.einsum('i->', niblock) - - rv[cc] += beta * np.dot(rblock, niblock) - - return rv - -@cython.boundscheck(False) -@cython.wraparound(False) -def cython_draw_ecor( \ - np.ndarray[np.double_t,ndim=1] r, \ - np.ndarray[np.double_t,ndim=1] Nvec, \ - np.ndarray[np.double_t,ndim=1] Jvec, \ - np.ndarray[np.int_t,ndim=2] Uinds): - """ - Given Jvec, draw new epoch-averaged residuals - - @param r: The timing residuals, array (n) - @param Nvec: The white noise amplitude, array (n) - @param Jvec: The jitter amplitude, array (k) - @param Uinds: The start/finish indices for the jitter blocks (k x 2) - - For this version, the residuals need to be sorted properly so that all the - blocks are continuous in memory. Here, there are n residuals, and k jitter - parameters. - - N = D + U*J*U.T - calculate: Norm(0, sqrt(J)) + (U^T * D^{-1} * U)^{-1}U.T D^{-1} r - """ - cdef unsigned int cc, ii, rows = len(r), cols = len(Jvec) - cdef double ji, nir, nisum - cdef np.ndarray[np.double_t,ndim=1] ni = np.empty(rows, 'd') - cdef np.ndarray[np.double_t,ndim=1] rv = np.random.randn(cols) - - for cc in range(cols): - rv[cc] *= sqrt(Jvec[cc]) - - ni = 1.0 / Nvec - - for cc in range(cols): - ji = 1.0 / Jvec[cc] - - nir = 0.0 - nisum = 0.0 - for ii in range(Uinds[cc,0],Uinds[cc,1]): - nisum += ni[ii] - nir += r[ii]*ni[ii] - - rv[cc] += nir / nisum - - return rv - - -@cython.boundscheck(False) -@cython.wraparound(False) -def cython_shermor_draw_ecor( \ - np.ndarray[np.double_t,ndim=1] r, \ - np.ndarray[np.double_t,ndim=1] Nvec, \ - np.ndarray[np.double_t,ndim=1] Jvec, \ - np.ndarray[np.int_t,ndim=2] Uinds): - """ - Do both the Sherman-Morrison block-inversion for Jitter, - and the draw of the ecor parameters together (Cythonized) - - @param r: The timing residuals, array (n) - @param Nvec: The white noise amplitude, array (n) - @param Jvec: The jitter amplitude, array (k) - @param Uinds: The start/finish indices for the jitter blocks (k x 2) - - For this version, the residuals need to be sorted properly so that all the - blocks are continuous in memory. Here, there are n residuals, and k jitter - parameters. - - N = D + U*J*U.T - calculate: r.T * N^-1 * r, log(det(N)), Norm(0, sqrt(J)) + (U^T * D^{-1} * U)^{-1}U.T D^{-1} r - """ - cdef unsigned int cc, ii, rows = len(r), cols = len(Jvec) - cdef double Jldet=0.0, ji, beta, xNx=0.0, nir, nisum - cdef np.ndarray[np.double_t,ndim=1] ni = np.empty(rows, 'd') - cdef np.ndarray[np.double_t,ndim=1] rv = np.random.randn(cols) - - ni = 1.0 / Nvec - - for cc in range(cols): - rv[cc] *= sqrt(Jvec[cc]) - - for cc in range(rows): - Jldet += log(Nvec[cc]) - xNx += r[cc]*r[cc]*ni[cc] - - for cc in range(cols): - nir = 0.0 - nisum = 0.0 - for ii in range(Uinds[cc,0],Uinds[cc,1]): - nisum += ni[ii] - nir += r[ii]*ni[ii] - - rv[cc] += nir / nisum - - if Jvec[cc] > 0.0: - ji = 1.0 / Jvec[cc] - - beta = 1.0 / (nisum + ji) - Jldet += log(Jvec[cc]) - log(beta) - xNx -= beta * nir * nir - - return Jldet, xNx, rv - - -@cython.boundscheck(False) -@cython.wraparound(False) -def cython_update_ea_residuals( \ - np.ndarray[np.double_t,ndim=1] gibbsresiduals, \ - np.ndarray[np.double_t,ndim=1] gibbssubresiduals, \ - np.ndarray[np.double_t,ndim=1] eat, \ - np.ndarray[np.int_t,ndim=2] Uinds): - """ - Given epoch-averaged residuals, update the residuals, and the subtracted - residuals, so that these can be further processed by the other conditional - probability density functions. - - @param gibbsresiduals: The timing residuals, array (n) - @param gibbssubresiduals: The white noise amplitude, array (n) - @param eat: epoch averaged residuals (k) - @param Uinds: The start/finish indices for the jitter blocks - (k x 2) - - """ - cdef unsigned int k = Uinds.shape[0], ii, cc - - for cc in range(Uinds.shape[0]): - for ii in range(Uinds[cc,0],Uinds[cc,1]): - gibbssubresiduals[ii] += eat[cc] - gibbsresiduals[ii] -= eat[cc] - - return gibbsresiduals, gibbssubresiduals - - -@cython.boundscheck(False) -@cython.wraparound(False) -def cython_Uj(np.ndarray[np.double_t,ndim=1] j, \ - np.ndarray[np.int_t,ndim=2] Uinds, nobs): - """ - Given epoch-averaged residuals (j), get the residuals. - Used in 'updateDetSources' - - @param j: epoch averaged residuals (k) - @param Uinds: The start/finish indices for the jitter blocks - (k x 2) - @param nobs: Number of observations (length return vector) - - """ - cdef unsigned int k = Uinds.shape[0], ii, cc - cdef np.ndarray[np.double_t,ndim=1] Uj = np.zeros(nobs, 'd') - - for cc in range(k): - for ii in range(Uinds[cc,0],Uinds[cc,1]): - Uj[ii] += j[cc] - - return Uj - -@cython.boundscheck(False) -@cython.wraparound(False) -def cython_UTx(np.ndarray[np.double_t,ndim=1] x, \ - np.ndarray[np.int_t,ndim=2] Uinds): - """ - Given residuals (x), get np.dot(U.T, x) - Used in 'updateDetSources' - - @param j: epoch averaged residuals (k) - @param Uinds: The start/finish indices for the jitter blocks - (k x 2) - - """ - cdef unsigned int k = Uinds.shape[0], ii, cc - cdef np.ndarray[np.double_t,ndim=1] UTx = np.zeros(k, 'd') - - for cc in range(k): - for ii in range(Uinds[cc,0],Uinds[cc,1]): - UTx[cc] += x[ii] - - return UTx - -@cython.boundscheck(False) -@cython.wraparound(False) -def cython_logdet_dN( \ - np.ndarray[np.double_t,ndim=1] Nvec, \ - np.ndarray[np.double_t,ndim=1] Jvec, \ - np.ndarray[np.double_t,ndim=1] dNvec, \ - np.ndarray[np.int_t,ndim=2] Uinds): - """ - Sherman-Morrison block-inversion for Jitter (Cythonized) - - Calculates Trace(N^{-1} dN/dNp), where: - - N^{-1} is the ecorr-include N inverse - - dN/dNp is the diagonal derivate of N wrt Np - - @param Nvec: The white noise amplitude, array (n) - @param Jvec: The jitter amplitude, array (k) - @param dNvec: The white noise derivative, array (n) - @param Uinds: The start/finish indices for the jitter blocks (k x 2) - - For this version, the residuals need to be sorted properly so that all the - blocks are continuous in memory. Here, there are n residuals, and k jitter - parameters. - """ - cdef unsigned int cc, ii, rows = len(Nvec), cols = len(Jvec) - cdef double tr=0.0, ji, nisum, Nnisum - cdef np.ndarray[np.double_t,ndim=1] ni = np.empty(rows, 'd') - cdef np.ndarray[np.double_t,ndim=1] Nni = np.empty(rows, 'd') - - ni = 1.0 / Nvec - Nni = dNvec / Nvec**2 - - for cc in range(rows): - tr += dNvec[cc] * ni[cc] - - for cc in range(cols): - if Jvec[cc] > 0.0: - ji = 1.0 / Jvec[cc] - - nisum = 0.0 - Nnisum = 0.0 - for ii in range(Uinds[cc,0],Uinds[cc,1]): - nisum += ni[ii] - Nnisum += Nni[ii] - - tr -= Nnisum / (nisum + ji) - - return tr - -@cython.boundscheck(False) -@cython.wraparound(False) -def cython_logdet_dJ( \ - np.ndarray[np.double_t,ndim=1] Nvec, \ - np.ndarray[np.double_t,ndim=1] Jvec, \ - np.ndarray[np.double_t,ndim=1] dJvec, \ - np.ndarray[np.int_t,ndim=2] Uinds): - """ - Sherman-Morrison block-inversion for Jitter (Cythonized) - - Calculates Trace(N^{-1} dN/dJp), where: - - N^{-1} is the ecorr-include N inverse - - dN/dJp = U dJ/dJp U^{T}, with dJ/dJp the diagnal derivative of J wrt - Jp - - @param Nvec: The white noise amplitude, array (n) - @param Jvec: The jitter amplitude, array (k) - @param dJvec: The jitter derivative, array (k) - @param Uinds: The start/finish indices for the jitter blocks (k x 2) - - For this version, the residuals need to be sorted properly so that all the - blocks are continuous in memory. Here, there are n residuals, and k jitter - parameters. - """ - cdef unsigned int cc, ii, rows = len(Nvec), cols = len(Jvec) - cdef double dJldet=0.0, ji, beta, nisum - cdef np.ndarray[np.double_t,ndim=1] ni = np.empty(rows, 'd') - - ni = 1.0 / Nvec - - for cc in range(cols): - if Jvec[cc] > 0.0: - ji = 1.0 / Jvec[cc] - - nisum = 0.0 - for ii in range(Uinds[cc,0],Uinds[cc,1]): - nisum += ni[ii] - - beta = 1.0 / (nisum + ji) - - dJldet += dJvec[cc]*(nisum - beta*nisum**2) - - return dJldet - -@cython.boundscheck(False) -@cython.wraparound(False) -def cython_logdet_dN_dN( \ - np.ndarray[np.double_t,ndim=1] Nvec, \ - np.ndarray[np.double_t,ndim=1] Jvec, \ - np.ndarray[np.double_t,ndim=1] dNvec1, \ - np.ndarray[np.double_t,ndim=1] dNvec2, \ - np.ndarray[np.int_t,ndim=2] Uinds): - """ - Sherman-Morrison block-inversion for Jitter (Cythonized) - - Calculates Trace(N^{-1} dN/dNp1 N^{-1} dN/dNp2), where: - - N^{-1} is the ecorr-include N inverse - - dN/dNpx is the diagonal derivate of N wrt Npx - - @param Nvec: The white noise amplitude, array (n) - @param Jvec: The jitter amplitude, array (k) - @param dNvec1: The white noise derivative, array (n) - @param dNvec2: The white noise derivative, array (n) - @param Uinds: The start/finish indices for the jitter blocks (k x 2) - - For this version, the residuals need to be sorted properly so that all the - blocks are continuous in memory. Here, there are n residuals, and k jitter - parameters. - """ - cdef unsigned int cc, ii, rows = len(Nvec), cols = len(Jvec) - cdef double tr=0.0, ji, nisum, Nnisum1, Nnisum2, NniNnisum, beta - cdef np.ndarray[np.double_t,ndim=1] ni = np.empty(rows, 'd') - cdef np.ndarray[np.double_t,ndim=1] Nni1 = np.empty(rows, 'd') - cdef np.ndarray[np.double_t,ndim=1] Nni2 = np.empty(rows, 'd') - - ni = 1.0 / Nvec - Nni1 = dNvec1 / Nvec**2 - Nni2 = dNvec2 / Nvec**2 - - for cc in range(rows): - tr += dNvec1[cc] * dNvec2[cc] * ni[cc]**2 - - for cc in range(cols): - if Jvec[cc] > 0.0: - ji = 1.0 / Jvec[cc] - - nisum = 0.0 - Nnisum1 = 0.0 - Nnisum2 = 0.0 - NniNnisum = 0.0 - for ii in range(Uinds[cc,0],Uinds[cc,1]): - nisum += ni[ii] - Nnisum1 += Nni1[ii] - Nnisum2 += Nni2[ii] - NniNnisum += Nni1[ii]*Nni2[ii]*Nvec[ii] - - beta = 1.0 / (nisum + ji) - - tr += Nnisum1 * Nnisum2 * beta**2 - tr -= 2 * NniNnisum * beta - - return tr - -@cython.boundscheck(False) -@cython.wraparound(False) -def cython_logdet_dN_dJ( \ - np.ndarray[np.double_t,ndim=1] Nvec, \ - np.ndarray[np.double_t,ndim=1] Jvec, \ - np.ndarray[np.double_t,ndim=1] dNvec, \ - np.ndarray[np.double_t,ndim=1] dJvec, \ - np.ndarray[np.int_t,ndim=2] Uinds): - """ - Sherman-Morrison block-inversion for Jitter (Cythonized) - - Calculates Trace(N^{-1} dN/dNp N^{-1} dN/dJp), where: - - N^{-1} is the ecorr-include N inverse - - dN/dNp is the diagonal derivate of N wrt Np - - dN/dJp = U dJ/dJp U^{T}, with dJ/dJp the diagnal derivative of J wrt - Jp - - @param Nvec: The white noise amplitude, array (n) - @param Jvec: The jitter amplitude, array (k) - @param dNvec: The white noise derivative, array (n) - @param dJvec: The white noise ecor derivative, array (k) - @param Uinds: The start/finish indices for the jitter blocks (k x 2) - - For this version, the residuals need to be sorted properly so that all the - blocks are continuous in memory. Here, there are n residuals, and k jitter - parameters. - """ - cdef unsigned int cc, ii, rows = len(Nvec), cols = len(Jvec) - cdef double tr=0.0, ji, nisum, Nnisum, beta - cdef np.ndarray[np.double_t,ndim=1] ni = np.empty(rows, 'd') - cdef np.ndarray[np.double_t,ndim=1] Nni = np.empty(rows, 'd') - - ni = 1.0 / Nvec - Nni = dNvec / Nvec**2 - - for cc in range(cols): - if Jvec[cc] > 0.0: - ji = 1.0 / Jvec[cc] - - nisum = 0.0 - Nnisum = 0.0 - for ii in range(Uinds[cc,0],Uinds[cc,1]): - nisum += ni[ii] - Nnisum += Nni[ii] - - beta = 1.0 / (nisum + ji) - - tr += Nnisum * dJvec[cc] - tr -= 2 * nisum * dJvec[cc] * Nnisum * beta - tr += Nnisum * nisum**2 * dJvec[cc] *beta**2 - - return tr - -@cython.boundscheck(False) -@cython.wraparound(False) -def cython_logdet_dJ_dJ( \ - np.ndarray[np.double_t,ndim=1] Nvec, \ - np.ndarray[np.double_t,ndim=1] Jvec, \ - np.ndarray[np.double_t,ndim=1] dJvec1, \ - np.ndarray[np.double_t,ndim=1] dJvec2, \ - np.ndarray[np.int_t,ndim=2] Uinds): - """ - Sherman-Morrison block-inversion for Jitter (Cythonized) - - Calculates Trace(N^{-1} dN/dJp1 N^{-1} dN/dJp2), where: - - N^{-1} is the ecorr-include N inverse - - dN/dJpx = U dJ/dJpx U^{T}, with dJ/dJpx the diagnal derivative of J wrt - Jpx - - @param Nvec: The white noise amplitude, array (n) - @param Jvec: The jitter amplitude, array (k) - @param dJvec1: The white noise derivative, array (k) - @param dJvec2: The white noise derivative, array (k) - @param Uinds: The start/finish indices for the jitter blocks (k x 2) - - For this version, the residuals need to be sorted properly so that all the - blocks are continuous in memory. Here, there are n residuals, and k jitter - parameters. - """ - cdef unsigned int cc, ii, rows = len(Nvec), cols = len(Jvec) - cdef double tr=0.0, ji, nisum, beta - cdef np.ndarray[np.double_t,ndim=1] ni = np.empty(rows, 'd') - - ni = 1.0 / Nvec - - for cc in range(cols): - if Jvec[cc] > 0.0: - ji = 1.0 / Jvec[cc] - - nisum = 0.0 - for ii in range(Uinds[cc,0],Uinds[cc,1]): - nisum += ni[ii] - - beta = 1.0 / (nisum + ji) - - tr += dJvec1[cc] * dJvec2[cc] * nisum**2 - tr -= 2 * dJvec1[cc] * dJvec2[cc] * beta * nisum**3 - tr += dJvec1[cc] * dJvec2[cc] * beta**2 * nisum**4 - - return tr - - -@cython.boundscheck(False) -@cython.wraparound(False) -cpdef double c_blas_block_shermor_2D_asymm( - np.ndarray[np.double_t,ndim=2] Z1, - np.ndarray[np.double_t,ndim=2] Z2, - np.ndarray[np.double_t,ndim=1] Nvec, - np.ndarray[np.double_t,ndim=1] Jvec, - Uinds, #Need to copy, because Cython can't always use numpy integer arrays. Bah - np.ndarray[np.double_t,ndim=2] ZNZ, - ): - cdef double d_Jldet - - cdef int n_Z1_rows = Z1.shape[0] - cdef int n_Z1_cols = Z1.shape[1] - cdef int n_Z2_cols = Z2.shape[1] - cdef int n_J_rows = len(Jvec) - cdef int n_Z1_row_major, n_Z2_row_major - - n_Z1_row_major = 0 if Z1.flags['F_CONTIGUOUS'] else 1 - n_Z2_row_major = 0 if Z2.flags['F_CONTIGUOUS'] else 1 - - # Hack, because somehow we can't pass integer arrays? - # This is a tiny bit of overhead right here - cdef int [:] Uinds_new - - # This makes it into a C-ordered (Row-Major) array we can pass - Uinds_new = np.ascontiguousarray(Uinds.flatten(), dtype=np.dtype("i")) - - assert Z1.shape[0] == Z2.shape[0] - assert Z1.shape[0] == len(Nvec) - assert Uinds.shape[0] == len(Jvec) - - blas_block_shermor_2D_asym( - n_Z1_rows, - n_Z1_cols, - n_Z1_row_major, - &Z1[0,0], - n_Z2_cols, - n_Z2_row_major, - &Z2[0,0], - &Nvec[0], - n_J_rows, - &Jvec[0], - &Uinds_new[0], - &ZNZ[0,0], - &d_Jldet) - - return d_Jldet - -def cython_blas_block_shermor_2D_asymm(Z1, Z2, Nvec, Jvec, Uinds): - """Wrapper for the C/Cython code""" - - ZNZ = np.zeros((Z1.shape[1], Z2.shape[1]), order='F') - Jldet = c_blas_block_shermor_2D_asymm(Z1, Z2, Nvec, Jvec, np.array(Uinds, order="C"), ZNZ) - - return Jldet, ZNZ - diff --git a/enterprise/fastshermanmorrison/fastshermanmorrison.c b/enterprise/fastshermanmorrison/fastshermanmorrison.c deleted file mode 100644 index f6f13c9d..00000000 --- a/enterprise/fastshermanmorrison/fastshermanmorrison.c +++ /dev/null @@ -1,233 +0,0 @@ -/* cython_fastshermor.c - * - * Rutger van Haasteren, April 19 2023, Hannover - * - */ - -#include -#include -#include - - -extern void dgemm_(char *transa, char *transb, int *m, int *n, int *k, - double *alpha, double *a, int *lda, double *b, - int *ldb, double *beta, double *c, int *ldc); -extern void dgemv_(char *trans, int *m, int *n, double *alpha, double *a, - int *lda, double *x, int *incx, double *beta, - double *y, int *incy); -extern void dger_(int *m, int *n, double *alpha, double *x, int *incx, - double *y, int *incy, double *a, int *lda); - -static void blas_block_shermor_2D_asym( - int n_Z1_rows, - int n_Z1_cols, - int n_Z1_row_major, - double *pd_Z1, - int n_Z2_cols, - int n_Z2_row_major, - double *pd_Z2, - double *pd_Nvec, - int n_J_rows, - double *pd_Jvec, - int *pn_Uinds, - double *pd_ZNZ, - double *pd_Jldet - ) { - /* C implementation of python_block_shermor_2D_asym, because the python - * overhead is large - * - * parameters - * ---------- - * - * :param n_Z1_rows: Number of rows of Z1 - * :param n_Z1_cols: Number of columns of Z1 - * :param n_Z1_row_major: 1 if Z1 is Row-Major, 0 if Column-Major - * :param pd_Z1: The Z1 matrix - * :param n_Z2_cols: Number of columns of Z2 - * :param n_Z2_row_major: 1 if Z2 is Row-Major, 0 if Column-Major - * :param pd_Z2: The Z2 matrix - * :param pd_Nvec: The Nvec vector - * :param n_J_rows: The number of Jvec elements - * :param pd_Jvec: The Jvec vector - * :param pn_Uinds: The matrix of quantization indices (Row-Major) - * :param pd_ZNZ: The return value of ZNZ (Column-Major) - * :param pd_Jldet: The return value of log(det(J)) - */ - - double d_galpha=1.0, d_gbeta=0.0, d_nisum=0.0, d_beta; - double *pd_Z1ni, *pd_ZNZ_add, *pd_ni, *pd_zn1, *pd_zn2; - int cc, i, j, m, n, k, lda, ldb, ldc, n_jblock, n_jblock_i, n_index; - char *transa, *transb; - - pd_Z1ni = malloc(n_Z1_rows*n_Z1_cols * sizeof(double)); - pd_ZNZ_add = calloc(n_Z1_rows*n_Z1_cols, sizeof(double)); - pd_ni = malloc(n_Z1_rows * sizeof(double)); - pd_zn1 = calloc(n_Z1_cols, sizeof(double)); - pd_zn2 = calloc(n_Z2_cols, sizeof(double)); - - /* openmp this? */ - for(i=0; i 0.0) { - - /* Note: pn_Uinds is row-major */ - d_nisum = 0.0; - n_jblock_i = pn_Uinds[2*cc]; - n_jblock = pn_Uinds[2*cc+1] - pn_Uinds[2*cc]; - for(i=pn_Uinds[2*cc]; i=2.7.0 coverage-conditional-plugin>=0.4.0 jupyter>=1.0.0 build==0.3.1.post1 -numpy>=1.16.3 -cython>=0.29.34 diff --git a/setup.py b/setup.py index fb2efe85..45970a9b 100644 --- a/setup.py +++ b/setup.py @@ -1,10 +1,7 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- -import numpy from setuptools import setup -from setuptools import Extension -from Cython.Build import cythonize with open("README.md", encoding="utf-8") as readme_file: readme = readme_file.read() @@ -17,14 +14,6 @@ "scikit-sparse>=0.4.5", "pint-pulsar>=0.8.3", "libstempo>=2.4.4", - "cython>=0.29.34", -] - -ext_modules=[ - Extension('enterprise.fastshermanmorrison.cython_fastshermanmorrison', - ['enterprise/fastshermanmorrison/cython_fastshermanmorrison.pyx'], - include_dirs = [numpy.get_include(), 'fastshermanmorrison/'], - extra_compile_args=["-O2", "-fno-wrapv"]) # 50% more efficient! ] test_requirements = [] @@ -62,5 +51,4 @@ ], test_suite="tests", tests_require=test_requirements, - ext_modules = cythonize(ext_modules) ) From fb7f533bd71d528c4d8c41c9e997243691259f9c Mon Sep 17 00:00:00 2001 From: Rutger van Haasteren Date: Mon, 24 Apr 2023 17:49:37 +0200 Subject: [PATCH 06/20] Removed fastshermanmorrison from setup.py --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index 45970a9b..2e5993a3 100644 --- a/setup.py +++ b/setup.py @@ -26,7 +26,7 @@ author="Justin A. Ellis", author_email="justin.ellis18@gmail.com", url="https://github.com/nanograv/enterprise", - packages=["enterprise", "enterprise.signals", "enterprise.fastshermanmorrison"], + packages=["enterprise", "enterprise.signals"], package_dir={"enterprise": "enterprise"}, include_package_data=True, package_data={"enterprise": ["datafiles/*", "datafiles/ephemeris/*", "datafiles/ng9/*", "datafiles/mdc_open1/*"]}, From 50cc052626824bd0d607006c0a68318409ed8de9 Mon Sep 17 00:00:00 2001 From: Rutger van Haasteren Date: Mon, 24 Apr 2023 18:05:01 +0200 Subject: [PATCH 07/20] Conform to black linter --- enterprise/signals/white_signals.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/enterprise/signals/white_signals.py b/enterprise/signals/white_signals.py index ac319820..cd740648 100644 --- a/enterprise/signals/white_signals.py +++ b/enterprise/signals/white_signals.py @@ -20,6 +20,7 @@ # logging.basicConfig(format="%(levelname)s: %(name)s: %(message)s", level=logging.INFO) logger = logging.getLogger(__name__) + def WhiteNoise(varianceFunction, selection=Selection(selections.no_selection), name=""): """Class factory for generic white noise signals.""" @@ -185,7 +186,7 @@ def EcorrKernelNoise( msg = "EcorrKernelNoise does not support method: {}".format(method) raise TypeError(msg) - if method == 'fast-sherman-morrison' and fastshermanmorrison is None: + if method == "fast-sherman-morrison" and fastshermanmorrison is None: msg = "Package `fastshermanmorrison` not installed. Fallback to sherman-morrison" logger.warning(msg) From 39b6b0e59393741b52bdac4b13be28ecda5d1381 Mon Sep 17 00:00:00 2001 From: Rutger van Haasteren Date: Wed, 26 Apr 2023 13:55:46 +0200 Subject: [PATCH 08/20] Added test for fast-sherman-morrison --- tests/test_white_signals.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/tests/test_white_signals.py b/tests/test_white_signals.py index 0736bf71..7f67d272 100644 --- a/tests/test_white_signals.py +++ b/tests/test_white_signals.py @@ -470,6 +470,10 @@ def test_ecorr_sherman_morrison(self): """Test of sherman-morrison ecorr signal and solve methods.""" self._ecorr_test(method="sherman-morrison") + def test_ecorr_fast_sherman_morrison(self): + """Test of fast-sherman-morrison ecorr signal and solve methods.""" + self._ecorr_test(method="fast-sherman-morrison") + def test_ecorr_block(self): """Test of block matrix ecorr signal and solve methods.""" self._ecorr_test(method="block") From 0e3435c370fcef368201a90818da986774b27b26 Mon Sep 17 00:00:00 2001 From: Rutger van Haasteren Date: Wed, 26 Apr 2023 13:59:21 +0200 Subject: [PATCH 09/20] Added requirements_tests.txt, so that the extensions can be loaded for testing. --- Makefile | 1 + requirements_tests.txt | 1 + 2 files changed, 2 insertions(+) create mode 100644 requirements_tests.txt diff --git a/Makefile b/Makefile index 70a5a1ed..43a249c3 100644 --- a/Makefile +++ b/Makefile @@ -31,6 +31,7 @@ init: @./.enterprise/bin/python3 -m pip install -U pip setuptools wheel @./.enterprise/bin/python3 -m pip install -r requirements.txt -U @./.enterprise/bin/python3 -m pip install -r requirements_dev.txt -U + @./.enterprise/bin/python3 -m pip install -r requirements_tests.txt -U @./.enterprise/bin/python3 -m pre_commit install --install-hooks --overwrite @./.enterprise/bin/python3 -m pip install -e . @echo "run source .enterprise/bin/activate to activate environment" diff --git a/requirements_tests.txt b/requirements_tests.txt new file mode 100644 index 00000000..7cfe1708 --- /dev/null +++ b/requirements_tests.txt @@ -0,0 +1 @@ +fastshermanmorrison-pulsar>=0.0.3 From 60a117ee2b3386e71524d39f5c43f7e5c9ee6314 Mon Sep 17 00:00:00 2001 From: Rutger van Haasteren Date: Wed, 26 Apr 2023 15:08:44 +0200 Subject: [PATCH 10/20] Added decorator to show fastshermanmorrison-pulsar warning only once --- enterprise/signals/utils.py | 12 ++++++++++++ enterprise/signals/white_signals.py | 6 ++++-- 2 files changed, 16 insertions(+), 2 deletions(-) diff --git a/enterprise/signals/utils.py b/enterprise/signals/utils.py index df949af0..9efe6d4c 100644 --- a/enterprise/signals/utils.py +++ b/enterprise/signals/utils.py @@ -31,6 +31,18 @@ logger = logging.getLogger(__name__) +def static_vars(**kwargs): + """This decorator allows a class factory to have a static variable, for + things such as checking whether a python package exists. If the package + does not exist, a warning is shown only once + """ + def decorate(func): + for k in kwargs: + setattr(func, k, kwargs[k]) + return func + return decorate + + class ConditionalGP: def __init__(self, pta, phiinv_method="cliques"): """This class allows the computation of conditional means and diff --git a/enterprise/signals/white_signals.py b/enterprise/signals/white_signals.py index cd740648..fe8c1317 100644 --- a/enterprise/signals/white_signals.py +++ b/enterprise/signals/white_signals.py @@ -119,7 +119,7 @@ def EquadNoise(*args, **kwargs): " or TNEquadNoise to obtain legacy enterprise definition for EQUAD only [tnequad^2]." ) - +@utils.static_vars(shown_fastshermanmorrison_warning=False) def EcorrKernelNoise( log10_ecorr=parameter.Uniform(-10, -5), selection=Selection(selections.no_selection), @@ -186,9 +186,11 @@ def EcorrKernelNoise( msg = "EcorrKernelNoise does not support method: {}".format(method) raise TypeError(msg) - if method == "fast-sherman-morrison" and fastshermanmorrison is None: + if method == "fast-sherman-morrison" and fastshermanmorrison is None \ + and not EcorrKernelNoise.shown_fastshermanmorrison_warning: msg = "Package `fastshermanmorrison` not installed. Fallback to sherman-morrison" logger.warning(msg) + EcorrKernelNoise.shown_fastshermanmorrison_warning = True class EcorrKernelNoise(signal_base.Signal): signal_type = "white noise" From 851b13ef1096693fe2b2ce25c56d0e200ddba2f7 Mon Sep 17 00:00:00 2001 From: Rutger van Haasteren Date: Wed, 26 Apr 2023 15:18:08 +0200 Subject: [PATCH 11/20] Proper linting corrections --- enterprise/signals/utils.py | 2 ++ enterprise/signals/white_signals.py | 8 ++++++-- 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/enterprise/signals/utils.py b/enterprise/signals/utils.py index 9efe6d4c..ee57f2f6 100644 --- a/enterprise/signals/utils.py +++ b/enterprise/signals/utils.py @@ -36,10 +36,12 @@ def static_vars(**kwargs): things such as checking whether a python package exists. If the package does not exist, a warning is shown only once """ + def decorate(func): for k in kwargs: setattr(func, k, kwargs[k]) return func + return decorate diff --git a/enterprise/signals/white_signals.py b/enterprise/signals/white_signals.py index fe8c1317..deb10c25 100644 --- a/enterprise/signals/white_signals.py +++ b/enterprise/signals/white_signals.py @@ -119,6 +119,7 @@ def EquadNoise(*args, **kwargs): " or TNEquadNoise to obtain legacy enterprise definition for EQUAD only [tnequad^2]." ) + @utils.static_vars(shown_fastshermanmorrison_warning=False) def EcorrKernelNoise( log10_ecorr=parameter.Uniform(-10, -5), @@ -186,8 +187,11 @@ def EcorrKernelNoise( msg = "EcorrKernelNoise does not support method: {}".format(method) raise TypeError(msg) - if method == "fast-sherman-morrison" and fastshermanmorrison is None \ - and not EcorrKernelNoise.shown_fastshermanmorrison_warning: + if ( + method == "fast-sherman-morrison" + and fastshermanmorrison is None + and not EcorrKernelNoise.shown_fastshermanmorrison_warning + ): msg = "Package `fastshermanmorrison` not installed. Fallback to sherman-morrison" logger.warning(msg) EcorrKernelNoise.shown_fastshermanmorrison_warning = True From bf1705e57042fae00fa62162f1b5fcfca22b2843 Mon Sep 17 00:00:00 2001 From: Rutger van Haasteren Date: Wed, 26 Apr 2023 15:48:48 +0200 Subject: [PATCH 12/20] Changed the EcorrKernelWarning machanism --- enterprise/signals/utils.py | 14 -------------- enterprise/signals/white_signals.py | 6 +++--- 2 files changed, 3 insertions(+), 17 deletions(-) diff --git a/enterprise/signals/utils.py b/enterprise/signals/utils.py index ee57f2f6..df949af0 100644 --- a/enterprise/signals/utils.py +++ b/enterprise/signals/utils.py @@ -31,20 +31,6 @@ logger = logging.getLogger(__name__) -def static_vars(**kwargs): - """This decorator allows a class factory to have a static variable, for - things such as checking whether a python package exists. If the package - does not exist, a warning is shown only once - """ - - def decorate(func): - for k in kwargs: - setattr(func, k, kwargs[k]) - return func - - return decorate - - class ConditionalGP: def __init__(self, pta, phiinv_method="cliques"): """This class allows the computation of conditional means and diff --git a/enterprise/signals/white_signals.py b/enterprise/signals/white_signals.py index deb10c25..5b9d8dd7 100644 --- a/enterprise/signals/white_signals.py +++ b/enterprise/signals/white_signals.py @@ -16,6 +16,7 @@ import fastshermanmorrison.fastshermanmorrison as fastshermanmorrison except ImportError: fastshermanmorrison = None + shown_fastshermanmorrison_warning = False # logging.basicConfig(format="%(levelname)s: %(name)s: %(message)s", level=logging.INFO) logger = logging.getLogger(__name__) @@ -120,7 +121,6 @@ def EquadNoise(*args, **kwargs): ) -@utils.static_vars(shown_fastshermanmorrison_warning=False) def EcorrKernelNoise( log10_ecorr=parameter.Uniform(-10, -5), selection=Selection(selections.no_selection), @@ -190,11 +190,11 @@ def EcorrKernelNoise( if ( method == "fast-sherman-morrison" and fastshermanmorrison is None - and not EcorrKernelNoise.shown_fastshermanmorrison_warning + and not shown_fastshermanmorrison_warning ): msg = "Package `fastshermanmorrison` not installed. Fallback to sherman-morrison" logger.warning(msg) - EcorrKernelNoise.shown_fastshermanmorrison_warning = True + shown_fastshermanmorrison_warning = True class EcorrKernelNoise(signal_base.Signal): signal_type = "white noise" From f8d71f5dd9c676ae7bc9943996bcf39de39d5423 Mon Sep 17 00:00:00 2001 From: Rutger van Haasteren Date: Wed, 26 Apr 2023 16:01:54 +0200 Subject: [PATCH 13/20] Blacking --- enterprise/signals/white_signals.py | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/enterprise/signals/white_signals.py b/enterprise/signals/white_signals.py index 5b9d8dd7..fd396840 100644 --- a/enterprise/signals/white_signals.py +++ b/enterprise/signals/white_signals.py @@ -187,11 +187,7 @@ def EcorrKernelNoise( msg = "EcorrKernelNoise does not support method: {}".format(method) raise TypeError(msg) - if ( - method == "fast-sherman-morrison" - and fastshermanmorrison is None - and not shown_fastshermanmorrison_warning - ): + if method == "fast-sherman-morrison" and fastshermanmorrison is None and not shown_fastshermanmorrison_warning: msg = "Package `fastshermanmorrison` not installed. Fallback to sherman-morrison" logger.warning(msg) shown_fastshermanmorrison_warning = True From 8c4bfa8a148725bd660a25ec444eb29ca71e56f4 Mon Sep 17 00:00:00 2001 From: Rutger van Haasteren Date: Wed, 26 Apr 2023 16:15:49 +0200 Subject: [PATCH 14/20] Lint update --- enterprise/signals/white_signals.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/enterprise/signals/white_signals.py b/enterprise/signals/white_signals.py index fd396840..ab28881b 100644 --- a/enterprise/signals/white_signals.py +++ b/enterprise/signals/white_signals.py @@ -16,7 +16,7 @@ import fastshermanmorrison.fastshermanmorrison as fastshermanmorrison except ImportError: fastshermanmorrison = None - shown_fastshermanmorrison_warning = False +shown_fastshermanmorrison_warning = False # logging.basicConfig(format="%(levelname)s: %(name)s: %(message)s", level=logging.INFO) logger = logging.getLogger(__name__) From fb01136445da5edd0c6e7cd1da51d0a98ab4c032 Mon Sep 17 00:00:00 2001 From: Rutger van Haasteren Date: Wed, 26 Apr 2023 18:35:32 +0200 Subject: [PATCH 15/20] Removed the have_shown stuff for now --- enterprise/signals/white_signals.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/enterprise/signals/white_signals.py b/enterprise/signals/white_signals.py index ab28881b..cd740648 100644 --- a/enterprise/signals/white_signals.py +++ b/enterprise/signals/white_signals.py @@ -16,7 +16,6 @@ import fastshermanmorrison.fastshermanmorrison as fastshermanmorrison except ImportError: fastshermanmorrison = None -shown_fastshermanmorrison_warning = False # logging.basicConfig(format="%(levelname)s: %(name)s: %(message)s", level=logging.INFO) logger = logging.getLogger(__name__) @@ -187,10 +186,9 @@ def EcorrKernelNoise( msg = "EcorrKernelNoise does not support method: {}".format(method) raise TypeError(msg) - if method == "fast-sherman-morrison" and fastshermanmorrison is None and not shown_fastshermanmorrison_warning: + if method == "fast-sherman-morrison" and fastshermanmorrison is None: msg = "Package `fastshermanmorrison` not installed. Fallback to sherman-morrison" logger.warning(msg) - shown_fastshermanmorrison_warning = True class EcorrKernelNoise(signal_base.Signal): signal_type = "white noise" From 3ff9921c8046b1a102dfceb72ae507d5571e3696 Mon Sep 17 00:00:00 2001 From: Rutger van Haasteren Date: Thu, 4 May 2023 17:21:25 +0200 Subject: [PATCH 16/20] Placed fastshermanmorrison into requirements_dev.txt --- Makefile | 1 - requirements_dev.txt | 1 + requirements_tests.txt | 1 - 3 files changed, 1 insertion(+), 2 deletions(-) delete mode 100644 requirements_tests.txt diff --git a/Makefile b/Makefile index 43a249c3..70a5a1ed 100644 --- a/Makefile +++ b/Makefile @@ -31,7 +31,6 @@ init: @./.enterprise/bin/python3 -m pip install -U pip setuptools wheel @./.enterprise/bin/python3 -m pip install -r requirements.txt -U @./.enterprise/bin/python3 -m pip install -r requirements_dev.txt -U - @./.enterprise/bin/python3 -m pip install -r requirements_tests.txt -U @./.enterprise/bin/python3 -m pre_commit install --install-hooks --overwrite @./.enterprise/bin/python3 -m pip install -e . @echo "run source .enterprise/bin/activate to activate environment" diff --git a/requirements_dev.txt b/requirements_dev.txt index 6d103434..b784df3d 100644 --- a/requirements_dev.txt +++ b/requirements_dev.txt @@ -15,3 +15,4 @@ pytest-cov>=2.7.0 coverage-conditional-plugin>=0.4.0 jupyter>=1.0.0 build==0.3.1.post1 +fastshermanmorrison-pulsar>=0.0.3 diff --git a/requirements_tests.txt b/requirements_tests.txt deleted file mode 100644 index 7cfe1708..00000000 --- a/requirements_tests.txt +++ /dev/null @@ -1 +0,0 @@ -fastshermanmorrison-pulsar>=0.0.3 From 37b4c9c2a8bf9020d8503d7ca74f7eee432a2381 Mon Sep 17 00:00:00 2001 From: Rutger van Haasteren Date: Mon, 25 Sep 2023 13:17:30 +0200 Subject: [PATCH 17/20] Provided mechanism for showing warning message about fastshermanmorrison, and updated fastshermanmorrison-pulsar in requirements_dev.txt --- enterprise/signals/white_signals.py | 8 +++++++- requirements_dev.txt | 2 +- 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/enterprise/signals/white_signals.py b/enterprise/signals/white_signals.py index cd740648..72d079d4 100644 --- a/enterprise/signals/white_signals.py +++ b/enterprise/signals/white_signals.py @@ -14,8 +14,11 @@ try: import fastshermanmorrison.fastshermanmorrison as fastshermanmorrison + + fsm_warning_issued = False except ImportError: fastshermanmorrison = None + fsm_warning_issued = False # logging.basicConfig(format="%(levelname)s: %(name)s: %(message)s", level=logging.INFO) logger = logging.getLogger(__name__) @@ -182,13 +185,16 @@ def EcorrKernelNoise( """ + global fsm_warning_issued + if method not in ["fast-sherman-morrison", "sherman-morrison", "block", "sparse"]: msg = "EcorrKernelNoise does not support method: {}".format(method) raise TypeError(msg) - if method == "fast-sherman-morrison" and fastshermanmorrison is None: + if method == "fast-sherman-morrison" and fastshermanmorrison is None and not fsm_warning_issued: msg = "Package `fastshermanmorrison` not installed. Fallback to sherman-morrison" logger.warning(msg) + fsm_warning_issued = True class EcorrKernelNoise(signal_base.Signal): signal_type = "white noise" diff --git a/requirements_dev.txt b/requirements_dev.txt index b784df3d..d173ce59 100644 --- a/requirements_dev.txt +++ b/requirements_dev.txt @@ -15,4 +15,4 @@ pytest-cov>=2.7.0 coverage-conditional-plugin>=0.4.0 jupyter>=1.0.0 build==0.3.1.post1 -fastshermanmorrison-pulsar>=0.0.3 +fastshermanmorrison-pulsar>=0.1.0 From 1df974af108703e76c06c40073896d1947a07c6b Mon Sep 17 00:00:00 2001 From: Rutger van Haasteren Date: Tue, 26 Sep 2023 18:34:10 +0200 Subject: [PATCH 18/20] Added pragmas for no coverage of exception results --- enterprise/signals/white_signals.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/enterprise/signals/white_signals.py b/enterprise/signals/white_signals.py index 72d079d4..f74f3d10 100644 --- a/enterprise/signals/white_signals.py +++ b/enterprise/signals/white_signals.py @@ -16,7 +16,7 @@ import fastshermanmorrison.fastshermanmorrison as fastshermanmorrison fsm_warning_issued = False -except ImportError: +except ImportError: #pragma: no cover fastshermanmorrison = None fsm_warning_issued = False @@ -191,7 +191,7 @@ def EcorrKernelNoise( msg = "EcorrKernelNoise does not support method: {}".format(method) raise TypeError(msg) - if method == "fast-sherman-morrison" and fastshermanmorrison is None and not fsm_warning_issued: + if method == "fast-sherman-morrison" and fastshermanmorrison is None and not fsm_warning_issued: # pragma: no cover msg = "Package `fastshermanmorrison` not installed. Fallback to sherman-morrison" logger.warning(msg) fsm_warning_issued = True @@ -239,7 +239,7 @@ def get_ndiag(self, params): elif method == "fast-sherman-morrison": if fastshermanmorrison: return self._get_ndiag_fast_sherman_morrison(params) - else: + else: # pragma: no cover return self._get_ndiag_sherman_morrison(params) elif method == "sparse": return self._get_ndiag_sparse(params) From 3698716cef980ffce1124149632e46b5f500a23d Mon Sep 17 00:00:00 2001 From: Rutger van Haasteren Date: Tue, 26 Sep 2023 20:45:13 +0200 Subject: [PATCH 19/20] Linting: forgot a space in a pragma --- enterprise/signals/white_signals.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/enterprise/signals/white_signals.py b/enterprise/signals/white_signals.py index f74f3d10..a1df6f86 100644 --- a/enterprise/signals/white_signals.py +++ b/enterprise/signals/white_signals.py @@ -16,7 +16,7 @@ import fastshermanmorrison.fastshermanmorrison as fastshermanmorrison fsm_warning_issued = False -except ImportError: #pragma: no cover +except ImportError: # pragma: no cover fastshermanmorrison = None fsm_warning_issued = False From 7f6795acac58a1e5340bcb219fbb4cc2a38fb25a Mon Sep 17 00:00:00 2001 From: Rutger van Haasteren Date: Wed, 8 Nov 2023 20:28:38 +0100 Subject: [PATCH 20/20] Updated version of fastshermanmorrison-pulsar to 0.4.0, which has binary builds --- requirements_dev.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/requirements_dev.txt b/requirements_dev.txt index d173ce59..e3ba1966 100644 --- a/requirements_dev.txt +++ b/requirements_dev.txt @@ -15,4 +15,4 @@ pytest-cov>=2.7.0 coverage-conditional-plugin>=0.4.0 jupyter>=1.0.0 build==0.3.1.post1 -fastshermanmorrison-pulsar>=0.1.0 +fastshermanmorrison-pulsar>=0.4.0