diff --git a/proseco/acq.py b/proseco/acq.py index ef2747b2..562b7bf3 100644 --- a/proseco/acq.py +++ b/proseco/acq.py @@ -1763,7 +1763,7 @@ def p_fid_id_spoiler(self, box_size, fid_id, acqs): f"not a candidate" ) else: - if fids.spoils(fid, self.acq, box_size): + if fids.spoils(fid, self.acq, box_size, fids.dither_acq): p_fid_id_spoiler = 0.0 self._p_fid_id_spoiler[box_size, fid_id] = p_fid_id_spoiler diff --git a/proseco/catalog.py b/proseco/catalog.py index 39bb171e..8ab21114 100644 --- a/proseco/catalog.py +++ b/proseco/catalog.py @@ -11,7 +11,6 @@ from . import __version__ as VERSION from . import characteristics as ACA -from . import characteristics_acq as ACQ from . import get_aca_catalog as _get_aca_catalog_package from .acq import AcqTable, get_acq_catalog, get_maxmag from .core import ( @@ -22,7 +21,7 @@ get_kwargs_from_starcheck_text, ) from .fid import FidTable, get_fid_catalog -from .guide import GuideTable, get_guide_catalog +from .guide import GuideTable, get_guide_candidates, get_guide_catalog, get_t_ccds_bonus from .monitor import BadMonitorError, get_mon_catalog # Colnames and types for final ACA catalog @@ -133,15 +132,29 @@ def _get_aca_catalog(**kwargs): if hasattr(aca.acqs, "dark_date"): aca.dark_date = aca.acqs.dark_date + # Get initial guide candidates + # This set of candidates is selected without considering fid lights, and + # this set of candidates is reused in fid, guide star selection, and fid + # optimization processing (if necessary). The "somewhat expensive" operation + # of reviewing star candidates for guide star criteria and nearby imposters + # is done only once per call to get_aca_catalog(). + # Use .copy() when passing to functions to prevent modification. + initial_guide_cands = get_guide_candidates(stars=aca.acqs.stars, **kwargs).copy() + # Note that aca.acqs.stars is a filtered version of aca.stars and includes # only stars that are in or near ACA FOV. Use this for fids and guides stars. aca.log("Starting get_fid_catalog") - aca.fids = get_fid_catalog(stars=aca.acqs.stars, acqs=aca.acqs, **kwargs) + aca.fids = get_fid_catalog( + stars=aca.acqs.stars, + acqs=aca.acqs, + guide_cands=initial_guide_cands.copy(), + **kwargs, + ) aca.acqs.fids = aca.fids if aca.optimize: aca.log("Starting optimize_acqs_fids") - aca.optimize_acqs_fids() + aca.optimize_acqs_fids(initial_guide_cands=initial_guide_cands.copy(), **kwargs) aca.acqs.fid_set = aca.fids["id"] @@ -153,6 +166,7 @@ def _get_aca_catalog(**kwargs): stars=aca.acqs.stars, fids=aca.fids, mons=aca.mons, + initial_guide_cands=initial_guide_cands.copy(), img_size=img_size_guide, **kwargs, ) @@ -409,15 +423,18 @@ def make_report(self, rootdir="."): self.acqs.make_report(rootdir=rootdir) self.guides.make_report(rootdir=rootdir) - def optimize_acqs_fids(self): + def optimize_acqs_fids(self, initial_guide_cands=None, **kwargs): """ Concurrently optimize acqs and fids in the case where there is not already a good (no spoilers) fid set available. - This updates the acqs and fids tables in place. + This updates the acqs and fids tables in place. The optimization considers + both acquisition star probability (P2) and guide star count to find the + best fid set. - :param acqs: AcqTable object - :param fids: FidTable object + :param initial_guide_cands: Table of initial guide candidates (typically from + get_guide_candidates()). If None, guide candidates are generated internally. + :param **kwargs: additional keyword arguments passed to get_guide_catalog() """ acqs = self.acqs fids = self.fids @@ -433,18 +450,39 @@ def optimize_acqs_fids(self): or len(self.fids.cand_fids) == 0 or len(self.fids.cand_fid_sets) == 0 ): + self.log("No acq-fid optimization required") return + # Only do these imports if it turns out we need optimization. + from chandra_aca.star_probs import guide_count + + from . import characteristics_acq as ACQ + + no_fid_guides = get_guide_catalog( + stars=acqs.stars, + initial_guide_cands=initial_guide_cands, + **kwargs, + ) + # Start with the no-fids optimum catalog and save required info to restore opt_P2 = -acqs.get_log_p_2_or_fewer() + t_ccds_bonus = get_t_ccds_bonus( + no_fid_guides["mag"], + no_fid_guides.t_ccd, + no_fid_guides.dyn_bgd_n_faint, + no_fid_guides.dyn_bgd_dt_ccd, + ) + opt_guide_count = guide_count(no_fid_guides["mag"], t_ccds_bonus) + orig_acq_idxs = acqs["idx"].tolist() orig_acq_halfws = acqs["halfw"].tolist() self.log( - f"Starting opt_P2={opt_P2:.2f}: ids={orig_acq_idxs} halfws={orig_acq_halfws}" + f"Starting conditions: P2={opt_P2:.2f} guide_count={opt_guide_count:.2f}" ) + self.log(f" acq ids={orig_acq_idxs} halfws={orig_acq_halfws}") + self.log(f" guide ids={no_fid_guides['idx'].tolist()}") - # If not at least 2 fids then punt on optimization. cand_fids = fids.cand_fids # Get list of fid_sets that are consistent with candidate fids. These @@ -455,26 +493,57 @@ def optimize_acqs_fids(self): spoiler_score = sum( cand_fids.get_id(fid_id)["spoiler_score"] for fid_id in fid_set ) - rows.append((fid_set, spoiler_score)) + fid_trap_score = sum( + cand_fids.get_id(fid_id)["fid_trap_spoiler"].astype(int) + for fid_id in fid_set + ) + rows.append((fid_set, spoiler_score, fid_trap_score)) # Make a table to keep track of candidate fid_sets along with the # ranking metric P2 and the acq catalog info halfws and star ids. - fid_sets = Table(rows=rows, names=("fid_ids", "spoiler_score")) + fid_sets = Table( + rows=rows, names=("fid_ids", "spoiler_score", "fid_trap_score") + ) fid_sets["P2"] = -99.0 # Marker for unfilled values + fid_sets["guide_count"] = -99.0 fid_sets["acq_halfws"] = None fid_sets["acq_idxs"] = None # Group the table into groups by spoiler score. This preserves the # original fid set ordering within a group. - fid_sets = fid_sets.group_by("spoiler_score") + fid_sets = fid_sets.group_by(("spoiler_score", "fid_trap_score")) # Iterate through each spoiler_score group and then within that iterate # over each fid set. + n_tries = 0 for fid_set_group in fid_sets.groups: spoiler_score = fid_set_group["spoiler_score"][0] - self.log(f"Checking fid sets with spoiler_score={spoiler_score}", level=1) + fid_trap_score = fid_set_group["fid_trap_score"][0] + self.log( + f"Checking fid sets with spoiler_score={spoiler_score} " + f"fid_trap_score={fid_trap_score}", + level=1, + ) for fid_set in fid_set_group: + # get the rows of candidate fids that correspond to the fid_ids + # in the current fid_set. + n_tries += 1 + fid_mask = [fid_id in fid_set["fid_ids"] for fid_id in cand_fids["id"]] + fids_for_set = cand_fids[fid_mask] + local_guides = get_guide_catalog( + stars=acqs.stars, + fids=fids_for_set, + initial_guide_cands=initial_guide_cands, + **kwargs, + ) + local_t_ccds_bonus = get_t_ccds_bonus( + local_guides["mag"], + local_guides.t_ccd, + local_guides.dyn_bgd_n_faint, + local_guides.dyn_bgd_dt_ccd, + ) + local_guide_count = guide_count(local_guides["mag"], local_t_ccds_bonus) # Set the internal acqs fid set. This does validation of the set # and also calls update_p_acq_column(). acqs.fid_set = fid_set["fid_ids"] @@ -487,10 +556,12 @@ def optimize_acqs_fids(self): # sets. Only optimize if P2 decreased by at least 0.001, # remembering P2 is the -log10(prob). fid_set_P2 = -acqs.get_log_p_2_or_fewer() - found_good_set = fid_set_P2 - opt_P2 > -0.001 + found_good_set = (fid_set_P2 - opt_P2 > -0.001) & ( + local_guide_count >= (opt_guide_count - 0.05) + ) if found_good_set: self.log( - f"No change in P2 for fid set {acqs.fid_set}, " + f"No change in P2 or guide_count for fid set {acqs.fid_set}, " f"skipping optimization" ) else: @@ -500,16 +571,19 @@ def optimize_acqs_fids(self): # Store optimization results fid_set["P2"] = -acqs.get_log_p_2_or_fewer() + fid_set["guide_count"] = local_guide_count fid_set["acq_idxs"] = acqs["idx"].tolist() fid_set["acq_halfws"] = acqs["halfw"].tolist() self.log( f"Fid set {fid_set['fid_ids']}: P2={fid_set['P2']:.2f} " + f" guide_count={fid_set['guide_count']:.2f} " f"acq_idxs={fid_set['acq_idxs']} halfws={fid_set['acq_halfws']}", level=2, ) if found_good_set: + self.log("Found acceptable fid set - stopping search") break else: # Put the catalog back to the original no-fid optimum and continue @@ -518,22 +592,43 @@ def optimize_acqs_fids(self): # Get the best fid set / acq catalog configuration so far. Fid sets not # yet considered have P2 = -99. - best_idx = np.argmax(fid_sets["P2"]) - best_P2 = fid_sets["P2"][best_idx] + + # Are there sets with P2 >= 2 and guide_count >= 4? + # P2 == 2 is considered the minimum operationally acceptable level + # and fractional guide_count == 4 is similarly considered the minimum + # acceptable level. In general, we have not put these tolerances into + # proseco, but have checked the catalogs with sparkles. However, here + # in optimization if there are cases that satisfy these minimums then + # we prefer those cases to ones that just maximize P2. + MIN_ACCEPTABLE_P2 = 2.0 + MIN_ACCEPTABLE_GUIDE_COUNT = 4.0 + passable = (fid_sets["P2"] >= MIN_ACCEPTABLE_P2) & ( + fid_sets["guide_count"] >= MIN_ACCEPTABLE_GUIDE_COUNT + ) + if np.any(passable): + fid_sets_passable = fid_sets[passable] + best_idx_passable = np.argmax(fid_sets_passable["P2"]) + best_P2 = fid_sets_passable["P2"][best_idx_passable] + best_idx = np.where(passable)[0][best_idx_passable] + + # If none passable then just get the best P2 + else: + best_idx = np.argmax(fid_sets["P2"]) + best_P2 = fid_sets["P2"][best_idx] # Get the row of the fid / acq stages table to determine the required minimum # P2 given the fid spoiler score. stage = ACQ.fid_acq_stages.loc[spoiler_score] stage_min_P2 = stage["min_P2"](opt_P2) - self.log( f"Best P2={best_P2:.2f} at idx={best_idx} vs. " - "stage_min_P2={stage_min_P2:.2f}", + f"stage_min_P2={stage_min_P2:.2f}", level=1, ) # If we have a winner then use that. - if best_P2 >= stage_min_P2: + if (best_P2 >= stage_min_P2) and np.any(passable): + self.log("Found acceptable acq-fid combination - stopping search") break # Set the acqs table to the best catalog @@ -547,7 +642,8 @@ def optimize_acqs_fids(self): self.log( f"Best acq-fid set: P2={best_P2:.2f} " - f"acq_idxs={best_acq_idxs} halfws={best_acq_halfws} fid_ids={acqs.fid_set}" + f"acq_idxs={best_acq_idxs} halfws={best_acq_halfws} " + f"fid_ids={acqs.fid_set} N opt runs={n_tries}" ) if best_P2 < stage_min_P2: diff --git a/proseco/fid.py b/proseco/fid.py index e2f4c109..cffef3cb 100644 --- a/proseco/fid.py +++ b/proseco/fid.py @@ -15,6 +15,7 @@ from . import characteristics as ACA from . import characteristics_acq as ACQ from . import characteristics_fid as FID +from . import guide from .core import ACACatalogTable, AliasAttribute, MetaAttribute from .jupiter import is_spoiled_by_jupiter @@ -36,6 +37,7 @@ def get_fid_catalog(obsid=0, **kwargs): :param focus_offset: SIM focus offset [steps] (default=0) :param sim_offset: SIM translation offset from nominal [steps] (default=0) :param acqs: AcqTable catalog. Optional but needed for actual fid selection. + :param guide_cands: Table of initial guide candidates (used for star spoilers) :param stars: stars table. Defaults to acqs.stars if available. :param dither_acq: acq dither size (2-element sequence (y, z), arcsec) :param dither_guide: guide dither size (2-element sequence (y, z), arcsec) @@ -72,14 +74,13 @@ def get_fid_catalog(obsid=0, **kwargs): # Add a `slot` column that makes sense fids.set_slot_column() - return fids class FidTable(ACACatalogTable): # Define base set of allowed keyword args to __init__. Subsequent MetaAttribute # or AliasAttribute properties will add to this. - allowed_kwargs = ACACatalogTable.allowed_kwargs | set(["acqs"]) + allowed_kwargs = ACACatalogTable.allowed_kwargs | set(["acqs", "guide_cands"]) # Catalog type when plotting (None | 'FID' | 'ACQ' | 'GUI') catalog_type = "FID" @@ -90,6 +91,7 @@ class FidTable(ACACatalogTable): cand_fids = MetaAttribute(is_kwarg=False) cand_fid_sets = MetaAttribute(is_kwarg=False) + guide_cands = MetaAttribute(default=None) required_attrs = ( "att", @@ -223,18 +225,23 @@ def set_initial_catalog(self): """Set initial fid catalog (fid set) if possible to the first set which is "perfect": - - No field stars spoiler any of the fid lights + - No field stars spoil any of the fid lights - Fid lights are not search spoilers for any of the current acq stars + - The fid set does not include any fid lights in the trap region if there + are guide candidates that would trigger the fid trap effect. If not possible then the table is still zero length and we will need to fall through to the optimization process. """ - # Start by getting the id of every fid that has a zero spoiler score, - # meaning no star spoils the fid as set previously in get_initial_candidates. + # Start by getting the id of every fid that has a zero spoiler score and + # is not affected by fid trap, meaning no star spoils the fid as set + # previously in get_initial_candidates. cand_fids = self.cand_fids unspoiled_fid_ids = set( - fid["id"] for fid in cand_fids if fid["spoiler_score"] == 0 + fid["id"] + for fid in cand_fids + if fid["spoiler_score"] == 0 and not fid["fid_trap_spoiler"] ) # Get list of fid_sets that are consistent with candidate fids. These @@ -248,32 +255,56 @@ def set_initial_catalog(self): if not ok_fid_sets: fid_set = () self.log("No acceptable fid sets (off-CCD or spoiled by field stars)") - - # If no acq stars then just pick the first allowed fid set. - elif self.acqs is None: + # If no stars then just pick the first allowed fid set. + elif self.acqs is None and self.guide_cands is None: fid_set = ok_fid_sets[0] - self.log(f"No acq stars available, using first OK fid set {fid_set}") + self.log(f"No acq/guide stars available, using first OK fid set {fid_set}") else: spoils_any_acq = {} - + spoils_any_guide_cand = {} for fid_set in ok_fid_sets: self.log(f"Checking fid set {fid_set} for acq star spoilers", level=1) for fid_id in fid_set: - if fid_id not in spoils_any_acq: + if (fid_id in spoils_any_acq) or (fid_id in spoils_any_guide_cand): + self.log(f"Fid {fid_id} known already spoiled", level=2) + break + if self.acqs is not None: fid = cand_fids.get_id(fid_id) - spoils_any_acq[fid_id] = any( - self.spoils(fid, acq, acq["halfw"]) for acq in self.acqs + if any( + self.spoils(fid, acq, acq["halfw"], self.dither_acq) + for acq in self.acqs + ): + spoils_any_acq[fid_id] = True + self.log(f"Fid {fid_id} spoils an acq star", level=2) + break + if self.guide_cands is not None: + fid = cand_fids.get_id(fid_id) + if any( + # For guide candidates use 25" box size. + self.spoils(fid, guide_cand, 25, self.dither_guide) + for guide_cand in self.guide_cands + ): + spoils_any_guide_cand[fid_id] = True + self.log(f"Fid {fid_id} spoils a guide candidate", level=2) + break + fid_trap, _ = guide.check_fid_trap( + self.guide_cands, [fid], self.dither_guide ) - if spoils_any_acq[fid_id]: - # Loser, don't bother with the rest. - self.log(f"Fid {fid_id} spoils an acq star", level=2) - break + if np.any(fid_trap): + spoils_any_guide_cand[fid_id] = True + self.log( + f"Fid {fid_id} spoils a guide candidate via trap", + level=2, + ) + break else: # We have a winner, none of the fid_ids in current fid set # will spoil any acquisition star. Break out of loop with # fid_set as the winner. - self.log(f"Fid set {fid_set} is fine for acq stars") + self.log( + f"Fid set {fid_set} is fine for acq stars and guide candidates" + ) break else: # Tried every set and none were acceptable. @@ -285,10 +316,9 @@ def set_initial_catalog(self): idxs = [cand_fids.get_id_idx(fid_id) for fid_id in sorted(fid_set)] for name, col in cand_fids.columns.items(): self[name] = col[idxs] - self.cand_fids = cand_fids - def spoils(self, fid, acq, box_size): + def spoils(self, fid, star, box_size, dither): """ Return true if ``fid`` could be within ``acq`` search box. @@ -300,14 +330,15 @@ def spoils(self, fid, acq, box_size): - Dither amplitude (since OBC adjusts search box for dither) :param fid: fid light (FidTable Row) - :param acq: acq star (AcqTable Row) + :param star: star (AcqTable Row or GuideTable Row) :param box_size: box size (arcsec) + :param dither: dither size (ACABox, arcsec) - :returns: True if ``fid`` could be within ``acq`` search box + :returns: True if ``fid`` could be within ``star`` search box """ - spoiler_margin = FID.spoiler_margin + self.dither_acq + box_size - dy = np.abs(fid["yang"] - acq["yang"]) - dz = np.abs(fid["zang"] - acq["zang"]) + spoiler_margin = FID.spoiler_margin + dither + box_size + dy = np.abs(fid["yang"] - star["yang"]) + dz = np.abs(fid["zang"] - star["zang"]) return dy < spoiler_margin.y and dz < spoiler_margin.z def get_fid_candidates(self): @@ -338,6 +369,7 @@ def get_fid_candidates(self): cand_fids["mag"] = np.full(shape, FID.fid_mag) # 7.000 cand_fids["spoilers"] = np.full(shape, None) # Filled in with Table of spoilers cand_fids["spoiler_score"] = np.full(shape, 0, dtype=np.int64) + cand_fids["fid_trap_spoiler"] = np.full(shape, False, dtype=bool) self.log(f"Initial candidate fid ids are {cand_fids['id'].tolist()}") @@ -371,8 +403,8 @@ def get_fid_candidates(self): cand_fids["idx"] = np.arange(len(cand_fids), dtype=np.int64) - # If stars are available then find stars that are bad for fid. - if self.stars: + # If stars or guide candidates are available then find stars that are bad for fid. + if self.stars or self.guide_cands is not None: for fid in cand_fids: self.set_spoilers_score(fid) @@ -432,21 +464,42 @@ def near_hot_or_bad_pixel(self, fid): return False def set_spoilers_score(self, fid): - """Get stars within FID.spoiler_margin (50 arcsec) + dither. Starcheck uses - 25" but this seems small: 20" (4 pix) positional err + 4 pixel readout - halfw + 2 pixel PSF width of spoiler star. + """ + Get stars within FID.spoiler_margin (50 arcsec) + dither and check for fid trap. + + Starcheck uses 25" but this seems small: 20" (4 pix) positional err + 4 pixel + readout halfw + 2 pixel PSF width of spoiler star. This sets the 'spoilers' column value to a table of spoilers stars (usually empty). - If also sets the 'spoiler_score' to 1 if there is a yellow spoiler - (4 <= star_mag - fid_mag < 5) or 4 for red spoiler (star_mag - fid_mag < 4). + It also sets the 'spoiler_score' based on: + - 1 for yellow spoiler (4 <= star_mag - fid_mag < 5) + - 4 for red spoiler (star_mag - fid_mag < 4) + + Additionally, if the fid is in the fid_trap region and there is a potential guide + candidate that would trigger the fid trap effect, the 'fid_trap_spoiler' flag is set. + + See https://cxc.cfa.harvard.edu/mta/ASPECT/aca_weird_pixels/ for more information on the + "fid trap" issue. + The spoiler score is used later to choose an acceptable set of fids and acq stars. :param fid: fid light (FidTable Row) """ + stars = self.stars[ACQ.spoiler_star_cols] dither = self.dither_guide + # Run guide star fid_trap checks if guide candidates are available + if self.guide_cands is not None: + fid_trap, _ = guide.check_fid_trap(self.guide_cands, [fid], dither) + if np.any(fid_trap): + fid["fid_trap_spoiler"] = True + self.log( + f"Fid {fid['id']} spoiled by fid trap potential from guide candidates", + level=1, + ) + # Potential spoiler by position spoil = ( np.abs(stars["yang"] - fid["yang"]) < FID.spoiler_margin + dither.y diff --git a/proseco/guide.py b/proseco/guide.py index 353fd212..eb2a20e2 100644 --- a/proseco/guide.py +++ b/proseco/guide.py @@ -14,7 +14,11 @@ from proseco.characteristics import MonFunc +if TYPE_CHECKING: + from astropy.table import Table + from . import characteristics as ACA +from . import characteristics_fid as FID from . import characteristics_guide as GUIDE from .core import ( ACACatalogTable, @@ -25,17 +29,94 @@ get_img_size, ) -if TYPE_CHECKING: - from astropy.table import Table - - CCD = ACA.CCD APL = AcaPsfLibrary() STAR_PAIR_DIST_CACHE = {} -def get_guide_catalog(obsid=0, **kwargs): +# Minimum number of "anchor stars" that are always evaluated *without* the bonus +# from dynamic background when dyn_bgd_n_faint > 0. This is mostly to avoid the +# situation where 4 stars are selected and 2 are faint bonus stars. In this case +# there would be only 2 anchor stars that ensure good tracking even without +# dyn bgd. +MIN_DYN_BGD_ANCHOR_STARS = 3 + + +def get_t_ccds_bonus(mags, t_ccd, dyn_bgd_n_faint, dyn_bgd_dt_ccd): + """Return array of t_ccds with dynamic background bonus applied. + + This adds ``dyn_bgd_dt_ccd`` to the effective CCD temperature for the + ``dyn_bgd_n_faint`` faintest stars, ensuring that at least MIN_DYN_BGD_ANCHOR_STARS + are evaluated without the bonus. See: + https://nbviewer.org/urls/cxc.harvard.edu/mta/ASPECT/ipynb/misc/guide-count-dyn-bgd.ipynb + + :param mags: array of star magnitudes + :param t_ccd: single T_ccd value (degC) + :param dyn_bgd_n_faint: number of faintest stars to apply bonus + :param dyn_bgd_dt_ccd: temperature offset to apply for dynamic background bonus (degC) + :returns: array of t_ccds (matching ``mags``) with dynamic background bonus applied + """ + t_ccds = np.full_like(mags, t_ccd) + + # If no bonus stars then just return the input t_ccd broadcast to all stars + if dyn_bgd_n_faint == 0: + return t_ccds + + idxs = np.argsort(mags) + n_faint = min(dyn_bgd_n_faint, len(t_ccds)) + idx_bonus = max(len(t_ccds) - n_faint, MIN_DYN_BGD_ANCHOR_STARS) + for idx in idxs[idx_bonus:]: + t_ccds[idx] += dyn_bgd_dt_ccd + + return t_ccds + + +def get_guide_candidates(obsid=0, **kwargs): + """ + Get initial guide star candidates. + + This function creates a GuideTable and performs the initial filtering + to generate candidate guide stars. The returned Table can be reused + across multiple calls to get_guide_catalog() during optimization by + passing a copy (via .copy()) to prevent modification. + + This explicitly ignores fid lights and monitor windows requests; these + considerations are handled later in get_guide_catalog(). + + :param obsid: obsid (default=0) + :param att: attitude (any object that can initialize Quat) + :param t_ccd: ACA CCD temperature (degC) + :param date: date of acquisition (any DateTime-compatible format) + :param dither: dither size (float or 2-element sequence (dither_y, dither_z), arcsec) + :param stars: astropy.Table of AGASC stars (will be fetched from agasc if None) + :param dark: ACAImage of dark map (fetched based on time and t_ccd if None) + :param **kwargs: additional keyword arguments passed to GuideTable + + :returns: Astropy Table of initial guide candidates + """ + guides = GuideTable() + kwargs.pop("fids", None) + kwargs.pop("mons", None) + guides.set_attrs_from_kwargs(obsid=obsid, **kwargs) + guides.set_stars() + + # Do a first cut of the stars to get a set of reasonable candidates. + cand_guides = guides.get_initial_guide_candidates() + + # For the candidates table, we do want it to have any manual include + # or exclude ids applied, so do that processing here. For example, if + # the user wants to manually include a guide star, that star should be + # explicitly included in the initial candidates list so that it can be + # "considered" by the fid selection and optimization later. + # These operations work in-place on cand_guides. + guides.process_include_ids(cand_guides, guides.stars) + guides.process_exclude_ids(cand_guides) + + return cand_guides + + +def get_guide_catalog(obsid=0, initial_guide_cands=None, **kwargs): """ Get a catalog of guide stars @@ -43,6 +124,9 @@ def get_guide_catalog(obsid=0, **kwargs): ``att``, ``t_ccd``, ``date``, and ``dither`` will be fetched via ``mica.starcheck`` if not explicitly provided here. + If ``initial_guide_cands`` is provided as an argument, it is assumed to be a + Table of initial guide candidates (typically from get_guide_candidates().copy()). + :param obsid: obsid (default=0) :param att: attitude (any object that can initialize Quat) :param t_ccd: ACA CCD temperature (degC) @@ -51,6 +135,8 @@ def get_guide_catalog(obsid=0, **kwargs): :param n_guide: number of guide stars to attempt to get :param fids: selected fids (used for guide star exclusion) :param stars: astropy.Table of AGASC stars (will be fetched from agasc if None) + :param initial_guide_cands: Astropy Table of initial guide candidates + (if None, a new one is created) :param include_ids: list of AGASC IDs of stars to include in guide catalog :param exclude_ids: list of AGASC IDs of stars to exclude from guide catalog :param dark: ACAImage of dark map (fetched based on time and t_ccd if None) @@ -58,24 +144,41 @@ def get_guide_catalog(obsid=0, **kwargs): :returns: GuideTable of acquisition stars """ - STAR_PAIR_DIST_CACHE.clear() guides = GuideTable() guides.set_attrs_from_kwargs(obsid=obsid, **kwargs) guides.set_stars() # Process monitor window requests, converting them into fake stars that - # are added to the include_ids list. + # are added to the include_ids list. This also modifies n_guide. guides.process_monitors_pre() # Do a first cut of the stars to get a set of reasonable candidates - guides.cand_guides = guides.get_initial_guide_candidates() - - # Process guide-from-monitor requests by finding corresponding star in - # cand_guides and adding to the include_ids list. - # guides.process_monitors_pre2() + if initial_guide_cands is None: + guides.cand_guides = guides.get_initial_guide_candidates() + else: + # Use the provided initial guide candidates table + guides.cand_guides = initial_guide_cands + + # Refilter candidates for the current fid set if needed + if guides.fids is not None and len(guides.fids) > 0: + guides.cand_guides = guides.filter_candidates_for_fids(guides.cand_guides) + + # Refilter candidates for monitor keep-out zones + guides.cand_guides = guides.filter_candidates_for_monitors(guides.cand_guides) + + # Put back any include/excludes + # Note explicitly that the monitor star processing in process_monitors_pre() + # can add stars to include_ids. In the standard flow in get_aca_catalog, + # that means that the initial guide candidates selection did not have those + # include_ids applied, so we need to do it here. The exclude_ids is probably + # not necessary to re-apply, but do it for symmetry. + # These operations work in-place on guides.cand_guides. + guides.process_include_ids(guides.cand_guides, guides.stars) + guides.process_exclude_ids(guides.cand_guides) # Run through search stages to select stars + STAR_PAIR_DIST_CACHE.clear() selected = guides.run_search_stages() # Transfer to table (which at this point is an empty table) @@ -186,6 +289,7 @@ class GuideTable(ACACatalogTable): cand_guides = MetaAttribute(is_kwarg=False) reject_info = MetaAttribute(default=[], is_kwarg=False) img_size = ImgSizeMetaAttribute() + fids = MetaAttribute(default=None) def reject(self, reject): """ @@ -206,29 +310,17 @@ def process_monitors_pre(self): monitors = self.mons.monitors is_mon = np.isin(monitors["function"], [MonFunc.MON_FIXED, MonFunc.MON_TRACK]) - if np.any(is_mon): - # For MON fixed and tracked, the dark cal map gets hacked to make a - # local blob of hot pixels. Make a copy to avoid modifying the global - # dark map which is shared by reference between acqs, guides and fids. - self.dark = self.dark.copy() - for monitor in monitors[is_mon]: - # Make a square blob of saturated pixels that will keep away any - # star selection. - is_track = monitor["function"] == MonFunc.MON_TRACK - dr = int(4 + np.ceil(self.dither.row if is_track else 0)) - dc = int(4 + np.ceil(self.dither.col if is_track else 0)) - row, col = int(monitor["row"]) + 512, int(monitor["col"]) + 512 - self.dark[row - dr : row + dr, col - dc : col + dc] = ( - ACA.bad_pixel_dark_current - ) + monitors = self.mons.monitors + is_mon = np.isin(monitors["function"], [MonFunc.MON_FIXED, MonFunc.MON_TRACK]) - # Reduce n_guide for each MON. On input the n_guide arg is the - # number of GUI + MON, but for guide selection we need to make the - # MON slots unavailable. - self.n_guide -= 1 - if self.n_guide < 2: - raise ValueError("too many MON requests leaving < 2 guide stars") + # Reduce n_guide for each MON. On input the n_guide arg is the + # number of GUI + MON, but for guide selection we need to make the + # MON slots unavailable. + n_mon = np.count_nonzero(is_mon) + self.n_guide -= n_mon + if self.n_guide < 2: + raise ValueError("too many MON requests leaving < 2 guide stars") is_guide = monitors["function"] == MonFunc.GUIDE for monitor in monitors[is_guide]: @@ -237,8 +329,6 @@ def process_monitors_pre(self): def process_monitors_post(self): """Post-processing of monitor windows. - - Clean up fake stars from stars and cand_guides - - Restore the dark current map - Set type, sz, res, dim columns for relevant entries to reflect status as monitor or guide-from-monitor in the guides catalog. """ @@ -246,11 +336,6 @@ def process_monitors_post(self): if self.mons is None or self.mons.monitors is None: return - # Mon window processing munges the dark cal to impose a keep-out zone. - # Put the dark current back to the standard one if possible - if self.acqs is not None: - self.dark = self.acqs.dark - # Find the guide stars that are actually from a MON request (converted # to guide) and set the size and type. for monitor in self.mons.monitors: @@ -259,6 +344,98 @@ def process_monitors_post(self): row["sz"] = "8x8" row["type"] = "GFM" # Guide From Monitor, changed in merge_catalog + def filter_candidates_for_fids(self, cand_guides): + """ + Filter candidate guide stars with regard to the selected fid lights. + + This method removes guide star candidates that would spoil or be spoiled by + fid lights + 1. Fid trap effect: candidates that trigger readout artifacts from fid lights + 2. Direct spoiling: candidates too close to fid lights (within spoiler margin) + + :param cand_guides: Table of candidate guide stars + :returns: filtered Table of candidate guide stars with spoiled stars removed + """ + if self.fids is None or len(self.fids) == 0: + return cand_guides + + # Collect all rows to exclude, then filter once at the end + all_spoiled = np.zeros(len(cand_guides), dtype=bool) + + # Handle the fid trap + fid_trap_spoilers, fid_rej = check_fid_trap( + cand_guides, fids=self.fids, dither=self.dither + ) + for rej in fid_rej: + rej["stage"] = 0 + self.reject(rej) + all_spoiled |= fid_trap_spoilers + + # Exclude guide stars directly spoiled by fid lights + spoiler_margin = FID.spoiler_margin + self.dither + 25 + for fid in self.fids: + dy = np.abs(fid["yang"] - cand_guides["yang"]) + dz = np.abs(fid["zang"] - cand_guides["zang"]) + spoiled = (dy < spoiler_margin.y) & (dz < spoiler_margin.z) + for star_id in cand_guides["id"][spoiled]: + self.reject( + { + "stage": 0, + "type": "spoiled by fid", + "id": star_id, + "text": f"Cand {star_id} rejected. Spoiled by fid {fid['id']}", + } + ) + all_spoiled |= spoiled + + # Filter once with combined mask + return cand_guides[~all_spoiled] + + def filter_candidates_for_monitors(self, cand_guides): + """Filter candidate guide stars that fall in monitor window keep-out zones. + + This directly excludes candidates near MON_FIXED and MON_TRACK windows, + replacing the old approach of modifying the dark map. + """ + if self.mons is None or self.mons.monitors is None: + return cand_guides + + monitors = self.mons.monitors + is_mon = np.isin(monitors["function"], [MonFunc.MON_FIXED, MonFunc.MON_TRACK]) + + # Collect all rows to exclude, then filter once at the end + all_in_keepout = np.zeros(len(cand_guides), dtype=bool) + + for monitor in monitors[is_mon]: + is_track = monitor["function"] == MonFunc.MON_TRACK + # Keep-out zone in pixels + dr = 4 + np.ceil(self.dither.row if is_track else 0) + dc = 4 + np.ceil(self.dither.col if is_track else 0) + + # Find candidates in the keep-out zone plus a pad + # The 7.5 pixel pad was selected to at least cover the previous dark-map + # modification approach + monitor_pad = 7.5 + drow = np.abs(cand_guides["row"] - monitor["row"]) + dcol = np.abs(cand_guides["col"] - monitor["col"]) + in_keepout = (drow < dr + monitor_pad) & (dcol < dc + monitor_pad) + + for star_id in cand_guides["id"][in_keepout]: + self.reject( + { + "stage": 0, + "type": "monitor keepout", + "id": star_id, + "monitor_id": monitor["id"], + "text": f"Cand {star_id} rejected. " + f"In monitor {monitor['id']} keep-out zone", + } + ) + all_in_keepout |= in_keepout + + # Filter once with combined mask + return cand_guides[~all_in_keepout] + def get_img_size(self, n_fids=None): """Get guide image readout size from ``img_size`` and ``n_fids``. @@ -614,7 +791,7 @@ def search_stage(self, stage): # Adopt the SAUSAGE convention of a bit array for errors # Not all items will be checked for each star (allow short circuit) - scol = "stat_{}".format(stage["Stage"]) + scol = f"stat_{stage['Stage']}" cand_guides[scol] = 0 n_sigma = stage["SigErrMultiplier"] @@ -760,12 +937,72 @@ def process_include_ids(self, cand_guides, stars): :param stars: stars table """ + if len(self.include_ids) == 0: + return + + # Check to see if the include_ids are already in cand_guides + missing_ids = [ + star_id for star_id in self.include_ids if star_id not in cand_guides["id"] + ] + if len(missing_ids) == 0: + self.log(f"All include_ids already in candidate list: {self.include_ids}") + return + row_max = CCD["row_max"] - CCD["row_pad"] - CCD["window_pad"] col_max = CCD["col_max"] - CCD["col_pad"] - CCD["window_pad"] ok = (np.abs(stars["row"]) < row_max) & (np.abs(stars["col"]) < col_max) - super().process_include_ids(cand_guides, stars[ok]) + # Let's get a smaller set of the possible stars in the field and + # populate the imposter mags for them and initialize other columns + isin_mask = np.isin(stars["id"], missing_ids, assume_unique=True) + new_stars = stars[isin_mask & ok].copy() + + # Now compute imposter mags for these stars + if len(new_stars) > 0: + imp_mag, imp_row, imp_col = get_imposter_mags( + new_stars, self.dark, self.dither + ) + new_stars["imp_mag"] = imp_mag + new_stars["imp_r"] = imp_row + new_stars["imp_c"] = imp_col + + # And add the other standard columns + # I'd like to add stage set to -1 and stat_1 to stat_6 set to 0 + new_stars["stage"] = -1 + for stage_num in range(1, len(GUIDE.stages) + 1): + scol = f"stat_{stage_num}" + new_stars[scol] = 0 + + # Reorder columns to match cand_guides + new_stars = new_stars[cand_guides.colnames] + + self.log(f"In process_include_ids trying to include {self.include_ids}") + super().process_include_ids(cand_guides, new_stars) + + def process_exclude_ids(self, cand_guides): + """ + Remove any stars from cand_guides that are in exclude_ids. + + This works in-place on cand_guides. + + :param cand_guides: candidate guide stars table + """ + # Deal with exclude_ids by removing rows from the candidate list + # This works in-place on cand_guides + for star_id in self.exclude_ids: + if star_id in cand_guides["id"]: + self.reject( + { + "stage": 0, + "type": "exclude_id", + "id": star_id, + "text": f"Cand {star_id} rejected. In exclude_ids", + } + ) + # Remove the row in-place + idx = np.where(cand_guides["id"] == star_id)[0] + cand_guides.remove_rows(idx) def get_candidates_mask(self, stars): """Get base filter for acceptable candidates. @@ -817,22 +1054,16 @@ def get_initial_guide_candidates(self): `has_spoiler_in_box` method and is the first of many spoiler checks. This spoiler check is not stage dependent. - 9. Filters the candidates to remova any that are spoiled by the fid trap (using - `check_fid_trap` method). - - 10. Puts any force include candidates from the include_ids parameter back in the candidate - list if they were filtered out by an earlier filter in this routine. + 9. If Jupiter is present, filters the candidates to remove any that are within + 15 columns of Jupiter. - 11. Filters/removes any candidates that are force excluded (in exclude_ids). - - 12. Uses the local dark current around each candidate to calculate an "imposter mag" + 10. Uses the local dark current around each candidate to calculate an "imposter mag" describing the brightest 2x2 in the region the star would dither over. This is saved to the candidate star table. """ stars = self.stars - dark = self.dark # Mark stars that are off chip offchip = (np.abs(stars["row"]) > CCD["row_max"]) | ( @@ -864,6 +1095,15 @@ def get_initial_guide_candidates(self): f"{len(cand_guides)} candidate guide stars" ) + # Add extra columns used by guide selection to the candidate table + cand_guides["stage"] = -1 + for stage_num in range(1, len(GUIDE.stages) + 1): + scol = f"stat_{stage_num}" + cand_guides[scol] = 0 + cand_guides["imp_mag"] = np.nan + cand_guides["imp_r"] = np.nan + cand_guides["imp_c"] = np.nan + bp, bp_rej = spoiled_by_bad_pixel(cand_guides, self.dither) for rej in bp_rej: rej["stage"] = 0 @@ -908,14 +1148,6 @@ def get_initial_guide_candidates(self): f"{len(cand_guides)} candidate guide stars" ) - fid_trap_spoilers, fid_rej = check_fid_trap( - cand_guides, fids=self.fids, dither=self.dither - ) - for rej in fid_rej: - rej["stage"] = 0 - self.reject(rej) - cand_guides = cand_guides[~fid_trap_spoilers] - if len(self.jupiter) > 0: from proseco.jupiter import check_spoiled_by_jupiter @@ -929,24 +1161,10 @@ def get_initial_guide_candidates(self): self.reject(rej) cand_guides = cand_guides[~spoiled_by_jupiter] - # Deal with include_ids by putting them back in candidate table if necessary - self.process_include_ids(cand_guides, stars) - - # Deal with exclude_ids by cutting from the candidate list - for star_id in self.exclude_ids: - if star_id in cand_guides["id"]: - self.reject( - { - "stage": 0, - "type": "exclude_id", - "id": star_id, - "text": f"Cand {star_id} rejected. In exclude_ids", - } - ) - cand_guides = cand_guides[cand_guides["id"] != star_id] - # Get the brightest 2x2 in the dark map for each candidate and save value and location - imp_mag, imp_row, imp_col = get_imposter_mags(cand_guides, dark, self.dither) + imp_mag, imp_row, imp_col = get_imposter_mags( + cand_guides, self.dark, self.dither + ) cand_guides["imp_mag"] = imp_mag cand_guides["imp_r"] = imp_row cand_guides["imp_c"] = imp_col @@ -961,7 +1179,8 @@ def in_bad_star_list(self, cand_guides): :param cand_guides: Table of candidate stars :returns: boolean mask where True means star is in bad star list """ - bad = np.isin(cand_guides["id"], list(ACA.bad_star_set)) + bad_star_ids = list(ACA.bad_star_set) + bad = np.isin(cand_guides["id"], bad_star_ids) # Set any matching bad stars as bad for plotting for bad_id in cand_guides["id"][bad]: @@ -975,6 +1194,8 @@ def check_fid_trap(cand_stars, fids, dither): """ Search for guide stars that would cause the fid trap issue and mark as spoilers. + The "fid trap" issue is documented at https://cxc.cfa.harvard.edu/mta/ASPECT/aca_weird_pixels/ . + :param cand_stars: candidate star Table :param fids: fid Table :param dither: dither ACABox diff --git a/proseco/tests/conftest.py b/proseco/tests/conftest.py index 2252606d..65ce8a58 100644 --- a/proseco/tests/conftest.py +++ b/proseco/tests/conftest.py @@ -49,3 +49,8 @@ def proseco_agasc_1p7(monkeypatch): def miniagasc_1p7(monkeypatch): agasc_file = get_agasc_filename("miniagasc_*", version="1p7") monkeypatch.setenv("AGASC_HDF5_FILE", agasc_file) + + +@pytest.fixture(autouse=True) +def use_kadi_flight_scenario(monkeypatch): + monkeypatch.setenv("KADI_SCENARIO", "flight") diff --git a/proseco/tests/test_catalog.py b/proseco/tests/test_catalog.py index 1531c60d..554af574 100644 --- a/proseco/tests/test_catalog.py +++ b/proseco/tests/test_catalog.py @@ -38,11 +38,7 @@ def test_allowed_kwargs(): } new_kwargs = FidTable.allowed_kwargs - ACACatalogTable.allowed_kwargs - assert new_kwargs == { - "acqs", - "include_ids", - "exclude_ids", - } + assert new_kwargs == {"acqs", "include_ids", "exclude_ids", "guide_cands"} @pytest.mark.skipif(not HAS_SC_ARCHIVE, reason="Test requires starcheck archive") @@ -701,7 +697,12 @@ def test_fid_trap_effect(): agasc_ids = [367148872, 367139768, 367144424, 367674552, 367657896] att = [10.659376, 40.980028, 181.012903] stars = StarsTable.from_agasc_ids(att, agasc_ids) - cat = get_aca_catalog(obsid=1576, stars=stars, raise_exc=True) + # Note that the fid lights need to be included, as guide/fid light optimization + # would otherwise unselect fid 6 and thus not apply the trap effect exclusion + # to star 367674552. + cat = get_aca_catalog( + obsid=1576, include_ids_fid=[1, 5, 6], stars=stars, raise_exc=True + ) assert 367674552 not in cat.guides["id"] # Obsid 2365 @@ -713,7 +714,11 @@ def test_fid_trap_effect(): stars = StarsTable.from_agasc_ids(att, agasc_ids) # Specify att kwarg explicitly to override what is found in mica.starcheck - cat = get_aca_catalog(obsid=2365, stars=stars, raise_exc=True, att=att) + # And include fid lights to ensure optimization applies the trap effect exclusion, + # otherwise fid 6 would be unselected and star 1184897704 would not be excluded. + cat = get_aca_catalog( + obsid=2365, stars=stars, raise_exc=True, att=att, include_ids_fid=[1, 5, 6] + ) assert 1184897704 not in cat.guides["id"] diff --git a/proseco/tests/test_fid.py b/proseco/tests/test_fid.py index 3774fbf0..0be05e5a 100644 --- a/proseco/tests/test_fid.py +++ b/proseco/tests/test_fid.py @@ -127,18 +127,31 @@ def test_fid_catalog_t_ccd(): def test_get_initial_catalog(FIDS): """Test basic catalog with no stars in field using standard 2-4-5 config.""" exp = [ - "", - " id yang zang row col mag spoiler_score idx slot", - "int64 float64 float64 float64 float64 float64 int64 int64 int64", - "----- -------- -------- ------- ------- ------- ------------- ----- -----", - " 1 918.09 -1741.51 -179.15 -344.83 7.00 0 0 -99", - " 2 -777.70 -1745.65 161.69 -346.09 7.00 0 1 0", - " 3 35.52 -1874.72 -1.77 -371.73 7.00 0 2 -99", - " 4 2135.73 163.01 -423.60 38.40 7.00 0 3 1", - " 5 -1830.77 156.55 373.88 35.74 7.00 0 4 2", - " 6 384.10 800.13 -70.59 165.37 7.00 0 5 -99", + " id yang zang row col mag spoiler_score fid_trap_spoiler idx slot", + "int64 float64 float64 float64 float64 float64 int64 bool int64 int64", + "----- -------- -------- ------- ------- ------- ------------- ---------------- ----- -----", + " 1 918.09 -1741.51 -179.15 -344.83 7.00 0 False 0 -99", + " 2 -777.70 -1745.65 161.69 -346.09 7.00 0 False 1 0", + " 3 35.52 -1874.72 -1.77 -371.73 7.00 0 False 2 -99", + " 4 2135.73 163.01 -423.60 38.40 7.00 0 False 3 1", + " 5 -1830.77 156.55 373.88 35.74 7.00 0 False 4 2", + " 6 384.10 800.13 -70.59 165.37 7.00 0 False 5 -99", ] - assert repr(FIDS.cand_fids).splitlines() == exp + cols = [ + "id", + "yang", + "zang", + "row", + "col", + "mag", + "spoiler_score", + "fid_trap_spoiler", + "idx", + "slot", + ] + assert ( + FIDS.cand_fids[cols].pformat(show_dtype=True, max_width=-1, max_lines=-1) == exp + ) assert np.all(FIDS["id"] == [2, 4, 5]) # Make catalogs with some fake stars (at exactly fid positions) that spoil @@ -152,18 +165,32 @@ def test_get_initial_catalog(FIDS): # Spoil fids 1, 2 fids2 = get_fid_catalog(stars=stars[:2], **STD_INFO) exp = [ - "", - " id yang zang row col mag spoiler_score idx slot", - "int64 float64 float64 float64 float64 float64 int64 int64 int64", - "----- -------- -------- ------- ------- ------- ------------- ----- -----", - " 1 918.09 -1741.51 -179.15 -344.83 7.00 4 0 -99", - " 2 -777.70 -1745.65 161.69 -346.09 7.00 4 1 -99", - " 3 35.52 -1874.72 -1.77 -371.73 7.00 0 2 0", - " 4 2135.73 163.01 -423.60 38.40 7.00 0 3 1", - " 5 -1830.77 156.55 373.88 35.74 7.00 0 4 2", - " 6 384.10 800.13 -70.59 165.37 7.00 0 5 -99", + " id yang zang row col mag spoiler_score fid_trap_spoiler idx slot", + "int64 float64 float64 float64 float64 float64 int64 bool int64 int64", + "----- -------- -------- ------- ------- ------- ------------- ---------------- ----- -----", + " 1 918.09 -1741.51 -179.15 -344.83 7.00 4 False 0 -99", + " 2 -777.70 -1745.65 161.69 -346.09 7.00 4 False 1 -99", + " 3 35.52 -1874.72 -1.77 -371.73 7.00 0 False 2 0", + " 4 2135.73 163.01 -423.60 38.40 7.00 0 False 3 1", + " 5 -1830.77 156.55 373.88 35.74 7.00 0 False 4 2", + " 6 384.10 800.13 -70.59 165.37 7.00 0 False 5 -99", + ] + cols = [ + "id", + "yang", + "zang", + "row", + "col", + "mag", + "spoiler_score", + "fid_trap_spoiler", + "idx", + "slot", ] - assert repr(fids2.cand_fids).splitlines() == exp + assert ( + fids2.cand_fids[cols].pformat(show_dtype=True, max_width=-1, max_lines=-1) + == exp + ) assert np.all(fids2["id"] == [3, 4, 5]) # Spoil fids 1, 2, 3 @@ -219,19 +246,33 @@ def test_fid_spoiling_acq(dither_z, FIDS): acqs["halfw"] = 100 fids5 = get_fid_catalog(acqs=acqs, **std_info) exp = [ - "", - " id yang zang row col mag spoiler_score idx slot", - "int64 float64 float64 float64 float64 float64 int64 int64 int64", - "----- -------- -------- ------- ------- ------- ------------- ----- -----", - " 1 918.09 -1741.51 -179.15 -344.83 7.00 0 0 0", - " 2 -777.70 -1745.65 161.69 -346.09 7.00 0 1 -99", - " 3 35.52 -1874.72 -1.77 -371.73 7.00 0 2 -99", - " 4 2135.73 163.01 -423.60 38.40 7.00 0 3 -99", - " 5 -1830.77 156.55 373.88 35.74 7.00 0 4 1", - " 6 384.10 800.13 -70.59 165.37 7.00 0 5 2", + " id yang zang row col mag spoiler_score fid_trap_spoiler idx slot", + "int64 float64 float64 float64 float64 float64 int64 bool int64 int64", + "----- -------- -------- ------- ------- ------- ------------- ---------------- ----- -----", + " 1 918.09 -1741.51 -179.15 -344.83 7.00 0 False 0 0", + " 2 -777.70 -1745.65 161.69 -346.09 7.00 0 False 1 -99", + " 3 35.52 -1874.72 -1.77 -371.73 7.00 0 False 2 -99", + " 4 2135.73 163.01 -423.60 38.40 7.00 0 False 3 -99", + " 5 -1830.77 156.55 373.88 35.74 7.00 0 False 4 1", + " 6 384.10 800.13 -70.59 165.37 7.00 0 False 5 2", + ] + cols = [ + "id", + "yang", + "zang", + "row", + "col", + "mag", + "spoiler_score", + "fid_trap_spoiler", + "idx", + "slot", ] - assert repr(fids5.cand_fids).splitlines() == exp + assert ( + fids5.cand_fids[cols].pformat(show_dtype=True, max_width=-1, max_lines=-1) + == exp + ) def test_fid_mult_spoilers(disable_fid_offsets, proseco_agasc_1p7): diff --git a/proseco/tests/test_guide.py b/proseco/tests/test_guide.py index bdb63f9e..c961c3ca 100644 --- a/proseco/tests/test_guide.py +++ b/proseco/tests/test_guide.py @@ -11,26 +11,78 @@ from chandra_aca.transform import count_rate_to_mag, mag_to_count_rate from Quaternion import Quat -from ..characteristics import CCD +from ..characteristics import CCD, MonCoord, MonFunc from ..characteristics_guide import mag_spoiler from ..core import StarsTable +from ..fid import get_fid_catalog from ..guide import ( GUIDE, + MIN_DYN_BGD_ANCHOR_STARS, GuideTable, check_column_spoilers, check_mag_spoilers, check_spoil_contrib, get_ax_range, + get_guide_candidates, get_guide_catalog, get_pixmag_for_offset, + get_t_ccds_bonus, run_cluster_checks, ) +from ..monitor import get_mon_catalog from ..report_guide import make_report from .test_common import DARK40, OBS_INFO, STD_INFO, mod_std_info HAS_SC_ARCHIVE = Path(mica.starcheck.starcheck.FILES["data_root"]).exists() +def test_get_t_ccds_bonus_1(): + mags = [1, 10, 2, 11, 3, 4] + t_ccd = 10 + + # Temps corresponding to two faintest stars are equal to t_ccd + dt_ccd + t_ccds = get_t_ccds_bonus(mags, t_ccd, dyn_bgd_n_faint=2, dyn_bgd_dt_ccd=-1) + assert np.all(t_ccds == [10, 9, 10, 9, 10, 10]) + + # Temps corresponding to three faintest stars are equal to t_ccd + dt_ccd + t_ccds = get_t_ccds_bonus(mags, t_ccd, dyn_bgd_n_faint=3, dyn_bgd_dt_ccd=-1) + assert np.all(t_ccds == [10, 9, 10, 9, 10, 9]) + + # Temps corresponding to just the three faintest stars are equal to t_ccd + + # dt_ccd because of the minimum number of anchor stars = 3. + t_ccds = get_t_ccds_bonus(mags, t_ccd, dyn_bgd_n_faint=4, dyn_bgd_dt_ccd=-1) + assert np.all(t_ccds == [10, 9, 10, 9, 10, 9]) + + t_ccds = get_t_ccds_bonus(mags, t_ccd, dyn_bgd_n_faint=0, dyn_bgd_dt_ccd=-1) + assert np.all(t_ccds == [10, 10, 10, 10, 10, 10]) + + +def test_get_t_ccds_bonus_min_anchor(): + mags = [1, 10, 2] + t_ccd = 10 + t_ccds = get_t_ccds_bonus(mags, t_ccd, dyn_bgd_n_faint=2, dyn_bgd_dt_ccd=-1) + # Assert that there are at least MIN_DYN_BGD_ANCHOR_STARS without bonus + assert np.count_nonzero(t_ccds == t_ccd) >= MIN_DYN_BGD_ANCHOR_STARS + + t_ccds = get_t_ccds_bonus(mags, t_ccd, dyn_bgd_n_faint=4, dyn_bgd_dt_ccd=-1) + assert np.count_nonzero(t_ccds == t_ccd) >= MIN_DYN_BGD_ANCHOR_STARS + + +def test_get_t_ccds_bonus_small_catalog(): + mags = [1] + t_ccd = 10 + t_ccds = get_t_ccds_bonus(mags, t_ccd, dyn_bgd_n_faint=2, dyn_bgd_dt_ccd=-1) + assert np.all(t_ccds == [10]) + + +def test_get_t_ccds_bonus_no_bonus_stars(): + """When dyn_bgd_n_faint=0, all t_ccds should equal input t_ccd.""" + mags = np.array([8.0, 9.0, 10.0, 10.5]) + t_ccd = -10.0 + t_ccds = get_t_ccds_bonus(mags, t_ccd, dyn_bgd_n_faint=0, dyn_bgd_dt_ccd=5.0) + assert np.allclose(t_ccds, -10.0) + + def test_select(proseco_agasc_1p7): """ Regression test that 5 expected agasc ids are selected at an arbitrary ra/dec/roll . @@ -844,6 +896,42 @@ def test_edge_star(dither): assert len(guides) == 4 +monitor_cases = [ + {"yang": 300, "zang": 500, "selected": False}, + {"yang": 233, "zang": 500, "selected": False}, + {"yang": 232, "zang": 500, "selected": True}, + {"yang": 367, "zang": 500, "selected": False}, + {"yang": 368, "zang": 500, "selected": True}, + {"yang": 300, "zang": 432, "selected": True}, + {"yang": 300, "zang": 433, "selected": False}, + {"yang": 300, "zang": 568, "selected": True}, + {"yang": 300, "zang": 567, "selected": False}, +] + + +@pytest.mark.parametrize("case", monitor_cases) +def test_monitor_star_keepout(case): + mon_y = 300 + mon_z = 500 + star_id = 999999 + stars = StarsTable.empty() + stars.add_fake_constellation(n_stars=5, mag=8) + stars.add_fake_star(id=star_id, yang=case["yang"], zang=case["zang"], mag=6.0) + kwargs = mod_std_info( + att=stars.att, + n_fid=0, + n_guide=5, + monitors=[[mon_y, mon_z, MonCoord.YAGZAG, 6.0, MonFunc.MON_TRACK]], + ) + mons = get_mon_catalog(**kwargs) + guide = get_guide_catalog( + **mod_std_info(att=stars.att, n_fid=0, n_guide=5), + stars=stars, + mons=mons, + ) + assert (star_id in guide["id"]) == case["selected"] + + def test_get_ax_range(): """ Confirm that the ranges from get_ax_range are reasonable for a variety of @@ -865,6 +953,239 @@ def test_get_ax_range(): assert n - extent - 2 < minus +def test_filter_candidates_for_fids(): + stars = StarsTable.empty() + stars.add_fake_constellation(n_stars=4) + + args = mod_std_info(detector="ACIS-I", include_ids_fid=[6]) + fids = get_fid_catalog(stars=stars, guide_cands=stars, **args) + + # For the nominal ACIS-I FID 6 position here of -69.87 345.27 + # a star will cause the fid trap if it has about the same distance + # to the register as the fid distance to the trap. + # Either positive or negative row works so this test is parameterized. + stars.add_fake_star(id=1001, mag=7.0, row=200, col=400) + + initial_guide_cands = get_guide_candidates(stars=stars, **args) + guide = get_guide_catalog( + stars=stars, initial_guide_cands=initial_guide_cands.copy(), **args + ) + # Stub in fids to guide for filtering test + guide.fids = fids + # Confirm the trap star is in the candidates before filtering. + assert 1001 in guide.cand_guides["id"] + # Run the filtering. + filtered = guide.filter_candidates_for_fids(guide.cand_guides) + assert 1001 not in filtered["id"] + + +def test_filter_candidates_for_monitors(): + stars = StarsTable.empty() + stars.add_fake_constellation(n_stars=4) + args = mod_std_info(detector="ACIS-I") + + monitors = [[50, -50, MonCoord.ROWCOL, 7.5, MonFunc.MON_TRACK]] + mons = get_mon_catalog(monitors=monitors, **args) + + stars.add_fake_star(id=1001, mag=7.5, row=50, col=-50) + + initial_guide_cands = get_guide_candidates(stars=stars, **args) + guide_with_mons = get_guide_catalog( + stars=stars, + initial_guide_cands=initial_guide_cands.copy(), + mons=mons, + **args, + ) + # The monitor window with MON_TRACK should exclude the star + assert 1001 not in guide_with_mons.cand_guides["id"] + assert 1001 not in guide_with_mons["id"] + # But it should still be in the initial candidates as the later filtering + # does not modify that table. + initial_guide_cands_table1 = initial_guide_cands.copy() + assert 1001 in initial_guide_cands_table1["id"] + + guide_without_mons = get_guide_catalog( + stars=stars, initial_guide_cands=initial_guide_cands.copy(), **args + ) + # Without monitors, the star should be present in the initial candidates + # and in the selected table. + assert 1001 in guide_without_mons.cand_guides["id"] + assert 1001 in guide_without_mons["id"] + initial_guide_cands_table2 = initial_guide_cands.copy() + assert 1001 in initial_guide_cands_table2["id"] + + +def test_initial_guide_candidate_reuse(): + stars = StarsTable.empty() + stars.add_fake_constellation(n_stars=5, mag=[8, 9, 9.5, 10, 10.5]) + kwargs = mod_std_info( + att=stars.att, + n_guide=5, + ) + # Get initial guide candidates + initial_guide_cands = get_guide_candidates(stars=stars, **kwargs) + + # Add 2 bright stars + stars.add_fake_star(id=1001, mag=6.0, yang=0, zang=0) + stars.add_fake_star(id=1002, mag=6.5, yang=100, zang=100) + + # Select guide stars using initial candidates + guides_with_reuse = get_guide_catalog( + stars=stars, initial_guide_cands=initial_guide_cands.copy(), **kwargs + ) + + # Confirm that the initial candidates were reused and the new bright stars + # were not selected even though they are brighter and in stars + assert 1001 not in guides_with_reuse["id"] + assert 1002 not in guides_with_reuse["id"] + + guides_without_reuse = get_guide_catalog(stars=stars, **kwargs) + # Confirm that the new bright stars were selected when not reusing initial candidates + assert 1001 in guides_without_reuse["id"] + assert 1002 in guides_without_reuse["id"] + + +def test_initial_guide_cands_do_not_change(): + """Test that .copy() prevents modifications to initial guide candidates.""" + + stars = StarsTable.empty() + stars.add_fake_constellation(n_stars=6, mag=[8, 9, 9.5, 10, 10.5, 11]) + kwargs = mod_std_info(att=stars.att, n_guide=5) + + # Get initial guide candidates + initial_guide_cands = get_guide_candidates(stars=stars, **kwargs) + + # Verify it returns an Astropy Table + assert isinstance(initial_guide_cands, Table) + + # Save original state for comparison + original_ids = set(initial_guide_cands["id"]) + original_len = len(initial_guide_cands) + original_colnames = list(initial_guide_cands.colnames) + + # Use it in multiple catalog selections with .copy() + for _ in range(3): + guides = get_guide_catalog( + stars=stars, initial_guide_cands=initial_guide_cands.copy(), **kwargs + ) + # Verify the guides table is mutable (as expected) + assert isinstance(guides, Table) + + # Verify initial_guide_cands remains unchanged after using .copy() + assert set(initial_guide_cands["id"]) == original_ids + assert len(initial_guide_cands) == original_len + assert list(initial_guide_cands.colnames) == original_colnames + + # Now test the in-place modification behavior directly by calling process_exclude_ids + # This demonstrates why .copy() is necessary when reusing candidate tables + from proseco.guide import GuideTable + + initial_guide_cands_nocopy = get_guide_candidates(stars=stars, **kwargs) + original_len_nocopy = len(initial_guide_cands_nocopy) + original_ids_nocopy = set(initial_guide_cands_nocopy["id"]) + + # Get one of the candidate IDs to exclude + exclude_id = initial_guide_cands_nocopy["id"][0] + + # Create a GuideTable instance and call process_exclude_ids directly + guides_test = GuideTable() + guides_test.set_attrs_from_kwargs( + stars=stars, exclude_ids_guide=[exclude_id], **kwargs + ) + + # Call process_exclude_ids which modifies the table in-place via remove_rows() + guides_test.process_exclude_ids(initial_guide_cands_nocopy) + + # Verify that initial_guide_cands_nocopy WAS modified when not using .copy() + # The excluded star should be removed from the original table + assert len(initial_guide_cands_nocopy) == original_len_nocopy - 1 + assert exclude_id not in initial_guide_cands_nocopy["id"] + assert set(initial_guide_cands_nocopy["id"]) != original_ids_nocopy + + # Verify the first test's table is STILL unchanged (wasn't affected by nocopy test) + assert set(initial_guide_cands["id"]) == original_ids + assert len(initial_guide_cands) == original_len + + +def test_initial_guide_candidate_reuse_with_fids(): + stars = StarsTable.empty() + kwargs = mod_std_info(detector="ACIS-I") + fids = get_fid_catalog(stars=stars, **kwargs) + + stars.add_fake_constellation(n_stars=3, mag=7) + # add one star that would spoil a fid light + stars.add_fake_star(id=1001, mag=7.0, row=fids[0]["row"], col=fids[0]["col"]) + + # Get initial guide candidates + initial_guide_cands = get_guide_candidates(stars=stars, **kwargs) + + guides_with_fids = get_guide_catalog( + stars=stars, + initial_guide_cands=initial_guide_cands.copy(), + fids=fids, + **kwargs, + ) + assert 1001 not in guides_with_fids["id"] + + # Select guide stars using initial candidates + guides_no_fids = get_guide_catalog( + stars=stars, + initial_guide_cands=initial_guide_cands.copy(), + **kwargs, + ) + # Confirm that the fid-spoiling star was selected when no fids are present + # and that it still exists in the initial candidates + assert 1001 in guides_no_fids["id"] + initial_guide_cands_table = initial_guide_cands.copy() + assert 1001 in initial_guide_cands_table["id"] + + +def test_process_include_ids_columns(): + """Test that process_include_ids properly initializes columns for include_ids stars. + + This verifies that when stars are added via include_ids, they have the imposter + magnitude columns properly initialized and not NaN or zero. They aren't used + for force-included stars anyway, but might as well go through the motions. + """ + stars = StarsTable.empty() + stars.add_fake_constellation(n_stars=5, mag=[8, 9, 9.5, 10, 10.5]) + # Add a star that would normally be filtered out (too faint or bad position) + # but we'll force include it + stars.add_fake_star(id=999999, mag=11.0, yang=500, zang=500) + kwargs = mod_std_info( + att=stars.att, + n_guide=5, + ) + guide_cands = get_guide_candidates(stars=stars, **kwargs) + + guide_cands_table = guide_cands.copy() + assert 999999 not in guide_cands_table["id"] + kwargs = mod_std_info( + att=stars.att, + n_guide=5, + include_ids_guide=[999999], # Force include the faint star + ) + + guides = get_guide_catalog( + stars=stars, initial_guide_cands=guide_cands.copy(), **kwargs + ) + + # Verify the star was included + assert 999999 in guides.cand_guides["id"] + + # Get the row for the included star + idx = np.where(guides.cand_guides["id"] == 999999)[0][0] + included_star = guides.cand_guides[idx] + + # Verify imposter magnitude columns are initialized and not NaN or zero + assert not np.isnan(included_star["imp_mag"]) + assert included_star["imp_mag"] != 0.0 + assert not np.isnan(included_star["imp_r"]) + assert included_star["imp_r"] != 0.0 + assert not np.isnan(included_star["imp_c"]) + assert included_star["imp_c"] != 0.0 + + def test_make_report_guide(tmpdir): """ Test making a guide report. Use a big-box dither here diff --git a/proseco/tests/test_guide_fid_optimization.py b/proseco/tests/test_guide_fid_optimization.py new file mode 100644 index 00000000..29338bd2 --- /dev/null +++ b/proseco/tests/test_guide_fid_optimization.py @@ -0,0 +1,166 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst +""" +Tests for guide candidate awareness in fid selection and acq-fid-guide optimization. +""" + +import numpy as np +import pytest + +from proseco import get_aca_catalog +from proseco.core import StarsTable +from proseco.fid import FidTable, get_fid_catalog +from proseco.guide import get_guide_catalog +from proseco.tests.test_common import mod_std_info + + +# Create a synthetic test for fid trap scoring and detection +@pytest.mark.parametrize("row", [-200, 200]) +def test_fid_trap_scoring(row): + stars = StarsTable.empty() + stars.add_fake_constellation(n_stars=4) + # For the nominal ACIS-I FID 6 position here of -69.87 345.27 + # a star will cause the fid trap if it has about the same distance + # to the register as the fid distance to the trap. + # Either positive or negative row works so this test is parameterized. + stars.add_fake_star(id=1001, mag=7.0, row=row, col=400) + args = mod_std_info(detector="ACIS-I") + fids = get_fid_catalog(stars=stars, guide_cands=stars, **args) + assert fids.cand_fids.get_id(6)["fid_trap_spoiler"] + # And for these cases, confirm 6 is not selected + assert 6 not in fids["id"] + + +trap_test_cases = [ + {"spoil_mags": [8, 10, 11, 8, 10, 11], "expected_fids": [2, 3, 5]}, + {"spoil_mags": [11, 8, 10, 10, 8, 10], "expected_fids": [1, 3, 4]}, + {"spoil_mags": [11, 7, 11, 7, 7, 11], "expected_fids": [1, 3, 6]}, +] + + +@pytest.mark.parametrize("case", trap_test_cases) +def test_fid_trap_optimize(case): + stars = StarsTable.empty() + stars.add_fake_constellation(n_stars=1, mag=[7.0]) + + # For the nominal ACIS-I FID 6 position here of -69.87 345.27 + # a star will cause the fid trap if it has about the same distance + # to the register as the fid distance to the trap. + # So this places a star to spoil fid 6 via the fid trap if used + # as a guide star. + stars.add_fake_star(id=1001, mag=7.0, row=200, col=400) + kwargs = mod_std_info(detector="ACIS-I", t_ccd=-5) + + initial_fids = get_fid_catalog(stars=stars, **kwargs) + + # Now, spoil all of the candidate fids on purpose so that optimization + # is required. These stars are all theoretically fine for acquisition and guide + # so optimization is going to need to decide to use the fid or the + # star in each case. + for fid, mag in zip(initial_fids.cand_fids, case["spoil_mags"]): + stars.add_fake_star( + id=2000 + fid["id"], + mag=mag, + row=fid["row"], + col=fid["col"], + ) + + fids = get_fid_catalog(stars=stars, guide_cands=stars, **kwargs) + # No fids should be selected as all are spoiled + assert len(fids) == 0 + + # However, the full catalog selection should select a fid catalog + # via optimization. + aca = get_aca_catalog(stars=stars, **kwargs, optimize=True) + assert len(aca.fids) == kwargs["n_fid"] + + # The selected fids should match expected + assert set(aca.fids["id"]) == set(case["expected_fids"]) + + # And none of the fid spoiler stars for the 3 selected fids should be + # used as guide stars in the the catalog + for fid in aca.fids: + dys = np.abs(fid["yang"] - aca.guides["yang"]) + dzs = np.abs(fid["zang"] - aca.guides["zang"]) + spoiled = (dys < 5) & (dzs < 5) + assert not np.any(spoiled) + + # And if fid 6 is used the trap star should not be used as a guide + if 6 in aca.fids["id"]: + assert 1001 not in aca.guides["id"] + + # And the log should show that optimization occurred + log_opt = [ + rec["data"] + for rec in aca.log_info["events"] + if "optimize_acqs_fids" == rec["func"] + ] + assert len(log_opt) > 2 + + +def test_fid_spoil_short_circuits_optimization(): + stars = StarsTable.empty() + stars.add_fake_constellation(n_stars=4, mag=[8, 9, 9.5, 10]) + + # For the nominal ACIS-I FID 6 position here of -69.87 345.27 + # a star will cause the fid trap if it has about the same distance + # to the register as the fid distance to the trap. + # Either positive or negative row works so this test is parameterized. + stars.add_fake_star(id=1001, mag=7.0, row=200, col=400) + kwargs = mod_std_info(detector="ACIS-I") + + aca = get_aca_catalog(stars=stars, **kwargs, optimize=True) + assert 6 not in aca.fids["id"] + + # And the log should show that the optimization function was called but not needed + log_opt = [ + rec["data"] + for rec in aca.log_info["events"] + if "optimize_acqs_fids" == rec["func"] + ] + assert "No acq-fid optimization required" == log_opt[-1] + + +guide_test_cases = [ + {"mags": [6, 6, 7, 7, 8, 8], "expected_fids": [3, 5, 6]}, + {"mags": [10, 6, 6, 10, 9, 9], "expected_fids": [1, 4, 5]}, + {"mags": [8, 8, 8, 8, 8, 8], "expected_fids": [1, 5, 6]}, + {"mags": [10, 10, 10, 6, 6, 6], "expected_fids": [1, 2, 3]}, +] + + +@pytest.mark.parametrize("case", guide_test_cases) +def test_guide_fid_optimization(case): + stars = StarsTable.empty() + kwargs = mod_std_info(detector="ACIS-I") + initial_fids = get_fid_catalog(stars=stars, **kwargs) + # Now, spoil all of the candidate fids on purpose with bright stars + for fid, mag in zip(initial_fids.cand_fids, case["mags"]): + stars.add_fake_star( + id=2000 + fid["id"], + mag=mag, + row=fid["row"], + col=fid["col"], + ) + aca = get_aca_catalog(stars=stars, **kwargs, optimize=True) + assert len(aca.fids) == kwargs["n_fid"] + assert set(aca.fids["id"]) == set(case["expected_fids"]) + + +def test_backward_compatibility_no_guide_cands(): + """Verify that fid selection works without guide_cands (backward compat).""" + kwargs = mod_std_info(n_fid=3) + + # Call without guide_cands argument (old behavior) + fids = get_fid_catalog(**kwargs) + + assert isinstance(fids, FidTable) + + +def test_guide_catalog_without_guides_arg(): + """Verify get_guide_catalog works without initial_guide_cands argument (backward compat).""" + kwargs = mod_std_info(n_guide=5) + + # Call without guides argument (old behavior) + guides = get_guide_catalog(**kwargs) + + assert len(guides) > 0 or len(guides.cand_guides) > 0 diff --git a/proseco/tests/test_mon_full_cat.py b/proseco/tests/test_mon_full_cat.py index 47d50cf8..d3cb7403 100644 --- a/proseco/tests/test_mon_full_cat.py +++ b/proseco/tests/test_mon_full_cat.py @@ -193,3 +193,79 @@ def test_mon_takes_guide(): assert aca[TEST_COLS + ["mag"]].pformat() == exp assert aca.n_guide == 4 # One of 5 available guide slots becomes a MON + + +def test_mon_takes_guide2(): + """Test that commanding a MON on top of a potential GUI star correctly + precludes use of that star and takes one of the available n_guide slots. + + Also checks that when a MON is identified as a star in the catalog, the + 'id' and 'mag' correspond to the catalog. This overrides the mag in + `monitors`. + + Unlike test_mon_takes_guide here the MON is AUTO rather than MON_TRACK. + Unlike test_full_catalog, here the MON position coincides with a star in + the catalog but that star is not a GUI candidate because of its CLASS. + """ + stars = StarsTable.empty() + stars.add_fake_constellation( + n_stars=4, mag=6.5 + np.arange(4) * 0.5, id=np.arange(100000, 100004) + ) + stars.add_fake_star(yang=750, zang=750, mag=8.0, id=100005, CLASS=3) + monitors = [[750, 750, MonCoord.YAGZAG, 8.123, MonFunc.AUTO]] + + aca = get_aca_catalog( + **mod_std_info(att=stars.att, n_fid=0, n_guide=5, n_acq=5), + stars=stars, + monitors=monitors, + ) + + # This has 4 + 1 stars in the "guide" catalog, where the MON has taken one slot. + exp = [ + "slot idx id type sz yang zang dim res halfw mag ", + "---- --- ------ ---- --- -------- -------- --- --- ----- -----", + " 0 1 100000 BOT 8x8 1500.00 0.00 28 1 160 6.50", + " 1 2 100001 BOT 8x8 0.00 1500.00 28 1 160 7.00", + " 2 3 100002 BOT 8x8 -1500.00 0.00 28 1 160 7.50", + " 3 4 100003 BOT 8x8 0.00 -1500.00 28 1 160 8.00", + " 7 5 100005 MON 8x8 750.00 750.00 0 0 20 8.00", + ] + + assert aca[TEST_COLS + ["mag"]].pformat() == exp + assert aca.n_guide == 4 # One of 5 available guide slots becomes a MON + + +def test_obsid_30843(): + """ + Test real case with obsid 30843 with an automatic monitor request that + ends up as a MON. This test most closely matches test_mon_takes_guide2 but + uses real data and was the basis for discovering/fixing the issue. + """ + args = { + "obsid": 30843.0, + "raise_exc": True, + "att": [-0.50655663, -0.04177606, 0.85826295, 0.07099192], + "date": "2025:074:22:39:37.088", + "detector": "ACIS-S", + "dither_acq": [15.9984, 15.9984], + "dither_guide": [15.9984, 15.9984], + "focus_offset": 0.0, + "man_angle": 124.90441066524069, + "n_acq": 8, + "n_fid": 3, + "n_guide": 5, + "sim_offset": -2442.0, + "t_ccd_acq": -6.367734460876512, + "t_ccd_guide": -3.9378722816242115, + "dyn_bgd_n_faint": 2, + "dyn_bgd_dt_ccd": -4.0, + "target_offset": [0.0, 0.0], + "target_name": "WR 25", + "duration": 39000.0, + "man_angle_next": 157.2577020420281, + "monitors": [[161.043292, -59.719753, 0.0, 8.8, 0.0]], + } + aca = get_aca_catalog(**args) + mon_id = 1130642984 + assert mon_id in aca["id"] + assert mon_id in aca.mons["id"] diff --git a/validate/make-regression-data.ipynb b/validate/make-regression-data.ipynb new file mode 100644 index 00000000..5971bddb --- /dev/null +++ b/validate/make-regression-data.ipynb @@ -0,0 +1,240 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Make regression data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "WARNING: IERSStaleWarning: leap-second file is expired. [astropy.utils.iers.iers]\n" + ] + } + ], + "source": [ + "import gzip\n", + "import pickle\n", + "import sys\n", + "from pathlib import Path\n", + "from astropy.table import Table\n", + "import astropy.units as u\n", + "from tqdm import tqdm\n", + "\n", + "sys.path.insert(0, str(Path.home() / \"git\" / \"proseco\"))\n", + "#sys.path.insert(0, str(Path.home() / \"git\" / \"sparkles\"))\n", + "\n", + "import kadi.commands as kc\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import parse_cm.paths\n", + "from astropy.table import Table\n", + "from cxotime import CxoTime\n", + "\n", + "from proseco import get_aca_catalog\n", + "import sparkles\n" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "5.17.1\n" + ] + } + ], + "source": [ + "import proseco\n", + "\n", + "print(proseco.__version__)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "4.31.0\n" + ] + } + ], + "source": [ + "print(sparkles.__version__)" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "# Get the last 2.5 years of approved schedules\n", + "obss = kc.get_observations(CxoTime.now() - 365 * 2.5 * u.day)\n", + "obss = [obs for obs in obss if obs[\"source\"] != \"CMD_EVT\"]\n", + "load_names = []\n", + "[load_names.append(obs[\"source\"]) for obs in obss if obs[\"source\"] not in load_names];\n" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "def print_acar(acar, load_name):\n", + " lines = []\n", + " lines.append(\"-\" * 40)\n", + " lines.append(f\"ACA Review for {load_name}\")\n", + " lines.append(f\"obsid: {acar.obsid}\")\n", + " lines.append(f\"target_name: {acar.target_name}\")\n", + " lines.append(f\"date: {acar.date}\")\n", + " lines.append(f\"P2: {-np.log10(acar.acqs.calc_p_safe())}\")\n", + " lines.append(f\"guide count: {acar.guide_count}\")\n", + " lines.extend(acar.pformat())\n", + " for msg in acar.messages:\n", + " lines.append(f\"MSG: {msg}\")\n", + " text = \"\\n\".join(lines)\n", + " return text\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# make an output dir that includes proseco version and sparkles version\n", + "outdir = Path(\"regression_data\") / f\"{proseco.__version__}_{sparkles.__version__}\"\n", + "outdir.mkdir(parents=True, exist_ok=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + " 91%|█████████ | 137/151 [21:01<01:55, 8.27s/it]/Users/jean/miniforge3/envs/ska3-latest/lib/python3.12/site-packages/chandra_aca/star_probs.py:438: UserWarning: \n", + "Model grid-local-quadratic-2025-10.fits.gz computed between -15 <= t_ccd <= 2, clipping input t_ccd(s) outside that range.\n", + " warnings.warn(\n", + "/Users/jean/miniforge3/envs/ska3-latest/lib/python3.12/site-packages/chandra_aca/star_probs.py:438: UserWarning: \n", + "Model grid-local-quadratic-2025-10.fits.gz computed between -15 <= t_ccd <= 2, clipping input t_ccd(s) outside that range.\n", + " warnings.warn(\n", + "/Users/jean/miniforge3/envs/ska3-latest/lib/python3.12/site-packages/chandra_aca/star_probs.py:438: UserWarning: \n", + "Model grid-local-quadratic-2025-10.fits.gz computed between -15 <= t_ccd <= 2, clipping input t_ccd(s) outside that range.\n", + " warnings.warn(\n", + "100%|██████████| 151/151 [23:06<00:00, 9.18s/it]\n" + ] + } + ], + "source": [ + "\n", + "out_text_file = outdir / \"catalogs.txt\"\n", + "# Remove existing text file\n", + "if out_text_file.exists():\n", + " out_text_file.unlink()\n", + "for load_name in tqdm(load_names):\n", + " path_pkl_out = outdir / f\"{load_name}_acas.pkl.gz\"\n", + " if path_pkl_out.exists():\n", + " new_acas = pickle.load(gzip.open(path_pkl_out, \"rb\"))\n", + " for obsid in new_acas:\n", + " acar = new_acas[obsid].get_review_table()\n", + " acar.run_aca_review()\n", + " text = print_acar(acar, load_name)\n", + " with open(out_text_file, \"a\") as f:\n", + " f.write(text + \"\\n\")\n", + " else:\n", + " path_pkl = parse_cm.paths.load_file_path(load_name, \"output/*.pkl.gz\")\n", + " acas = pickle.load(gzip.open(path_pkl, \"rb\"))\n", + " # Use a little table to keep everything in time order\n", + " obsids = Table([{\"date\": acas[obsid].date, \"obsid\": obsid} for obsid in acas])\n", + " obsids.sort(\"date\")\n", + " new_acas = {}\n", + " for obsid in obsids:\n", + " args = acas[obsid[\"obsid\"]].call_args.copy()\n", + " args[\"raise_exc\"] = True\n", + " has_edits = False\n", + " for kwarg in [\"include_ids_guide\", \"include_ids_acq\", \"include_halfws_acq\",\n", + " \"exclude_ids_guide\", \"exclude_ids_acq\"]:\n", + " if kwarg in args:\n", + " args.pop(kwarg)\n", + " has_edits = True\n", + " if has_edits:\n", + " # Do it once without the edits and again with the as-run edits\n", + " new_aca = get_aca_catalog(**args)\n", + " new_acas[obsid[\"obsid\"] + .1] = new_aca\n", + " acar = new_aca.get_review_table()\n", + " acar.run_aca_review()\n", + " text = print_acar(acar, load_name)\n", + " with open(out_text_file, \"a\") as f:\n", + " f.write(text + \"\\n\")\n", + " # With edits it looks like there's more of a chance with a force-include\n", + " # not working, so catch exceptions here\n", + " try:\n", + " new_aca = get_aca_catalog(**acas[obsid[\"obsid\"]].call_args)\n", + " except ValueError as e:\n", + " text = f\"ERROR for obsid {obsid['obsid']} load {load_name}: {e}\"\n", + " with open(out_text_file, \"a\") as f:\n", + " f.write(text + \"\\n\")\n", + " continue\n", + " new_acas[obsid[\"obsid\"]] = new_aca\n", + " acar = new_aca.get_review_table()\n", + " acar.run_aca_review()\n", + " text = print_acar(acar, load_name)\n", + " with open(out_text_file, \"a\") as f:\n", + " f.write(text + \"\\n\")\n", + "\n", + " with gzip.open(path_pkl_out, \"wb\") as f:\n", + " pickle.dump(new_acas, f)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "ska3-latest", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.8" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/validate/pr410-regression-optimization.ipynb b/validate/pr410-regression-optimization.ipynb new file mode 100644 index 00000000..1ba3594d --- /dev/null +++ b/validate/pr410-regression-optimization.ipynb @@ -0,0 +1,597 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "959bd284", + "metadata": {}, + "source": [ + "## Review regression outputs for PR410\n", + "\n", + "This notebook uses the reselected pkl files made with make-regression-data.ipynb.\n", + "\n", + "The pickle files over a date range for dir1 (master and last release) and dir2 (current PR output) are reviewed for diffs.\n", + "\n", + "For each diff, the logs and other metadata are reviewed to infer if the catalog change is likely reasonable based on new guide candidate / fid light optimization." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "58f4bdb6", + "metadata": {}, + "outputs": [], + "source": [ + "import kadi.commands as kc\n", + "from cxotime import CxoTime\n", + "from astropy.table import Table\n", + "import astropy.units as u\n", + "from pathlib import Path\n", + "import pickle\n", + "import gzip\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "55121928", + "metadata": {}, + "outputs": [], + "source": [ + "dir1 = \"regression_data/5.17.2.dev1+g6fa410d_4.31.0\" # generated from the opt-log branch\n", + "dir2 = \"regression_data/5.17.2.dev51+ge5ee6d8_4.31.0\" # from guide-opt-fid branch" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "c307d58c", + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "# Get the last 2.5 years of approved schedules\n", + "obss = kc.get_observations(\"2023-07-04\", \"2025-12-30\")\n", + "obss = [obs for obs in obss if obs[\"source\"] != \"CMD_EVT\"]\n", + "load_names = []\n", + "[load_names.append(obs[\"source\"]) for obs in obss if obs[\"source\"] not in load_names];\n" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "1b0e8c48", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "'JUL0323A'" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "load_names[0]" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "236ba6f2", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "'DEC2925B'" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "load_names[-1]" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "edb0f2a7", + "metadata": {}, + "outputs": [], + "source": [ + "from proseco.diff import catalog_diff" + ] + }, + { + "cell_type": "markdown", + "id": "b039604e", + "metadata": {}, + "source": [ + "This method checks the new \"fid_trap_spoiler\" flag in the fid candidate table to see if one or more of the candidate fids could have been spoiled by a guide star and its relationship to the fid trap.\n", + "\n", + "This method also checks the selection logs to see if any fid lights were removed from candidacy because they would spoil guide star candidate(s) (not via fid trap issue).\n", + "\n", + "The intent with the method was just to scan all of the new pkl files and their proseco catalogs to confirm that each change corresponded to at least one of those kinds of changes and throw an exception otherwise. The method did not throw any exceptions." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "b3dc113a", + "metadata": {}, + "outputs": [], + "source": [ + "def review_cat(cat1, cat2):\n", + " cdiff = catalog_diff(cat1, cat2, style=\"unified\")\n", + " if len(cat1.fids) == 0 and len(cat2.fids) == 0:\n", + " if cdiff.text != \"\":\n", + " raise ValueError(\"No fids but catalog diff\")\n", + " return None\n", + " trap_spoiler = np.any(cat2.fids.cand_fids[\"fid_trap_spoiler\"])\n", + " direct_guide_spoiler = np.any([entry[\"data\"].endswith(\"spoils a guide candidate\") for entry in cat2.fids.log_info[\"events\"]])\n", + " if cdiff.text == \"\":\n", + " return None\n", + " elif trap_spoiler and direct_guide_spoiler:\n", + " return \"Both trap and direct guide spoilers\"\n", + " elif trap_spoiler:\n", + " return \"Trap spoiler\"\n", + " elif direct_guide_spoiler:\n", + " return \"Direct guide spoiler\"\n", + " else:\n", + " raise ValueError(\"Diff text but no obvious cause\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "57a35e8a", + "metadata": {}, + "outputs": [], + "source": [ + "from collections import defaultdict" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "9d896fd1", + "metadata": {}, + "outputs": [], + "source": [ + "counts = defaultdict(int)\n", + "for load_name in load_names:\n", + " pkl1 = Path(f\"{dir1}/{load_name}_acas.pkl.gz\")\n", + " pkl2 = Path(f\"{dir2}/{load_name}_acas.pkl.gz\")\n", + " if not (pkl1.exists() and pkl2.exists()):\n", + " continue\n", + " acas1 = pickle.load(gzip.open(pkl1, \"rb\"))\n", + " acas2 = pickle.load(gzip.open(pkl2, \"rb\"))\n", + " for obsid in acas1.keys():\n", + " cmp = review_cat(acas1[obsid], acas2[obsid])\n", + " if cmp is not None:\n", + " counts[cmp] = counts[cmp] + 1\n", + " else:\n", + " counts[\"no_diff\"] = counts[\"no_diff\"] + 1\n" + ] + }, + { + "cell_type": "markdown", + "id": "427571ad", + "metadata": {}, + "source": [ + "The review of the data shows that in most cases (~90%) there is no change in the catalog, but in the other 10% of cases, the use of the guide star information does change the fid selection and the resulting catalog (either because of the fid trap effect or because the fid would be within a guide star search/track region). " + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "d9fe6ef6", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "defaultdict(int,\n", + " {'no_diff': 6122,\n", + " 'Trap spoiler': 557,\n", + " 'Direct guide spoiler': 66,\n", + " 'Both trap and direct guide spoilers': 39})" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "counts" + ] + }, + { + "cell_type": "markdown", + "id": "eee471a0", + "metadata": {}, + "source": [ + "In the next chunk of code, we review how the P2 and guide count behave in the catalogs that have changed." + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "666a0c7a", + "metadata": {}, + "outputs": [], + "source": [ + "compare = []\n", + "for load_name in load_names:\n", + " pkl1 = Path(f\"{dir1}/{load_name}_acas.pkl.gz\")\n", + " pkl2 = Path(f\"{dir2}/{load_name}_acas.pkl.gz\")\n", + " if not (pkl1.exists() and pkl2.exists()):\n", + " continue\n", + " acas1 = pickle.load(gzip.open(pkl1, \"rb\"))\n", + " acas2 = pickle.load(gzip.open(pkl2, \"rb\"))\n", + " for obsid in acas1.keys():\n", + " cmp = review_cat(acas1[obsid], acas2[obsid])\n", + " if cmp is not None:\n", + " # What are the p2 and guide_count of the original and new catalog?\n", + " acar1 = acas1[obsid].get_review_table()\n", + " guide_count_1 = acar1.guide_count\n", + " acar2 = acas2[obsid].get_review_table()\n", + " guide_count_2 = acar2.guide_count\n", + " p2_1 = -np.log10(acar1.acqs.calc_p_safe())\n", + " p2_2 = -np.log10(acar2.acqs.calc_p_safe())\n", + " compare.append((load_name, obsid, cmp, guide_count_1, guide_count_2, p2_1, p2_2))\n" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "id": "e6cc3ce8", + "metadata": {}, + "outputs": [], + "source": [ + "# Oops, should made that easier to Table\n", + "compare_list_of_dicts = []\n", + "for row in compare:\n", + " compare_list_of_dicts.append({\n", + " \"load_name\": row[0],\n", + " \"obsid\": row[1],\n", + " \"diff_type\": row[2],\n", + " \"guide_count_1\": row[3],\n", + " \"guide_count_2\": row[4],\n", + " \"p2_1\": row[5],\n", + " \"p2_2\": row[6],\n", + " })" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4d4026c6", + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "dat = Table(compare_list_of_dicts)" + ] + }, + { + "cell_type": "markdown", + "id": "7b7f1789", + "metadata": {}, + "source": [ + "The 4 plots in the cells below show that the P2 is overwhelmingly unchanged, but in a few cases is a bit worse. The guide_count is generally the same or better. This is not unexpected as the PR now includes the guide candidates in the fid selection process (biasing an improvement in guide_count) and it is not surprising that, with the current optimization strategy, a P2 value lower than the previous estimate could be picked, especially if that P2 value is > 2.0 (the passable threshold in the optimization logic)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e5519d50", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAioAAAGwCAYAAACHJU4LAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAARHVJREFUeJzt3Qt8znX/x/HPzgdszBzmPOfDmkO4c8xpkkPUXWi5k1O5nZOKEEJU/7qJolQo90ipW5hzOaUiNjMrpxQ3Qwsbxo7X//H9dU8bw3b9rmvX77p+r+fjcbV2Xde+19d31+G979HNYrFYBAAAwIDcHV0BAACA2yGoAAAAwyKoAAAAwyKoAAAAwyKoAAAAwyKoAAAAwyKoAAAAw/IUJ5adnS1nzpyREiVKiJubm6OrAwAACkBt4Xb58mWpUKGCuLu7u25QUSGlcuXKjq4GAACwwqlTp6RSpUquG1RUT0rOPzQgIMBuj5ORkSGbNm2Szp07i5eXl5gRbUA78FzgNcH7Au+PtvqMSElJ0Toacj7HXTao5Az3qJBi76Di7++vPYaZg4rZ20ChHWgDngu8HnhfsN17Y0GmbTCZFgAAGBZBBQAAGJZDg4qa8TtmzBipWrWq+Pn5ScuWLWXv3r2OrBIAADAQhwaVwYMHy+bNm+WTTz6RgwcPapNwOnXqJKdPn3ZktQAAgNmDyrVr12TVqlXy+uuvS9u2baVmzZoydepUCQ0NlQULFjiqWgAAwEActuonMzNTsrKyxNfXN8/1agho165d+f5MWlqadsm9vClntrG62EtO2fZ8DKOjDWgHngu8Jnhf4P3RVp8Rhfk8dbOo7eEcRM1J8fb2lqioKClXrpwsX75cnnzySalVq5YcPnz4lvurHpdp06bdcr36ebUsCgAAGF9qaqpERkZKcnLyXbcXcWhQOX78uAwcOFB27NghHh4e0qRJE6ldu7bs379fEhISCtSjojaMSUpKsvs+KmouTUREhGn3EKENaAeeC7wmeF/g/dFWnxHq8zs4OLhAQcWhG77VqFFDtm/fLlevXtUqHRISIn369NHmqeTHx8dHu9xMNUpRBIiiehwjow1oB54LvCZ4X+D9Ue9nRGE+Sw2xM22xYsW0y8WLF2Xjxo3aBFsAAOA4WdkW2XPigpy/fF3KlvCV5qFB4uFe9AcAOzSoqFCiRp7q1Kkjx44dk+eff177/wEDBjiyWgAAmNqG+ESZtiZBEpOv37guJNBXpvSoL13CQsyzj4oamxo+fLjUrVtXm0TbunVr7VAjsw+vAADgiB6Ub48lyT+X7ZOhy/bnCSnK2eTr8s9l+7UQY5oeld69e2sXAADgOBviE2X8FwflUurtlw2rlTdq4Ef1tETUL19kw0Cc9QMAgMlDytBl++8YUnKHFdXTouauFBWCCgAAJnUtPUvGfXag0D+nJtgWFUOs+gEAAEVrVnSCvL/zhFizm5paBVRUCCoAAJhs0uyo5TGy7mDhJ8WqWSnlA/9cqlxUCCoAAJgkoMzZclgWfHNcMnXsSa+WKBflfioEFQAATDBhdtSKWEnPzLa6DEfto0JQAQDABKt69Hi2Uy0Z0aGW+XamBQAA9hvq2XbonAz7t/UhpaS/l8x+5J4i70XJjaACAICLWRt7WhvqsXagR/WbLOnfTFrXKeOQXpTcCCoAALiQIR/vlc0J53WV8XTbULm/XlkxAjZ8AwDARUxfe0h3SBnSJlQmdK0vRkGPCgAALjAfZe6WI/Lhrl91lTO/b2Pp3qiCGAlBBQAAJ7Y29oyMW3VArmdYv/TYx9Nd5vZt5NBJs7dDUAEAwEl7Ufq8t1t+/O2SrnK63VNe3n68icMnzd4OQQUAACcTHZcoY1fGynUdG7g937mWDGlbU7w9jT1dlaACAIAT9aKMXhEja+MKf05PbuGVAmR4h9riDAgqAAA4SS/K858fkKvpWbpDylcj2oizIKgAAGBwM9clyKKdJ3SVUbmkr6wfc78U93Wuj37nqi0AACYb6hkZtU+i48/pKmfe442lR0NjLTsuKIIKAAAGXXb87MpYyci2WF1GST9Pmf33cEMuOy4oggoAAAbz1Ec/yLYjSbrKaFq1pHz6TEvDLjsuKIIKAAAGEvbyermSbv2yY6Vrg3Ly7j+aiisgqAAAYBA1X1onOrZGueGBe5x3qOdmxt7lBQAAE7hwJV2qj7dNSFHKlvAVV0GPCgAADtR0+iZJupphk7LcRKR8oK80Dw0SV0FQAQDAQepOWq9rG/zccqbMTulR3+kn0OZGUAEAwAH7o7R7fYvNQooSVMxbejaqIIF+3lr5rhJWCCoAABShDfGJMuzf+0XH9ihSroS3bHu+g8SeuiRbEs7Kl7Gn5Y+r6fLRt79ql5BAX61nxZn3T8nBZFoAAIqA6uWYs/mwDF2mL6SUKe4tP0yMED9vD0m+9mc4uXDTHJezydfln8v2a6HI2dGjAgCAna368ZQ893mc7nI61AmWjwb87UbwmbYmQfLLPOo6NfCjbo+oX96ph4EIKgAA2NG9r2ySP1L1rerx93KTfZMf0HpRcuw5cUESk6/f9mdUWFG3q/u1qFFanBVBBQAAO0jPzJY6k9bn2+NRGHMebSi9mla65frzl28fUqy5n1ERVAAAsLEpX8XL0t2/6SrD28NNfpr+4G2HbQq6qZuzb/5GUAEAwIbunb5J/tC5gVvdcv6y4dn2d7yP2tRNre5RE2ctLrz5G6t+AACwgeTUDKk2fp3ukBJWocRdQ4qielrUEmTFzYU3fyOoAACg0yv73KTprG9s0o4jOtQq8H27hIXIgn5NtJ6T3NT36npX2EeFoR8AAHT0ojR8ZZPN/u63Zklxl7AQ7f5qdY+aOKvmpKjhHmfvSclBUAEAwAptX98qJy/krKixTSiwdkmxh7ubUy9BNuzQT2ZmpkyaNElCQ0PFz89PqlevLq+88opkZ9vu7AMAAGyt6YzNuUKK7Tn7kmKX6VF57bXXZOHChbJ06VJp0KCB/PjjjzJgwAAJDAyU0aNHO7JqAADcQu0G+8g7OyXpSrpdW8fZlxS7TFD57rvvpGfPntKtWzft+2rVqsny5cu1wAIAgJGoc3Oe/TRWrmXo6/X3dHfTAo8rLyl2maDSunVrrUflyJEjUrt2bTlw4IDs2rVL5syZk+/909LStEuOlJQU7WtGRoZ2sZecsu35GEZHG9AOPBd4TZj5fWHjoXMyYsUB3eW89eg94u3pLiNXHNBCSe6wkjPLZeKDdSQ7K1Oys3Q/nGGfC4V5zrhZLBa9u/taTT30Sy+9pA0BeXh4SFZWlsycOVMmTJiQ7/2nTp0q06ZNu+X6qKgo8ff3L4IaAwDMRJ1ynHDRTRYddtcxadYiPpIts++zSM5CnAN/uMkXv7rLpfS/yivpbZFHqmVLw9IO+1guMqmpqRIZGSnJyckSEBBg3KCyYsUKef755+WNN97Q5qjExsbKmDFj5K233pL+/fsXqEelcuXKkpSUdNd/qB4q+W3evFkiIiLEy8tLzIg2oB14LvCaMNv7wurYM/L8qnjdZ/U8eV8lmdztz43ZclPDPz/+dlHOX06TsiV8pGnVUk67pDijkM8F9fkdHBxcoKDi0KEfFVLGjx8vffv21b6/55575LfffpNZs2blG1R8fHy0y81UoxTFi6SoHsfIaAPagecCrwkzvC+0mb1FTl366w9ja6jIcXjGg9pQT35Ui7WuXU7M+FzwKsTzxdPRXT/u7nl/gWoIiOXJAABHUL0cDV5eL9cz9fWjBPl7yv6XH7BZvczMoUGlR48e2pyUKlWqaEM/MTEx2rDPwIEDHVktAIAJrTlwRkYuj9FdTrtaQbJkUAub1AkODirz5s2TyZMny7Bhw+T8+fNSoUIFeeaZZ+Tll1/mdwMAKDKDl+6VLT+d113OkDahMjGf+Shw0qBSokQJbSny7ZYjAwBgT1euZ8r9b3yt+8RjL7csiX35ASnmd+s8SujDWT8AAFN6aN5OiTv9535cegQX95bJ96TedtIs9CGoAABMfqCgdfy93WXXCx2lhI+bREdH26xuyIugAgAwlac+/F53SOlYt6x8+FQz0+zM60gEFQCAKaRnZku3t3fI0fNXdZUzoFUVmdLjHpvVC3dGUAEAuPzeKMP+vU87r0evZ9qGyoSurOopSgQVAIBLn3g8dNl+3eV4uokkTL/9LrOwH1ocAOCSouPO2CSktKtdWo7N6kZIcRB6VAAALmdt7BkZsSJG91/ycx9vLD0aVrBZvVB4BBUAgEuZvvaQfLjrV11l3FslUFYObeW0pxm7EoIKAMAlXEvPko5vfiNnkvWdejy/byPp3qiizeoFfQgqAACnX9XT+73dsu+3S7rKaVKlpHw2tCW9KAZDUAEAOK21sadl9KexkmXRV87cPo2kZ2N6UYyIoAIAcEqDl+6RLT/9rrucdyObSNfwEJvUCbZHUAEAOJ2Bi/fI14f1h5SF/ZpIlzBCipERVAAATmXgEv0hJcDHQ2KmPMB8FCdAUAEAOM9ZPXO3y9HfU3WVE1YhQNaOamOzesG+CCoAAMObseaQfPCtvr1R1I4oc3s3koeaMGnWmRBUAACGXnoc8a9t8ovOXpTQYH/ZMrYdQz1OiKACADCk6LhEGfNpjKTrXHs8uFWoTOrBicfOiqACADCc6WsT5MNdJ3SV0atRBXn90YYcJujkCCoAAINtg/+1nElO11XOu5GNpWs4hwm6AoIKAMAQBi7+Qb4+nKS7HDZwcy0EFQCAwyfMNpq2US6nZekqx8fTXeb2bcQGbi6GoAIAcJgN8Ynyz2X7RedRPdI1rJzMi7yXVT0uiKACAHCIr/b/V0atPKCrDA93kXl9OavHlRFUAABFPtTzyDs75cDpy7rK6VK/jLzTrxm9KC6OoAIAKNKhnqHL9usuZ1DrUJncnb1RzICgAgAoEmtjT8uIFbG6y3mmbahM6EpIMQuCCgDA7l5Zc0g+0nlWT+2y/rJ21P1s4GYyBBUAgF03cGs5a4tcvJapq5x5jzeWHg3ZwM2MCCoAALsYsGSPfPPz77rKCPD1lJiXOzNh1sQIKgAAm2s6Y5MkXcnQVUb7umVk8VPNbVYnOCeCCgDAprrN3a4rpISU8Javn+8gft4eNq0XnBNBBQBgs/1RRkTtk0OJV6wuo7iPh3w3MYLfCG4gqAAAdEnPzJZxn/0oq+PO6SqntL+n7Hv5AX4byIOgAgCw2pe/uMnoaVt0t2C7OsGyZMDf+E3gFgQVAIBVer27Ww6dc9fVeu5uInP7svQYt0dQAQAUei5Ki1e3yPkr6SLiZnXrVQvyla3jOrD0GHekLwrrVK1aNXFzc7vlMnz4cEdWCwBwG9FxiVJnUvT/Qor17qkYINte6EhIgbF7VPbu3StZWVk3vo+Pj5eIiAh57LHHHFktAEA+vSijV8TI2rhE3W0zqHU1mdy9AW0M4weVMmXK5Pl+9uzZUqNGDbn//vsdVicAwK0nHr+4Kk6SdW6D37xaSVk2uAVn9cA556ikp6fLsmXLZOzYsdrwT37S0tK0S46UlBTta0ZGhnaxl5yy7fkYRkcb0A48F8z5mlgXlyhjPjuoq4xi3h4y6+EG8mBYeRFLlmRk/NWT7grM8lywZRsUpq3cLBaLRQxg5cqVEhkZKSdPnpQKFfI/eGrq1Kkybdq0W66PiooSf3//IqglAJjH58fdZOd5dx0TZi3i654ts5pbtNU9QI7U1FTtMz85OVkCAgLEKYLKAw88IN7e3rJmzZrb3ie/HpXKlStLUlLSXf+heqjkt3nzZm3+jJeXl5gRbUA78Fwwz2tCbeDW8rVvJPm6vp6PKqX8ZOvYNuLqXPm5YK82UJ/fwcHBBQoqhhj6+e2332TLli3yxRdf3PF+Pj4+2uVmqlGK4slRVI9jZLQB7cBzwbVfEzPXHZJFO3/VXc6bj4bL35tWFjNxteeCPdugMO3k0OXJORYvXixly5aVbt26OboqAGDaVT2PvLvLJiFFjfL83+YjWpmAXg4PKtnZ2VpQ6d+/v3h6GqKDBwBMtzdKg5c3yP6TyTYpT8WTxOTrsufEBZuUB3NzeDJQQz5qAu3AgQMdXRUAMJ1Z0Qny3o4Tdin7/OXrdikX5uLwoNK5c2cxyHxeADANNSwzZ8thu4UUpWwJX7uVDfNweFABABT9Bm4vfHZAUtLss5+JmqNSPtBXmocG2aV8mIvD56gAAIo2pAxdtl93SImoX/Z/u6vk7RHP2S5lSo/6nOMDmyCoAICJhntGrYjRXc67kU1k0ZPNZF7fhlLSO+9tqidlQb8m0iUsRPfjAApDPwDg4tQGbkt3n5CoH36T9Ezr5wQ2rhQonw9rdaOn5IEG5STj1ywpU/8++SM1U5uTooZ7cm4HbIGgAgAubOa6BFm0U9+E2XrlissXw1uLn7fHLbepTPK30CDTb3QG+yGoAICLDvP0Xvit7NOxN4qvp7u81buRdA1nGAeOQ1ABABcTHXdGRq+IkYxs68sY1aGGjO5Uh2EcOBxBBQBcyNSv4mXJ7t90lTG3TyPp2biizeoE6EFQAQAXGeppOWuLnLucrqscH0938fFiQSiMg2cjALjA3ih1JkXrDik5K4T+uWy/ViZgBAQVAHBia2NPaxu4ZeqYj5JbzuLlaWsSOP0YhkBQAQAnNX1tgoxYEWvzcjn9GEbCHBUAcML5KI8u+FZiTlm/9LggOP0YRkBQAQAnsubAn0uPs4vg0HlOP4YRMPQDAE5i8NK9MnK5vpBSu2xx+emVLhIS6HvjAMGbqevV7Zx+DCMgqACAExi0ZI9s+em8rjKGtKkmm8ber22Fr043Vm4OK5x+DKMhqACAwU356qBs/fl3q38+JMBHjsx4UCZ2a3DjOnW6sTrlWJ12nBunH8NomKMCAAY+8fj97cfl96sZVpURVqGErHi6pRT3zf+tXoWViPrlZc+JC9rEWU4/hhERVADAYGZFJ8h7O/SdeDywVTV5ucdfPSi34+HuJi1qlNb1WIA9MfQDAAYyfW287pCirI8/y+6ycAkEFQAwyN4oz3y8Vz7cpe9AwRxnk6+zFT5cAkEFABwsOi5Rak2Mlo0J+lb15MZW+HAVBBUAcPB8lGFR++2ygRtb4cMVEFQAwIEHCtpiPsrdsBU+nBlBBQAcMB9lzuYjug8UbFIlsED3Yyt8ODOCCgAU8XyUhtM2ypytR60uQ+0eO6RNqHw2tBVb4cPlsY8KABRRL4o6THBtXKKucu6rHiQfD/ybeHv++Xem2gr/n8v2a+El9zQXtsKHq6BHBQCKoBclfOpGXSHFy0Pk3cjGsuLpFjdCisJW+HB19KgAgB3NXJcgi3bqmzBbLsBbdo/vpO0imx+2wocrI6gAgJ3MXHdIFu38VVcZHeuWkQ+fan7X+7EVPlwVQQUA7LGqZ9NhXSGlTY0geb9/c/Hz9rBp3QBnQ1ABABvPRxn7Waxcz8i2uox3I5tI1/AQfi8AQQUAjDMfRc1AeSeyMSEFyIUeFQCwwVDPyKh9Eh1/Tlc579CTAtyCoAIAOqw5kCgvfHFQMq0f6ZGSfp4y++/h2uodAHkRVADASv93wE1OfXfQ6vZT26Esfqq5tKwZfNulx/n13uw5cUE7v0dtjd88NKjAPws4I4IKABRSema2tHljm5xL1bdn5vzIJtKmdpkC339DfKJMW5MgicnXb1wXEuir7U5LbwxcFTvTAkAhTF97SGpPWi/nUtJzbVRfOAE+7rKwX5NChQsVUtRW+blDinI2+bp2vbodcEUODyqnT5+Wfv36SenSpcXf318aNWok+/btc3S1AOAWD83fKR/u0reBW/d7QiRmSpdChRQ13KN6UnKf5ZMj5zp1u7of4GocOvRz8eJFadWqlbRv317Wr18vZcuWlePHj0vJkiUdWS0AyEMFgOFR+yTuvym6WuaZtqEyoWv9Qv+cmpNyc09KbiqeqNvV/VrUKK2rjoDRODSovPbaa1K5cmVZvHjxjeuqVat22/unpaVplxwpKX++aWRkZGgXe8kp256PYXS0Ae1g1ufCxkPnZOJ/4iX5epbVZTSrWlKWPNVUO0zQmnZLvHS1wPfLyAiQomK258Lt0A5S6DYozHPGzWKxOKyvsH79+vLAAw/If//7X9m+fbtUrFhRhg0bJkOGDMn3/lOnTpVp06bdcn1UVJQ2bAQAthST5CZLjuaMkFszH8Ui7ctnS69QfW+zR5PdZH7C3bfSH1E/S2oFMvwD40tNTZXIyEhJTk6WgIAA4wYVX19f7evYsWPlsccekz179siYMWPkvffekyeffLJAPSqqRyYpKemu/1A9VPLbvHmzREREiJeXl5gRbUA7mOm5oIZ63tl2XOZ/80u+80IKQq0YnvNYuDx4T3mb1KfdmzvkXEpavvVREap8oI98M7ZtkS5VNsNzoSBoByl0G6jP7+Dg4AIFFYcO/WRnZ0vTpk3l1Vdf1b5v3LixHDp0SBYsWJBvUPHx8dEuN1ONUhQvkqJ6HCOjDWgHV38uRMedkRdWxcmVNOuHekKD/WXL2HY2Cw2qlac+1EBb3aNKzB1Wch5hSo8G4uvjbZPHK3T9XPS5UFi0gxS4DQrzfHHoqp+QkBBt+Ce3evXqycmTJx1WJwDm3Rul7/vfybCoGKtDioebyNw+jeSbce1t3rOhVgkt6NdEygf+2ROdQ32vrmcfFbgqh/aoqBU/hw8fznPdkSNHpGrVqg6rEwDzmRWdIO/vOGH1MI9SvYy/bH7Wdr0o+VFhJKJ+eXamhak4NKg8++yz0rJlS23op3fv3toclffff1+7AEBRhZT3dlh/4rEyqGVVmfxQmBQFFYRYggwzcWhQadasmXz55ZcyYcIEeeWVVyQ0NFTmzJkjTzzxhCOrBcAE1ATV3ceStJ4U61m0CbO97q1iw5oBsDqoJCYmytatWyUoKEg6deok3t5/Tdy6evWqvPnmm/Lyyy8Xpkjp3r27dgGAohIdlyiTVsfLhatqG3zrBPp5yt8rp0m3cE48BuypwJNp9+7dq018HT58uDz66KMSFhamrdDJceXKlXz3OAEAI5m+Nl6GRe3XFVLGdKwlP4xvLw1Ls2cJYJig8tJLL8kjjzyibXt/7tw5ba30/fffLzExMfatIQDYyICP9siHu36z+udL+XtphwmOiahdpPuVAGZW4KEfdVDgO++8I+7u7lKiRAnt/9XqnI4dO8rGjRulShXGaAEYdz5Ki1lb5Pxl63tReoSXlzl9mxBQACPPUbl+Pe+hWC+88IIWXDp37iwfffSRresGADbZwG3k8hjJsnKUxtfTXd7q3Ui6MhcFMHZQUXNSdu/eLeHh4XmuHzdunKhd+B9//HF71A8ArDZzXYIs2mndqh41sKMmys7t25heFMAZgora0n7btm0ydOjQW257/vnntbCitr4HAGcOKT6ebjKuc13p37KadtoxAMcq8Ktw8ODBsmzZstveroaBTpzQt2kSANhquMfanpR/9W4sQ9pWJ6QAzjhH5YcffpCvvvpKOyVR7aOi5qYAgNEmzqo9UqwxqHVV5qIAzhpU1A6yjz32mPj6+oqnp6e2uZu6jBkzxr41BIACHCi4dPevsvfXC3I1LVMuXM0odJu1r11GJncvmm3wAdghqKjzeJ566ilZuHChFlRmzJihXQgqABzZezJmxX5ZE3dWVzlVg/xk8cDmNqsXAAfMUVGnHKt5KCqk5EygvXTpkiQlJdmwOgBQ8G3wG7y8QXdI6Vi3jGx/oQPNDjh7j4raIr9kyZI3vvfx8RE/Pz9JSUmR4OBge9UPAOxy4nG98iXki2GtxM/bgxYGXGUyrdqBNjAw8Mb32dnZ2iGF8fF/TVx76KGHbFtDAMg11PP21qO6Q8qQNqEysVt92hVwtaDSv3//W6575plnbvy/m5ubZGVl2aZmAJDLhvhEmbI6Xs7p2AY/qJiXzOgZJl3DK9C2gKsFFdV7AgCO6kWZu/Wo1WWMaF9DWtUsI81Dg9hlFnDlHhUAKOpelPGr4uTStUyryyhXwluejahDQAGcFEEFgCGtjjktoz+N1V3OtJ5hhBTAiRFUABjO4KV7ZMtPv+sqw8/LXf7Vp5F0CQuxWb0AFD2CCgCXCyndwsrL25FN6EkBXABBBYBhJs3O2XxYd0h5pm2oTOjK0mPAtEFl4sSJ0q5dO2nVqpX4+/vbp1YATBVQ5m09Igu2H5e0TItVZXh5uMm4zrVlQCtOPQbE7EFl3759Mm/ePElLS5MmTZpooeX++++X1q1bS/Hixe1TSwAuuw3+2JWxcj1T3/YH8x5vzFwUwOxn/eTYsGGDXLx4UbZt2yY9e/aUmJgY6dOnjwQFBcl9991nn1oCcDkz1yXIsKj9ukJKSX8vWdivCSEFcGFWzVHx8PCQFi1aaOGkVKlSUqJECfnPf/4jx48ft30NAbjcUM+o5ftl3UHrDxN0E5FRHWtpFw939R0AV1XooLJgwQLZvn27dlHb5bdp00Yb+pk8ebKEh4fbp5YAXGYDtxc+OyApafqO2ngnsol0DWfZMWAGhQ4qw4cPlzJlyshzzz0nQ4cOlYCAAPvUDIDLhZShy/brKsPH013m9mVvFMBMCh1UvvjiC9mxY4esWLFCXn75ZWnYsKE2oVZdVO8KE2oB3Cw9M1tGRMXoahh1Xo+rbIWvhr/2nLgg5y9fl7IlfDmDCLBlUOnVq5d2UZKTk2Xnzp3y+eefaxNr1enJajUQAOTuSRn32QHJzLZu6XFJP0+Z/fdwl5kwq9pj2poESUy+fuO6kEBfmdKjvsv8GwGHT6a9cOGCNkdFrfxRl/j4eCldurQ2VwUAcn8o/3PZfrEuooh0DSsn8yLvdYlelDu1x9nk69r1C1jBBOgPKmrCbEJCgrbip23btjJkyBBt2CcsLKywRQFwUWpoY/fRJK0nxdqQ0qleWXm3X1NxpTZRPSn5tYe6TkUxdXtE/fIuE8wAhwSVp59+mmAC4LYfxvO/Pibv7TguqenWr+wZ0iZUJnZzrW3w1ZyU3MM9+YUVdbu6X4sapYu0boBLBZURI0ZoX9PT0+XEiRNSo0YN8fTkyCDA7NbGnpFxqw7I9QzrN3B7uGEFee2xhuLtWei9KA1PTZy15f0Asyj0u8G1a9dk0KBB2jk/DRo0kJMnT2rXjxo1SmbPnm2POgIwuAFL9siIFTG6Qsqg1tXkX483dsmQoqjVPba8H2AWhX5HGD9+vBw4cECbROvr+9cLqlOnTvLpp5/aun4ADOxaepY0mLxevvlZ34nHEfXLyuTuDcSVNQ8N0lb33G72ibpe3a7uB+AvhR6zUVvlq0CizvVRy5Fz1K9fny30ARMZ8sk+2XbkD11l+Hq5y//9PVy6N6oork5NkFVLkNXqHvXOmXtSbc47qbqdibSAzqDy+++/S9myZW+5/urVq3mCCwDXnTA7/gc3uZZtfUgp7uMhQ9rUkBEdaprqg1ntk6KWIN+8j0p59lEBbBdUmjVrJuvWrZORI0dq3+eEk0WLFmkHFQJwXWofkJHLYyQj28PqMoKKecn3Ezq57FyUgoQVtQSZnWkBOwWVWbNmSZcuXbS9VDIzM2Xu3Lly6NAh+e6777RN4AC47qoeNWFWr1cfvse0ISWH6kViCTJQMIV+t2jZsqV8++23kpqaqi1N3rRpk5QrV04LKvfee2+hypo6darWI5P7Ur58+cJWCYCdTVsdrzuklPL3koXsvAqgkKzaAOWee+6RpUuXii2oJc5btmy58b2Hh/VdygBsf5hg8xmb5dL1TKvLqFraT159OFzuq17aVPNRANiGw3dqU5vFFbQXRR14mPvQw5SUFO1rRkaGdrGXnLLt+RhGRxuYrx1mRv8sS777c58ka1Uu6StbxrTR/j87K1Oyrd+s1nDM9Fy4HdqAdrD2uVCY142bxWIp0FEc7u7ud13Vo25X81YKM/TzxhtvSGBgoPj4+Mjf/vY3efXVV6V69eq3vf+0adNuuT4qKkrbgA6AfpnZIpN/dJPULDUybG0PiEWCvLJlSlNrT/oB4MpSU1MlMjJSkpOTJSAgwDZBZfXq1be9bffu3TJv3jxRRamdawtq/fr1WmVr164t586dkxkzZsjPP/+sTc5VpzEXpEelcuXKkpSUdNd/qB4q+W3evFkiIiLEy8tLzIg2MEc7vLbhsHzw7W+6y3nqvsoysVs9cWWu/lwoCNqAdrD2uaA+v4ODgwsUVAo89NOzZ89brlOhYsKECbJmzRp54oknZPr06VIYDz74YJ55L2p5s5qgq+a/jB079pb7q14XdbmZapSieKMoqscxMtrAddthVnSC7pBSI9hf1o+531SrelzxuVBYtAHtUNjnQmFeM1a9m5w5c0aGDBki4eHh2lBPbGysFi6qVKkiehQrVkwLLEePHtVVDoDCuXI9U97bcUJXs73du5FsHdfeVCEFgP0V6h1FddG8+OKLUrNmTW14ZuvWrVpvSlhYmE0qo4Z1fvrpJwkJCbFJeQDuvsvsiH/vl7CpG61uKj8vd23Z8UNNXH8bfABFr8BDP6+//rq89tpr2gqd5cuX5zsUVFjjxo2THj16aD0x58+f1+aoqHGr/v376y4bwJ0Dyvyvj8rbXx+VLOsPPJayxb3lu5c6sewYgOODijo12c/PT+tNUcM8t9tH5Ysvvijwg//3v/+Vxx9/XJsMW6ZMGe2gw++//16qVq1a4DIAFE50XKI8//kBuZquZ62wRYL8vWTPpAiaH4AxgsqTTz5p80MHV6xYYdPyANzZ9LUJ8uEufXNRlPqB2bJ6XAeaG4BxgsqSJUvsWxMAdjVoyQ+y9eckXWX4e7nJ9+M7yDdbrJ/TAgBOtTMtAPsb8NH38s2RP/SV0aqaTOnRwNQ7sQIoegQVwMXP6mnz2tdy7vJfGyUW1vMP1JYhbWqw7BiAQxBUABelNnDTuzfKkDahMrx9LZvVCQAKi52ZABdki5DSqV5Zmditvs3qBADWoEcFcMHhnkU79YWUQa2ryuTuttnIEQD0IKgALrKB2/fH/5DvfkmS2FOXJNvKQ4t9Pd3lrd6NpGs4u0MDMAaCCuDkNsQnyvhVcXLpWqaucrqHh8jcvo3ZZRaAoRBUACcPKUOX7ddVRo1gP1k/ph2regAYEkEFcOK5KGNXxuoqIywkQNaObmOzOgGArRFUACe0Oua0vLAqTtIyrT9RMLi4FyEFgOERVAAn89C8nRJ3OkVXGR3qBMtHA/5mszoBgL0QVAAnGuq5//WvJTGl8LvMVi7lK54e7nJf9SB5uXuY+Hl72KWOAGBrBBXACQLKkx9+L9+fuGjVzwcV85Ztz3dgNQ8Ap0RQAQxs5roE3Zu3zegZRkgB4LQIKoBBDV66V7b8dF5XGc+0DWXzNgBOjaACGHCop9+i72TPb5esLqOUv5fM7BUmXcMr2LRuAFDUCCqAix0mOKJdDXm2cx2GewC4BIIKYBDT1sTL4m9/01VGeMUAGdelrs3qBACORlABDGDgkj3y9c+/6w4pX41kl1kAroWgAjjYQ/N3Stx/rd/AzV1E5vRuJA81qWjTegGAERBUAAfJyrbIvzb/rCuk3Fu1pKx8piXzUQC4LIIK4IBVPeNXHZA1B85IhpVH9Xi6i/yrT2Pp0ZBVPQBcG0EFcLIN3Lo2KCfznriXXhQApkBQAYpomOfRBd9KzKlkXeUMblVNJvVoYLN6AYDREVQAO1NDPKNXxEi2RV85Q9qEysRu9W1VLQBwCgQVwI4GLdkjW3UuO/Z0d5O3+zZil1kApkRQAeykx7wdcvD0ZV1l1CjjL5uebcd8FACmRVAB7LCq5x8ffq8rpPh6ucnrD4fLQ00q2bRuAOBsCCqAwc7qGdWhpozuVJteFAAgqAC2W9UzanmMrDuYqKuciPplZWznOvxaAOB/6FEBdIqOOyPPrYyRa5n6ymFVDwDciqAC6DBjTYJ88K2+oZ7yAd6y44WO4q22mwUA5EFQARy49LhTvTLyQf/m/A4A4DYIKoAVnv5kv3xzJMnqtivl5ym7J3QSP28P2h8A7oCgAhTCleuZ8so+N/kj3bqQ4qZ6YlqFyqQe7DALAAVBUAEKqPvc7RKfeEVErOsFaVG9lCwdeB9zUQCgEAwze2/WrFni5uYmY8aMcXRVgFs2cKs9af3/Qor1K3qWP92SkAIAztijsnfvXnn//fclPDzc0VUB8pi2Jl4Wf/ubrlaZ82hD6dWUHWYBwCl7VK5cuSJPPPGELFq0SEqVKuXo6gA3tHlti+6QojZwI6QAgBP3qAwfPly6desmnTp1khkzZtzxvmlpadolR0pKivY1IyNDu9hLTtn2fAyjM1MbXEvPkhazv5GrGdm6yhnUqqqM71LH5drMTM+FO6EdaAOeC9a/Hgrz/uFmsVgs4iArVqyQmTNnakM/vr6+0q5dO2nUqJHMmTMn3/tPnTpVpk2bdsv1UVFR4u/vXwQ1hqtbeMhNfkpx/9/6HGtYpJxvtrzQ0CLs3wYA+UtNTZXIyEhJTk6WgIAAMWRQOXXqlDRt2lQ2bdokDRs21K67W1DJr0elcuXKkpSUdNd/qB4q+W3evFkiIiLEy8tLzMgMbdBkxtdyOU3fPvgDWlSRl7rWFVdmhudCQdAOtAHPBetfD+rzOzg4uEBBxWFDP/v27ZPz58/Lvffee+O6rKws2bFjh8yfP18LJB4eeZeB+vj4aJebqUYpijfMonocI3PVNmg9e4uukKL6X96JbCxdwyuIWbjqc6GwaAfagOdC4V8PhXnvcFhQ6dixoxw8eDDPdQMGDJC6devKiy++eEtIAewhOTVD2r+xRS5cs34+ip+Xu8RP6yIe7tYOFwEADBdUSpQoIWFhYXmuK1asmJQuXfqW6wFby8q2yH0zN8vvV/VNCA0LKSFrR7e1Wb0AAAZb9QMUtbWxZ2TEihjd5bzdu6E81IT9UQDANEFl27Ztjq4CXNyQj/fK5oTzusoI8MySPZO7iK+Pt83qBQBwgqAC2HOoZ2TUPt0hpX3tYOlV+izzUQCgiBBU4PKi4xLlhVUH5Epalq5y5vdtJA80KCvR0dE2qxsA4M4IKnBpM9clyKKdJ3SVUaaYl3w/MULrRTH7bqwAUNQIKnBZtjhQ8B9/qyzTH+awTABwFIIKXHI+ymMLd8v+k5d0ldOpXllCCgA4GEEFLhVQ5m09KvO/OSaZ2fpOhhjSpppM7NbAZnUDAFiHoAKXsObAGRn7aYzoPPBYmlYtKVFDWog3JwoCgCEQVOD0bLE3Ss6qnu6NKtqkTgAA2yCowOknzOoNKYG+nvLao+HSJSzEZvUCANgGQQVOOx9lxL/3yfpD53SV0z08ROb2bcwGbgBgUAQVOJ3ouDMyakWsrgmz3u4ic/o2lq7hFWxaNwCAbRFU4FS9KKOW75d1B8/qKufeKiVl5dCW9KIAgBMgqMBptsF/dmWMpGVa34viJiKD24TKxG71bVo3AID9EFRgeDPXHZJFO3/VVUbzaqVk2eD7WHYMAE6GoAKDn3i8X6Lj9Q31hFUorg31AACcD0EFhrQ29ow8uzJWMnTuMBteMUC+GtnGZvUCABQtggoMZ/DSvbLlJ317o3i6u8mbvRtKTzZwAwCnRlCBoQxc/IN8fThJVxn3Vi0pK59hVQ8AuAKCCgwhPTNbus7dLsd+T7W6DC8PN/nXY2obfPZGAQBXQVCBS6zq6RpWTuZF3sveKADgYggqcOiqnj7v7ZYff7ukq5wh7I0CAC6LoAKH2BCfKONXxcmla5lWl+Hr5S5vPdZIuoZzmCAAuCqCChyyy+ywqP26yuh2T3l5+/EmDPUAgIsjqKBIh3rmbj4ib39zTFc58x5vLD0aMmEWAMyAoIIiG+p5buUBuZqepaucdyM58RgAzISgAqcY6vH2EG2op0sY81EAwEwIKrDr3igTVsXJqpjTusppUrmkfPZPNnADADMiqMCwe6Mog1pXk8ndG9ikTgAA50NQgc0NXrpHtvz0u64yint7yuuPhrP0GABMjqACm67qGRm1X3dIGd2xpozqWJulxwAAggqMs4FboK+HvPZoQybMAgBuoEcFukXHnZFhUTG6ymADNwBAfggq0GV1zGkZ/WmsrjI4qwcAcDsEFVg9H6X3e7tln44DBYt5e8gbjzZkwiwA4LYIKnDILrNjOtaUkUyYBQDcBUEFBZZtEZn/zXGZ+/VxXa32bmQTelEAAAVCUEGBrI8/K5P2esjVLOtDir+3u7zVuxGregAABUZQwV3Nik6Q93acEBE3q1rL18tdhratISM71mJvFABAobiLAy1YsEDCw8MlICBAu7Ro0ULWr1/vyCohn6XHf4YU6+eiHJrWRcZEsIEbAMDJelQqVaoks2fPlpo1a2rfL126VHr27CkxMTHSoAHnuzh6Vc/3v/whz648YHUZc/s0kp6NK9q0XgAAc3FoUOnRo0ee72fOnKn1snz//fcEFQeKjkuUSavj5cLVdKvLiKhflpACAHCdOSpZWVny2WefydWrV7UhoPykpaVplxwpKSna14yMDO1iLzll2/MxjODClXTpNGeHXE7L1lXOoJZVZfyDdVyyvczyXLgT2oB24LnAa0Lv+0Jh3kPdLBaLRRzo4MGDWjC5fv26FC9eXKKioqRr16753nfq1Kkybdq0W65XP+Pv718EtXVdE/e6yZVMd6snzHq7W6RR6WzpU90ing6d+QQAMLrU1FSJjIyU5ORkbY6qoYNKenq6nDx5Ui5duiSrVq2SDz74QLZv3y7169cvUI9K5cqVJSkp6a7/UD1U8tu8ebNERESIl5eXuJp6UzZJpo5OlIcbhcish8NMsaLH1Z8LBUEb0A48F3hN6H1fUJ/fwcHBBQoqDh/68fb2vjGZtmnTprJ3716ZO3euvPfee7fc18fHR7vcTDVKUXxoFNXjFOWE2eYz9IUUtQ3+//VubIqQ4srPBWvQBrQDzwVeE9a+LxTm/dPhQeVmqoMnd68J7LfseHhUjOjtTnuzd0PThRQAQNFxaFB56aWX5MEHH9SGby5fviwrVqyQbdu2yYYNGxxZLRNt4Ga9kEBfmdKjPrvMAgBcN6icO3dO/vGPf0hiYqIEBgZqm7+pkKLGuGAfaw7o2cDNIv3vqyJd7qkozUOD6EkBALh2UPnwww8d+fCm8ntKmnSes00upmZaXUbbstkyqVs908/NAAAUHcPNUYHtJ8zeM2WDpGbo2xulUkkf+XuNqzarFwAABcGOFy6+w2yNl6J1h5QAX0/55rn7bVYvAAAKih4VFzVzXYIs2qlvwqzSvnawLB74N1PvxAoAcByCigsO9Yxavl/WHTyrq5xSfp6ye0In8fP2sFndAAAoLIKKi+2N8tzKWLmWqW93lEGtq8nk7pxeDQBwPIKKi7DFUM8jjSrK7EfDxZvDegAABkFQcYGhnuH/3icbDp3TVc67kY2la3gFm9ULAABbIKg4+aqeYVH7dZezsF8TdpgFABgSQcXE2+D7e7vLwald2GEWAGBYBBUn9NX+/+oKKcW83WXbuA5SJuDWk6gBADASgoqTzUcZGbVfouOtX3o8pE01mdiNFT0AAOdAUHGipcejVsRKZrb1S4/nPNpQejWtZNN6AQBgTwQVk8xH6VSvLCEFAOB0OOvHCVb26A0pEfXLygf9m9msTgAAFBV6VAw+J2XS6nirfz7Qx12+n9iZbfABAE6LoGJA19Kz5NXoBIk9dUkuXE23qozSxbxk3+TONq8bAABFiaBiMEM+3iubE87rKqNykJ/sfKGDzeoEAICjMEfFxULKoJZVCSkAAJdBj4pB5qLsOvK7rpBSu1wxWTuyLQcKAgBcCkHFwTbEJ8q0NQmSmHzdqp/38nCTfz3WULo3qmjzugEA4GgEFQeHlH8u2y/WbOFW3MdDBreuLiM71uKsHgCAyyKoOHC4R/WkFDaktK5ZWoa3ryXNQ4MIKAAAl0dQcZA9Jy5YNdyz6Mlm7IsCADANVv04yPnL163aYdbP28Mu9QEAwIgIKg5StoRvoUOK6k0BAMBMGPopgrkoaphH9aCocJIzt0R9DQn0lbPJ1287T8Xbw016N60kE7s1oCcFAGBKBJUiXnqswsmUHvWlS1iI9lWt+nETyRNW1PfK24831u4HAIBZMfRj56XHN0+YVT0o6np1uwohC/o1kfKBeYeB1PfqekIKAMDs6FEp4qXHlv/1mKjbI+qX18KI+prf8BAAAGZHUHHA0mMVVtTt6n4tapTWQon6CgAA8mLox4FLj61ZogwAgJkQVBy49LiwS5QBADAbhn5svOxYudvSY7f/TZhV9wMAALdHULFCdFyiTFodLxeupue77FgFlrstPVa3M2EWAIA7Y+inkGZFJ8iwqP15Qor8b3JszrJjhaXHAADoR49KIaw/eFbe23Hitrdbci07Vr0lLD0GAEAfgkoBZVtEXln7013vl3vZscLSYwAArMfQTwEdT3GTi6kZBbovy44BAHCBoDJr1ixp1qyZlChRQsqWLSu9evWSw4cPixGlFCyjaFh2DACACwSV7du3y/Dhw+X777+XzZs3S2ZmpnTu3FmuXr3q8KXH3x3/Q1bHnta+qu8DvAr2s6WLebPsGAAAV5ijsmHDhjzfL168WOtZ2bdvn7Rt2/aW+6elpWmXHCkpKdrXjIwM7WILGw+dkxnRP8vZlL8ep1yAjzxY7s+v53Jdn58p3etKdlamZGeJS8lpX1u1s7OiHWgDngu8Hnhf0P/eWJjPEjeLxZLfnmQOcezYMalVq5YcPHhQwsLCbrl96tSpMm3atFuuj4qKEn9/f92Pf+APN/noSE4nU+5DAf9sog4h2fJ1Yn63/3kfdXvPaoZpTgAADCk1NVUiIyMlOTlZAgICnCOoqGr07NlTLl68KDt37sz3Pvn1qFSuXFmSkpLu+g+9GzW80+7NHXl6Um6qoZQP8JWXHqwjr64/nOd+QcW8ZGq3evLgPeXFVan0q4bnIiIixMurgONgLoh2oA14LvB64H1B/3uj+vwODg4uUFAxzPLkESNGSFxcnOzateu29/Hx8dEuN1ONovfD88fjf9whpChu2u1lAvzl2/Edb7t9vquzRVu7AtqBNuC5wOuB9wXr3xsL8zliiKAycuRI+eqrr2THjh1SqVIlw594zN4oAAAUDU9HD/eokPLll1/Ktm3bJDQ01GF14cRjAACMx6FBRS1NVhNhV69ere2lcvbsWe36wMBA8fPzK9K63O3EYzVHRd3OiccAAJhkH5UFCxZoE2natWsnISEhNy6ffvppkdcl58Rj5ebZJjnfT3ywrmnmogAAYAQOH/oxkpwTj9XBgurMnhzlA9U+KqnyQINyDq0fAABmY4jJtEaS34nHjSuVkI0b1ju6agAAmA5BJR83r+ox+26sAAA4CqcnAwAAwyKoAAAAwyKoAAAAwyKoAAAAwyKoAAAAwyKoAAAAwyKoAAAAwyKoAAAAwyKoAAAAw3LqnWlzzgpKSUmx6+OonWlTU1O1x/Hy8hIzog1oB54LvCZ4X+D90VafETmf2wU588+pg8rly5e1r5UrV3Z0VQAAgBWf44GBgXe8j5vFaEcYF0J2dracOXNGSpQoIW5ubnZ7HJX8VBg6deqUBAQEiBnRBrQDzwVeE7wv8P5oq88IFT1USKlQoYK4u7u7bo+K+sdVqlSpyB5PNb5Zg0oO2oB24LnAa4L3Bd4fbfEZcbeelBxMpgUAAIZFUAEAAIZFUCkAHx8fmTJlivbVrGgD2oHnAq8J3hd4f3TEZ4RTT6YFAACujR4VAABgWAQVAABgWAQVAABgWAQVAABgWASVO5g1a5Y0a9ZM2/m2bNmy0qtXLzl8+LCYyYIFCyQ8PPzGJj4tWrSQ9evXi9mfF2on5DFjxoiZTJ06Vft3576UL19ezOb06dPSr18/KV26tPj7+0ujRo1k3759YibVqlW75bmgLsOHDxezyMzMlEmTJkloaKj4+flJ9erV5ZVXXtF2TDeby5cva++HVatW1dqiZcuWsnfvXpuV79Q709rb9u3btReeCivqSTlx4kTp3LmzJCQkSLFixcQM1M6/s2fPlpo1a2rfL126VHr27CkxMTHSoEEDMRv14nv//fe18GZG6ne+ZcuWG997eHiImVy8eFFatWol7du31wK7+gPm+PHjUrJkSTHb6yArK+vG9/Hx8RIRESGPPfaYmMVrr70mCxcu1N4T1evixx9/lAEDBmi7rY4ePVrMZPDgwdpz4JNPPtG2xF+2bJl06tRJ+6ysWLGi/gdQy5NRMOfPn1dLuS3bt283dZOVKlXK8sEHH1jM5vLly5ZatWpZNm/ebLn//vsto0ePtpjJlClTLA0bNrSY2Ysvvmhp3bq1o6thOOq1UKNGDUt2drbFLLp162YZOHBgnuseeeQRS79+/SxmkpqaavHw8LCsXbs2z/XqvWLixIk2eQyGfgohOTlZ+xoUFCRmpP6CWrFihVy9elUbAjIb1bvWrVs37S8Fszp69Kj2F5Pq7u7bt6/88ssvYiZfffWVNG3aVOs5UL0pjRs3lkWLFomZpaena39BDxw40K6HwxpN69atZevWrXLkyBHt+wMHDsiuXbuka9euYiaZmZnaZ4Ovr2+e69UQkGoPm7BJ3DEB9ZdCjx49TPnXVFxcnKVYsWJaag4MDLSsW7fOYjbLly+3hIWFWa5du6Z9b8YelejoaMvnn3+uPR9yepXKlStnSUpKspiFj4+PdpkwYYJl//79loULF1p8fX0tS5cutZjVp59+qr03nD592mK2z4Tx48db3NzcLJ6entrXV1991WJGLVq00N4P1HMgMzPT8sknn2jtUbt2bZuUT1ApoGHDhlmqVq1qOXXqlMVs0tLSLEePHrXs3btXe2EGBwdbDh06ZDGLkydPWsqWLWuJjY29cZ0Zg8rNrly5ogWVN99802IWXl5e2ptybiNHjrTcd999FrPq3LmzpXv37hYz/vFSqVIl7asK7x9//LElKCjIsmTJEovZHDt2zNK2bVttaoQKrc2aNbM88cQTlnr16tmkfIJKAYwYMUJ7Qv7yyy82aXRn17FjR8vTTz9tMYsvv/zyxgsw56K+V38xqP9Xf0GYVadOnSxDhw61mEWVKlUsgwYNynPdu+++a6lQoYLFjH799VeLu7u75T//+Y/FbNRnwvz58/NcN336dEudOnUsZv7j5cyZM9r/9+7d29K1a1eblMuqnzsPi8nIkSPlyy+/lG3btmnj8vizXdLS0kzTFB07dpSDBw/muU7N7q9bt668+OKLplv5kkM9B3766Sdp06aNmIVa8XPzFgVqjoJalmlGixcv1ubqqLlbZpOamiru7nmnear3AjMuT86hVsOqi1odt3HjRnn99dfFFggqd5k8GRUVJatXr9b2Ujl79qx2vVp+piYKmcFLL70kDz74oFSuXFlbK68m06rQtmHDBjEL9bsPCwvLc516Map9NG6+3pWNGzdOevToIVWqVJHz58/LjBkzJCUlRfr37y9m8eyzz2p7RLz66qvSu3dv2bNnj7ZcXV3MRn0gq6Cifv+enub7KFGvhZkzZ2qvB7U8WW3Z8NZbb2mTis1GhRL1B2ydOnXk2LFj8vzzz2v/r/6gswmb9Mu4KNU8+V0WL15sMQu1/E7NzfH29raUKVNGG/bZtGmTxezMOEelT58+lpCQEG2ehhrqUEsxzTRXKceaNWu0idVqUm3dunUt77//vsWMNm7cqL0fHj582GJGKSkp2nuAGg5UE6qrV6+uLcdVc/rMOKG6evXq2udE+fLlLcOHD7dcunTJZuW7qf/YJvIAAADYFvuoAAAAwyKoAAAAwyKoAAAAwyKoAAAAwyKoAAAAwyKoAAAAwyKoAAAAwyKoAAAAwyKoALCJX3/9Vdzc3CQ2NrbAP7NkyRIpWbKkw+sBwLgIKgBuOHXqlAwaNEgqVKgg3t7e2mF7o0ePlj/++OOuraTOg0pMTCzU+Ud9+vTRDvUrau3atdPCjLr4+PhI7dq1tfN7srKytNvVeVY9e/aUkJAQ7VynRo0ayb///e8irycAggqA//nll1+kadOmWnBYvny5drjYwoULZevWrdKiRQu5cOHCbdsqPT1dOzm2fPnyhTqgTh3uqU7fdYQhQ4ZowUqdhjxq1CiZNGmS/N///Z922+7duyU8PFxWrVolcXFx2kFzTz75pKxZs8YhdQVMzWanBgFwal26dLFUqlTJkpqamuf6xMREi7+/v2Xo0KE3rlMHVU6fPt3Sv39/S0BAgOXJJ5+0nDhxQjukLiYm5sb9Vq9ebalZs6Z2aFu7du0sS5Ys0e5z8eJF7XZ1wGdgYOCN+0+ZMsXSsGFDy8cff6w9hipbHYaoDoDLsX79ekurVq20nwsKCrJ069bNcuzYsRu351ePghwq2alTJ8t9991325/p2rWrZcCAAQVoSQC2xNAPAK23RB3VPmzYMK2XIzfVS/LEE0/Ip59+qh3lnuONN97Qhnn27dsnkydPzneuyKOPPiq9evXS5os888wzMnHixLu29vHjx+U///mPrF27Vrts375dZs+efeP2q1evytixY2Xv3r1ab4+7u7s8/PDDkp2dres3qf7dGRkZt709OTlZgoKCdD0GgMIreB8tAJd19OhRLYTUq1cv39vV9RcvXpTff//9xlBNhw4dZNy4cXmCSW5q2KhOnTpaoFHU/8fHx8vMmTPvWBcVONQk2xIlSmjf/+Mf/9ACSc7P/f3vf89z/w8//FCrU0JCQqHmx+R+vE2bNmlBbcyYMfne5/PPP9eC0XvvvVfo8gHoQ48KgLvK6UlRk09zqPksd6LmfjRr1izPdc2bN7/rY1WrVu1GSFHUhNbz58/n6XGJjIyU6tWrS0BAgISGhmrXnzx5slC/yXfffVeKFy8uvr6+8tBDD0m/fv1kypQpt9xPTax96qmnZNGiRdKgQYNCPQYA/QgqAKRmzZpaCFG9Evn5+eefpVSpUhIcHHzjOrUa5m7hJnewybnubry8vPJ8r8rIPazTo0cPbRWSCg4//PCDdsmZ0FsYajhLDUmp4HPt2jWtZ8bf3z/PfdSwk3q8t956S5tMC6DoEVQASOnSpSUiIkLrZVAf2rmdPXtWW5qrlhLfHDzupG7dutpwSW4//vijrtZWAeWnn37SVuh07NjxxpCUNQIDA7WAppZVqxVL+fWkdOvWTZsf8/TTT+uqNwDrEVQAaObPny9paWnywAMPyI4dO7Q9VTZs2KAFmIoVK951bsnN1ORZ1RPz4osvakueV65cqc09UQoTeHJTvToqVL3//vva8umvv/5am1hrazkhRS1bVnNiVFhTlzst0QZgHwQVAJpatWppPR41atTQek/UV9WT0L59e/nuu+8KveJFzR1Rk1C/+OILbU+SBQsW3Fj1ozZZs+oNy91dVqxYoa00UhNnn3322RuTdW1JBarU1FSZNWuWNkcm5/LII4/Y/LEA3JmbWqN8l/sAgE2oXhm1Gkj11gBAQbA8GYDdqDkvauWPGq759ttvtd6PESNG0OIACoygAsCu+7PMmDFDm9tRpUoVee6552TChAm0OIACY+gHAAAYFpNpAQCAYRFUAACAYRFUAACAYRFUAACAYRFUAACAYRFUAACAYRFUAACAYRFUAACAGNX/AxyKwdD11dxyAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "plt.plot(dat['p2_1'], dat['p2_2'], 'o')\n", + "plt.xlabel('Original P2')\n", + "plt.ylabel('New P2')\n", + "plt.grid(True)" + ] + }, + { + "cell_type": "code", + "execution_count": 41, + "id": "378c62ba", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjoAAAGwCAYAAACgi8/jAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAKXVJREFUeJzt3Qd8FGX+x/FfSEIQJKEJEghFFCT0EpDiCSpVQcCCooionEiU5qGgqICn2EAsodlQT0/0BM4TTgSpd6CGdqgBDxRFmkjvSUjm//o9/9t9pZOQbCb7zOf9eg3Jzu7OPvvssPPNU2ZCHMdxBAAAwEKl3C4AAABAoBB0AACAtQg6AADAWgQdAABgLYIOAACwFkEHAABYi6ADAACsFSYel56eLnv27JHy5ctLSEiI28UBAAD5oKcBPH78uERHR0upUrm323g26CQkJJglJSVFfvzxR7eLAwAAzsOvv/4qNWvWzPX+EK+fGfno0aNSoUIFU1GRkZFuF0dSU1Pliy++kK5du0p4eLjbxbEKdUv9Biv2Xeo3WKUG8Jh27NgxiYmJkSNHjkhUVFSuj/Nsi46Pr7tKQ05JCTply5Y1ZSHoULfBhH2Xug1W7LvBXbfnGnbCYGQAAGAtgg4AALAWQQcAAFjLs0FHZ1zFxsZKXFyc20UBAAAB4tmgEx8fL0lJSZKYmOh2UQAAQIB4NugAAAD7EXQAAIC1CDoAAMBaBB0AAGAtgg4AALCWZ4MO08sBALCfZ4MO08sBALCfZ4MOAACwH0EHAABYK8ztAgAAAHfVGbswINuNCHXk+TbiKlp0AACAtQg6AADAWgQdAABgLc8GHc6jAwCA/TwbdDiPDgAA9vNs0AEAAPYj6AAAAGsRdAAAgLUIOgAAwFoEHQAAYC2CDgAAsBZBBwAAWIugAwAArOXZoMOZkQEAsJ9ngw5nRgYAwH6eDToAAMB+BB0AAGAtgg4AALAWQQcAAFiLoAMAAKxF0AEAANYi6AAAAGsRdAAAgLUIOgAAwFqeDTpcAgIAAPt5NuhwCQgAAOzn2aADAADsR9ABAADWIugAAABrEXQAAIC1CDoAAMBaBB0AAGAtgg4AALAWQQcAAFiLoAMAAKxF0AEAANYi6AAAAGsRdAAAgLUIOgAAwFqeDToJCQkSGxsrcXFxbhcFAAAEiGeDTnx8vCQlJUliYqLbRQEAAAHi2aADAADsR9ABAADWIugAAABrEXQAAIC1CDoAAMBaBB0AAGAtgg4AALAWQQcAAFiLoAMAAKxF0AEAANYi6AAAAGsRdAAAgLUIOgAAwFoEHQAAYC2CDgAAsBZBBwAAWIugAwAArEXQAQAA1iLoAAAAa3k26CQkJEhsbKzExcW5XRQAABAgng068fHxkpSUJImJiW4XBQAABIhngw4AALAfQQcAAFiLoAMAAKxF0AEAANYi6AAAAGsRdAAAgLUIOgAAwFoEHQAAYC2CDgAAsBZBBwAAWIugAwAArEXQAQAA1iLoAAAAaxF0AACAtQg6AADAWgQdAABgLYIOAACwFkEHAABYi6ADAACsRdABAADWIugAAABrEXQAAIC1CDoAAMBaBB0AAGAtgg4AALAWQQcAAFiLoAMAAKxF0AEAANYi6AAAAGsFfdD59ddfpVOnThIbGytNmzaVjz/+2O0iAQCAEiJMglxYWJhMmzZNmjdvLvv375eWLVtKz549pVy5cm4XDQAAuCzog0716tXNoqpWrSqVKlWSQ4cOEXQAAID7XVerVq2SXr16SXR0tISEhMiCBQuyPWb69OlSt25dKVOmjLRq1UpWr16d47bWrVsn6enpEhMTUwwlBwAAJZ3rLTonT56UZs2ayeDBg+XGG2/Mdv/cuXNl5MiRJux06NBBZs2aJT169JCkpCSpVauW/3EHDx6UO++8U9544408Xy85OdksPseOHTM/U1NTzeI2XxlKQllsQ91Sv8GKfZf6DbSIUCcw2y3lBOyYlt9thjiOE5h3dx60RWf+/PnSp08f/7q2bduacTczZszwr2vYsKF5zOTJk81tDS5dunSRIUOGyMCBA/N8jQkTJsjEiROzrf/ggw+kbNmyRfp+AABAYJw6dUoGDBggR48elcjIyJLbopOXlJQUWb9+vYwdOzbT+q5du8qaNWvM75rT7rrrLrn66qvPGXLUuHHjZPTo0ZladLSrS7eZV0UVF02oS5YsMcEtPDzc7eJYhbqlfoMV+y71G2iNJywOWIvOU63TA3JM8/XInEuJDjoHDhyQtLQ0qVatWqb1envfvn3m93//+9+me0unlvvG97z33nvSpEmTHLcZERFhlqz0AyhJwaKklccm1C31G6zYd6nfQElOCwm6fTe/2yvRQSdjl1ZG2orjW9exY0czABkAAKDEzbrKS5UqVSQ0NNTfeuOj58vJ2soDAAAQVEGndOnSZjq5jlnJSG+3b9/etXIBAIDg4HrX1YkTJ2T79u3+2zt27JBNmzaZE//p9HEdOKyDjFu3bi3t2rWT2bNny86dO2Xo0KGFet2EhASz6BggAABgJ9eDjp7kr3Pnzv7bvhlRgwYNkjlz5kj//v3NOXImTZoke/fulcaNG8uiRYukdu3ahXrd+Ph4s+io7aioqEK/DwAAUPK4HnT0gpznOpXPsGHDzAIAAGDNGB0AAIDCIOgAAABrEXQAAIC1PBt0dMZVbGysxMXFuV0UAAAQIJ4NOjrjSq+AnpiY6HZRAABAgHg26AAAAPsRdAAAgLUIOgAAwFoEHQAAYC2CDgAAsJZngw7TywEAsJ9ngw7TywEAsJ9ngw4AALAfQQcAAFiLoAMAAKxF0AEAANYi6AAAAGsRdAAAgLU8G3Q4jw4AAPbzbNDhPDoAANjPs0EHAADYj6ADAACsRdABAADWIugAAABrEXQAAIC1CDoAAMBaBB0AAGAtzwYdThgIAID9PBt0OGEgAAD282zQAQAA9iPoAAAAaxF0AACAtQg6AADAWgQdAABgLYIOAACwFkEHAABYi6ADAACsRdABAADWIugAAABreTbocK0rAADs59mgw7WuAACwn2eDDgAAsB9BBwAAWIugAwAArHVeQeeSSy6RgwcPZlt/5MgRcx8AAEDQBp2ff/5Z0tLSsq1PTk6W3bt3F0W5AAAACi2sIA/+9NNP/b8vXrxYoqKi/Lc1+Hz55ZdSp06dwpcKAACguINOnz59zM+QkBAZNGhQpvvCw8NNyJkyZUpRlAsAAKB4g056err5WbduXUlMTJQqVaoUvgQAAAAlIej47Nixo+hLAgAAUBKCjtLxOLrs37/f39Lj89ZbbxVF2QAAAIo/6EycOFEmTZokrVu3lurVq5sxOwAAAFYEnZkzZ8qcOXNk4MCBRV8iAAAAN8+jk5KSIu3bty+qMgAAAJScoHPvvffKBx98IMEsISFBYmNjJS4uzu2iAACAktR1debMGZk9e7YsXbpUmjZtas6hk9HUqVOlpIuPjzfLsWPHMp34EAAAeDzobN68WZo3b25+/+677zLdx8BkAAAQ1EFn+fLlRV8SAACAkjBGBwAAwNoWnc6dO+fZRbVs2bLClAkAAMC9oOMbn+OTmpoqmzZtMuN1sl7sEwAAIKiCzksvvZTj+gkTJsiJEycKWyYAAICSN0bnjjvu4DpXAADAzqCzdu1aKVOmTFFuEgAAoHi7rvr165fptuM4snfvXlm3bp08/vjj518aAAAAt4NO1jMJlypVSho0aGCuaN61a9eiKhsAAEDxB5233367cK8KAABQUoOOz/r162XLli3mnDp6gcwWLVoUXckAAADcCDr79++XW2+9VVasWCEVKlQwY3SOHj1qTiT44YcfykUXXVTYcgEAALgz6+rBBx80V/3+/vvv5dChQ3L48GFzskBdN3z48MKXCgAAwK0Wnc8//1yWLl0qDRs29K/TrquEhAQGIwMAgOBu0UlPT5fw8PBs63Wd3gcAABC0Qefqq6+WESNGyJ49e/zrdu/eLaNGjZJrrrmmKMsHAABQvEHntddek+PHj0udOnWkXr16cumll0rdunXNuldfffX8SwMAAOD2GJ2YmBjZsGGDLFmyRLZu3WpmXekYnWuvvbYoywYAAFB8LTrLli0zgUZnV6kuXbqYGVg60youLk4aNWokq1evlmCgA6f1vWi5AQCAnQoUdKZNmyZDhgyRyMjIHC8Lcd9998nUqVMlGMTHx0tSUpIkJia6XRQAAFASgs5//vMf6d69e67363Wu9GzJAAAAQRd0fvvttxynlfuEhYXJ77//XhTlAgAAKN6gU6NGDfn2229zvX/z5s1SvXr1wpcKAACguINOz5495YknnpAzZ85ku+/06dPy5JNPyvXXX18U5QIAACje6eXjx4+XefPmSf369eWBBx6QBg0amCuX6xXMdRZTWlqaPPbYY4UvFQAAQHEHnWrVqsmaNWvk/vvvl3Hjxpnz5ygNO926dZPp06ebxwAAAATlCQNr164tixYtMlcs3759uwk7l112mVSsWDEwJQQAACjOMyMrDTacbA8AAFh3rSsAAIBgQNABAADWIugAAABrEXQAAIC1CDoAAMBaBB0AAGAtgg4AALAWQQcAAFiLoAMAAKxF0AEAANYi6AAAAGsRdAAAgLUIOgAAwFoEHQAAYC2CDgAAsBZBBwAAWIugAwAArEXQAQAA1iLoAAAAaxF0AACAtawIOn379pWKFSvKTTfd5HZRAABACWJF0Bk+fLi8++67bhcDAACUMFYEnc6dO0v58uXdLgYAAChhXA86q1atkl69ekl0dLSEhITIggULsj1m+vTpUrduXSlTpoy0atVKVq9e7UpZAQBAcHE96Jw8eVKaNWsmr732Wo73z507V0aOHCmPPfaYbNy4Ua688krp0aOH7Ny5s9jLCgAAgkuY2wXQ0KJLbqZOnSr33HOP3Hvvveb2tGnTZPHixTJjxgyZPHlygV8vOTnZLD7Hjh0zP1NTU83iNl8ZSkJZbEPdUr/Bin2X+g20iFAnMNst5QTsmJbfbboedPKSkpIi69evl7Fjx2Za37VrV1mzZs15bVPD0cSJE7Ot/+KLL6Rs2bJSUixZssTtIliLuqV+gxX7LvUbKM+3Cb66PXXqVPAHnQMHDkhaWppUq1Yt03q9vW/fPv/tbt26yYYNG0w3WM2aNWX+/PkSFxeX4zbHjRsno0ePztSiExMTY8JTZGSkuE0Tqu4QXbp0kfDwcLeLYxXqlvoNVuy71G+gNZ6wOGAtOk+1Tg/IMc3XIxPUQcdHByln5DhOpnXalZVfERERZslKP4CSFCxKWnlsQt1Sv8GKfZf6DZTktMzH2WDYd/O7PdcHI+elSpUqEhoamqn1Ru3fvz9bKw8AAEBQBZ3SpUub6eRZ+/b0dvv27V0rFwAACA6ud12dOHFCtm/f7r+9Y8cO2bRpk1SqVElq1aplxtMMHDhQWrduLe3atZPZs2ebqeVDhw4t1OsmJCSYRccAAQAAO7kedNatW2fObOzjGyg8aNAgmTNnjvTv318OHjwokyZNkr1790rjxo1l0aJFUrt27UK9bnx8vFl0MFNUVFSh3wcAACh5XA86nTp1MoOL8zJs2DCzAAAAWDNGBwAAoDAIOgAAwFoEHQAAYC3PBh2dcRUbG5vrGZQBAEDw82zQ0RlXSUlJkpiY6HZRAABAgHg26AAAAPsRdAAAgLUIOgAAwFoEHQAAYC2CDgAAsJZngw7TywEAsJ9ngw7TywEAsJ9ngw4AALAfQQcAAFiLoAMAAKxF0AEAANYi6AAAAGsRdAAAgLU8G3Q4jw4AAPbzbNDhPDoAANjPs0EHAADYj6ADAACsRdABAADWIugAAABrEXQAAIC1CDoAAMBaBB0AAGAtgg4AALCWZ4MOZ0YGAMB+ng06nBkZAAD7eTboAAAA+xF0AACAtQg6AADAWgQdAABgLYIOAACwFkEHAABYi6ADAACsRdABAADWIugAAABreTbocAkIAADs59mgwyUgAACwn2eDDgAAsB9BBwAAWIugAwAArEXQAQAA1iLoAAAAaxF0AACAtQg6AADAWgQdAABgLYIOAACwFkEHAABYi6ADAACsRdABAADWIugAAABrhYlHJSQkmCUtLc3togAAXFBn7ELzMyLUkefbiDSesFiS00IKvd2fn72uCEqHouLZFp34+HhJSkqSxMREt4sCAAACxLNBBwAA2I+gAwAArEXQAQAA1iLoAAAAaxF0AACAtQg6AADAWgQdAABgLYIOAACwFkEHAABYi6ADAACsRdABAADWIugAAABrEXQAAIC1CDoAAMBaBB0AAGAtgg4AALAWQQcAAFiLoAMAAKxF0AEAANbybNBJSEiQ2NhYiYuLc7soAAAgQDwbdOLj4yUpKUkSExPdLgoAAAgQzwYdAABgP4IOAACwFkEHAABYi6ADAACsRdABAADWIugAAABrEXQAAIC1CDoAAMBaBB0AAGAtgg4AALAWQQcAAFiLoAMAAKxF0AEAANYi6AAAAGsRdAAAgLUIOgAAwFoEHQAAYC2CDgAAsBZBBwAAWIugAwAArEXQAQAA1iLoAAAAaxF0AACAtQg6AADAWgQdAABgLYIOAACwFkEHAABYi6ADAACsRdABAADWsiLofPbZZ9KgQQO57LLL5I033nC7OAAAoIQIkyB39uxZGT16tCxfvlwiIyOlZcuW0q9fP6lUqZLbRQMAAC4L+hadb775Rho1aiQ1atSQ8uXLS8+ePWXx4sVuFwsAAJQArgedVatWSa9evSQ6OlpCQkJkwYIF2R4zffp0qVu3rpQpU0ZatWolq1ev9t+3Z88eE3J8atasKbt37y628gMAgJLL9a6rkydPSrNmzWTw4MFy4403Zrt/7ty5MnLkSBN2OnToILNmzZIePXpIUlKS1KpVSxzHyfYcDUy5SU5ONovPsWPHzM/U1FSzuM1XhpJQFttQt9RvsGLfDYyI0P8/fkSUyvyzsILx+zvif3VR5Nv9X50Gok7yu80QJ6ek4BINKPPnz5c+ffr417Vt29aMu5kxY4Z/XcOGDc1jJk+eLGvWrJEXXnjBPE+NGDHCPGfAgAE5vsaECRNk4sSJ2dZ/8MEHUrZs2YC8LwAAULROnTpljvVHjx41Y3SDMuikpKSY8PHxxx9L3759/Y/TMLNp0yZZuXKlGYyswWfFihX+wchfffWVVK5cOd8tOjExMXLgwIE8K+p8NJ6w+LzS71Ot0+XxdaUkOT33lqnvJnQrZOm8R9P/kiVLpEuXLhIeHu52cazjtfo9n//f+ZX1/3cw1G2g6iOQ33W+Muf3excF56vbQOy7evyuUqXKOYOO611XedHwkZaWJtWqVcu0Xm/v27fP/B4WFiZTpkyRzp07S3p6ujz88MO5hhwVERFhlqz0AyjqDyE57fz/w+h/tryeX1K/7IJBID5reK9+C/P/+1xyq7+SXLeBqo9Avt+sZT7X9y7OXyD23fxur0QHndzG3GgjVMZ1vXv3NgsAAECJmnWVF22SCg0N9bfe+Ozfvz9bKw8AAEBQBZ3SpUub6eTaN52R3m7fvn2htp2QkCCxsbESFxdXyFICAICSyvWuqxMnTsj27dv9t3fs2GEGGuuZjXX6uJ71eODAgdK6dWtp166dzJ49W3bu3ClDhw4t1OvGx8ebRQczRUVFFcE7AQAAJY3rQWfdunVmILGPBhs1aNAgmTNnjvTv318OHjwokyZNkr1790rjxo1l0aJFUrt2bRdLDQAAgoHrQadTp045nvQvo2HDhpkFAADAmjE6AAAAhUHQAQAA1vJs0GHWFQAA9vNs0NEZV3ph0MTERLeLAgAAAsSzQQcAANiPoAMAAKxF0AEAANZy/Tw6bvOdw0fPkFzU0pNPFfg5aaGOnDqVJmnJoZKex1V0A1Fe26WmpsqpU6dM3ZXUK0AHM6/V7/n8/86vrP+/g6FuA1Ufgfyu85U5v9+7KDhf3QZi3/XtG+c6F1+Ic65HWG7Xrl0SExPjdjEAAMB5+PXXX6VmzZq53u/5oJOeni579uyR8uXLS0iI+0leE6oGL/3gIiMj3S6OVahb6jdYse9Sv8HqWACPadpOc/z4cYmOjpZSpXIfieP5riutnLySoFt0hyDoULfBiH2Xug1W7LvBV7f5uSg3g5EBAIC1CDoAAMBaBJ0SJiIiQp588knzE9RtMGHfpW6DFfuu3XXr+cHIAADAXrToAAAAaxF0AACAtQg6AADAWgQdAABgLYJOCfD0009L+/btpWzZslKhQoV8nxFywoQJ5oyQF1xwgXTq1Em+//77gJc12Bw+fFgGDhxoTiqli/5+5MiRPJ9z4sQJeeCBB8yJJLVuGzZsKDNmzCi2Mttct2rLli3Su3dv8xw9I/kVV1whO3fuLJYye6F+fe677z5ztvdp06YFtJxeqFu91tgjjzwiTZo0kXLlypnv3TvvvNOcVR8i06dPl7p160qZMmWkVatWsnr16jyrZeXKleZx+vhLLrlEZs6cGdhq1GtdwV1PPPGEM3XqVGf06NFOVFRUvp7z7LPPOuXLl3c++eQT59tvv3X69+/vVK9e3Tl27FjAyxtMunfv7jRu3NhZs2aNWfT366+/Ps/n3HvvvU69evWc5cuXOzt27HBmzZrlhIaGOgsWLCi2cttat9u3b3cqVarkjBkzxtmwYYPz448/Op999pnz22+/FVu5ba5fn/nz5zvNmjVzoqOjnZdeeingZbW9bo8cOeJce+21zty5c52tW7c6a9euddq2beu0atXK8boPP/zQCQ8Pd15//XUnKSnJGTFihFOuXDnnl19+yfHxP/30k1O2bFnzOH28Pk+f/7e//S1gZSTolCBvv/12voJOenq6c/HFF5uw43PmzBnz3JkzZwa4lMFD/xNplv/qq6/86/QLStfpl1VuGjVq5EyaNCnTupYtWzrjx48PaHm9ULcayO+4445iKqX36lft2rXLqVGjhvPdd985tWvXJugUYd1m9M0335jn5HZA94o2bdo4Q4cOzbTu8ssvd8aOHZvj4x9++GFzf0b33Xefc8UVVwSsjHRdBaEdO3bIvn37pGvXrv51ejKmq666StasWeNq2UqStWvXmmbptm3b+tdpN4muy6ueOnbsKJ9++qns3r3bdBEuX75c/vvf/0q3bt2KqeR21q1eQHfhwoVSv359U5dVq1Y1z1+wYEExltzufVfrWLthxowZI40aNSqm0nqjbrM6evSo6RrM73ADG6WkpMj69eszHYuU3s6tLrX+sz5evw/WrVtnuggDgaAThDTkqGrVqmVar7d99+H/60kPplnpurzq6ZVXXpHY2FgzRqd06dLSvXt30wetAQjnX7f79+8345+effZZU6dffPGF9O3bV/r162f67FH4ffe5556TsLAwGT58ONVZxHWb0ZkzZ2Ts2LEyYMAAT198+cCBA5KWllagY5Guz+nxZ8+eNdsLBIJOgOhAYU37eS2aYAtDt5GRtj5kXef1us2pPs5VTxp0vvrqK9Oqo3+tTJkyRYYNGyZLly4V2wWybrW1Qd1www0yatQoad68uTlYXH/99YEfjOiB+tV99eWXX5Y5c+Z44nuguL8XfLTV4dZbbzX7s/4BBCnwsSinx+f2uRSFMD6kwNBZO/qfIS916tQ5r21ffPHF/mRcvXr1TH8xZ03KXq7bzZs3y2+//Zbtvt9//z3Xejp9+rQ8+uijMn/+fLnuuuvMuqZNm8qmTZvkxRdflGuvvVZsFsi6rVKlimlt0NayjHRW27/+9S/xgkDWr8500e+AWrVq+dfpX9sPPfSQmXn1888/i80CWbcZQ84tt9xihg8sW7bM0605vv/ToaGh2Vpv8joW6fErp8frd0PlypUlEAg6AdwBdAkEncanO8uSJUukRYsW/r5Sbf7Xpmvb5bdu27VrZ/rRv/nmG2nTpo1Z9/XXX5t1Op0/ty8yXUqVytzYqf+ZfS0SNgtk3Wo3YFxcnPzwww+Z1uv4p9q1a4sXBLJ+dWxO1iCuYx90/eDBg8V2gazbjCFn27ZtZtxeoA7KwaR06dJmmrgei7Qb2kdva8ttbvX/j3/8I9M67cZu3bq1hIeHB6agARvmjHzTUfsbN250Jk6c6Fx44YXmd12OHz/uf0yDBg2cefPm+W/rjCudZaXrdHr5bbfdxvTyXKaRNm3a1Myq0KVJkybZppFmrdurrrrKzLzS6eU6FVJnw5UpU8aZPn06e3Uh61Z/16mks2fPdrZt2+a8+uqrZur+6tWrqdsi2HezYtZV0dRtamqq07t3b6dmzZrOpk2bnL179/qX5ORkT++7H/5vevmbb75pZrSNHDnSTC//+eefzf06+2rgwIHZppePGjXKPF6fx/RyDxg0aJCZpph10QOtj97WA27GKeZPPvmkmWYeERHh/OEPfzCBB5kdPHjQuf322805h3TR3w8fPpzpMVnrVr+87rrrLnMOEg04+oU3ZcoUU+coXN0q/WK79NJLTd3quV44P1HR1m9GBJ2iqVs9n1ZO39FZv6e9KiEhwexrpUuXNqfiWLlyZabjm/7xmNGKFSucFi1amMfXqVPHmTFjRkDLF6L/BKatCAAAwF3MugIAANYi6AAAAGsRdAAAgLUIOgAAwFoEHQAAYC2CDgAAsBZBBwAAWIugAwAArEXQAVyiV+pdsGBBUNa/XiBSy68XOy2p9BpPzzzzjHiFXrW8QoUKBXrOXXfdJX369AloOV577TXp3bt3kb4GUBAEHSAA9Oq8Dz74oFxyySUSEREhMTEx0qtXL/nyyy+tqG99P3v37pXGjRsXajsalnxL+fLlzYX95s2b57//9ddflyuvvFIqVqxoFr1opV6M8Vz0CtULFy40n4FX9O/f31wgtSBefvllE0wCaciQIZKYmOiZK9Sj5CHoAAFo7dAr+i5btkyef/55+fbbb+Xzzz+Xzp07S3x8vBX1rVdzv/jiiyUsLKzQ23r77bdNaNKDYbNmzeTmm2+WtWvXmvtWrFght912m7latK6rVauWdO3aVXbv3p3nNrUVQbej4ckL9MraF1xwgVStWrVAz4uKiipwK1BBadAfMGCAvPrqqwF9HSBXAb2SFuBBPXr0cGrUqOGcOHEi230ZLxyo//1ef/11p0+fPs4FF1xgLnT597//3X//2bNnnbvvvttc9E4vgFm/fn1n2rRpmbanF8y74YYbnBdeeMFc4LVSpUrOsGHDnJSUFP9j9uzZ4/Ts2dNsQ7f1/vvvZ7vY45EjR5whQ4Y4F110kbnIYefOnc1VmnPju8jhxo0bzW29sKHeXrp0qdOqVSvzftq1a+ds3bo1z7rS58yfP99/W8utVzbWKx7nROtEy/fOO+/kus20tDSnQoUKzmeffZZpvb7np59+2hk8eLBz4YUXOjExMc6sWbMyPWbXrl3OLbfcYp6vdalXrNb3qjZv3uyEhIQ4v//+u7l96NAhc/umm27yP/+ZZ55xrrjiCqcwtPwTJ040+5Be9FAvfPrPf/4zW93PnTvXXCxRL+r71ltvmQtQRkVFZdrWU089ZT5Tfb/33HOP88gjj5jtZd1/fHR7Dz74oDNmzBinYsWKTrVq1czFgzPSC9w2btzYfE56Ne/777/fOX78uP/+nMqhF3HU93Lq1KlC1Q1wPmjRAYrQoUOHTOuNttyUK1cu2/1Z/3qeOHGi3HLLLaarpWfPnnL77bebbaj09HSpWbOmfPTRR5KUlCRPPPGEPProo+Z2Rtra8eOPP5qf77zzjumKyNgdceedd8qePXtM68gnn3wis2fPlv3792f8Y0euu+460922aNEiWb9+vbRs2VKuueYaf1ny67HHHpMpU6bIunXrTGvP3XffXaDnh4eHm+dpC0VOTp06Ze6rVKlSrtvQujxy5IjpBstKy6brN27cKMOGDZP7779ftm7d6t+2trpdeOGFsmrVKtPVor93795dUlJSTDdd5cqVZeXKlebx+hi9rT99tI6vuuoqKQztTtJyvvjii+a9dOvWzYxx2bZtW6bHPfLIIzJ8+HDZsmWLeUxW77//vjz99NPy3HPPmc9UW8NmzJhxztfXfUj33a+//tq0SE6aNEmWLFniv79UqVLyyiuvyHfffWceqy2XDz/8cJ7b1DrXzy0/3Y5AkTuveAQgR19//bX5a3vevHnnrCF93Pjx4/23tQVIWwgy/vWelbbW3HjjjZn+IteWCm3p8Ln55pud/v37m9+3bNliXicxMdF//7Zt28w6X4vOl19+6URGRjpnzpzJ9Fr16tXL1uKRnxYdn4ULF5p1p0+fzleLjr6+tkDoukWLFuX6/rVceW1TtxcaGuqkp6dnWq/1dMcdd/hv6/1Vq1Z1ZsyYYW6/+eabToMGDTI9Lzk52bROLV682Nzu16+f88ADD5jfR44c6Tz00ENOlSpVnO+//95JTU01LSd5fX75ER0dbVqeMoqLizPvPWPdZ23dy9qS0rZtWyc+Pj7TYzp06HDOFp2OHTtme21tCcrNRx995FSuXDnXcvhoC9GcOXPyfO9AINCiAxTtHw7mpw6uzY+mTZv6f9e/onVMScbWlpkzZ5q/hi+66CLTuqCDc3fu3JlpG40aNTJjZnyqV6/u38YPP/xgWki0hcbn0ksvNQN7ffSv/RMnTpjWCX0N37Jjxw7TUlQQGd+PlkNlfD850TE4+nply5aVqVOnmpaMHj16ZHucti789a9/NYOVy5Qpk+v2Tp8+bcaF5PQZZCyf3q/jjHzl03rYvn27+Qx8daAtR2fOnPHXQ6dOnUyrjdKWHW0B+sMf/mB+1zFG+todOnTIsVw6Ayxj/Wb9HNWxY8dM61vWbehtbbnJKKcWq4z0s2/Tpk2mdVlv5yRjHWXdn5S2HHbp0kVq1Khh6kpbDA8ePCgnT57Mc7s6hkhbzYDiVviRhAD8LrvsMnMA1YNSfqbtaldNRvpc7bJS2kU1atQo043Rrl07c1B54YUXTJdCfrfhC15ZZVyvj9WDme8AnlFBB6pmLIsvaPjKkpuXXnrJzKaKjIzMdTCthh8NCkuXLs12IM6qSpUq5oCq3U2lS5fOtXy+MvrKpz91ELl2+WSlQdMXdEaMGGECkXbd6IwwDUEadLS7TJ+f2wDooUOHmm5Kn+jo6FzfQ9aQpp9X1nU5dY3mZzvnklcd/fLLL6aLVd/LU089ZYKgdvHdc889uXY3+mg3qK8egeJE0AGKkH7x63iJhIQEM34i68FID4b5DQ+rV6+W9u3bm7EkPgVtYbn88svl7NmzZkyKHoSVHqS1HD7a2qPjc7Tlp06dOlLctFVFW5lyo+Huz3/+syxevPicrRiqefPm5qeOa/L9nh9aD3PnzjVhS0NXTnzjdLQ8OkNMH6djciZPniyHDx/Oc3yO7ht5jS1Suj0NQBoetKXIZ82aNflqjcmoQYMGZkyMnk/IR8dOFYY+X/cnDd86VkdlHTOWE91vtWWsRYsWhXp94HzQdQUUsenTp0taWpo5MOngXx1Eqi08OoBTW2bySw/+emDRA7yeH+Xxxx833SMFDTraWvLHP/7RHPQ08Ojv2o3g+2tf79dyaQuUvpZOj9cD6/jx4wt9YCws7a7Scrz11lsmhGkg00W72nKjrQYaWgp63hYdCK6tQTfccIMJmdp1py012oKza9cu8xitMw0gf/nLX0zrjtIWJm090nMk+dYVxpgxY8wAYg1d2v00duxYc2JGLUdB6DmE3nzzTTNgWPdBDWc6uDm/3ao5qVevngk6OlX8p59+kvfee890r56L1qeeU0qfDxQ3gg5QxOrWrSsbNmww4zceeugh0wqgYxr0QJifWS8+2j3Qr18/cyK4tm3bmnEQGVt38uvdd9+VatWqmQN03759zQnctHvFN85FD3w620rv11lS9evXl1tvvdUEHn2e26FRQ8RNN91kutd8i3Zl5UXDXE5dUHnRMUI6g0pnJ2m9N2zY0NSHjrvJ2MKjn6sGWV+o0frTLizVsWNHKSxtCdT9RpcmTZqYWXyffvqp6RYtaHAbN26c/OlPfzLBT4Obngk5r/FN56ItZDqOSoOY7tdax9qadS46tkr3O8ANIToi2ZVXBuAKbZ3QMxvreBedQm4j7SbRrpsPP/ywQK1ottPArV2F2hJTXHQsk+5n2iqpJygEihtjdADL6XlOtKtHWwf0DMR6zhPtBso4BsQ22mqhLVkHDhwQr9IB2dqtpGPGdFaetqpouM14TpzioLPI9LMg5MAttOgAltNxN9oNomMqtMtKBzhPmzZNateu7XbREEDa5abXV9Nu1OTkZNPCpeOdtFsO8BKCDgAAsBaDkQEAgLUIOgAAwFoEHQAAYC2CDgAAsBZBBwAAWIugAwAArEXQAQAA1iLoAAAAsdX/ATX/+k66V0d9AAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "plt.hist(dat[\"p2_2\"] - dat[\"p2_1\"], bins=20, log=True);\n", + "plt.xlabel('Change in P2 (new - original)')\n", + "plt.ylabel('Count')\n", + "plt.grid(True)" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "id": "b7e50cc7", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAGwCAYAAABVdURTAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAQm1JREFUeJzt3Qd4VGXaxvEnhBAIEkogJEgv0puAgEoRacIirLifLipY1pUVG4ggWIAFBV1XEUWwoLRV3BVYcSkLrhQLKE1BmkgXEyItAUJCgPmu510nO0lmkpnJ1DP/33WNYc6cOTnzzpnM7VujbDabTQAAACyiRLBPAAAAwJcINwAAwFIINwAAwFIINwAAwFIINwAAwFIINwAAwFIINwAAwFJKSoS5fPmy/Pzzz1KuXDmJiooK9ukAAAA36LR8Z86ckWrVqkmJEoXXzURcuNFgU6NGjWCfBgAA8MKRI0ekevXqhe4TceFGa2zshRMfH+/TY+fk5MjKlSulZ8+eEhMT49NjWw1lRVlxXfEZDBf8vQqN8srIyDCVE/bv8cJEXLixN0VpsPFHuImLizPHJdxQVlxXgcdnkLLiugo+f38O3elSQodiAABgKYQbAABgKYQbAABgKYQbAABgKYQbAABgKYQbAABgKYQbAABgKYQbAABgKYQbAABgKRE3QzEAAP5y4eJlmbf+oBw6mSm1KsXJoPa15NsjpyXtTJYklist19SpJNEl8s6we+myTTbsOyHr9x/X+XelY70E6VA3Ic9+us83B04WOI79uV/tOy4/nz4v1SqWkWvrVpYO9RLM8zbsPyGf7/1Ftv+ULnGlouWaOgky5NraUqpkCTmbdVEe+3CL7Ek9K+XLxMjwGxtIqZho+frASbNIZXyZkpKRdVFsl0UysnLkstj+u3jl+YvyQ9oZiY+NkZ5NkuTu6+uY4+m51Bu77NczLiGPrl9p/nVwSt/ICjfjx4+XCRMm5NlWtWpVSU1NdfmctWvXyogRI2THjh1mZdBRo0bJ0KFDA3C2AAC4NnnZTnn78wNy2fa/bROX7sqzT3L50jKuXxPp3SzZ3F/xfYo8uWi7nM7Myd3n9dU/SoW4GJlyS3Ozn+4z4ZOdkpKelec4N7dMlg83/ZTnuWr66n0myGg0OnfhUp7HVu1Kk+eX75LKZUvJL2cv5G4/cuq83Ddvs8dv76bDp2Xyit1yY+NE+XRXmtOGodpPLg14wAl6zU3Tpk3l008/zb0fHR3tct8DBw5Inz595P7775f58+fLl19+KQ8++KBUqVJFBg4cGKAzBoDgc/V/8ghesHlz3YEi90tNz5I/zd8iM+682twfOn+L0/00sOhjD3SuI2+tOyAOecnQoFPY78vMF2oc2WySJ9gUl55b3mBTUKADTtDDTcmSJSUpKcmtfWfOnCk1a9aUqVOnmvuNGzeWTZs2yUsvvUS4ARAxXP2fvGONAALbFKU1Nu4GAY2g+v5dvny5yP31uPmDTbiqHcCAE/Rws3fvXtO8FBsbK+3bt5fnn39e6tat63Tf9evXmyXUHfXq1UtmzZplViF1tvpodna2uTkuma50f735kv14vj6uFVFWlBXXlXf+veOYPLzguwJfePYagddubym9mlblMxjAv1ezvzqYpymqKLqrYzAtjCfHDQc5xfh+9OS5UTbtHRQky5cvl8zMTLnqqqvk2LFjMmnSJNm9e7fpT5OQ8N/OUI50v7vvvlvGjh2bu+2rr76S6667Tn7++WdJTk52q1+Pev/9982S7ADgT/rltC8jSjJyRK4o+d//az9zUSQ+RqRevE08aUnSY03YEi2nTYuCsyfapEIpkXFXX/LouCiejw6UkM9TGXxctMvyaseia6tc0bwwaNAgSU9Pl/j4+NCtubnpppty/928eXPp2LGj1KtXT+bMmWM6DTsTFZX3E2vPZvm3240ZMybPsbTmpkaNGqYGqKjC8SZVrlq1Snr06OG0FgmUFdeVf4XCZ1D7wmw6dErSzmTLweOZ8uGmI3LsjPP+DUnxsfJ0n0Zu1bQoHcVyesOmQvaIMsGnSpMO0r5OpZAvq3BRVFkd++qgfL78h6CcW3gpIX369Pb62faWl7BolnJUtmxZE3K0qcoZ7ZuTfyRVWlqa6bfjrKZHaXOX3vLTC9RfH2h/HttqKCvKygrXlX3477q9v8iWw6flTNZFt553LCPbNDFp51J3+sqcyHTvuLqfu68/nD6Dwe5E7aqs7r6unkxZ8YPbTUh6xknlS5s+N66Cr52+PCs1TcUU41rz5LkhFW60b8yuXbukU6dOTh/Xmp1PPvkkz7aVK1dK27Ztw+bDCcD6w3/d5di5tEeTpCK/qPUL3R3u7hdOQrkTtc7xcn+nOm6NlrK/w3rehY2WstPj6mgpFe4Z52AAR0sFtZFw5MiRZt4aHeL99ddfy6233mqqnYYMGZLbpDR48ODc/XU+m0OHDplmJg1B7777rulMrMcBgGAN/y3O/1nbO5dqjURRtKZCv9BdRSDdro/rflYLNtpZOn8nXHsnan082Mb0aWKGbRdVkaQ1NvaaOr3NvPNqM6dNfhXjYsxjelzdX5/nSN9n/X3OnqtKlwytPkARNc/NTz/9JL///e/l+PHjZq6aDh06yIYNG6RWrVrm8ZSUFDl8+HDu/nXq1JFly5bJ8OHDZfr06WaU1bRp0xgGDiCkh/+6Q5taiqI1O/p//PqFrt+hNhc1Alaa70aborTGxuaDmi9/0yDyeM9GHs1QrAFHz93MMrz/uPx86rxUq1BGrq33v1mG7fs4a5Ib1bux0xmK085kyfC/f1fkOddOKCMHT5z3a7lE3AzFCxYsKPTx2bNnF9jWpUsX2bKl8Go8APA3/QLzZV8Id5uS9ItO/08+fxNNUog00fiafqEXNmzaseZLly0INm2iuq9T3ulMijovDSlnsnNk0Zajua91+pp9eZrddB9nx9Ht1zWobG6O1u874db5PjeghYz86DtTC+ary7ln4yrSQFLl4dt7S+nYUhIMIdXnBgDChf6fuS/YO5d60pRU2P/JW407NVqe7BfKzW6u5i5yt8O5sybMVBehxX7dae2QvTbQF+67tpY8eVND08oSzOsxtBrlACBMaJNDcRWnKcn+f/L9W11pflox2IRqJ2ptKtOakY+/PWp+6v3iHGvcx9+7bHZTWkvn6e+I/rUJU0UVcd1pcPpj5zpSXC2ujJdnbm4moYCaGwDwwl0da8tzy3YVq2nKqk1JvuRuDUSgOlH7atSWfVj7iyt2FzocvDjNbr1dNGGWj4uRe66tY2r/7Ofy8bfF65TdvXEVeWfINRIqCDcA4Mfhv/rFd3u7mlK7cpxULhtrvo2Pn822dFOSL4VSJ+qilr5wp/lIg8Rr/9kr73yxX85mu17c0hfNbpcu26R8mVIyqncj+XLvL7Jq5zFJz7poFuV85dMf5INvDsnvr6kpOZcuS2qGd816jarGyeJhnaVMKdeLXgcD4QYAijE6RuWf50YnTO/WsIr8oVM9AowPhEInan1/Jy/bXaxRW1rrM+Lv3xW6Yrevmt1WOKlhyi81I1te+dT5pLnuuPe62vJsv6YSigg3AODj4b/aZKU1O/CdYHei1vXBNAx42nxkb35atTNV3v3yoFe/u1LZGI+a3Va46KDsSz2aJIZssFGEGwDww/Bf+J6r4dCBoAufetp85E7tiTt+2+pKt0NcYfMC+cIVsdEy5bfN5TetrpRQRrgBAIStQK03pau4e9J85Mvak+6/dvz1xbxAxTG8+1XyULf6YdFPjHADAAhLgVxvql68zaziroudFjVqy5e1J54up+GP+X4qxMXIlFuah9WoPsINAIThCtXFPfevD5yUzcejJEH7iNRPDMlzz1/GbWpVlM2HTpn7B4+fc9oZtjgT3xVWVltPRMn/taku01bvK3LUls5946vaE09Hgvlqvp+2tSpIh7oJ0rHuf5eBCMXrozCEGwCw0ArVnp17tMzduykkz91ZGev3a1HzCvlyvan8ZSV79+UuVKnDqV2N2vJF7UnZ2Gj56+9amtegYcndEF3UvEDuGt69YYElHcIJ4QYAgjxVfqCEy7m7Ok93J0z0ZOI7VzVwrs4h/ddQo/1PdO4iZ4HD29qT2pXKSOuaFeWWq6vLtfUry7+/T5V2z30qJ8/9b5K/ooKofV6gocVYTqF8mZK5i3aGK8INAFhwhepwPXdf9lcpqgbFVQ3cM32byMSlhZfVgo2H5YvR3ZyWlae1JzoC6cWBLaRPi2q52yYv2+l0gkg9Vw0uOseMvlf5g5WuVn/kZKYZwaf/9sYLA1uE3PXrKcINAFh0hepwPHdfjvYprAalsFqsB98vvNajqLIqbFZluz7NqkrdKuXM87Vvi2OYWLbt5yJnvtY5c/Rmr8nRoPPoB5vlX9uPibdKRUfJtN+3Donau+Ii3ABABKxQHS7n7ovfX9R6U0XVYvniXF3NqlxUs5Ke29Mff+/2OaT+WpNTskSUXCzGQmd9myfJtN9fHfY1NnaEGwAI4xWqrXbuxf397qw35avaoaLO1ZtZlXXfk+fcnDHQIYx5G2w61qkkc+5rb7kZtQk3ABCmK1R7c+6FfaknlC0lqennzeicYA1tL+5oH3fWmypu7ZAn77OnsyoHquYsvnRJ+Xps95Bb8NJXCDcAEIYrVHtKz+nmlsmF9uU4ce6CDP/7d+bfwRoe7k5/FTt7KT9WyMil4tYOBfp99nfNWdSvP1+8tYVlg42yVj0UAPiZvS+F/p+7I70fKkOpXfXlWPJditv724eHa8fbUCnj/FnCXuaPdm8g/VtdaWpI3Akc9tohV3vqdn38jUGtA/4+28/NX5JC/Dr1FWpuACDMVqj2hqf9TIIxPDz/nDNrn7ghd0bi/DMUF6fM3a2B0/e5V7NkWf9jmqz8/Gvp2am932dzdjw3Xy9+Oe/ea8z8OaF8nfoK4QYAwmyFam9405cjkMPDC5v1WWtl7Hx1Hq5GM+Xvs6Pvc/s6leTELpv5GYhg4OrciuP+TrWl01VVJFIQbgAgAhSnL4e/O7kGa+bkUK6B03O7fNkmTyzcJueyLxXrWD2aJMpTfZtKJCHcAEAEKM4oJH92cnVnzpmxi7fL+ZzLkhTv+/ARqjVwOpHfg+9vLdYxysREm47D/Vr+b+bjSEG4AYAI4MkopEAObXenL5DO+zL8w2/Nv0NxkU9f+9e3P8vDC4oXbB67sb48fONVIVELFQyMlgKACOFqFJIzgRra7mmTVyBHcWmtks7588m2FNmbHmXu+9tzS3fKQwu2et2ZWFcTn3nn1fJYj4YRG2wUNTcAEEHs/UwcRwBlZF2SiUt3Fdqx1tcjoewjn/YeO+vRcQI1iqtgB+do+eiv62T8zU39VmukwebtzwtfU8qV6+snyAOd60XMaKiiEG4AIMLkHwEUExNjhjz7q2Ots5FQemhvK0L8PYrLVQfnYxnZfuvgrH1svA02SfGxMufe9oQaB4QbAIDfOta6Cgq+aOHxxyiuojo4+6PWyNPFMvPT2iRqa/Kizw0AwC8KCwq+4I9RXEV1cHasNfLl7/RksUy7CnExpn+NlTtXe4uaGwCAX3i7+vaDXevKgo0/yalzFwK+QKm7tUHFqTXK3/8oNcOzY/VuWlXu6lBbOri53EQkItwAAPzC2wDQMClenv9ts6AsUOpubZC3tUbO+h9VKhvj9vNfv721/KZV5M1b4ymapQCEJPsw3I+/PWp+BmIYLnz7XngbAPR5wVqg1N1FNb2pNbL3P8pfm+Vuk9T9neoQbNxEzQ2AkFPYOkP0Lwif98LTWZHzNzcFY3kEdxfV9PQcitv/SIPNU32bePnsyEPNDYCQ4ur/bgM5eRt8817Yg4IqKgq4Cg72UVy6eKb+DOTClQVrjWK9rjVyt/9RpbKl8t2PkTcGtSbYeIiaGwAhIxjDcCONvTNryulzsv/XWXdj/PheuFrhOv88N76eNLC4HGuNTFnt+FYeuq2zlI7NGz583f/omb6NJal8mZBbyDPcEG4AhAxPhuGG4mKHoc6TWXd9+V44a16yz1Acyl/i9lqjnJx4WfbTVrfPL/9oKH1t7vY/0mDDtV18hBsAISMQw3Ajlaez7vr6vXA2SaAVv8Rd9VF6pm+TQvsfBWKR0khCnxsAIcPfw3AjVVFNTEofdxwFxXvh2z5Kw97fIje3/G94jArSIqWRhHADIGT4cxhuJPNm1l3eC98HyCXfpcj0Qa0DPrw9EtEsBSBk+GsYbqTzpomJ98I/AbJi2Vj5YnS3gA5vj0TU3AAIKcGavM3KvG1i4r3wT4AMxvD2SEPNDYCQE4zJ26ysqMn0CuvMynvhHvoohRbCDYCQ5Gx0DYLT3Md74d8ACd+jWQoAIoA/Zt2Fe7Mx018s8Ki5AYAI4etZd+HebMyhNvtyJAiZmpvJkydLVFSUPPbYYy73WbNmjdkn/2337t0BPVcACFf2JqZ+LZKlQXkb/Zh8TAOMjob64P4O8urtrcxPvU+wicCam40bN8pbb70lLVq0cGv/PXv2SHx8fO79KlWq+PHsAACRxL58QmpGlpw8my3ly0TLwULW4cqPPkrBF/Rwc/bsWbnjjjvk7bfflkmTJrn1nMTERKlQoYLfzw0AEFmcLZ9Q1DpcCD1BDzfDhg2Tvn37Svfu3d0ON61bt5asrCxp0qSJPP3003LDDTe43Dc7O9vc7DIyMszPnJwcc/Ml+/F8fVwroqwoK64rPoOh5t87jsnDC75zOtpJpf66Dtdrt7eUXk2rBvjswkeOn74LPTleUMPNggULZMuWLaZZyh3Jycmm+apNmzYmsMybN09uvPFG0xenc+fOLvvyTJgwocD2lStXSlxcnPjDqlWr/HJcK6KsKCuuKz6DoUCX1ZqwJfrXYON6PiWb2OTpRd9KzsFLwrRLgf37npmZ6fa+UTabzVVI9asjR45I27ZtTcho2bKl2da1a1dp1aqVTJ061e3j9OvXz3QqXrJkids1NzVq1JDjx4/n6bfjq1Spb2aPHj0kJsadltnIRVlRVlxXfAZDydcHTsqd725ye//597aV9sxZE9C/7/r9XblyZUlPTy/y+ztoNTebN2+WtLQ0Uwtjd+nSJVm3bp28/vrrJpBER0cXeZwOHTrI/PnzXT4eGxtrbvlpgfsrgPjz2FZDWVFWXFd8BkPBicyLHu/P3/nA/n335FhBCzfanLR9+/Y82+655x5p1KiRjB492q1go7Zu3WqaqwAA8PfyCd7uj8AKWrgpV66cNGvWLM+2smXLSkJCQu72MWPGyNGjR2Xu3LnmvjZX1a5dW5o2bSoXLlwwNTYLFy40NwAA/LV8giPdj2UUQlvQR0sVJiUlRQ4fPpx7XwPNyJEjTeApU6aMCTlLly6VPn36BPU8AQDhM3+Ns8VYHdffKkxUEetwITSEVLjRUU+OZs+enef+qFGjzA0AAE8s2/azPP3x93LyXE6eGhjHZRFcLZ/wv/1jZVw/5rkJByEVbgAA8LXJy3bKm+sOFNiuAUZrahwXDnVcfyvPDMU7v2MdrjBCuAEAWNaybSlOg42d9q/RmhoNNI5NVLr+luPQ5mVHv6UpKoyEzMKZAAD4uo+NNkUVRWtwtKYG1kG4AQBYkgaWk+cuuLWvdjKGddAsBQCw5Gio1PTzbj+XeWushXADALDkat6VypZy67mVysYwb43FEG4AAGFdW/P6Zz/KK5/+UOCxU242SU3q34zOwhZDuAEAhG1tzfglOyQ143+LIztyZ1XoBzrXkT4tqvn83BBchBsAQFgGG52jxp0AY296cpzAL6FsKZnYv5n0acHahFZEuAEAhF1TlPavcTfYqGd+01SS4ks7XXoB1kO4AQCEFR0R5Wx5hMJosHGcmA/WRrgBAITVMO+9x866/Tytm0liFe+IQ7gBAITdMG9PsIp35CHcAAAs03HYUf5VvxE5CDcAAMt0HLYb3r2BPNStAZ2GIxThBgBgmY7D1NZAEW4AACHJ3cUsH7qhnjSoWo4h3shFuAEAhCR3F7O8rn4VhnkjjxJ57wIAEBp0oj1tZnI11Z5u18d1P8AR4QYAEJJ0BmEd7aTyBxz7fYZ5wxnCDQAgZOkw7hl3Xm0m4nOk93U7w7zhDH1uAAAhTQNMjyZJuTMUszYUikK4AQCERRMVa0PBXYQbAEDA1oWi1gWBQLgBAAR0XSgm2oO/0aEYAODXdaHyzzKcmp5ltuvjgD8QbgAAAV0Xyr5NH9f9AF8j3AAAAr4ulEYafVz3A3yNcAMACNq6UO7uB3iCDsUAAJ+PhnJ3XSh39wM8QbgBAHht2bYUefrj7+XkuQt5RkM907eJ+amdh531qon6dZZh1oWCP9AsBQDwyuRlO+XB97fkCTbya1+aYe9vkZtbJpv7rAuFQCPcAAA8tmzbz/LmugMuH9famiXfpcj0Qa1ZFwoBR7MUAMDjPjbaFFUUrcGpWDZWvhjdjRmKEVCEGwCAR7Tz8MlzOW7tq52MWRcKgUazFADAI54M32Y0FIKBcAMA8Ii7gSWhbClGQyEoCDcAAI/o8G0d5l2Uif2bmSYpINAINwAAj2hgGdevSYEh3o4e6FxH+rT471BwINAINwAAj/Vuliwz7ry6QA1OpbIx8sag1jKmTxNKFUHDaCkAgNcBp0eTJIZ5I+QQbgAAXmOYN0IRzVIAAMBSCDcAAMBSCDcAAMBSQibcTJ48WaKiouSxxx4rdL+1a9dKmzZtpHTp0lK3bl2ZOXNmwM4RAABYMNzce++9cubMmQLbz507Zx7zxsaNG+Wtt96SFi1aFLrfgQMHpE+fPtKpUyfZunWrjB07Vh555BFZuHChV78XAABYj8fhZs6cOXL+/PkC23Xb3LlzPT6Bs2fPyh133CFvv/22VKxYsdB9tZamZs2aMnXqVGncuLH84Q9/MIHqpZde8vj3AgCACB8KnpGRITabzdy05kabhewuXboky5Ytk8TERI9PYNiwYdK3b1/p3r27TJo0qdB9169fLz179syzrVevXjJr1izJycmRmJiYAs/Jzs42N8fXoXR/vfmS/Xi+Pq4VUVaUFdcVn8Fwwd+r0CgvT47ndripUKGC6ROjt6uuuqrA47p9woQJ7p+liCxYsEC2bNlimqXckZqaKlWrVs2zTe9fvHhRjh8/LsnJyU778jg7r5UrV0pcXJz4w6pVq/xyXCuirCgrris+g+GCv1fBLa/MzEzfh5vVq1ebWptu3bqZPi6VKlXKfaxUqVJSq1YtqVatmtu/+MiRI/Loo4+akOFYC1QUDVGO9JycbbcbM2aMjBgxIk/NTY0aNUwNUHx8vPg6Veqb2aNHD6e1SKCsuK78i88gZcV1FXz++hzaW158Gm66dOmS26lXw0GJEsUbaLV582ZJS0szI58cm7fWrVsnr7/+umlKio6OzvOcpKQkU3vjSI9RsmRJSUhIcPp7YmNjzS0/LXB/BRB/HttqKCvKiuuKz2C44O9VcMvLk2N5vPyC1tCcPn1avvnmGxMsLl++nOfxwYMHu3WcG2+8UbZv355n2z333CONGjWS0aNHFwg2qmPHjvLJJ5/k2aY1P23btiVMAAAA78KNhgsd3aRDv8uVK5enOUj/7W640ec2a9Ysz7ayZcuaGhj7dm1SOnr0aO4orKFDh5paHW1muv/++00HY+1M/MEHH3j6MgAAgEV53Lb0+OOP5851ozU4p06dyr2dPHnSpyeXkpIihw8fzr1fp04dMyprzZo10qpVK5k4caJMmzZNBg4c6NPfCwAAIqjmRmtSdOI8f4w00tDiaPbs2U77/ugIKwAAAJ/U3Oi8Mps2bfL0aQAAAKFZc6MT7j3xxBOyc+dOad68eYGOvDfffLMvzw8AAMC/4UY78qo///nPBR7TDsU6nBsAACBswk3+od8AAAChpHgz8QEAAIR7zY2z5ihHzz77bHHOBwAAILDhZvHixQXWkNAlGXQJhHr16hFuAABAeIWbrVu3Ol3M6u6775bf/va3vjovAACA4PW50dW1tbnqmWee8cXhAAAAgt+hWJdiSE9P99XhAAAAAtMspWs5ObLZbGYNqHnz5knv3r29OwsAAIBghZtXXnklz/0SJUpIlSpVZMiQIWYVbwAAgLAKNzoyCgAAwJJ9bn766SezSjgAAEDYhhtdfkFHRpUvX15q1aolNWvWlAoVKsjEiRNZmgEAAIRfs9RTTz0ls2bNkilTpsh1111nOhR/+eWXMn78eMnKypLnnnvOP2cKAADgj3AzZ84ceeedd+Tmm2/O3dayZUu58sor5cEHHyTcAACA8GqWOnnypDRq1KjAdt2mjwEAAIRVuNFamtdff73Adt2mjwEAAIRVs9SLL74offv2lU8//VQ6duwoUVFR8tVXX8mRI0dk2bJl/jlLAAAAf9XcdOnSRX744QezSKYuuaBNUbfccovs2bNHOnXq5OnhAAAAgltzo6pVq0bHYQAAEN41N3v37pXf//73kpGRUeAxXTBz0KBBsn//fl+fHwAAgH/CzV/+8hepUaOGxMfHF3hMJ/TTx3QfAIgUly7bZP2+E/Lxt0fNT70PIIyapdatW2dW/nbl//7v/0ztDQBEghXfp8iET3ZKSnpW7rbk8qVlXL8m0rtZclDPDYh0btfcHDp0SBITE10+XrlyZTNiCgAiIdj8af6WPMFGpaZnme36OIAwCDfa9LRv3z6Xj//4449Om6wAwEq06UlrbJw1QNm36eM0UQFhEG46d+4sr732msvHp02bxlBwAJb3zYGTBWps8gccfVz3AxDi4WbMmDGyfPlyufXWW+Wbb74xI6T09vXXX8vAgQPl3//+t9kHAKws7UyWT/cDEMQOxa1bt5aPPvpI7r33Xlm8eHGexxISEuTvf/+7XH311X44RQAIHYnlSvt0PwBBnsTvN7/5jelYvGLFCtPHxmazyVVXXSU9e/aUuLg4P5weAISWa+pUMqOitPOws343USKSVL602Q9AmMxQXKZMGbP0AgBEougSUWa4t46K0iDjGHD0vtLHdT8AYbK2FABEOp3HZsadV5saGkd6X7czzw0QhmtLAUCk0wDTo0mSGRWlnYe1j402RVFjAwQf4QYAvKRBpmO9BMoPCDE0SwEAAEvxKtzoTMVPP/20WSU8LS3NbNMRVDt27PD1+QEAAPg33Kxdu1aaN29uJu9btGiRnD171mzftm2bjBs3ztPDAQAABDfcPPnkkzJp0iRZtWqVlCpVKnf7DTfcIOvXr/ft2QEAAPg73Gzfvt3pPDdVqlSREydOeHo4AACA4IabChUqSEpKSoHtW7dulSuvvNJX5wUAABCYcDNo0CAZPXq0pKamSlRUlFy+fFm+/PJLGTlypAwePNi7swAAAAhWuHnuueekZs2appZGOxM3adJEOnfuLNdee60ZQQUAABBWk/jFxMTI3/72N/nzn/9smqK05kZXDG/QoIF/zhAAACAQMxTXq1fP3AAAAMIu3IwYMcLtA7788stu7ztjxgxzO3jwoLnftGlTefbZZ+Wmm25yuv+aNWvMkPP8du3aJY0aNXL79wIAgAgPN9r85Gjz5s1y6dIladiwobn/ww8/SHR0tLRp08ajX169enWZMmWK1K9f39yfM2eO9O/f3/w+DTqu7NmzR+Lj4/MMQwcAAHA73KxevTpPzUy5cuVMEKlYsaLZdurUKbnnnnukU6dOHpVqv379CnRW1pqcDRs2FBpuEhMTzZB0AACAYve5+etf/yorV67MDTZK/62zFvfs2VMef/xx8YbWBP3jH/+Qc+fOSceOHQvdVzswZ2VlmZFaOkLLWVOVXXZ2trnZZWRkmJ85OTnm5kv24/n6uFZEWVFWXFd8BsMFf69Co7w8OV6UzWazeXJwrbX5+OOPpVu3bnm2f/bZZ6ZJ6cyZMx7PeKxhRsPKFVdcIe+//7706dPHZXPUunXrTPOXBpZ58+bJzJkzTV8cHY7uzPjx42XChAkFtuvviYuL8+hcAQBAcGRmZpq59tLT0/N0TfFJuNGJ+nTxTK3B6dChg9mmzUhPPPGECRjaXOWJCxcuyOHDh+X06dOycOFCeeedd8zxtVbG3aYtnUxwyZIlbtfc1KhRQ44fP15k4XiTKnXNrR49epgh86CsuK4Ci88gZcV1FXz++hzq93flypXdCjceN0tpTYnORnznnXfmVhGVLFlS7rvvPvnLX/7i8cnq4pv2DsVt27aVjRs3yquvvipvvvmmW8/XgDV//nyXj8fGxppbflrg/gog/jy21VBWlBXXFZ/BcMHfq+CWlyfH8jjcaFPOG2+8YYLMvn37RCt+NJyULVtWfEGP51jTUhQdWZWcnOyT3w0AACJ4Ej8NMy1atCjWLx87dqyZ00abibSvzoIFC0z/mRUrVpjHx4wZI0ePHpW5c+ea+1OnTpXatWubkVTanKU1NtqUpTcAAAC3w80tt9wis2fPNm1c+u/CLFq0yO2SPXbsmNx1111mlfHy5cubsKTBRtvplG7X/jh2Gmi0SUwDT5kyZUzIWbp0qcsOyAAAIPK4FW40eGinXfu/fWXWrFmFPq6BytGoUaPMDQAAoFjh5r333nP6bwAAgFBTItgnAAAAENQOxXXq1MltonJm//79xT0nAACAwIWbxx57LM99netGh2NrR2CdyA8AACCsws2jjz7qdPv06dNl06ZNvjgnAACA4Pe50flqmG8GAABYJtx89NFHUqlSJV8dDgAAIDDNUq1bt87ToViXS0hNTZVffvnFLMsAAAAQVuFmwIABee6XKFFCqlSpIl27dpVGjRr58twAAAD8H27GjRvn+W8BAAAI9YUzd+zYIZcuXcq9Hx0dbdZ6AgAACIsOxZ9//rm0a9cu936HDh1M/5tWrVqZmy56+emnn/rrPAEAAHwbbrSzsK7g7Wj16tVy4MABMyuxzn8zY8YMdw8HAAAQ3HCzceNGueaaa/Jsq169utSqVUtq165tgs/69ev9cY4AAAC+DzdHjx6V5OTk3Ptz5syRpKSk3Ps6x82JEyfc/80AAADBDDflypUzTVB2t9xyi8TFxeXe18fi4+N9f4YAAAD+CDft27eXuXPnunx89uzZZh8AAICwGAo+YsQI6d69uyQkJJjVvxMTE832tLQ0eeGFF2T+/PmycuVKf54rAACA78LNDTfcIK+99poMHz5cXn75ZdMEpcswpKenS8mSJWXq1KnSrVs3dw8HAAAQ/En8HnzwQenXr59ZJHPv3r1mW4MGDeTWW2+VGjVq+OcMAQAA/DlDsYYYrb0BAAAI6w7FAAAA4YBwAwAALIVwAwAALIVwAwAAIjvcPPXUU7Jq1SrJzMz0zxkBAAAEMtxs3rxZBg4cKBUrVpSOHTvKmDFjZMWKFXL27NninAcAAEBwwo0GmVOnTsmaNWukf//+snXrVrntttvMwpkdOnTwzVkBAAAEap4bFR0dbWptNNBoDY4uqvnPf/5T9u3b5+15AAAABKfmZsaMGXL77bdLcnKydOrUyawnpT+1ueqXX37xzVkBAAAEquZm2LBhUqVKFXn88cdl6NChZo0pAACAsK25WbRokdxxxx2yYMECszJ4+/btZfTo0bJ8+XI6FQMAgPCruRkwYIC5KV0R/PPPPzcLaWrnYl0lPDs72x/nCQAA4L8OxSdPnpS1a9eaEVN6+/777yUhIUG6dOnizeEAAACCF25atGghO3fuNCOlOnfuLPfff7907dpVmjVr5ruzAgAACFS4+eMf/0iYAQAA1gk3Dz30kPl54cIFOXDggNSrV09KlvSqdQsAACD4o6XOnz8v9913n8TFxUnTpk3l8OHDZvsjjzwiU6ZM8f0ZAgAA+DPcPPnkk/Ldd9+ZjsSlS5fO3d69e3f58MMPPT0cAACAT3ncnqTLLGiI0XWkdOi3XZMmTVh+AQAAhF/NjS6xoJP35Xfu3Lk8YQcAACAswk27du1k6dKlufftgebtt982i2kCAACEVbPU5MmTpXfv3maum4sXL8qrr74qO3bskPXr15uJ/QAAAMKq5ubaa6+VL7/8UjIzM80wcF0VvGrVqibctGnTxj9nCQAA4CavJqhp3ry5zJkzx5unAgAAhFbNjS/NmDHDLOcQHx9vbtpnR1cXL4w2fWkNkQ5Dr1u3rsycOTNg5wsAACwUbkqUKCHR0dGF3jydqbh69epm4r9NmzaZW7du3czq4tqHxxmdEblPnz7SqVMn2bp1q4wdO9ZMHrhw4UKPfi8AALAut9PI4sWLXT721VdfyWuvvSY2m82jX96vX78895977jlTm7NhwwYz+3F+WktTs2ZNmTp1qrnfuHFjE4peeuklGThwoEe/GwAARHi40RqV/Hbv3i1jxoyRTz75RO644w6ZOHGi1ydy6dIl+cc//mHmy3E1pFw7Lffs2TPPtl69esmsWbMkJydHYmJiCjwnOzvb3OwyMjLMT91fb75kP56vj2tFlBVlxXXFZzBc8PcqNMrLk+N51aH4559/lnHjxplOxRouvv32W2nWrJk3h5Lt27ebMJOVlSVXXHGFqSHS2Y6dSU1NNSOzHOl9HZJ+/PhxSU5Odjp0fcKECQW26ygvXR/LH1atWuWX41oRZUVZcV3xGQwX/L0KbnnpKG2/hJv09HR5/vnnTRNUq1at5D//+Y/p/1IcDRs2NOHo9OnTpu/MkCFDTKdhVwEn/yzI9qYwV7Mja83SiBEj8tTc1KhRw9QAaSdmX6dKfTN79OjhtBYJlBXXlX/xGaSsuK6Cz1+fQ3vLi0/DzYsvvigvvPCCJCUlyQcffOC0mcobpUqVkvr165t/t23bVjZu3GgmBnzzzTcL7Ku/W2tvHKWlpZmOzAkJCU6PHxsba275aYH7K4D489hWQ1lRVlxXfAbDBX+vgltenhyrpCergZcpU8YEEW2OcjXPzaJFi6Q4tCbGsY+MI22+0v49+ZuXNBQRJgAAgEfhZvDgwT5fGFOHct90002mmejMmTOyYMECWbNmjaxYsSK3Seno0aMyd+5cc3/o0KHy+uuvm2am+++/33Qw1s7EWpMEAADgUbiZPXu2z0vs2LFjctddd0lKSoqUL1/eTOinwUbb6ZRuP3z4cO7+derUkWXLlsnw4cNl+vTpUq1aNZk2bRrDwAEAQPFGS/mK1rp4Gqi6dOkiW7Zs8eNZAQCAcBbU5RcAAAB8jXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAsJajhZvLkydKuXTspV66cJCYmyoABA2TPnj2FPmfNmjUSFRVV4LZ79+6AnTcAAAhdQQ03a9eulWHDhsmGDRtk1apVcvHiRenZs6ecO3euyOdqCEpJScm9NWjQICDnDAAAQlvJYP7yFStW5Ln/3nvvmRqczZs3S+fOnQt9ru5XoUIFP58hAAAIN0ENN/mlp6ebn5UqVSpy39atW0tWVpY0adJEnn76abnhhhuc7pednW1udhkZGeZnTk6OufmS/Xi+Pq4VUVaUFdcVn8Fwwd+r0CgvT44XZbPZbBIC9DT69+8vp06dks8//7zQ5qh169ZJmzZtTGiZN2+ezJw50/TFcVbbM378eJkwYUKB7e+//77ExcX5/HUAAADfy8zMlEGDBpmKkPj4+PAIN9r3ZunSpfLFF19I9erVPXpuv379TKfiJUuWuFVzU6NGDTl+/HiRheNNqtS+Qz169JCYmBifHttqKCvKiuuKz2C44O9VaJSXfn9XrlzZrXATEs1SDz/8sAkmWiPjabBRHTp0kPnz5zt9LDY21tzy0wL3VwDx57GthrKirLiu+AyGC/5eBbe8PDlWUMONVhppsFm8eLFpVqpTp45Xx9m6daskJyf7/PwAAED4KRnspijt+/Lxxx+buW5SU1PN9vLly0uZMmXMv8eMGSNHjx6VuXPnmvtTp06V2rVrS9OmTeXChQumxmbhwoXmBgAAENRwM2PGDPOza9euBYaE33333ebfOofN4cOHcx/TQDNy5EgTeDQAacjRvjp9+vQJ8NkDAIBQFPRmqaLMnj07z/1Ro0aZGwAAgDOsLQUAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcOMjly7b5OsDJ2Xz8SjzU+8DAIAICzeTJ0+Wdu3aSbly5SQxMVEGDBgge/bsKfJ5a9eulTZt2kjp0qWlbt26MnPmTAmmFd+nyPUvfCZ3vrtJ5u6NNj/1vm4HAAARFG40pAwbNkw2bNggq1atkosXL0rPnj3l3LlzLp9z4MAB6dOnj3Tq1Em2bt0qY8eOlUceeUQWLlwowaAB5k/zt0hKelae7anpWWY7AQcAgMAqKUG0YsWKPPffe+89U4OzefNm6dy5s9PnaC1NzZo1ZerUqeZ+48aNZdOmTfLSSy/JwIEDJZC06WnCJzvFWQOUbosSMY/3aJIk0SX0HgAAsHS4yS89Pd38rFSpkst91q9fb2p3HPXq1UtmzZolOTk5EhMTk+ex7Oxsc7PLyMgwP3VfvRWH9q3JX2OTP+Do4+t/TJP2dVy/pkhkL/vivgeRgLKirLiu+AyGkxw//X335HghE25sNpuMGDFCrr/+emnWrJnL/VJTU6Vq1ap5tul9bdI6fvy4JCcnF+jXM2HChALHWblypcTFxRXrnLXzsEh0kfut/PxrObGLDsbOaHMk3ENZuY+yoqz8gesquOWVmZkZfuHmoYcekm3btskXX3xR5L5RUVEFgpGz7WrMmDEmNDnW3NSoUcPU/sTHxxfrnBMOnJS5ezcVuV/PTu2puXGSwPXC79GjR4HaNlBW3uK6oqz8gesqNMrL3vISNuHm4YcfliVLlsi6deukevXqhe6blJRkam8cpaWlScmSJSUhIaHA/rGxseaWnxZ4cQu9Y/1ESS5f2nQedlYvo1ErqXxpsx99bpzzxfsQKSgryorris9gJP/NivHgWEEdLaU1Llpjs2jRIvnss8+kTp06RT6nY8eOBaq6tImpbdu2Af+S1MAyrl8T8+/8dUb2+/o4wQYAgMAJarjRYeDz58+X999/38x1ozUyejt//nyeZqXBgwfn3h86dKgcOnTINDXt2rVL3n33XdOZeOTIkUF5Db2bJcuMO682NTSO9L5u18cBAEDgBLVZasaMGeZn165dCwwJv/vuu82/U1JS5PDhw7mPae3OsmXLZPjw4TJ9+nSpVq2aTJs2LeDDwB1pgNHh3joqSjsPax8bmqIAAIjAcGPvCFyY2bNnF9jWpUsX2bJli4QSbXrS4d46Kkp/0hQFAEBwsLYUAACwFMINAACwFMINAACwFMINAACwFMINAACwFMINAACwFMINAACwFMINAACwFMINAACwlJBYFTwYsyJ7snS6J8u8Z2ZmmmOz0jVlxXUVeHwGKSuuq+Dz1+fQ/r3tzuoGERduzpw5Y37WqFEj2KcCAAC8+B4vX758oftE2dyJQBZy+fJl+fnnn80q5FFRUT49tqZKDU1HjhyR+Ph4nx7baigryorris9guODvVWiUl8YVDTa6YHaJEoX3qom4mhstkOrVq/v1d+ibSbihrLiugofPIGXFdWXNz2FRNTZ2dCgGAACWQrgBAACWQrjxodjYWBk3bpz5CcqK6yrw+AxSVlxXwRcKn8OI61AMAACsjZobAABgKYQbAABgKYQbAABgKYQbAABgKYQbN02ePFnatWtnZjZOTEyUAQMGyJ49e4p83tq1a6VNmzZSunRpqVu3rsycOVOszpuyWrNmjZkxOv9t9+7dYmUzZsyQFi1a5E521bFjR1m+fHmhz4nEa8qbsorUa8rVZ1Jf+2OPPVbofpF6bXlaVpF8bY0fP77A605KSgq564pw4yZ9c4YNGyYbNmyQVatWycWLF6Vnz55y7tw5l885cOCA9OnTRzp16iRbt26VsWPHyiOPPCILFy4UK/OmrOw0BKWkpOTeGjRoIFams2VPmTJFNm3aZG7dunWT/v37y44dO5zuH6nXlDdlFanXVH4bN26Ut956ywTDwkTyteVpWUX6tdW0adM8r3v79u2hd13pUHB4Li0tTYfQ29auXetyn1GjRtkaNWqUZ9sDDzxg69ChQ0QVuTtltXr1arPPqVOnbJGuYsWKtnfeecfpY1xT7pcV15TNdubMGVuDBg1sq1atsnXp0sX26KOPurzuIv3a8qSsIvnaGjdunK1ly5Zu7x+s64qaGy+lp6ebn5UqVXK5z/r1602NhaNevXqZ/+vUJeEjhTtlZde6dWtJTk6WG2+8UVavXi2R5NKlS7JgwQJTw6VNLs5wTblfVnaRfE1pDWrfvn2le/fuRe4b6deWJ2UV6dfW3r17zeKVderUkdtvv132798fctdVxC2c6Qs67+GIESPk+uuvl2bNmrncLzU1VapWrZpnm97XZprjx4+bD4XVuVtWWhZaHaztstnZ2TJv3jzzB0Pbtjt37ixWplW6+gWdlZUlV1xxhSxevFiaNGnidN9Iv6Y8KatIvqaUhr8tW7aYphZ3RPK15WlZRfK11b59e5k7d65cddVVcuzYMZk0aZJce+21pnk4ISEhZK4rwo0XHnroIdm2bZt88cUXRe6rna0c2SeEzr890suqYcOG5manX2BHjhyRl156yfJ/LPR1f/vtt3L69GnTDj1kyBDTb8nVl3YkX1OelFUkX1P6Oh999FFZuXKl6cTprki8trwpq0i+tm666abcfzdv3ty89nr16smcOXPM/8iGynVFs5SHHn74YVmyZImpgtQOjoXRHuSaWh2lpaVJyZIlnSbcSC4rZzp06GCqP62uVKlSUr9+fWnbtq0ZqdGyZUt59dVXne4b6deUJ2UVydfU5s2bzXWhNQt6behNQ+C0adPMv7VZL79Ivba8KatIvrbyK1u2rAk5rl57sK4ram7cpElTv6y1GlyrHrWtsSiaaD/55JM82/T/DvQPc0xMjFiVN2XljPast3JVeGHlp1XdzkTqNeVNWUXyNaVNJPlHsNxzzz3SqFEjGT16tERHRxd4TqReW96UVSRfW/np52/Xrl1mNJQzQbuu/Npd2UL+9Kc/2cqXL29bs2aNLSUlJfeWmZmZu8+TTz5pu+uuu3Lv79+/3xYXF2cbPny4befOnbZZs2bZYmJibB999JHNyrwpq1deecW2ePFi2w8//GD7/vvvzeN6eS5cuNBmZWPGjLGtW7fOduDAAdu2bdtsY8eOtZUoUcK2cuVK8zjXlPdlFanXlCv5RwBxbXlfVpF8bT3++OPmb7t+v23YsMH2m9/8xlauXDnbwYMHQ+q6Ity4W1AiTm/vvfde7j5DhgwxHwpHehG0bt3aVqpUKVvt2rVtM2bMsFmdN2X1wgsv2OrVq2crXbq0Gd57/fXX25YuXWqzunvvvddWq1Ytc31UqVLFduONN+Z+WSuuKe/LKlKvKXe/sLm2vC+rSL62brvtNltycrIJKNWqVbPdcsstth07doTcdRWl//FfvRAAAEBg0aEYAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGiEAHDx40K/LqCtvumj17tlSoUCHo51Ec7ryG8ePHS6tWrQJyPgD8g3ADhKkjR47IfffdJ9WqVTOrZdeqVUseffRROXHiRJHPrVGjhqSkpEizZs3c/n233Xab/PDDDxIMP/74o9x7771Ss2ZNiY2NlSuvvNIsePi3v/1NLl68GJKvQVdC1gVk69ata85Zy7xfv37yn//8RwJNA+Q///nPgP9eIFhYFRwIQ/v37zer7V511VXywQcfmJXXd+zYIU888YQsX75cNmzYIJUqVXL63AsXLpgwlJSU5NHvLFOmjLkF2jfffCPdu3eXpk2byvTp081qzWfPnpWdO3fKzJkzTUBr2bJlSL0GrZG67rrrTC3Riy++KC1atJCcnBz597//LcOGDZPdu3f7/RyAiOb31asA+Fzv3r1t1atXz7PSutLV13UF3qFDh+Zu08UmJ06caBa0i4+Ptw0ePNisrK0f/61bt+bu9/HHH9vq169vFgPs2rWrbfbs2WafU6dOmcd14VNd7d1u3LhxtpYtW9rmzp1rfoceWxfVy8jIyN1n+fLltuuuu848r1KlSra+ffvafvzxx9zHnZ2Ho8uXL9saN25sa9Omje3SpUsu91GrV6/Oc75Kj6vb9Pc4ew1q8uTJtsTERNsVV1xhFuccPXq0eV2O3n33XVujRo1ssbGxtoYNG9qmT59uK8xNN91ku/LKK21nz54t8Jjj+R06dMh2880328qWLWtWVv7d735nS01NzX1c37P+/fvneb4u6Oi4MKH+++GHH7Y98cQTZhHHqlWrmvfGTt8bxwVs9T5gdTRLAWHm5MmTpgbgwQcfLFALobUxd9xxh3z44Yf6Py652//yl7+YGo7NmzfLM88847Sm4dZbb5UBAwaY/i8PPPCAPPXUU0Wey759+0xzx7/+9S9zW7t2rUyZMiX38XPnzsmIESNk48aNpjmmRIkS8tvf/lYuX77s1mvVc9m1a5eMHDnSPNdVk4u3/v73v8u4cePkueeek02bNklycrK88cYbefZ5++23TVnoPnouzz//vCnDOXPmuHx/VqxYYWpoypYtW+Bxe58ffX+0vHV/LbdVq1aZ8tSmM0/puejv+vrrr01N0Z///GdzPKVlr9577z3TFGm/D1gZzVJAmNm7d6/5YmzcuLHTx3X7qVOn5JdffpHExESzrVu3biYgOIYZR9q807BhQxOClP77+++/N1/ohdGQop10y5UrZ+7fddddJsTYnzdw4MA8+8+aNcuckzYpudPfx94/Rs/HLi0tzfRjsdMvcw163pg6darpy/OHP/zB3J80aZJ8+umnkpWVlbvPxIkT5a9//avccsst5r42Aer5v/nmmzJkyBCn/YP0/dHms8Lo79m2bZscOHDA9MdR8+bNM81vGkDatWvn9uvQZi8NaapBgwby+uuvm/ehR48eUqVKldxQ5WlTJBCuqLkBLMZeY+NYo9G2bdtCn7Nnz54CX6bXXHNNkb+rdu3aucFGac2Hhg87rYkYNGiQCSPx8fEmGKjDhw978IryvpaEhARTo6M3/cLWPkTe0poY7bvkyPG+BkR7x+0rrrgi96YhSF+bu+Xv6ndrqLEHG9WkSRPzmvQxT2i4cZT/fQAiDTU3QJipX7+++eLU2gNt1shPO6tWrFhRKleunLvNWfNI/i/k/F/Gjs1arsTExOS5r8dwbHLS0UH65a1NOzqqSx/TGht3A4nWQthfk314dnR0tCkDVbLk//6E2ZutHM9bO/EWh/216Pm3b98+z2N6Hq7OWctBA4qz96ewMs+/XV9T/vfB2Wsq6n0AIg01N0CY0ZoLbW7QviHnz58vMPxYh0drvw1P+qJoE0r+vhjaB6U4dEi6fsE//fTTZti2vbnME61btzbn9tJLLxX5ZW1vftF+JXZFzZ+j56Qjyxw53q9ataoZdq6j0zRQOd7stVD56Si1Xr16mZFd2ucov9OnT+fW0mgNltYM2WlgTU9Pz21y1Nfk+HrceU3OaPi5dOmSx88DwhXhBghD2qciOzvbfImuW7fOfEFqJ1YNPfplXFRfmfy0A7HWjowePdr0c9GOttqXpjgddrX2SIPYW2+9ZfqhfPbZZ6ZzsSf0d2tHWG0206HVS5YsMX2O7MPAtdnIXoOigUNriXQSPn0NS5cuNX1lCqPzAr377rvmps/Rfis6pN6RHm/y5Mny6quvmn22b99uzunll192eVwNnhomtGlv4cKF5pw16E2bNi232UuHt2tzknYA37JlixnyPnjwYOnSpUtuM6L2ldKQOXfuXHMMPT/tC+UpbT7UPjgafj0NmEA4ItwAYUibPvRLr169eqaWRn/+8Y9/lBtuuEHWr1/vco4bV7QW4qOPPpJFixaZL9wZM2bkjpbSCei8oU0qCxYsMCO0tClq+PDhuR2WPdGhQwdzDO1UrCOQtMbj2muvNfP7vPLKK/KnP/0pt3ZCt2lI03lvXnjhBdM3pjBads8++6wJdW3atJFDhw7lHs9OOxu/8847Juw1b97chA/9t6uaG6WPaWDR9+Pxxx83r1+DpwYMLVvHifU0BHbu3NmEHe2bpCPd7DS86sisUaNGmT5RZ86cMQHIUxrydPSUhj+tDQOsLkrHgwf7JACEHq390doRx2YTAAgHdCgGkNuUorUD2pT05ZdfmlqWhx56iNIBEHYINwAM7dOhzTg6qZyu4aTNKWPGjKF0AIQdmqUAAICl0KEYAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAACIlfw/CW1AZAnVYpQAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "\n", + "plt.plot(dat[\"guide_count_1\"], dat[\"guide_count_2\"], 'o')\n", + "plt.xlabel('Original Guide Count')\n", + "plt.ylabel('New Guide Count')\n", + "plt.grid(True)" + ] + }, + { + "cell_type": "code", + "execution_count": 42, + "id": "9a28eb0e", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjoAAAGwCAYAAACgi8/jAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAK81JREFUeJzt3Ql4E9X+//FvoaUIArLLJpuAFAQEKgJ6FWQRFERQ8aKIuFyRoixuoKgUF64iiEtBcIGrVwRFQB9FEZRNcCkoohavokV2kaKsUko6/+d7/r/kSfemaZvkzPv1PFEymU7OnExmPjnnzEyU4ziOAAAAWKhMqAsAAABQUgg6AADAWgQdAABgLYIOAACwFkEHAABYi6ADAACsRdABAADWihaXy8zMlD179kilSpUkKioq1MUBAACFoJcBPHLkiNStW1fKlMm73cb1QUdDToMGDQpTpwAAIMzs3LlT6tevn+frrg862pLjrajKlSsXW8VnZGTIxx9/LL169ZKYmJhiW65NqCPqiO2I71q4YH8UeXV0+PBh01DhPY7nxbVBJykpyTw8Ho95riGnuINOhQoVzDLDYYMIR9QRdcR2xHctXLA/itw6KmjYiWsHIyckJEhKSookJyeHuigAAKCEuDboAAAA+xF0AACAtQg6AADAWgQdAABgLYIOAACwlmuDjp5aHhcXJ/Hx8aEuCgAAKCGuDTqcXg4AgP1cG3QAAID9CDoAAMBaBB0AAGAtgg4AALAWQQcAAFiLoAMAAKwVLS6+jo4+PB5Pib5P60nLJd2T/y3ki2L7vy8v9mUCAGAb17bocB0dAADs59qgAwAA7EfQAQAA1iLoAAAAaxF0AACAtQg6AADAWgQdAABgLYIOAACwFkEHAABYy7VBR6+KHBcXJ/Hx8aEuCgAAKCGuDTpcGRkAAPu5NugAAAD7EXQAAIC1CDoAAMBaBB0AAGAtgg4AALAWQQcAAFiLoAMAAKxF0AEAANYi6AAAAGu5NuhwCwgAAOzn2qDDLSAAALCfa4MOAACwH0EHAABYi6ADAACsRdABAADWIugAAABrEXQAAIC1CDoAAMBaBB0AAGAtgg4AALAWQQcAAFiLoAMAAKxF0AEAANYi6AAAAGu5NugkJSVJXFycxMfHh7ooAACghLg26CQkJEhKSookJyeHuigAAKCEuDboAAAA+xF0AACAtQg6AADAWgQdAABgLYIOAACwFkEHAABYi6ADAACsRdABAADWIugAAABrEXQAAIC1CDoAAMBaBB0AAGAtgg4AALAWQQcAAFiLoAMAAKxF0AEAANYi6AAAAGsRdAAAgLUIOgAAwFoEHQAAYC3XBp2kpCSJi4uT+Pj4UBcFAACUENcGnYSEBElJSZHk5ORQFwUAAJQQ1wYdAABgP4IOAACwFkEHAABYi6ADAACsRdABAADWIugAAABrEXQAAIC1CDoAAMBaBB0AAGAtgg4AALAWQQcAAFiLoAMAAKxF0AEAANYi6AAAAGsRdAAAgLUIOgAAwFoEHQAAYC2CDgAAsBZBBwAAWIugAwAArEXQAQAA1iLoAAAAaxF0AACAtQg6AADAWgQdAABgLYIOAACwFkEHAABYi6ADAACsRdABAADWivigs3PnTrnkkkskLi5O2rRpI2+//XaoiwQAAMJEtES46OhomTFjhrRr1072798v7du3l759+0rFihVDXTQAABBiER906tSpYx6qVq1aUq1aNTl48CBBBwAAhL7rau3atdKvXz+pW7euREVFydKlS3PMM3PmTGncuLGUL19eOnToIOvWrct1WRs3bpTMzExp0KBBKZQcAACEu5C36Bw7dkzatm0rw4cPl0GDBuV4feHChTJmzBgTdrp27SqzZ8+WPn36SEpKipx11lm++dLS0uTGG2+Ul19+Od/3S09PNw+vw4cPm/9nZGSYR3HxLiu2jFNsy8xt+ZHMuw42rEtJoY6oI7YjvmvhIiPM9tmFLUeU4zglcyQuAm3RWbJkiQwYMMA3rVOnTmbczaxZs3zTWrZsaeaZMmWKea7BpWfPnnLbbbfJ0KFD832PSZMmSWJiYo7p8+fPlwoVKhTr+gAAgJJx/PhxGTJkiBw6dEgqV64cvi06+Tl58qRs2rRJxo8fn2V6r169ZMOGDebfmtNuuukm6d69e4EhR02YMEHGjRuXpUVHu7p0mflVVFGS5ooVK+ShjWUkPTNKitv3k3pLpPPWkYbUmJiYUBcnLFFH1BHbEd+1cJERZvtsb49MQcI66Bw4cEA8Ho/Url07y3R9vm/fPvPv9evXm+4tPbXcO77n9ddfl3PPPTfXZcbGxppHdvqhlcQHpyEn3VP8QSccNrLiUlJ1bxPqiDpiO+K7Fi5iwmSfXdgyhHXQ8e/S8qetON5pF154oRmADAAAEHZnXeWnRo0aUrZsWV/rjZdeLyd7Kw8AAEBEBZ1y5cqZ08m1T9CfPu/SpUvIygUAACJDyLuujh49Ktu2bfM9T01Nlc2bN5sL/+np4zpwWAcZd+zYUTp37ixz5syRHTt2yIgRI4J636SkJPPQMUAAAMBOIQ86epG/bt26+Z57z4gaNmyYzJs3TwYPHmyukTN58mTZu3evtG7dWpYtWyYNGzYM6n0TEhLMQ0dtV6lSJej1AAAA4SfkQUdvyFnQpXxGjhxpHgAAANaM0QEAAAgGQQcAAFiLoAMAAKzl2qCjZ1zFxcVJfHx8qIsCAABKiGuDjp5xpXdAT05ODnVRAABACXFt0AEAAPYj6AAAAGsRdAAAgLUIOgAAwFoEHQAAYC3XBh1OLwcAwH6uDTqcXg4AgP1cG3QAAID9CDoAAMBaBB0AAGAtgg4AALAWQQcAAFiLoAMAAKzl2qDDdXQAALCfa4MO19EBAMB+rg06AADAfgQdAABgLYIOAACwFkEHAABYi6ADAACsRdABAADWIugAAABruTbocMFAAADs59qgwwUDAQCwn2uDDgAAsB9BBwAAWIugAwAArEXQAQAA1iLoAAAAaxF0AACAtQg6AADAWgQdAABgLYIOAACwFkEHAABYy7VBh3tdAQBgP9cGHe51BQCA/VwbdAAAgP0IOgAAwFoEHQAAYC2CDgAAsBZBBwAAWIugAwAArFWkoNOkSRNJS0vLMf2vv/4yrwEAAERs0Nm+fbt4PJ4c09PT02X37t3FUS4AAICgRQcy83vvvef79/Lly6VKlSq+5xp8PvnkE2nUqFHwpQIAACjtoDNgwADz/6ioKBk2bFiW12JiYkzImTZtWnGUCwAAoHSDTmZmpvl/48aNJTk5WWrUqBF8CQAAAMIh6HilpqYWf0kAAADCIegoHY+jj/379/taerxeffXV4igbAABA6QedxMREmTx5snTs2FHq1KljxuxEmqSkJPPI7ewxAADg4qDz4osvyrx582To0KESqRISEszj8OHDWc4eAwAALr+OzsmTJ6VLly7FXxoAAIBQB51bb71V5s+fX5zlAAAACI+uqxMnTsicOXNk5cqV0qZNG3MNHX/Tp08vrvIBAACUbtDZsmWLtGvXzvz7+++/z/JaJA5MBgAAdipS0Fm1alXxlwQAACAcxugAAABY26LTrVu3fLuoPv3002DKBAAAELqg4x2f45WRkSGbN28243Wy3+wTAAAgooLOM888k+v0SZMmydGjR4MtEwAAQPiN0bnhhhu4zxUAALAz6Hz++edSvnz54lwkAABA6XZdDRw4MMtzx3Fk7969snHjRnnooYeKXhoAAIBQB53sN8EsU6aMtGjRwtzRvFevXsVVNgAAgNIPOnPnzg3uXQEAAMI16Hht2rRJtm7daq6pExcXJ+edd17xlQwAACAUQWf//v1y3XXXyerVq+WMM84wY3QOHTpkLiS4YMECqVmzZrDlAgAACM1ZV3feeaccPnxYfvjhBzl48KD8+eef5mKBOu2uu+4KvlQAAAChatH56KOPZOXKldKyZUvfNO26SkpKYjAyAACI7BadzMxMiYmJyTFdp+lrAAAAERt0unfvLqNHj5Y9e/b4pu3evVvGjh0rl156aXGWDwAAoHSDzgsvvCBHjhyRRo0aSdOmTeXss8+Wxo0bm2nPP/+8RALtZtPutvj4+FAXBQAAhNMYnQYNGsjXX38tK1askB9//NGcdaWhoUePHhIpEhISzEMHUGe/ACIAAHBhi86nn35qAo2GA9WzZ09zBpaeaaUtI61atZJ169aVVFkBAABKLujMmDFDbrvtNqlcuXKO17RV5Pbbb5fp06cHVgIAAIBwCDrffvutXHbZZXm+rve50qslAwAARFzQ+f3333M9rdwrOjpa/vjjj+IoFwAAQOkGnXr16sl3332X5+tbtmyROnXqBF8qAACA0g46ffv2lYcfflhOnDiR47W///5bHnnkEbniiiuKo1wAAACle3r5xIkTZfHixdK8eXMZNWqUtGjRwty5XO9grtel8Xg88uCDDwZfKgAAgNIOOrVr15YNGzbIHXfcIRMmTDDXz1Eadnr37i0zZ8408wAAAETkBQMbNmwoy5YtM3cs37Ztmwk7zZo1k6pVq5ZMCQEAAErzyshKgw23TwAAANbd6woAACASEHQAAIC1CDoAAMBaBB0AAGAtgg4AALAWQQcAAFiLoAMAAKxF0AEAANYi6AAAAGsRdAAAgLUIOgAAwFoEHQAAYC2CDgAAsBZBBwAAWCs61AWAezQa/0GW57FlHXnqfJHWk5ZLuieqyMvd/u/Li6F0AAAb0aIDAACsRdABAADWIugAAABrEXQAAIC1rAg6V111lVStWlWuvvrqUBcFAACEESuCzl133SWvvfZaqIsBAADCjBVBp1u3blKpUqVQFwMAAISZkAedtWvXSr9+/aRu3boSFRUlS5cuzTHPzJkzpXHjxlK+fHnp0KGDrFu3LiRlBQAAkSXkFww8duyYtG3bVoYPHy6DBg3K8frChQtlzJgxJux07dpVZs+eLX369JGUlBQ566yzAn6/9PR08/A6fPiw+X9GRoZ5FBfvsmLLOMW2zNyWH0n0AoFZnv9f3QRbR5FYF4Gum83rGCzqiDpiO3Lndy2jkOWIchynZI7ERaAtOkuWLJEBAwb4pnXq1Enat28vs2bN8k1r2bKlmWfKlCm+aatXr5YXXnhBFi1alO97TJo0SRITE3NMnz9/vlSoUKHY1gUAAJSc48ePy5AhQ+TQoUNSuXLl8G3Ryc/Jkydl06ZNMn78+CzTe/XqJRs2bCjSMidMmCDjxo3L0qLToEEDs8z8KqooSXPFihXy0MYykp5Z9Nsb5OX7Sb0l0uitHvxpS86jHTODrqNIrItAt6OePXtKTExMqIsTlqgj6ojtyJ3ftcP/1yNTkLAOOgcOHBCPxyO1a9fOMl2f79u3z/e8d+/e8vXXX5tusPr165tWofj4+FyXGRsbax7Z6YdWEh+cHsCDuY9TXsJhIwtUXvUQbB1FYl0EqqS2T5tQR9QR25G7vmsxhSxDWAcd/y4tf9rb5j9t+fKsLQUAAABhcdZVfmrUqCFly5bN0nqj9u/fn6OVBwAAIKKCTrly5czp5Non6E+fd+nSJWTlAgAAkSHkXVdHjx6Vbdu2+Z6npqbK5s2bpVq1aub0cR04PHToUOnYsaN07txZ5syZIzt27JARI0YE9b5JSUnmoWOAIlGj8R+EuggAAIS9kAedjRs3misbe3nPiBo2bJjMmzdPBg8eLGlpaTJ58mTZu3evtG7dWpYtWyYNGzYM6n0TEhLMQ0dtV6lSJej1AAAA4SfkQeeSSy4xg4vzM3LkSPMAAACwZowOAABAMAg6AADAWgQdAABgLdcGHT3jKi4uLs8rKAMAgMjn2qCjZ1zpHdCTk5NDXRQAAFBCXBt0AACA/Qg6AADAWgQdAABgLYIOAACwFkEHAABYy7VBh9PLAQCwn2uDDqeXAwBgP9cGHQAAYD+CDgAAsBZBBwAAWIugAwAArEXQAQAA1iLoAAAAa0WLi6+jow+PxxPqogAAEFKNxn9Q4DyxZR156nyR1pOWS7onqtDL3v7vyyWUXNuiw3V0AACwn2uDDgAAsB9BBwAAWIugAwAArEXQAQAA1iLoAAAAaxF0AACAtQg6AADAWlwwkAsGIsiLaBVFqC+gBQBu4doWHS4YCACA/VwbdAAAgP0IOgAAwFoEHQAAYC2CDgAAsBZBBwAAWIugAwAArEXQAQAA1iLoAAAAaxF0AACAtQg6AADAWtzrintdIYzvoRVb1pGnzhdpPWm5pHuiJJS4PxeASOTaFh3udQUAgP1cG3QAAID9CDoAAMBaBB0AAGAtgg4AALAWQQcAAFiLoAMAAKxF0AEAANYi6AAAAGsRdAAAgLUIOgAAwFoEHQAAYC2CDgAAsBZBBwAAWCtaXCopKck8PB5PqIsCRIRG4z8osWVv//flJbZsAO7m2hadhIQESUlJkeTk5FAXBQAAlBDXBh0AAGA/gg4AALAWQQcAAFiLoAMAAKxF0AEAANYi6AAAAGsRdAAAgLUIOgAAwFoEHQAAYC2CDgAAsBZBBwAAWIugAwAArEXQAQAA1iLoAAAAaxF0AACAtQg6AADAWgQdAABgLYIOAACwFkEHAABYi6ADAACsFS0ulZSUZB4ejyfURUGQGo3/gDoEAOTKtS06CQkJkpKSIsnJyaEuCgAAKCGuDToAAMB+BB0AAGAtgg4AALAWQQcAAFiLoAMAAKxF0AEAANYi6AAAAGsRdAAAgLUIOgAAwFoEHQAAYC2CDgAAsBZBBwAAWIugAwAArEXQAQAA1iLoAAAAaxF0AACAtQg6AADAWgQdAABgLYIOAACwFkEHAABYi6ADAACsRdABAADWIugAAABrEXQAAIC1CDoAAMBaBB0AAGAtgg4AALAWQQcAAFiLoAMAAKxlRdB5//33pUWLFtKsWTN5+eWXQ10cAAAQJqIlwp06dUrGjRsnq1atksqVK0v79u1l4MCBUq1atVAXDQAAhFjEt+h89dVX0qpVK6lXr55UqlRJ+vbtK8uXLw91sQAAQBgIedBZu3at9OvXT+rWrStRUVGydOnSHPPMnDlTGjduLOXLl5cOHTrIunXrfK/t2bPHhByv+vXry+7du0ut/AAAIHyFvOvq2LFj0rZtWxk+fLgMGjQox+sLFy6UMWPGmLDTtWtXmT17tvTp00dSUlLkrLPOEsdxcvyNBqa8pKenm4fX4cOHzf8zMjLMo7h4lxVbJmf5IFnqhjrKm1vqKJjvnvdvi/P7axvqiDoqSGxZp8T2RyX13SzscqOc3JJCiGhAWbJkiQwYMMA3rVOnTmbczaxZs3zTWrZsaeaZMmWKbNiwQaZOnWr+To0ePdr8zZAhQ3J9j0mTJkliYmKO6fPnz5cKFSqUyHoBAIDidfz4cXOsP3TokBmjG5FB5+TJkyZ8vP3223LVVVf55tMws3nzZlmzZo0ZjKzBZ/Xq1b7ByF988YVUr1690C06DRo0kAMHDuRbUUVJmitWrJCHNpaR9My8W5jcTH8VPNoxkzqijtiOLPiufT+pt0Qy7z67Z8+eEhMTE9SyWk+yc5xobBG3o5LaNvT4XaNGjQKDTsi7rvKj4cPj8Ujt2rWzTNfn+/btM/+Ojo6WadOmSbdu3SQzM1Puu+++PEOOio2NNY/sdMMOduPOjW4M6R6CDnXEdlTS+K6Fto5KYv8ZCsVxLLB9n58e4HZUUttGYZcb1kEnrzE32gjlP61///7mAQAAEFZnXeVHm6TKli3ra73x2r9/f45WHgAAgIgKOuXKlTOnk2u/qT993qVLl6CWnZSUJHFxcRIfHx9kKQEAQLgKedfV0aNHZdu2bb7nqampZqCxXtlYTx/Xqx4PHTpUOnbsKJ07d5Y5c+bIjh07ZMSIEUG9b0JCgnnoYKYqVaoUw5oAAIBwE/Kgs3HjRjOQ2EuDjRo2bJjMmzdPBg8eLGlpaTJ58mTZu3evtG7dWpYtWyYNGzYMYakBAEAkCHnQueSSS3K96J+/kSNHmgcAAIA1Y3QAAACCQdABAADWIugAAABruTbocHo5AAD2c23Q0VPL9Q7oycnJoS4KAAAoIa4NOgAAwH4EHQAAYK2QX0cn1LzX8NErJBenjIwMOX78uHjSy0qm5XeyLSpPWUeOH/dQR9QR25EF37Xi3oeWNu8+W9cj2LttZ6YfFxt5irgdldS24V1uQdfii3IKmsNyu3btkgYNGoS6GAAAoAh27twp9evXz/N11wedzMxM2bNnj1SqVEmioqKKNWlqgNIPoHLlysW2XJtQR9QR2xHftXDB/ijy6kjbaY4cOSJ169aVMmXyHonj+q4rrZz8kmCwdGMIhw0inFFH1BHbEd+1cMH+KLLqqDA35WYwMgAAsBZBBwAAWIugU0JiY2PlkUceMf8HdcR2VHL4rlFHbEelIzZCj2uuH4wMAADsRYsOAACwFkEHAABYi6ADAACsRdABAADWIugEYebMmdK4cWMpX768dOjQQdatW5fv/GvWrDHz6fxNmjSRF198UWwXSB0tXrxYevbsKTVr1jQXo+rcubMsX75cbBfoduS1fv16iY6Olnbt2ontAq2j9PR0efDBB6Vhw4bmDJGmTZvKq6++KjYLtI7eeOMNadu2rVSoUEHq1Kkjw4cPl7S0NLHR2rVrpV+/fuYKunoF/KVLlxb4N27bX68NsI4ian+t97pC4BYsWODExMQ4L730kpOSkuKMHj3aqVixovPbb7/lOv+vv/7qVKhQwcyn8+vf6d8vWrTI2uoPtI709SeffNL56quvnJ9++smZMGGC+fuvv/7asVWgdeT1119/OU2aNHF69erltG3b1rFZUeqof//+TqdOnZwVK1Y4qampzpdffumsX7/esVWgdbRu3TqnTJkyzrPPPmv2Tfq8VatWzoABAxwbLVu2zHnwwQedd955R+/t6CxZsiTf+d24v14WYB1F0v6aoFNE559/vjNixIgs08455xxn/Pjxuc5/3333mdf93X777c4FF1zg2CrQOspNXFyck5iY6NiqqHU0ePBgZ+LEic4jjzxifdAJtI4+/PBDp0qVKk5aWprjFoHW0dSpU01Q9vfcc8859evXd2xXmIO4G/fXgdZRJO2v6boqgpMnT8qmTZukV69eWabr8w0bNuT6N59//nmO+Xv37i0bN26UjIwMsU1R6ii3G67qDduqVasmNipqHc2dO1d++eUXc+Eu2xWljt577z3p2LGjPPXUU1KvXj1p3ry53HPPPfL333+LjYpSR126dJFdu3bJsmXLzI0Rf//9d1m0aJFcfvnlpVTq8Oa2/XVxCOf9tetv6lkUBw4cEI/HI7Vr184yXZ/v27cv17/R6bnNf+rUKbM87SN3ex1lN23aNDl27Jhce+21YqOi1NHPP/8s48ePN+MvdHyO7YpSR7/++qt89tlnZmzFkiVLzDJGjhwpBw8etHKcTlHqSIOOjtEZPHiwnDhxwuyH+vfvL88//3wplTq8uW1/XRzCeX9Ni04QdMCWP/1llH1aQfPnNt3NdeT15ptvyqRJk2ThwoVSq1YtsVlh60gPZkOGDJHExETTSuEmgWxH+stSX9MD+fnnny99+/aV6dOny7x586xt1Qm0jlJSUuSuu+6Shx9+2LQGffTRR5KamiojRowopdKGPzfur4sq3PfX9v8kLAE1atSQsmXL5vi1tH///hy/ArzOPPPMXOfXX+XVq1cX2xSljrz0y3LLLbfI22+/LT169BBbBVpH2iysTefffPONjBo1yndQ1x2wbkcff/yxdO/eXdy+Hemvbe2yqlKlim9ay5YtTT1pd02zZs3E7XU0ZcoU6dq1q9x7773meZs2baRixYpy0UUXyWOPPeb6Fgu37a+DEQn7a1p0iqBcuXLmtMMVK1Zkma7PtUk4N3rqXfb59cCkYwliYmLENkWpI+8vg5tuuknmz59v/XiBQOtIT+H87rvvZPPmzb6H/gJv0aKF+XenTp3ENkXZjvQAvmfPHjl69Khv2k8//SRlypSR+vXri22KUkfHjx839eFPw5J/y4WbuW1/XVQRs78O9WjoSD+d85VXXjGnH44ZM8aczrl9+3bzup7tMHTo0BynK44dO9bMr39n++mKgdbR/PnznejoaCcpKcnZu3ev76GnUtsq0DrKzg1nXQVaR0eOHDFnD1199dXODz/84KxZs8Zp1qyZc+uttzq2CrSO5s6da75rM2fOdH755Rfns88+czp27GjO3rKRbhPffPONeehhb/r06ebf3tPv2V87AddRJO2vCTpB0A+4YcOGTrly5Zz27dubHarXsGHDnIsvvjjL/KtXr3bOO+88M3+jRo2cWbNmObYLpI703/oFy/7Q+WwW6HbktqBTlDraunWr06NHD+e0004zoWfcuHHO8ePHHZsFWkd6OrmeDqx1VKdOHef66693du3a5dho1apV+e5b2F87AddRJO2vo/Q/oW5VAgAAKAmM0QEAANYi6AAAAGsRdAAAgLUIOgAAwFoEHQAAYC2CDgAAsBZBBwAAWIugAwAArEXQgavonYeXLl0qkWj79u2m/Hpfq3Ch97kZMGBAvvNccsklMmbMGHGroUOHyhNPPCFuoXeJP+OMM4p9Owq2HC+88IL079+/WN8DkYGgA2vo3YbvvPNOadKkicTGxkqDBg2kX79+8sknn4gNdH327t0rrVu3DnpZq1atkiuuuEJq1qwp5cuXl6ZNm8rgwYNl7dq1AS3n2WefNQeU0qB3bb/mmmvMHbm1zM2bN5fbbrvN3LCzNK1evdoEzr/++qvAebds2SIffPCB2S7dQrejQD+T0tiOdFtJTk6Wzz77rETfB+GHoAMraGuH3sH5008/laeeesrc5fujjz6Sbt26SUJCgthA7y595plnSnR0dFDLmTlzplx66aVSvXp1WbhwoWzdulVef/11c6frsWPHBrSsKlWqBPzrvSjef/99ueCCCyQ9PV3eeOMNX5n1/R966CEJV9qKoOGsUqVK4gYZGRly2mmnSa1atcJuO9IfP0OGDJHnn3++RN8HYSjUN9sCikOfPn2cevXqOUePHs3x2p9//un7t27yL730kjNgwABzM8Ozzz7beffdd32vnzp1yrn55pvNTVfLly/vNG/e3JkxY0aW5elN66688kpn6tSpzplnnulUq1bNGTlypHPy5EnfPHv27HH69u1rlqHLeuONN8wNF5955hnfPHqX39tuu82pWbOmU6lSJadbt27O5s2b81zH1NRUU369o7D/TfhWrlzpdOjQwaxP586dnR9//DHPZeidiPUu12PHjs319czMzHxvGKrl1/XIXhdeWv96h2O9c7bWzdNPP21u/jd69GjfPOnp6c69997r1K1b16lQoYK5Y7auS16OHTvm1KhRw3xmufH/fPXGufHx8ebGlvr+999/v5ORkeF7PftnoHQddV0Ls414P4PC3MTQ4/E4Z5xxhvP+++9nma5lePzxx53hw4c7p59+utOgQQNn9uzZWebRm2tee+215u91++rfv795b7VlyxYnKirK+eOPP8zzgwcPmud6t3avJ554wrngggucYGj5ExMTzfdK61Pr6cMPP/S97q2LhQsXms84NjbWefXVV82d0atUqZJlWY8++qjZznV9b7nlFvO5+G9b2bcjXd6dd95ptpOqVas6tWvXzvIZqWnTpjmtW7c225DeuPWOO+4wd+D2yq0cun3outh+g1dkRYsOIt7BgwdN64223FSsWDHH69l/KSYmJsq1115ruhX69u0r119/vVmGyszMlPr168tbb70lKSkp8vDDD8sDDzxgnmfv+vnll1/M///zn/+YZnf/pvcbb7xR9uzZY7o53nnnHZkzZ47s37/f97oeTy+//HLT3bZs2TLZtGmTtG/f3rS0eMtSWA8++KBMmzZNNm7caFp7br755jzn1bLor+777rsv19e1SyYY9957r6mTJUuWyMcff2zWX9fN3/Dhw2X9+vWyYMEC8xloi8dll10mP//8c67LXL58uRw4cCDPMns/3927d5vPMz4+Xr799luZNWuWvPLKK/LYY48FvB55bSPafah1qP73v/+ZrkTtdsmN/q12b3Xs2DHHa/p56XTtjhs5cqTccccd8uOPP5rXjh8/bloiTz/9dNOVqF0t+m+to5MnT5quS22NW7NmjZlf59Hn/t2OWu8XX3yxBEPXS8v59NNPm3Xp3bu3GeOS/XO6//775a677jKtbDpPdtoC9/jjj8uTTz5ptoWzzjrLfDYF0e+Vfp+//PJL00o7efJkWbFihe/1MmXKyHPPPSfff/+9mVdbc/PaRry0znX7/+qrrwKqC0S4bMEHiDhffvml+WW5ePHiAufV+SZOnJilBUJ/Dfv/Us1OW2sGDRqU5den/irX1h+va665xhk8eLD599atW837JCcn+17/+eefzTRva8Inn3ziVK5c2Tlx4kSW92ratGmOX/eFadHx+uCDD8y0v//+O9dljBgxwryvv0WLFpkWGO9DWwyK0qKjv6b11/KCBQt8r6elpZlWEW+LzrZt20x97969O8tyL730UmfChAm5lvnJJ58066QtF/l54IEHnBYtWmRplUpKSjKtCNo6EUiLTn7biLfe/VuScrNkyRKnbNmyWcrjLcMNN9zge66v16pVy5k1a5Z5/sorr+RYD20F03pcvny5eT5w4EBn1KhR5t9jxoxx7r77btPq9cMPP5gWLF3n/LbpwtAWN2158qetZfp98N8es7d4Zm9J6dSpk5OQkJBlnq5duxbYonPhhRfmeG9tCcrLW2+95VSvXj3PcnhpC9G8efPyXXfYhRYdRLz/f2wqfGtEmzZtfP/WX4w6fsK/teXFF180v/x0oK7+kn7ppZdkx44dWZbRqlUrM2bGq06dOr5l6C99bVnRFhqvs88+W6pWrep7rr9sjx49an6J63t4H6mpqaalKBD+66PlUP7rk132etJf4Xomlw6aPXbsmHg8HikKLbe2OHTu3Nk3rVq1atKiRQvf86+//tp8XjqQ2H+9tXUir/X2fr4F0RYFfW//9evataup5127dgW0LgVtI4Xx999/m3EhuW2X/svX13XslXf5um1s27bNvKe3frQeT5w44asjPZNNW22U1p22AP3jH/8w/9YBt/reuu650TPA/Os++7atDh8+bFoksy9Dn2s9+8utxcqffh/OP//8LNOyP8+Nfx1l/44pbTns2bOn1KtXz9SVtqKmpaWZbTg/OoZIW83gHsGNagTCQLNmzczBQnfAhTlFNSYmJstz/VvtslLaRaUDcrXJXg+augOdOnWqaT4v7DLyOjD7T9d5dcftPVj5C3RQpn9ZvAdVb1lyq6tDhw6ZLjM9uCo92GkQyz7IWbsGsq+LNvvnpTCBRMulAVEP5v5B0VuO3GgoUtq14x+icnv/7KEiewgu7Drl9/kWVo0aNcwBVcNfuXLlCr18/b8OrNcun+w0fHuDzujRo00g0q6biy66yIQgDTraXaZ/n9cA6BEjRphuOa+6devmuQ651Wf2abl1FxdmOQXJr45+++0306Wo6/Loo4+aIKhdfLfccku+26jSLkhvPcIdaNFBxNOdnLZKJCUl5fprrjCnAXutW7fOnH2k4ybOO+88EwACbWE555xz5NSpU2b8hZcekPzLoa09GjY0XOh7+D/0AFlSrr76anMA0fESBdGDgZbR/6CU3zV8tOy67C+++MI37c8//8xyqrHWqbYY6S/z7OvtDV7Z9erVy9SJjtPIjbde4+LiZMOGDVnKq8/1gK+/+r3rpONq/FsutBUtEN7QUlDLV7t27cz/daxXIHTb0HEweuZS9jrSs5OUd5yOjj9q27atVK5c2YzJ0aBT0Pgc/b74LzO3s/h0eRqAsp+KrfXZsmXLgNZHW/Syj4nR8WTB0L/X75j+INGz8TQMawtUQfS7rC1juh3CPQg6sIKeMq0HHm0S18GieqDQFh4drJhfK0B2uuPXnagOgNUDtJ66rF0BgQadHj16yL/+9S+zg9fAo//WJnPvL1t9XculLVD6Xnp6vB5EJk6cGPRBID86EFQPDjrQdNiwYab5X99bu5S0rpS3pUVbDf744w8TMPQAoUHyww8/zHPZ2iKjv6h1QLJeu0hbGvRCcNqK4qUHJB3Yq90MixcvNiFD61eDlw7Kzo22GLz88suma00Hw65cudKUWetJB5/qr3ql4XTnzp3mmjXa+vPuu+/KI488IuPGjfOVoXv37ua0dA20Wj6tg+wtSwVp2LCh+Rz1lHetH+0ay42GKg0tgV63RetHg92VV15pyql1pAFGW3C8XXD6/tpV9d///td8Tt6uHm090rr3TguGfo76ueglCLT7afz48SboajkCoZ+HDgrXAcP6vdRwpoObgxn4rtd90qCjp4r/+uuv5jPVLueCaH3qdbb07+EeBB1YoXHjxuZgrWMV7r77bvOLV/vvdadfmDM8vPSgOXDgQHPRs06dOpk+fz2ABuq1114zF7bTg9FVV11lLlamLQt6oTulO3k9sOvrepaUBoDrrrvOHMD170qSHnj0jCg9SGsLj3ZnaTeAHlD17LVzzz3XzKe/3DVAasDRVgMNbffcc0++y9ZuPl0nDSQa5i688ELTjeJv7ty5Jujo56S/9nVe7RrUM5ryogd9DYLaYqTXQtEw+c9//tN0w3nPqtJWG61TLaeWVz9LDV4aHr0mTJhgyqcXS9R11qAZ6EFP30fPytIDv35Wo0aNynNeDbi5dUHlp0KFCuYMKg2lui3q56DbiI670ZYWL93WNdx7Q41uU9qFpbTeg6VnUulnpA/dJnTbeO+998z2Emhw03rXbUeDn25nGoC934Wi0Nay6dOnmyCm33Wt4ylTphT4d2+++ab5LsJdonREcqgLAdhOf4nrgVxbI/QUcriDdpNomNNT6QNpWbSd/gjRrkptiSkt2oKn3z1tqfV2AcIdGIwMlAC9pod2aegvYR0Tol0sjRo1Mq0JcA9ttdDWPb0OkFvpgGztVtJxdNpNqK0qGvj9r4lTGnQMj34WhBz3oUUHKAE67kab/HX8gHZZ6QDnGTNmmPEdgJtol5vec067lvUWHtrCpd2J2i0HlAaCDgAAsBaDkQEAgLUIOgAAwFoEHQAAYC2CDgAAsBZBBwAAWIugAwAArEXQAQAA1iLoAAAAsdX/AwBM2JZvsifsAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "plt.hist(dat[\"guide_count_2\"] - dat[\"guide_count_1\"], bins=20, log=True);\n", + "plt.xlabel('Change in Guide Count (new - original)')\n", + "plt.ylabel('Count')\n", + "plt.grid(True)" + ] + }, + { + "cell_type": "code", + "execution_count": 43, + "id": "e42d1627", + "metadata": {}, + "outputs": [], + "source": [ + "import re\n", + "def get_n_optimizations(aca):\n", + " for event in aca.log_info[\"events\"]:\n", + " m = re.search(r'.*N opt runs=(\\d+)', event['data'])\n", + " if m:\n", + " return int(m.group(1))\n", + " return 0\n" + ] + }, + { + "cell_type": "markdown", + "id": "433e1580", + "metadata": {}, + "source": [ + "Here we've again scanned all of the new pkl catalogs and checked to see which actually got to catalog \"optimization\" *caused* by this PR. While doing that review, this also counts up the catalogs with optimization that is not new and looks for weird cases (optimized without this PR and not optimized with this PR code)." + ] + }, + { + "cell_type": "code", + "execution_count": 44, + "id": "1b27c96a", + "metadata": {}, + "outputs": [], + "source": [ + "guide_opt_counts = []\n", + "opt_counts = []\n", + "weird_counts = []\n", + "for load_name in load_names:\n", + " pkl1 = Path(f\"{dir1}/{load_name}_acas.pkl.gz\")\n", + " pkl2 = Path(f\"{dir2}/{load_name}_acas.pkl.gz\")\n", + " if not pkl1.exists():\n", + " print(f\"Skipping {load_name} missing pkl1\")\n", + " continue\n", + " if not pkl2.exists():\n", + " print(f\"Skipping {load_name} missing pkl\")\n", + " continue\n", + " acas2 = pickle.load(gzip.open(pkl2, \"rb\"))\n", + " acas1 = pickle.load(gzip.open(pkl1, \"rb\"))\n", + " for obsid in acas2.keys():\n", + " if obsid > 38000:\n", + " continue\n", + " new_aca1 = acas1[obsid]\n", + " new_aca2 = acas2[obsid]\n", + " n_opt1 = get_n_optimizations(new_aca1)\n", + " n_opt2 = get_n_optimizations(new_aca2)\n", + " if n_opt2 > 1 and n_opt1 <= 1:\n", + " guide_opt_counts.append((load_name, obsid, n_opt2))\n", + " if n_opt2 > 1 and n_opt1 > 1:\n", + " opt_counts.append((load_name, obsid, n_opt1, n_opt2))\n", + " if n_opt1 > 1 and n_opt2 <= 1:\n", + " weird_counts.append((load_name, obsid, n_opt1))\n" + ] + }, + { + "cell_type": "markdown", + "id": "869054ed", + "metadata": {}, + "source": [ + "It looks like there are no weird cases" + ] + }, + { + "cell_type": "code", + "execution_count": 45, + "id": "7b9e882a", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[]" + ] + }, + "execution_count": 45, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "weird_counts" + ] + }, + { + "cell_type": "markdown", + "id": "ed9d265c", + "metadata": {}, + "source": [ + "It looks like the \"worst\" case was in obsid 28280 where with the new guide count and fid spoiling information the optimization needed to run 10 different fid sets to get to one that either satisfied the requirements or had the closest P2 to the no-fid P2. Given the frequency of needing to dig that deep into optimization, this seems reasonable.\n", + "\n", + "The list of tuples is week, obsid, number of optimization attempts." + ] + }, + { + "cell_type": "code", + "execution_count": 46, + "id": "a92025a6", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[('JUN1024A', np.float64(26571.0), 2),\n", + " ('JUN1124A', np.float64(26571.0), 2),\n", + " ('AUG1224A', np.float64(28280.0), 10),\n", + " ('OCT0724A', np.float64(28382.0), 2),\n", + " ('DEC0924A', np.float64(30674.0), 2),\n", + " ('DEC1024A', np.float64(30674.0), 2),\n", + " ('JUL2125A', np.float64(31009.0), 4)]" + ] + }, + "execution_count": 46, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "guide_opt_counts" + ] + }, + { + "cell_type": "markdown", + "id": "7930b8e9", + "metadata": {}, + "source": [ + "The fraction of *extra* optimization due to the new information seems reasonable. Below are just the optimization counts for obsids before the guide star information was used. The data structure there is just week, obsid, number of previous optimization attempts, number of current optimization attempts." + ] + }, + { + "cell_type": "code", + "execution_count": 47, + "id": "17abd780", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[('AUG2123A', np.float64(24782.0), 4, 4),\n", + " ('AUG2123A', np.float64(28488.0), 4, 4),\n", + " ('FEB2624A', np.float64(26650.0), 4, 4),\n", + " ('JUL0124A', np.float64(29472.0), 4, 4),\n", + " ('SEP2324A', np.float64(29899.0), 2, 2),\n", + " ('DEC0224A', np.float64(27983.0), 4, 4),\n", + " ('DEC0224A', np.float64(27979.0), 4, 4),\n", + " ('FEB0325A', np.float64(25501.0), 2, 2),\n", + " ('FEB0325A', np.float64(25505.0), 4, 4),\n", + " ('FEB2425A', np.float64(30284.0), 4, 3),\n", + " ('JUN2325A', np.float64(28777.0), 4, 4),\n", + " ('JUL1425A', np.float64(29577.0), 2, 2),\n", + " ('AUG0425A', np.float64(30446.0), 4, 4),\n", + " ('SEP1525A', np.float64(30136.0), 3, 3)]" + ] + }, + "execution_count": 47, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "opt_counts" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "af889fd1", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f4384977", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "ska3-flight-2026.1rc1", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.13.10" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}