Merge pull request #54 from jzuhone/cea

acis_thermal_check/apps/acisfp_check.py

@@ -221,7 +221,7 @@ def make_prediction_plots(self, outdir, states, temps, load_start):
 
         # Now plot any perigee passages that occur between xmin and xmax
         # for eachpassage in perigee_passages:
-        paint_perigee(self.perigee_passages, states, plots)
+        paint_perigee(self.perigee_passages, plots)
 
         plots['default'] = plots[f"{self.name}_3"]
 

acis_thermal_check/apps/cea_check.py

@@ -0,0 +1,138 @@
+#!/usr/bin/env python
+
+"""
+========================
+cea_check
+========================
+
+This code generates backstop load review outputs for checking the HRC
+CEA temperature 2CEAHVPT. It also generates CEA model validation
+plots comparing predicted values to telemetry for the previous three
+weeks.
+"""
+
+# Matplotlib setup
+# Use Agg backend for command-line (non-interactive) operation
+import matplotlib
+matplotlib.use('Agg')
+
+import sys
+from acis_thermal_check import \
+    ACISThermalCheck, \
+    get_options, \
+    mylog
+from acis_thermal_check.utils import PredictPlot
+from Ska.Matplotlib import pointpair
+
+
+class CEACheck(ACISThermalCheck):
+    def __init__(self):
+        valid_limits = {'2CEAHVPT': [(1, 2.0), (50, 1.0), (99, 2.0)],
+                        'PITCH': [(1, 3.0), (99, 3.0)],
+                        'TSCPOS': [(1, 2.5), (99, 2.5)]
+                        }
+        hist_limit = [5.0]
+        limits_map = {}
+        other_telem = ["2imonst", "2sponst", "2s2onst"]
+        super(CEACheck, self).__init__("2ceahvpt", "cea", valid_limits,
+                                       hist_limit, limits_map=limits_map, 
+                                       other_telem=other_telem)
+
+    def make_prediction_viols(self, temps, states, load_start):
+        """
+        Find limit violations where predicted temperature is above the
+        specified limits.
+
+        Parameters
+        ----------
+        temps : dict of NumPy arrays
+            NumPy arrays corresponding to the modeled temperatures
+        states : NumPy record array
+            Commanded states
+        load_start : float
+            The start time of the load, used so that we only report
+            violations for times later than this time for the model
+            run.
+        """
+        mylog.info('Checking for limit violations')
+
+        temp = temps[self.name]
+        times = self.predict_model.times
+
+        # Only check this violation when HRC is on
+        mask = self.predict_model.comp['2imonst_on'].dvals
+        mask |= self.predict_model.comp['2sponst_on'].dvals
+        hi_viols = self._make_prediction_viols(
+            times, temp, load_start, self.limits["planning_hi"].value,
+            "planning", "max", mask=mask)
+        viols = {"hi":
+                     {"name": f"Hot ({self.limits['planning_hi'].value} C)",
+                      "type": "Max",
+                      "values": hi_viols}
+                 }
+        return viols
+
+    def _calc_model_supp(self, model, state_times, states, ephem, state0):
+        """
+        Update to initialize the cea0 pseudo-node. If 2ceahvpt
+        has an initial value (T_cea) - which it does at
+        prediction time (gets it from state0), then T_cea0 
+        is set to that.  If we are running the validation,
+        T_cea is set to None so we use the dvals in model.comp
+        """
+        for node in ["cea0", "cea1"]:
+            if state0 is None:
+                T_cea = model.comp["2ceahvpt"].dvals
+            else:
+                T_cea = state0["2ceahvpt"]
+            model.comp[node].set_data(T_cea, model.times)
+        model.comp["2ps5aon_on"].set_data(True)
+        model.comp["2ps5bon_on"].set_data(False)
+        model.comp["2imonst_on"].set_data(states["hrc_i"] == "ON", state_times)
+        model.comp["2sponst_on"].set_data(states["hrc_s"] == "ON", state_times)
+        model.comp["2s2onst_on"].set_data(states["hrc_15v"] == "ON", state_times)
+        model.comp["224pcast_off"].set_data(states["hrc_15v"] == "ON", state_times)
+        model.comp["215pcast_off"].set_data(states["hrc_15v"] == "ON", state_times)
+
+    def _make_state_plots(self, plots, num_figs, w1, plot_start,
+                          states, load_start):
+        # Make a plot of ACIS HRC states
+        plots['hrc'] = PredictPlot(
+            fig_id=num_figs+1,
+            title='HRC States',
+            xlabel='Date',
+            x=pointpair(states['tstart'], states['tstop']),
+            y=pointpair(states['hrc_i']),
+            yy=pointpair(states['hrc_s']),
+            ylabel='HRC-I/S',
+            x2=pointpair(states['tstart'], states['tstop']),
+            y2=pointpair(states['hrc_15v']),
+            ylabel2='HRC 15 V',
+            linewidth2=4.0,
+            xmin=plot_start,
+            width=w1, load_start=load_start)
+        plots['hrc'].ax.lines[0].set_label('HRC-I')
+        plots['hrc'].ax.lines[1].set_label('HRC-S')
+        plots['hrc'].ax.legend(fancybox=True, framealpha=0.5, loc=2)
+        plots['hrc'].filename = 'hrc.png'
+
+        num_figs += 1
+        super()._make_state_plots(plots, num_figs, w1, plot_start,
+                                  states, load_start)
+
+
+def main():
+    args = get_options(use_acis_opts=False)
+    cea_check = CEACheck()
+    try:
+        cea_check.run(args)
+    except Exception as msg:
+        if args.traceback:
+            raise
+        else:
+            print("ERROR:", msg)
+            sys.exit(1)
+
+
+if __name__ == '__main__':
+    main()

acis_thermal_check/main.py

@@ -125,7 +125,7 @@ def __init__(self, msid, name, validation_limits, hist_limit,
         if hist_ops is None:
             hist_ops = ["greater_equal"]*len(hist_limit)
         self.hist_ops = hist_ops
-        self.perigee_passages = []
+        self.perigee_passages = defaultdict(list)
         self.write_pickle = False
         self.limits = {}
 
@@ -163,7 +163,13 @@ def run(self, args, override_limits=None):
         proc = self._setup_proc_and_logger(args, model_spec)
 
         # Record the selected state builder in the class attributes
-        self.state_builder = make_state_builder(args.state_builder, args)
+        # If there is no "state_builder" command line argument assume
+        # kadi
+        hrc_states = self.name in ["cea"]
+        state_builder = getattr(args, "state_builder", "kadi")
+        mylog.info(f"ACISThermalCheck is using the '{state_builder}' state builder.")
+        self.state_builder = make_state_builder(state_builder, args, 
+                                                hrc_states=hrc_states)
 
         # If args.run_start is not none, write validation and prediction
         # data to a pickle later
@@ -615,59 +621,26 @@ def write_temps(self, outdir, times, temps):
         temp_table.write(outfile, format='ascii', delimiter='\t', overwrite=True)
 
     def _gather_perigee(self, plot_start, load_start):
-        import glob
-        # The first step is to build a list of all the perigee passages.
-
         # Gather the perigee passages that occur from the
         # beginning of the model run up to the start of the load
         # from kadi
         rzs = events.rad_zones.filter(plot_start, load_start)
         for rz in rzs:
-            self.perigee_passages.append([rz.start, rz.perigee])
-        for rz in rzs:
-            self.perigee_passages.append([rz.stop, rz.perigee])
-
-        # We will get the load passages from the relevant CRM pad time file
-        # (e.g. DO12143_CRM_Pad.txt) inside the bsdir directory
-        # Each line is either an inbound or outbound ECS
-        #
-        # The reason we are doing this is because we want to draw vertical
-        # lines denoting each perigee passage on the plots
-        #
-        # Open the file
-        crm_file_path = glob.glob(f"{self.bsdir}/*CRM*")[0]
-        crm_file = open(crm_file_path, 'r')
-
-        alines = crm_file.readlines()
-
-        idx = None
-        # Keep reading until you hit the last header line which is all "*"'s
-        for i, aline in enumerate(alines):
-            if len(aline) > 0 and aline[0] == "*":
-                idx = i+1
-                break
-
-        if idx is None:
-            raise RuntimeError("Couldn't find the end of the CRM Pad Time file header!")
-
-        # Found the last line of the header. Start processing Perigee Passages
-
-        # While there are still lines to be read
-        for aline in alines[idx:]:
-            # create an empty Peri. Passage instance location
-            passage = []
-
-            # split the CRM Pad Time file line read in and extract the
-            # relevant information
-            splitline = aline.split()
-            passage.append(splitline[6])  # Radzone entry/exit
-            passage.append(splitline[9])  # Perigee Passage time
-
-            # append this passage to the passages list
-            self.perigee_passages.append(passage)
-
-        # Done with the CRM Pad Time file - close it
-        crm_file.close()
+            self.perigee_passages["entry"].append(rz.start)
+            self.perigee_passages["perigee"].append(rz.perigee)
+            self.perigee_passages["exit"].append(rz.stop)
+
+        # obtain the rest from the backstop file
+        for cmd in self.state_builder.bs_cmds:
+            if cmd["tlmsid"] == "OORMPDS":
+                key = "entry"
+            elif cmd["tlmsid"] == "OORMPEN":
+                key = "exit"
+            elif cmd["type"] == "ORBPOINT" and cmd["event_type"] == "EPERIGEE":
+                key = "perigee"
+            else:
+                continue
+            self.perigee_passages[key].append(cmd["time"])
 
     def _make_state_plots(self, plots, num_figs, w1, plot_start,
                           states, load_start):
@@ -796,7 +769,7 @@ def make_prediction_plots(self, outdir, states, temps, load_start):
 
         # Now plot any perigee passages that occur between xmin and xmax
         # for eachpassage in perigee_passages:
-        paint_perigee(self.perigee_passages, states, plots)
+        paint_perigee(self.perigee_passages, plots)
 
         # Now write all of the plots after possible
         # customizations have been made
@@ -1069,6 +1042,35 @@ def make_validation_plots(self, tlm, model_spec, outdir):
 
         fig_id += 1
 
+        if self.name == "cea":
+            for msid in ["2imonst", "2sponst", "2s2onst"]:
+                fig = plt.figure(10 + fig_id, figsize=(12, 6))
+                fig.clf()
+                comp_hrc = model.comp[f"{msid}_on"].dvals
+                tlm_hrc = np.char.strip(tlm[msid]) == "ON"
+                ticklocs, fig, ax = plot_cxctime(model.times, comp_hrc, label="Model",
+                                                 fig=fig, ls='-', lw=4, color=thermal_red, zorder=9)
+                ticklocs, fig, ax = plot_cxctime(model.times, tlm_hrc, label="Data",
+                                                 fig=fig, ls='-', lw=2, color=thermal_blue, zorder=10)
+                ax.grid()
+                ax.set_xlim(xmin, xmax)
+                ax.set_yticks([0, 1])
+                ax.set_yticklabels(["OFF", "ON"])
+                ax.set_ylim([-0.1, 1.1])
+                # add lines for perigee passages
+                for rz in rzs:
+                    ptimes = cxctime2plotdate([rz.tstart, rz.tstop])
+                    for ptime in ptimes:
+                        ax.axvline(ptime, ls='--', color='C2',
+                                   linewidth=2, zorder=2)
+                ax.legend(fancybox=True, framealpha=0.5, loc=2)
+                plots[msid] = {
+                    "lines": {"fig": fig,
+                              "ax": ax,
+                              "filename": f'{msid}_valid.png'}
+                }
+                fig_id += 1 
+
         if 'earthheat__fptemp' in model.comp:
 
             fig = plt.figure(10 + fig_id, figsize=(12, 6))
@@ -1249,8 +1251,6 @@ def _setup_proc_and_logger(self, args, model_spec):
         mylog.info('# model_spec file MD5sum = %s' % md5sum)
         mylog.info('Command line options:\n%s\n' % pformat(args.__dict__))
 
-        mylog.info("ACISThermalCheck is using the '%s' state builder." % args.state_builder)
-
         if args.backstop_file is None:
             self.bsdir = None
         else:

acis_thermal_check/state_builder.py

@@ -25,6 +25,7 @@ def __init__(self, logger=None):
             # Make a logger but with no output
             logger = logging.getLogger('statebuilder-no-logger')
         self.logger = logger
+        self.state_keys = STATE_KEYS.copy()
 
     def get_prediction_states(self, tlm):
         """
@@ -54,7 +55,7 @@ def get_validation_states(self, datestart, datestop):
                          (start.date, stop.date))
 
         with kadi_states.disable_grating_move_duration():
-            states = kadi_states.get_states(start, stop, state_keys=STATE_KEYS,
+            states = kadi_states.get_states(start, stop, state_keys=self.state_keys,
                                             merge_identical=True)
 
         # Set start and end state date/times to match telemetry span.  Extend the
@@ -79,7 +80,7 @@ class KadiStateBuilder(StateBuilder):
     be used for validation only.
     """
     def __init__(self, interrupt=False, backstop_file=None,
-                 logger=None):
+                 logger=None, hrc_states=False):
         """
         Give the KadiStateBuilder arguments that were passed in
         from the command line, and set up the connection to the
@@ -94,8 +95,12 @@ def __init__(self, interrupt=False, backstop_file=None,
             file will be searched for within this directory.
         logger : Logger object, optional
             The Python Logger object to be used when logging.
+        hrc_states : boolean, optional
+            Whether or not to add HRC-specific states. Default: False
         """
         super().__init__(logger=logger)
+        if hrc_states:
+            self.state_keys += ["hrc_15v", "hrc_i", "hrc_s"]
 
         # Note: `interrupt` is ignored in this class. This concept is not needed
         # since backstop 6.9, which provides the RUNNING_LOAD_TERMINATION_TIME
@@ -181,7 +186,7 @@ def get_prediction_states(self, tbegin):
         # but this could probably be set to True.
         with kadi_states.disable_grating_move_duration():
             states = kadi_states.get_states(cmds=cmds, start=tbegin, stop=sched_stop,
-                                            state_keys=STATE_KEYS,
+                                            state_keys=self.state_keys,
                                             merge_identical=False)
 
         # Make the column order match legacy Chandra.cmd_states.
@@ -311,7 +316,7 @@ def get_prediction_states(self, tbegin):
         with kadi_states.disable_grating_move_duration():
             states = kadi_states.get_states(cmds=bs_cmds, start=tbegin,
                                             stop=sched_stop,
-                                            state_keys=STATE_KEYS)
+                                            state_keys=self.state_keys)
 
         # Make the column order match legacy Chandra.cmd_states.
         states = states[sorted(states.colnames)]

acis_thermal_check/templates/index_template.rst

@@ -74,6 +74,9 @@ No {{proc.msid}} {{viols[key]["name"]}} Violations
 {% endfor %}
 
 .. image:: {{plots.default.filename}}
+{% if proc.msid == "2CEAHVPT" %}
+.. image:: {{plots.hrc.filename}}
+{% endif %}
 .. image:: {{plots.pow_sim.filename}}
 {% if proc.msid == "FPTEMP" %}
 .. image:: {{plots.roll_taco.filename}}
@@ -131,8 +134,15 @@ No Validation Violations
 
 {% if msid == "ccd_count" %}
 
-CCD/FEP Count
--------------
+ACIS CCD/FEP Count
+------------------
+
+.. image:: {{plot.lines.filename}}
+
+{% elif msid in ["2imonst", "2sponst", "2s2onst"] %}
+
+{{ msid.upper() }}
+---------------------
 
 .. image:: {{plot.lines.filename}}