---
title: ASTR 597 W23 HW 4
author: Andy Tzanidakis and friends
tags: LSST Class
---

[TOC]

# Homework 4 - ASTR 597 W23
- Author: Anastasios (Andy) Tzanidakis
- Course ID: ASTR 597 W23

# Prelude & Setup

Let's confirm that rubin_sim is working and that we can use MAF on the current baseline cadence.

Let's find a basic metric to run. `CoaddM5` is provides the final coadded depth of the ten-year survey. Let's make plots of this metric for the `r` band as a spatial (healpix) map as well as a 1-d histogram, and then compute some summary statistics.

# choose the spatial slicer
slicer = maf.slicers.HealpixSlicer(nside=64)

# choose the metric
metric =  maf.metrics.Coaddm5Metric()

# only r-band
constraint = "filter = 'r'"

# choose how to plot"
plot_funcs = [maf.plots.HealpixSkyMap(), maf.plots.HealpixHistogram()]
plot_dict = {'nside': 64, 'colorMin': 0}

# define the summary metrics
summary_metrics = [maf.metrics.MinMetric(), maf.metrics.MedianMetric(), maf.metrics.MaxMetric(), 
                 maf.metrics.PercentileMetric(percentile=25), maf.metrics.PercentileMetric(percentile=75)]

# wrap it up
bundle = maf.metricBundles.MetricBundle(metric, slicer, constraint=constraint, run_name=run_name,
                                        plot_dict=plot_dict, plot_funcs=plot_funcs,
                                        summary_metrics=summary_metrics)

out_dir = 'example'

bdict = {'example':bundle}
bgroup = maf.metricBundles.MetricBundleGroup(bdict, opsim_fname, out_dir=out_dir)
bgroup.run_all()
bgroup.plot_all(closefigs=False)

bdict['example'].summary_values


# Problem 1

## Problem 1a-Counts

Let's get a sense of how the LSST divides its exposures. Use `UniSlicer`, `CountMetric`, and appropriate SQL constraints to provide counts of the total exposures by filter. The OpSim schema is described [here](https://rubin-sim.lsst.io/rs_scheduler/output_schema.html).

We will use the `UniSlicer` slicer and `CountMetric` to count the number of exposures per filter:

example1_slicer = maf.UniSlicer()
example1_metric = maf.CountMetric(col="numExposures")

n = []
for flts in tqdm(list('ugrizy')):

    example1_constraint = f"filter = '{flts}'"


    example1_bundle = maf.MetricBundle(
        example1_metric,
        example1_slicer,
        example1_constraint,
        run_name=run_name,
    )

    example1_bg = maf.MetricBundleGroup(
        [example1_bundle], opsim_fname, 
        out_dir='.'
    )

    example1_bg.run_all()

    n_exp = example1_bundle.metric_values.data[0]
    n.append(n_exp)

Below, we can see the number of exposures across each filter - with the u-band having the least amount of exposures. Overall, when we sum this up we find that LSST will deliver over 2 million exposures across all filters.


# Problem 2

## Problem 2a - Cloud Coverage

What is the mean cloud cover for exposures taken during the simulated baseline survey? Plot a histogram of the cloud cover.

# choose the spatial slicer
slicer = maf.slicers.UniSlicer()

# choose the metric
bins = np.arange(0, 1, step=0.1) # cloud fraction cover
metric =  maf.metrics.HistogramMetric(bins=bins, bin_col="cloud", col="cloud", statistic="count")

out_dir = 'rubin_sim_outputs'

# wrap it up
bundle = maf.metricBundles.MetricBundle(metric, slicer, constraint=None, run_name=run_name,
                                        summary_metrics=[maf.metrics.MedianMetric()])

bdict = {'q2':bundle}
bgroup = maf.metricBundles.MetricBundleGroup(bdict, opsim_fname, out_dir=out_dir)
bgroup.run_all()

bc = [0.5*(bins[i]+bins[i+1]) for i in range(0, len(bins)-1)] # bin centers

# hacky way to evaluate the median given the number of exposures
np.median(np.concatenate([np.zeros(shape=1648774)+0.05, 
                         np.zeros(shape=261497)+0.15, 
                         np.zeros(shape=171478)+0.25]))

Below we can see that overall the median sky coverage per visit. Overall, we find that the median cloud coverage fraction to be about 5% for any given exposure. Not bad -- looks like the weather overall in Cerro Pachón Ridge is pretty good! :+1:


# Problem 3

## Problem 3a - Metric Mesh

MAF was built to compare evaluate the impact of different observing strategies on various science cases. You can always write new metrics and run them yourself on multiple OpSim simulations. This can be time consuming, though. For many standard metrics, the metric values for various runs are stored in queryable tables--see the [04_Getting_Data](https://github.com/lsst/rubin_sim_notebooks/blob/main/maf/tutorial/04_Getting_Data.ipynb) tutorial notebook.

In this exercise we'll get a small taste of the tradeoffs confronted by the Survey Cadence Optimization Committee.

Make a mesh plot of the metric summary values for the metric sets `"DESC WFD"` (Cosmology), `"TVS XRB"` (Galactic Transients), and `"SSO discovery"` (Solar System science) on run families `"baseline"`, `"galactic plane footprint"`, amd `"rolling"`. Set `baseline_run='baseline_v2.0_10yrs'` to compare to one of the earlier baseline cadences. Qualitatively, how does the new (v3.0) baseline cadence look for these metrics? What other trends do you see across these run families?

We can use the code directly given to us from the 04_Getting_data tutorial. We will use the `baseline`, `galactic_plane footprint`, and `rolling` families, and metric set: `["DESC WFD", "TVS XRB", "SSO discovery"]`. In practice this is how we generate the following diagrams:

def plot_summary_mesh(summary, metric_sets, run_label_map, metric_label_map, baseline_run="baseline_v3.0_10yrs"):
    """Plot the summary metric mesh given summary, baseline_run, metric_sets, run_labels, and metric label maps."""
    
    maf.plot_run_metric_mesh(summary,
                             baseline_run=baseline_run, 
                             metric_set=metric_sets.loc[this_metric_set], run_label_map=run_label_map,
                             metric_label_map=metric_label_map)
# family runs & metric sets
family_runs = maf.get_family_runs()
metric_sets = maf.get_metric_sets()

# metric list
metric_list = ['SSO discovery', 'TVS XRB', 'DESC WFD']

# family list
fam_list = ['baseline', 'galactic plane footprint', 'rolling']

summary = maf.get_metric_summaries(this_family, this_metric_set)

# Plotting 
label_map = family_runs.loc[this_family, ["run", "brief"]].set_index("run")["brief"]

metric_map = metric_sets.loc[this_metric_set, "short_name"].droplevel(
    "metric set")

# v3.0
plot_summary_mesh(summary, metric_sets, family_runs,
                  label_map, metric_map, baseline_run="baseline_v3.0_10yrs")

# v2.0
plot_summary_mesh(summary, metric_sets, family_runs, 
                  label_map, metric_map, baseline_run="baseline_v2.0_10yrs")

Below we show the mesh summary plot of baselines v.3.0 and v.2.0:

- Overall `v.3.0` is better for transients (i.e SNe rates) and cosmological metrics (i.e coadded depths), while the solar system metrics seem to be slightly impacted compared to `v.2.0`

- We see that the the largest change between the two baselines is affected mostly by the number of SN Ia rates and super-luminous SNe (or broadly the transients)


# Problem 4

## Problem 4a - Kilonovae Light Curves

Some metrics take simulated events and generate the data points as they would be observed by LSST. Typically these lightcurves are only used internally to compute higher-level metrics, but some metrics make it possible to return the generated lightcurves.

Use the `KNePopMetric` with `output_lc=True` to generate random kilonova lightcurves as observed by LSST. Plot an event with at least ten data points.

The [kilonova](https://github.com/lsst/rubin_sim_notebooks/blob/main/maf/science/KNeMetric.ipynb) and [XRB](https://github.com/lsst/rubin_sim_notebooks/blob/main/maf/science/XRB_Metric.ipynb) science notebooks will be helpful.

We follow the tutorial from the XRB and KNe metrics to start with this problem. We will inject 10,000 simulated KNe events across the sky. We will also have a simple function to read all the light curve photometry and removed the masked detections:

baseline_file = get_baseline()
opsim = os.path.basename(baseline_file).replace('.db','')

summaryMetrics = [maf.SumMetric(metric_name='Total detected'), 
                  maf.CountMetric(metric_name='Total lightcurves in footprint'),
                  maf.CountMetric(metric_name='Total lightcurves on sky', mask_val=0),
                  maf.MeanMetric(metric_name='Fraction detected in footprint'), 
                  maf.MeanMetric(mask_val=0, metric_name='Fraction detected of total')]

n_events = 10_000

# Kilonova parameters - set up to run on a particular subset of the KNe models
# All models
filename = maf.get_kne_filename(None)
slicer_all = maf.generate_kn_pop_slicer(n_events=n_events, n_files=len(filename), d_min=10, d_max=600)
metric_all = maf.KNePopMetric(output_lc=True, file_list=filename, metric_name='KNePopMetric_all')

bundle_all = maf.MetricBundle(metric_all, slicer_all, None, 
                              info_label='all models', run_name=opsim, summary_metrics=summaryMetrics)


bundle_all = maf.MetricBundle(metric_all, slicer_all, None, 
                              info_label='all models', run_name=opsim, summary_metrics=summaryMetrics)

g.run_all()

# Mask out the masked indicies
bad = ~bdict['all'].metric_values.mask # masked out values
all_lcs = bdict['all'].metric_values[bad]

# Fetch all the light curves
lcs = [get_lc(v) for v in all_lcs.tolist()]


def fetch_lc(lc_index, lc_lib=all_lcs):
    """Fetch the light curve given the light curve index. 
    
    Parameters
    ----------
    lc_index (int): light curve index
    lc_lib (list): list of all light curves
    
    Returns 
    -------
    light curve (time, magnitude, magitude error, filters)
  
    """
    lc = lc_lib[lc_index]
    t, mags, flts = lc['lc'][0], lc['lc'][1], lc['lc'][-1]
    mag_err = lc['lc'][2]
    rm = mags<99
    t, mags, flts = t[rm], mags[rm], flts[rm]
    mag_err = mag_err[rm]
    
    return t, mags, mag_err, flts

Here I will also search through the light curve indexes to find the ones that match the conditions we're looking for in our simulated KNe light curve:
- Number of detections >10
- Brightest magnitude (in any filter) < 21 mag
- Overall duration of event <30 days

We can also imagine (if we wanted to do a more serious filtering criterion) to estimate roughly the timescale of the KNe by measuring the dm/dt or fitting a slope from the brightest detection to the fading section of the light curve.

# Let's sort through and look through some good light curves
good_lcs_index = [] # list of good quality KNe light curves. 
for j in tqdm(range(len(lcs))):
    # select
    _lc_ = all_lcs[j]
    p_det = all_lcs[j]['multi_detect'] #
    
    if p_det>0:
        t, mags, flts = _lc_['lc'][0], _lc_['lc'][1], _lc_['lc'][-1]
        rm = mags<99
        t, mags, flts = t[rm], mags[rm], flts[rm]
        
        if len(flts)>15 and min(mags)<21 and (t[np.argmax(mags)]-t[np.argmin(mags)])<30:
            good_lcs_index.append(j)

Voila! Below we have a simulated KNe light curve that was caught early enough that has both a rise and a fading detections. Just as we expected, the KNe will be very fast fading, for this particular model within 5 days of peak it will likely no longer visible with single visit exposures with LSST :female_mage:


# Acknowledgements
Many concepts and ideas for the workflow of this project were inspired by conversations and code I shared with Tom Wagg, David Wang, Tobin Wainer, and Jake Kurlander.