LSST Applications  21.0.0-172-gfb10e10a+18fedfabac,22.0.0+297cba6710,22.0.0+80564b0ff1,22.0.0+8d77f4f51a,22.0.0+a28f4c53b1,22.0.0+dcf3732eb2,22.0.1-1-g7d6de66+2a20fdde0d,22.0.1-1-g8e32f31+297cba6710,22.0.1-1-geca5380+7fa3b7d9b6,22.0.1-12-g44dc1dc+2a20fdde0d,22.0.1-15-g6a90155+515f58c32b,22.0.1-16-g9282f48+790f5f2caa,22.0.1-2-g92698f7+dcf3732eb2,22.0.1-2-ga9b0f51+7fa3b7d9b6,22.0.1-2-gd1925c9+bf4f0e694f,22.0.1-24-g1ad7a390+a9625a72a8,22.0.1-25-g5bf6245+3ad8ecd50b,22.0.1-25-gb120d7b+8b5510f75f,22.0.1-27-g97737f7+2a20fdde0d,22.0.1-32-gf62ce7b1+aa4237961e,22.0.1-4-g0b3f228+2a20fdde0d,22.0.1-4-g243d05b+871c1b8305,22.0.1-4-g3a563be+32dcf1063f,22.0.1-4-g44f2e3d+9e4ab0f4fa,22.0.1-42-gca6935d93+ba5e5ca3eb,22.0.1-5-g15c806e+85460ae5f3,22.0.1-5-g58711c4+611d128589,22.0.1-5-g75bb458+99c117b92f,22.0.1-6-g1c63a23+7fa3b7d9b6,22.0.1-6-g50866e6+84ff5a128b,22.0.1-6-g8d3140d+720564cf76,22.0.1-6-gd805d02+cc5644f571,22.0.1-8-ge5750ce+85460ae5f3,master-g6e05de7fdc+babf819c66,master-g99da0e417a+8d77f4f51a,w.2021.48
LSST Data Management Base Package
check_logged_chi2.py
Go to the documentation of this file.
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 """
22 Extract chi2 and degrees of freedom values logged by one or more jointcal runs,
23 print warnings about oddities, and make plots.
24 """
25 
26 import argparse
27 import dataclasses
28 import itertools
29 import os.path
30 import re
31 
32 import numpy as np
33 import matplotlib
34 matplotlib.use("Agg")
35 import matplotlib.pyplot as plt # noqa: E402
36 import seaborn as sns # noqa: E402
37 sns.set_style("ticks", {"legend.frameon": True})
38 sns.set_context("talk")
39 
40 
41 @dataclasses.dataclass
42 class 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)
50 
51  def __post_init__(self):
52  # ensure the array types are correct
53  self.rawraw = np.array(self.rawraw, dtype=np.float64)
54  self.ndofndof = np.array(self.ndofndof, dtype=np.int)
55  self.reducedreduced = np.array(self.reducedreduced, dtype=np.float64)
56  self.init_countinit_count = self._find_init_find_init()
57 
58  def _find_init(self):
59  """Return the index of the first "fit step", after initialization.
60 
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}")
74 
75 
76 class LogParser:
77  """Parse a jointcal logfile to extract chi2 values and plot them.
78 
79  Call the instance with the path to a file to check it for anamolous chi2
80  and output plots to your current directory.
81 
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.matchermatcher = re.compile(chi2_re)
96  self.plotplot = plot
97  self.verboseverbose = verbose
98 
99  # Reuse the Figure to speed up plotting and save memory.
100  self.figfig = plt.figure(figsize=(15, 8))
101 
102  # How to find the beginning and end of the relevant parts of the log
103  # to scan for chi2 values.
104  self.section_startsection_start = {"astrometry": "Starting astrometric fitting...",
105  "photometry": "Starting photometric fitting..."}
106  self.section_endsection_end = {"astrometry": "Updating WCS for visit:",
107  "photometry": "Updating PhotoCalib for visit:"}
108 
109  def __call__(self, logfile):
110  """Parse logfile to extract chi2 values and generate and save plots.
111 
112  The plot output is written to the current directory, with the name
113  derived from the basename of ``logfile``.
114 
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.verboseverbose:
122  print("Processing:", title)
123 
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_extract_chi2(opened_log, "astrometry")
128  increased = self._find_chi2_increase_find_chi2_increase(astrometry, title, "astrometry")
129  photometry = self._extract_chi2_extract_chi2(opened_log, "photometry")
130  increased |= self._find_chi2_increase_find_chi2_increase(photometry, title, "photometry")
131 
132  if astrometry is None and photometry is None and self.verboseverbose:
133  print(f"WARNING: No chi2 values found in {logfile}.")
134 
135  if increased or self.plotplot:
136  self._plot_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)
140 
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
157 
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_startsection_start[section]
162  end = self.section_endsection_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
171 
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.matchermatcher.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"))
183 
184  # No chi2 values were found (e.g., photometry wasn't run).
185  if len(kind) == 0:
186  return None
187 
188  return Chi2Data(kind, np.array(chi2, dtype=np.float64),
189  np.array(ndof, dtype=int), np.array(reduced, dtype=np.float64))
190 
191  def _plot(self, astrometry, photometry, title):
192  """Generate plots of chi2 values.
193 
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())
204 
205  self.figfig.clf()
206  ax0, ax1 = self.figfig.subplots(ncols=2, gridspec_kw={"wspace": 0.05})
207 
208  self.figfig.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()
214 
215  if astrometry is not None:
216  patch1, patch2 = self._plot_axes_plot_axes(ax0, ax1, astrometry, palette, label="astrometry")
217 
218  if photometry is not None:
219  patch3, patch4 = self._plot_axes_plot_axes(ax0, ax1, photometry, palette, label="photometry")
220 
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)
225 
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)
229 
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)
243 
244  ax0.set_xlabel("Iteration #", fontsize=20)
245  ax0.set_ylabel(r"$\chi ^2$", fontsize=20)
246 
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)
253 
254  ax1.set_xlabel("Iteration #", fontsize=20)
255  ax1.set_ylabel("# degrees of freedom", fontsize=20)
256 
257  return patch1[0], patch2[0]
258 
259 
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()
270 
271 
272 def 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)
int max
def _plot(self, astrometry, photometry, title)
def _extract_chi2(self, opened_log, section)
def __init__(self, plot=True, verbose=True)
def _find_chi2_increase(self, chi2Data, title, label, threshold=1)
def _plot_axes(self, ax0, ax1, chi2Data, palette, label="")