diff --git a/src/toast/templates/offset/offset.py b/src/toast/templates/offset/offset.py index ed8430b81..13a5f24e8 100644 --- a/src/toast/templates/offset/offset.py +++ b/src/toast/templates/offset/offset.py @@ -487,27 +487,48 @@ def _initialize(self, new_data): # not invert "M" and solve M x = b using the Cholesky # decomposition of M (*not* M^{-1}). icenter = noisefilter.size // 2 - wband = min(self.precond_width, icenter) - precond_width = max( - wband, min(self.precond_width, n_amp_view) - ) - preconditioner = np.zeros( - [precond_width, n_amp_view], dtype=np.float64 - ) - if detnoise != 0: - preconditioner[0, :] = 1.0 / offsetvar_slice - preconditioner[:wband, :] += np.repeat( - noisefilter[icenter : icenter + wband, np.newaxis], - n_amp_view, - 1, - ) - lower = True - preconditioner = scipy.linalg.cholesky_banded( - preconditioner, - overwrite_ab=True, - lower=lower, - check_finite=True, - ) + try_width = self.precond_width + + good_cholesky = False + while not good_cholesky: + wband = min(try_width, icenter) + precond_width = max( + wband, min(try_width, n_amp_view) + ) + preconditioner = np.zeros( + [precond_width, n_amp_view], dtype=np.float64 + ) + if detnoise != 0: + preconditioner[0, :] = 1.0 / offsetvar_slice + preconditioner[:wband, :] += np.repeat( + noisefilter[icenter : icenter + wband, np.newaxis], + n_amp_view, + 1, + ) + lower = True + try: + preconditioner = scipy.linalg.cholesky_banded( + preconditioner, + overwrite_ab=True, + lower=lower, + check_finite=True, + ) + good_cholesky = True + except scipy.linalg.LinAlgError: + if try_width < icenter and try_width < n_amp_view: + # Still room to increase the band width + try_width *= 2 + msg = f"obs {ob.name}, det {det}, view {ivw} " + msg += " scipy cholesky_banded fail, increasing" + msg += f" width to {try_width}" + log.debug(msg) + else: + # Width is now the size of the data interval + msg = f"obs {ob.name}, det {det}, view {ivw} " + msg += " scipy cholesky_banded fail with max " + msg += "width of {try_width}" + log.error(msg) + raise RuntimeError(msg) if self.debug_plots is not None: axprec.plot( np.arange(len(preconditioner)), diff --git a/src/toast/tests/ops_mapmaker.py b/src/toast/tests/ops_mapmaker.py index 3468aeb76..a4b4b17de 100644 --- a/src/toast/tests/ops_mapmaker.py +++ b/src/toast/tests/ops_mapmaker.py @@ -595,7 +595,7 @@ def test_compare_madam_diagpre(self): step_time=step_seconds * u.second, use_noise_prior=True, precond_width=1, - # debug_plots=testdir, + debug_plots=testdir, ) tmatrix = ops.TemplateMatrix(templates=[tmpl]) @@ -613,6 +613,7 @@ def test_compare_madam_diagpre(self): solve_rcond_threshold=1.0e-4, map_rcond_threshold=1.0e-4, iter_max=50, + output_dir=testdir, write_hits=False, write_map=False, write_cov=False, @@ -867,8 +868,8 @@ def test_compare_madam_bandpre(self): noise_model=default_model.noise_model, step_time=step_seconds * u.second, use_noise_prior=True, - precond_width=10, - # debug_plots=testdir, + precond_width=2, + debug_plots=testdir, ) tmatrix = ops.TemplateMatrix(templates=[tmpl]) @@ -886,6 +887,7 @@ def test_compare_madam_bandpre(self): solve_rcond_threshold=1.0e-4, map_rcond_threshold=1.0e-4, iter_max=50, + output_dir=testdir, write_hits=False, write_map=False, write_cov=False,