diff --git a/scripts/covmatrix_plot.py b/scripts/covmatrix_plot.py index 79d7ff5658a0792b862b7fe27e01d3781e853fa8..50d6eb62faeacaa6d0cf287669af672a88985894 100755 --- a/scripts/covmatrix_plot.py +++ b/scripts/covmatrix_plot.py @@ -8,6 +8,7 @@ from typing import TYPE_CHECKING from h5py import File, Group from matplotlib import pyplot as plt from matplotlib.backends.backend_pdf import PdfPages +from matplotlib.colors import SymLogNorm from numpy import arange, diagonal, isnan, quantile from numpy.typing import NDArray @@ -26,7 +27,6 @@ def main(opts: Namespace) -> None: model = group["model"][:] edges = group["edges"][:] - blocksize = edges.size - 1 sfile, sgroup, smodel = None, None, None if opts.subtract: @@ -48,6 +48,11 @@ def main(opts: Namespace) -> None: except KeyError: elements = None + blocksize = edges.size - 1 + block_selected = 1 + block_selection = slice(blocksize * block_selected, blocksize * (block_selected + 1)) + block_name = elements[block_selected] + if opts.output is not None: ofile = opts.output and opts.output[0] or opts.input.replace(".hdf5", ".pdf") assert ofile != opts.input and ofile.endswith(".pdf") @@ -156,8 +161,8 @@ def main(opts: Namespace) -> None: ylabel="bin", title=f"Covariance matrix {name}{title_suffix}", ) - pcolor_with_blocks(matrix_cov, blocks=elements, cmap=cmap) - if pdf: + hm = pcolor_with_blocks(matrix_cov, blocks=elements, cmap=cmap) + if pdf and hm is not None: pdf.savefig() if vmax is not None: @@ -166,12 +171,25 @@ def main(opts: Namespace) -> None: 111, xlabel="", ylabel="bin", - title=f"Covariance matrix (truncated) {name}{title_suffix}", + title=f"Covariance matrix {name} (truncated){title_suffix}", ) - pcolor_with_blocks(matrix_cov, blocks=elements, cmap=cmap, vmax=vmax) - if pdf: + hm = pcolor_with_blocks(matrix_cov, blocks=elements, cmap=cmap, vmax=vmax) + if pdf and hm is not None: pdf.savefig() + plt.figure(figsize=figsize_2d) + ax = plt.subplot( + 111, + xlabel="", + ylabel="bin", + title=f"Covariance matrix {name} (symlog){title_suffix}", + ) + hm = pcolor_with_blocks( + matrix_cov, blocks=elements, cmap=cmap, symlog=True + ) + if pdf and hm is not None: + pdf.savefig() + plt.figure(figsize=figsize_2d) ax = plt.subplot( 111, @@ -179,8 +197,8 @@ def main(opts: Namespace) -> None: ylabel="bin", title=f"Covariance matrix {name} (blocks){title_suffix}", ) - pcolor_with_blocks(bmatrix_cov, blocks=elements, cmap=cmap) - if pdf: + hm = pcolor_with_blocks(bmatrix_cov, blocks=elements, cmap=cmap) + if pdf and hm is not None: pdf.savefig() plt.figure(figsize=figsize_2d) @@ -188,12 +206,39 @@ def main(opts: Namespace) -> None: 111, xlabel="", ylabel="bin", - title=f"Covariance matrix {name} ({elements[0]}){title_suffix}", + title=f"Covariance matrix {name} (blocks, symlog){title_suffix}", ) - pcolor_with_blocks( - matrix_cov[:blocksize, :blocksize], blocks=elements[:1], cmap=cmap + hm = pcolor_with_blocks(bmatrix_cov, blocks=elements, cmap=cmap, symlog=True) + if pdf and hm is not None: + pdf.savefig() + + plt.figure(figsize=figsize_2d) + ax = plt.subplot( + 111, + xlabel="", + ylabel="bin", + title=f"Covariance matrix {name} (block {block_name}){title_suffix}", ) - if pdf: + hm = pcolor_with_blocks( + matrix_cov[block_selection,block_selection], blocks=elements[:1], cmap=cmap + ) + if pdf and hm is not None: + pdf.savefig() + + plt.figure(figsize=figsize_2d) + ax = plt.subplot( + 111, + xlabel="", + ylabel="bin", + title=f"Covariance matrix {name} (block {block_name}, symlog){title_suffix}", + ) + hm = pcolor_with_blocks( + matrix_cov[block_selection,block_selection], + blocks=elements[:1], + cmap=cmap, + symlog=True, + ) + if pdf and hm is not None: pdf.savefig() plt.figure(figsize=figsize_2d) @@ -203,8 +248,8 @@ def main(opts: Namespace) -> None: ylabel="bin", title=rf"Relative covariance matrix {name}, %虏{title_suffix}", ) - pcolor_with_blocks(matrix_cov_rel, blocks=elements, cmap=cmap) - if pdf: + hm = pcolor_with_blocks(matrix_cov_rel, blocks=elements, cmap=cmap) + if pdf and hm is not None: pdf.savefig() if vmax_rel is not None: @@ -213,12 +258,12 @@ def main(opts: Namespace) -> None: 111, xlabel="", ylabel="bin", - title=rf"Relative covariance matrix {name}, %虏{title_suffix}", + title=rf"Relative covariance matrix {name} (truncated), %虏{title_suffix}", ) - pcolor_with_blocks( + hm = pcolor_with_blocks( matrix_cov_rel, blocks=elements, cmap=cmap, vmax=vmax_rel ) - if pdf: + if pdf and hm is not None: pdf.savefig() plt.figure(figsize=figsize_2d) @@ -226,12 +271,41 @@ def main(opts: Namespace) -> None: 111, xlabel="", ylabel="bin", - title=f"Relative covariance matrix {name} ({elements[0]}){title_suffix}", + title=rf"Relative covariance matrix {name} (symlog), %虏{title_suffix}", ) - pcolor_with_blocks( - matrix_cov_rel[:blocksize, :blocksize], blocks=elements[:1], cmap=cmap + hm = pcolor_with_blocks( + matrix_cov_rel, blocks=elements, cmap=cmap, symlog=True ) - if pdf: + if pdf and hm is not None: + pdf.savefig() + + plt.figure(figsize=figsize_2d) + ax = plt.subplot( + 111, + xlabel="", + ylabel="bin", + title=f"Relative covariance matrix {name} (block {block_name}){title_suffix}", + ) + hm = pcolor_with_blocks( + matrix_cov_rel[block_selection,block_selection], blocks=elements[:1], cmap=cmap + ) + if pdf and hm is not None: + pdf.savefig() + + plt.figure(figsize=figsize_2d) + ax = plt.subplot( + 111, + xlabel="", + ylabel="bin", + title=f"Relative covariance matrix {name} (block {block_name}, symlog){title_suffix}", + ) + hm = pcolor_with_blocks( + matrix_cov_rel[block_selection,block_selection], + blocks=elements[:1], + cmap=cmap, + symlog=True, + ) + if pdf and hm is not None: pdf.savefig() plt.figure(figsize=figsize_2d) @@ -241,8 +315,8 @@ def main(opts: Namespace) -> None: ylabel="bin", title=f"Correlation matrix {name}{title_suffix}", ) - pcolor_with_blocks(matrix_cor, blocks=elements, cmap=cmap) - if pdf: + hm = pcolor_with_blocks(matrix_cor, blocks=elements, cmap=cmap) + if pdf and hm is not None: pdf.savefig() plt.figure(figsize=figsize_2d) @@ -250,16 +324,16 @@ def main(opts: Namespace) -> None: 111, xlabel="", ylabel="bin", - title=f"Correlation matrix {name} ({elements[0]}){title_suffix}", + title=f"Correlation matrix {name} (block {block_name}){title_suffix}", ) hm = pcolor_with_blocks( - matrix_cor[:blocksize, :blocksize], + matrix_cor[block_selection,block_selection], blocks=elements[:1], pcolormesh=True, cmap=cmap, ) # heatmap_show_values(hm, lower_triangle=True) - if pdf: + if pdf and hm is not None: pdf.savefig() plt.figure(figsize=figsize_2d) @@ -276,7 +350,7 @@ def main(opts: Namespace) -> None: logger.info(f"Plot {name}") - if pdf: + if pdf and hm is not None: pdf.savefig() if opts.show: @@ -285,7 +359,7 @@ def main(opts: Namespace) -> None: if not opts.show: plt.close("all") - if pdf: + if pdf and hm is not None: pdf.__exit__(None, None, None) logger.info(f"Write output file: {ofile}") @@ -394,28 +468,41 @@ def pcolor_with_blocks( pcolormesh: bool = False, sep_kwargs: Mapping = {}, vmax: float | None = None, + symlog: bool = False, **kwargs, ): from numpy import fabs - from numpy.ma import array dmin = data.min() dmax = data.max() + if dmin==dmax==0: + return None + bound = max(fabs(dmin), dmax) if vmax is not None: vmin, vmax = -vmax, vmax else: vmin, vmax = -bound, bound + if symlog: + pcolor_kwargs = dict( + norm=SymLogNorm( + linthresh=bound * 0.01, linscale=0.05, vmin=vmin, vmax=vmax + ), + **kwargs, + ) + else: + pcolor_kwargs = dict(vmin=vmin, vmax=vmax, **kwargs) + # fdata = fabs(data) # data = array(data, mask=(fdata < 1.0e-9)) ax = plt.gca() ax.set_aspect("equal") if pcolormesh: - hm = ax.pcolormesh(data, *args, vmin=vmin, vmax=vmax, **kwargs) + hm = ax.pcolormesh(data, *args, **pcolor_kwargs) else: - hm = ax.pcolorfast(data, *args, vmin=vmin, vmax=vmax, **kwargs) + hm = ax.pcolorfast(data, *args, **pcolor_kwargs) hm.set_rasterized(True) if colorbar: add_colorbar(hm, rasterized=True)