Coverage for python / lsst / jointcal / check_logged_chi2.py: 0%
133 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-05-05 08:42 +0000
1# This file is part of jointcal.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
21"""
22Extract chi2 and degrees of freedom values logged by one or more jointcal runs,
23print warnings about oddities, and make plots.
24"""
26import argparse
27import dataclasses
28import itertools
29import os.path
30import re
32import numpy as np
33import matplotlib
34matplotlib.use("Agg")
35import matplotlib.pyplot as plt # noqa: E402
36import seaborn as sns # noqa: E402
37sns.set_style("ticks", {"legend.frameon": True})
38sns.set_context("talk")
41@dataclasses.dataclass
42class Chi2Data:
43 """Store the chi2 values read in from a jointcal log file.
44 """
45 kind: list()
46 raw: np.ndarray
47 ndof: np.ndarray
48 reduced: np.ndarray
49 init_count: int = dataclasses.field(init=False)
51 def __post_init__(self):
52 # ensure the array types are correct
53 self.raw = np.array(self.raw, dtype=np.float64)
54 self.ndof = np.array(self.ndof, dtype=np.int)
55 self.reduced = np.array(self.reduced, dtype=np.float64)
56 self.init_count = self._find_init()
58 def _find_init(self):
59 """Return the index of the first "fit step", after initialization.
61 NOTE
62 ----
63 There are never more than ~25 items in the list, so search optimization
64 is not worth the trouble.
65 """
66 # Logs pre-DM-25779
67 if "Fit prepared" in self.kind:
68 return self.kind.index("Fit prepared") + 1
69 # Logs post-DM-25779
70 elif "Fit iteration 0" in self.kind:
71 return self.kind.index("Fit iteration 0")
72 else:
73 raise RuntimeError(f"Cannot find end of initialization sequence in {self.kind}")
76class LogParser:
77 """Parse a jointcal logfile to extract chi2 values and plot them.
79 Call the instance with the path to a file to check it for anamolous chi2
80 and output plots to your current directory.
82 Parameters
83 ----------
84 plot : `bool`
85 Make plots for each file (saved to the current working directory)?
86 verbose : `bool`
87 Print extra updates during processing?
88 """
89 def __init__(self, plot=True, verbose=True):
90 # This regular expression extracts the chi2 values, and the "kind" of
91 # chi2 (e.g. "Initial", "Fit iteration").
92 # Chi2 values in the log look like this, for example:
93 # jointcal INFO: Initial chi2/ndof : 2.50373e+16/532674=4.7003e+10
94 chi2_re = "jointcal INFO: (?P<kind>.+) chi2/ndof : (?P<chi2>.+)/(?P<ndof>.+)=(?P<reduced_chi2>.+)"
95 self.matcher = re.compile(chi2_re)
96 self.plot = plot
97 self.verbose = verbose
99 # Reuse the Figure to speed up plotting and save memory.
100 self.fig = plt.figure(figsize=(15, 8))
102 # How to find the beginning and end of the relevant parts of the log
103 # to scan for chi2 values.
104 self.section_start = {"astrometry": "Starting astrometric fitting...",
105 "photometry": "Starting photometric fitting..."}
106 self.section_end = {"astrometry": "Updating WCS for visit:",
107 "photometry": "Updating PhotoCalib for visit:"}
109 def __call__(self, logfile):
110 """Parse logfile to extract chi2 values and generate and save plots.
112 The plot output is written to the current directory, with the name
113 derived from the basename of ``logfile``.
115 Parameters
116 ----------
117 logfile : `str`
118 The filename of the jointcal log to process.
119 """
120 title = os.path.basename(logfile)
121 if self.verbose:
122 print("Processing:", title)
124 with open(logfile) as opened_log:
125 # Astrometry is always run first, so we can scan for that until the
126 # end of that section, and then continue scanning for photometry.
127 astrometry = self._extract_chi2(opened_log, "astrometry")
128 increased = self._find_chi2_increase(astrometry, title, "astrometry")
129 photometry = self._extract_chi2(opened_log, "photometry")
130 increased |= self._find_chi2_increase(photometry, title, "photometry")
132 if astrometry is None and photometry is None and self.verbose:
133 print(f"WARNING: No chi2 values found in {logfile}.")
135 if increased or self.plot:
136 self._plot(astrometry, photometry, title)
137 plotfile = f"{os.path.splitext(title)[0]}.png"
138 plt.savefig(plotfile, bbox_inches="tight")
139 print("Saved plot:", plotfile)
141 def _find_chi2_increase(self, chi2Data, title, label, threshold=1):
142 """Return True and print a message if the raw chi2 increases
143 markedly.
144 """
145 if chi2Data is None:
146 return False
147 diff = np.diff(chi2Data.raw)
148 ratio = diff/chi2Data.raw[:-1]
149 if np.any(ratio > threshold):
150 increased = np.where(ratio > threshold)[0]
151 print(f"{title} has increasing {label} chi2:")
152 for x in zip(chi2Data.raw[increased], chi2Data.raw[increased + 1],
153 ratio[increased], diff[increased]):
154 print(f"{x[0]:.6} -> {x[1]:.6} (ratio: {x[2]:.6}, diff: {x[3]:.6})")
155 return True
156 return False
158 def _extract_chi2(self, opened_log, section):
159 """Return the values extracted from the chi2 statements in the logfile.
160 """
161 start = self.section_start[section]
162 end = self.section_end[section]
163 kind = []
164 chi2 = []
165 ndof = []
166 reduced = []
167 # Skip over lines until we get to the section start line.
168 for line in opened_log:
169 if start in line:
170 break
172 for line in opened_log:
173 # Stop parsing at the section end line.
174 if end in line:
175 break
176 if "chi2" in line:
177 match = self.matcher.search(line)
178 if match is not None:
179 kind.append(match.group("kind"))
180 chi2.append(match.group("chi2"))
181 ndof.append(match.group("ndof"))
182 reduced.append(match.group("reduced_chi2"))
184 # No chi2 values were found (e.g., photometry wasn't run).
185 if len(kind) == 0:
186 return None
188 return Chi2Data(kind, np.array(chi2, dtype=np.float64),
189 np.array(ndof, dtype=int), np.array(reduced, dtype=np.float64))
191 def _plot(self, astrometry, photometry, title):
192 """Generate plots of chi2 values.
194 Parameters
195 ----------
196 astrometry : `Chi2Data` or None
197 The as-read astrometry data, or None if there is none to plot.
198 photometry : `Chi2Data` or None
199 The as-read photometry data, or None if there is none to plot.
200 title : `str`
201 Title for the whole plot.
202 """
203 palette = itertools.cycle(sns.color_palette())
205 self.fig.clf()
206 ax0, ax1 = self.fig.subplots(ncols=2, gridspec_kw={"wspace": 0.05})
208 self.fig.suptitle(title)
209 # Use a log scale if any of the chi2 values are very large.
210 if max(getattr(astrometry, "raw", [0])) > 100 or max(getattr(photometry, "raw", [0])) > 100:
211 ax0.set_yscale("log")
212 ax1.yaxis.set_label_position("right")
213 ax1.yaxis.tick_right()
215 if astrometry is not None:
216 patch1, patch2 = self._plot_axes(ax0, ax1, astrometry, palette, label="astrometry")
218 if photometry is not None:
219 patch3, patch4 = self._plot_axes(ax0, ax1, photometry, palette, label="photometry")
221 # Let matplotlib figure out the best legend location: if there is data
222 # in the "upper right", we definitely want to see it.
223 handles, labels = ax0.get_legend_handles_labels()
224 ax1.legend(handles, labels)
226 def _plot_axes(self, ax0, ax1, chi2Data, palette, label=""):
227 """Make the chi2 and degrees of freedom subplots."""
228 xrange = np.arange(0, len(chi2Data.raw), dtype=float)
230 # mark chi2=1
231 ax0.axhline(1, color='grey', ls='--')
232 # mark the separation between initialization and iteration
233 ax0.axvline(chi2Data.init_count-0.5, color='grey', lw=0.9)
234 color = next(palette)
235 patch1 = ax0.plot(xrange[:chi2Data.init_count], chi2Data.raw[:chi2Data.init_count], '*', ms=10,
236 label=f"{label} pre-init", color=color)
237 patch2 = ax0.plot(xrange[chi2Data.init_count:], chi2Data.raw[chi2Data.init_count:], 'o', ms=10,
238 label=f"{label} post-init", color=color)
239 patch1 = ax0.plot(xrange[:chi2Data.init_count], chi2Data.reduced[:chi2Data.init_count], '*',
240 markerfacecolor="none", ms=10, color=color)
241 patch2 = ax0.plot(xrange[chi2Data.init_count:], chi2Data.reduced[chi2Data.init_count:], 'o',
242 markerfacecolor="none", ms=10, label=f"{label} reduced", color=color)
244 ax0.set_xlabel("Iteration #", fontsize=20)
245 ax0.set_ylabel(r"$\chi ^2$", fontsize=20)
247 # mark the separation between initialization and iteration
248 ax1.axvline(chi2Data.init_count-0.5, color='grey', lw=0.9)
249 ax1.plot(xrange[:chi2Data.init_count], chi2Data.ndof[:chi2Data.init_count], '*', ms=10,
250 label="pre-init", color=color)
251 ax1.plot(xrange[chi2Data.init_count:], chi2Data.ndof[chi2Data.init_count:], 'o', ms=10,
252 label="post-init", color=color)
254 ax1.set_xlabel("Iteration #", fontsize=20)
255 ax1.set_ylabel("# degrees of freedom", fontsize=20)
257 return patch1[0], patch2[0]
260def parse_args():
261 parser = argparse.ArgumentParser(description=__doc__,
262 formatter_class=argparse.RawDescriptionHelpFormatter)
263 parser.add_argument("files", metavar="files", nargs="+",
264 help="Log file(s) to extract chi2 values from.")
265 parser.add_argument("--plot", action="store_true",
266 help="Generate a plot PNG for each log file, otherwise just for questionable ones.")
267 parser.add_argument("-v", "--verbose", action="store_true",
268 help="Print extra information during processing.")
269 return parser.parse_args()
272def main():
273 args = parse_args()
274 log_parser = LogParser(plot=args.plot, verbose=args.verbose)
275 for file in args.files:
276 log_parser(file)