<jupyter_start><jupyter_text>Develop Log_dt Metric<NAME> (Drexel), <NAME> (Drexel), with input from <NAME> (Yale); August 2020 This notebook is an exploration of the distribution of pairs of observations.  We are comparing DDF observations to WFD observations in the same part of the sky and comparing 4 different opSims created with the DDFs specifically in mind (or that serve as a baseline for them).<jupyter_code>#Please enter your SciServer username between the single quotes below!
your_username = ''<jupyter_output><empty_output><jupyter_text>Now we can start computing metrics to use for comparing different opSims. --> 0. Software Setup (repeat each session)<jupyter_code># import matplotlib to show plots inline.
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import matplotlib as mpl
import os
pgf_with_latex = {                      # setup matplotlib to use latex for output
    "axes.linewidth":1.2,               # width of box, 2 is too wide, 1 is too narrow
    "pgf.texsystem": "pdflatex",        # change this if using xetex or lautex
#     "text.usetex": True,                # use LaTeX to write all text
    "font.family": "serif",
    "font.serif": [],                   # blank entries should cause plots to inherit fonts from the document
    "font.sans-serif": [],
    "font.monospace": [],
    "axes.labelsize": 20,               # LaTeX default is 10pt font, font size of axis text label
    "axes.labelpad" : 6,                # Distance between label and axis    
    "axes.formatter.limits":[-99,99],   # use sci notation if log10 of axis range is smaller than first or larger than second.
                                        # GTR: Actually *don't* -- should change the axis label instead.  E.g., "Flux Density (10^-17 ergs/s/cm^2)" 
                                        # This is a hack b/c there doesn't seem to be an rcParams version of
                                        # axes.ticklabel_format(style='plain')
    #"axes.formatter.style":"plain",    # Turn off multiplicative offsets (sci notation) to the axes [GTR: Doesn't work]
    "axes.formatter.useoffset":False,   # Turn off additive offsets to the axes
    "font.size": 20,
    "legend.fontsize": 20,              # Make the legend/label fonts a little smaller
    "xtick.labelsize": 20,              # Font size of numbers 
    "ytick.labelsize": 20,
    "xtick.direction": "in",
    "ytick.direction": "in",
    "xtick.minor.visible": True,
    "ytick.minor.visible": True,
    "xtick.major.top":True,
    "xtick.minor.top":True,
    "ytick.major.right":True,
    "ytick.minor.right":True,
    'xtick.major.width':1.0, 
    'xtick.minor.width':1.0, 
    'ytick.major.width':1.0, 
    'ytick.minor.width':1.0, 
    'xtick.major.size':8,             # size of tickmarks in points
    'xtick.minor.size':4, 
    'ytick.major.size':8, 
    'ytick.minor.size':4,
    'xtick.major.pad':8,               # distance between box and numbers
    'ytick.major.pad':8,
#     "figure.figsize": figsze(1,1),     # default fig size of 0.9 textwidth
#     "pgf.preamble": [
#         r"\usepackage[utf8x]{inputenc}",    # use utf8 fonts because your computer can handle it
#         r"\usepackage[T1]{fontenc}",        # plots will be generated using this preamble
#         ]
    }

mpl.rcParams.update(pgf_with_latex)<jupyter_output><empty_output><jupyter_text>Import the sims_maf modules needed.<jupyter_code># import our python modules
import lsst.sims.maf.db as db
import lsst.sims.maf.metrics as metrics
import lsst.sims.maf.slicers as slicers
import lsst.sims.maf.stackers as stackers
import lsst.sims.maf.plots as plots
import lsst.sims.maf.metricBundles as metricBundles
# add opsimUtils module path to search
import sys
sys.path.insert(0, '../Scripts_NBs/')
from opsimUtils import *<jupyter_output><empty_output><jupyter_text>1. Define `LogTgaps` Metric<jupyter_code># the base 10 log of 30sec and 10 years in dasy are
# but the bins really aren't needed for the metric, the metric only return Tgaps. 
np.log10([30/3600/24, 3650])
from lsst.sims.maf.metrics import BaseMetric

class LogTgapsMetric(BaseMetric):
    """Histogram the times of the gaps between observations.


    Measure the gaps between observations.  By default, only gaps
    between neighboring visits are computed.  If allGaps is set to true, all gaps are
    computed (i.e., if there are observations at 10, 20, 30 and 40 the default will
    return a histogram of [10,10,10] while allGaps returns a histogram of [10,10,10,20,20,30])

    Parameters
    ----------
    timesCol : str, opt
        The column name for the exposure times.  Values assumed to be in days.
        Default observationStartMJD.
    allGaps : bool, opt
        Histogram the gaps between all observations (True) or just successive observations (False)?
        Default is False. If all gaps are used, this metric can become significantly slower.
    bins : np.ndarray, opt
        The bins to use for the histogram of time gaps (in days, or same units as timesCol).
        Default values are bins from 0 to 2 hours, in 5 minute intervals.

    Returns a histogram at each slice point; these histograms can be combined and plotted using the
    'SummaryHistogram plotter'.
     """

    def __init__(self, timesCol='observationStartMJD', allGaps=False, bins=np.logspace(-3.46, 3.56, 100),
                 units='days', **kwargs):
        # Pass the same bins to the plotter.
        self.bins = bins
        self.timesCol = timesCol
        super(LogTgapsMetric, self).__init__(col=[self.timesCol], metricDtype='object', units=units, **kwargs)
        self.allGaps = allGaps

    def run(self, dataSlice, slicePoint=None):
        if dataSlice.size < 2:
            return self.badval
        times = np.sort(dataSlice[self.timesCol])
        if self.allGaps:
            allDiffs = []
            for i in np.arange(1,times.size,1):
                allDiffs.append((times-np.roll(times,i))[i:])
            dts = np.concatenate(allDiffs)
        else:
            dts = np.diff(times)
#         result, bins = np.histogram(dts, self.bins)
        return dts<jupyter_output><empty_output><jupyter_text>Load up the 4 opSims that we want to compare.  The AGN SC's proposal, the DESC proposal, the "Daily" DDF opSim, and "baseline" (for comparison).  With just these 4 I haven't bothered to set things up as a loop over them.<jupyter_code>#Read in an opSim (in this case "baseline_v1.4_10yrs.db")
opsdb1 = db.OpsimDatabase('/home/idies/workspace/lsst_cadence/FBS_1.5/baseline_v1.5_10yrs.db')
opsdb2 = db.OpsimDatabase('/home/idies/workspace/lsst_cadence/FBS_1.5/daily_ddf_v1.5_10yrs.db')
opsdb3 = db.OpsimDatabase('/home/idies/workspace/lsst_cadence/FBS_1.5/agnddf_v1.5_10yrs.db')
opsdb4 = db.OpsimDatabase('/home/idies/workspace/lsst_cadence/FBS_1.5/descddf_v1.5_10yrs.db')

if your_username == '':
    raise Exception('Please provide your username!  See the top of this notebook.')
    
outDir = '/home/idies/scr/tmp/MAFOutput/LogTGaps/ResultDbs/'
if not os.path.exists(outDir):
    os.makedirs(outDir)
    
resultsDb = db.ResultsDb(outDir=outDir)
#The DDFs aren't recorded in the same way for each opSim, so we need to figure out which "proposal" they correspond to.

print(ddfInfo(opsdb1, 'XMM-LSS'))
print(ddfInfo(opsdb2, 'XMM-LSS'))
print(ddfInfo(opsdb3, 'XMM-LSS'))
print(ddfInfo(opsdb4, 'XMM-LSS'))<jupyter_output>{'proposalId': 6, 'Coord': (35.707, -4.72)}
{'proposalId': 5, 'Coord': (35.707, -4.72)}
{'proposalId': 5, 'Coord': (35.707, -4.72)}
{'proposalId': 4, 'Coord': (35.707, -4.72)}<jupyter_text>If we wanted to look at the COSMOS field instead of XMM-LSS, then run this cell instead.print(ddfInfo(opsdb1, 'COSMOS'))print(ddfInfo(opsdb2, 'COSMOS'))print(ddfInfo(opsdb3, 'COSMOS'))print(ddfInfo(opsdb4, 'COSMOS'))  2. Run Metric at the center of the DDF using pointslicer<jupyter_code># metric = is the number of pairs of observations, distributed in log space, thus "log t-gaps", which is based on one of the project-defined metrics.
metric1 = LogTgapsMetric(allGaps=True)

# slicer = a grouping or subdivision of visits for the simulated survey
# based on their position on the sky (using a Healpix grid)
#COSMOS
#slicer1 = slicers.UserPointsSlicer(ra=[150.11], dec=[2.14])
#XMMLSS
slicer1 = slicers.UserPointsSlicer(ra=[35.707], dec=[-4.72])

# constraint = the sql query (or 'select') that selects all visits in r band
constraint1 = 'filter = "r"'<jupyter_output><empty_output><jupyter_text>Now set up the metric bundles.  Note that I made one for each opSim.  It is notclear that I need to do this way as I'm not changing the metric, but this is also where we extract the results from, so I think that I have to?<jupyter_code># MetricBundle = combination of the metric, slicer, and sqlconstraint
ddf_nvisit1 = metricBundles.MetricBundle(metric1, slicer1, constraint1)
bundleDict1 = {'logtgaps': ddf_nvisit1}
# MetricBundle = combination of the metric, slicer, and sqlconstraint
ddf_nvisit2 = metricBundles.MetricBundle(metric1, slicer1, constraint1)
bundleDict2 = {'logtgaps': ddf_nvisit2}
# MetricBundle = combination of the metric, slicer, and sqlconstraint
ddf_nvisit3 = metricBundles.MetricBundle(metric1, slicer1, constraint1)
bundleDict3 = {'logtgaps': ddf_nvisit3}
# MetricBundle = combination of the metric, slicer, and sqlconstraint
ddf_nvisit4 = metricBundles.MetricBundle(metric1, slicer1, constraint1)
bundleDict4 = {'logtgaps': ddf_nvisit4}<jupyter_output><empty_output><jupyter_text>Now run the metric on each of the opSims defined above.<jupyter_code>group1 = metricBundles.MetricBundleGroup(bundleDict1, opsdb1, outDir=outDir, resultsDb=resultsDb)
group1.runAll()
group2 = metricBundles.MetricBundleGroup(bundleDict2, opsdb2, outDir=outDir, resultsDb=resultsDb)
group2.runAll()
group3 = metricBundles.MetricBundleGroup(bundleDict3, opsdb3, outDir=outDir, resultsDb=resultsDb)
group3.runAll()
group4 = metricBundles.MetricBundleGroup(bundleDict4, opsdb4, outDir=outDir, resultsDb=resultsDb)
group4.runAll()<jupyter_output>Querying database SummaryAllProps with constraint filter = "r" for columns ['observationStartMJD', 'fieldDec', 'fieldRA']
Found 473116 visits
Running:  ['logtgaps']
Completed metric generation.
Running reduce methods.
Running summary statistics.
Completed.<jupyter_text>Now make a plot comparing the pair distribution in log delta T for each of those opSims.<jupyter_code>mask1 = ddf_nvisit1.metricValues.mask
data1 = ddf_nvisit1.metricValues.data[~mask1]
mask2 = ddf_nvisit2.metricValues.mask
data2 = ddf_nvisit2.metricValues.data[~mask2]
mask3 = ddf_nvisit3.metricValues.mask
data3 = ddf_nvisit3.metricValues.data[~mask3]
mask4 = ddf_nvisit4.metricValues.mask
data4 = ddf_nvisit4.metricValues.data[~mask4]
fig = plt.figure(figsize=(10,8))
plt.hist(np.log10(data1[0]), bins=50, log=True, histtype='step', label='baseline')
plt.hist(np.log10(data2[0]), bins=50, log=True, histtype='step', label='daily')
plt.hist(np.log10(data3[0]), bins=50, log=True, histtype='step', label='agn', linewidth=2)
plt.hist(np.log10(data4[0]), bins=50, log=True, histtype='step', label='desc')
plt.xlim(-3, 3.5)
plt.ylim(bottom=1)
plt.xlabel('log10(days)')
plt.ylabel('Number of pairs')
plt.legend(loc='upper left')
#plt.savefig('XMMLSS_r_manydbs.png')<jupyter_output><empty_output><jupyter_text>OK, that was the result where we took all the observations in the DDF, regardless of whether they were "DDF" or "WFD" observations.  Now let's split them up.The point is this: I'd like to define a metric in the WFD area.  However, no matter what we do there are going to be gaps and some timescales with inherently fewer (but non-zero) observations.  So, I'm thinking that it would be a good idea to scale the result from the WFD by the DDFs (as they represent the "best" that we could hope to achieve).  That would remove the upward-right sloping aspect of the log delta t plot above and maybe make it easier to compare different WFD results.  I was thinking of averaging over all of the DDFs for the baseline opSim as the vector to normalize by.  A further question is how to normalize the WFD and the DDF results since clearly there are many more observation pairs in the DDFs. Modify the above contraint to keep only the **WFD** observations ("not like DDF").<jupyter_code># metric = the "maximum" of the "airmass" for each group of visits in the slicer
#metric1b = LogTgapsMetric(allGaps=True)

# slicer = a grouping or subdivision of visits for the simulated survey
# based on their position on the sky (using a Healpix grid)
#slicer1b = slicers.UserPointsSlicer(ra=[150.11], dec=[2.14])

# constraint = the sql query (or 'select') that selects all visits in r band
constraint1b = 'filter = "r"'
constraint1b += ' and note not like "DD%"'<jupyter_output><empty_output><jupyter_text>The next 3 cells define the metric bundle (and dict), run the metric on just the first opSim and extracts the data for plotting.<jupyter_code># MetricBundle = combination of the metric, slicer, and sqlconstraint
ddf_nvisit1b = metricBundles.MetricBundle(metric1, slicer1, constraint1b)
bundleDict1b = {'logtgaps': ddf_nvisit1b}
group1b = metricBundles.MetricBundleGroup(bundleDict1b, opsdb1, outDir=outDir, resultsDb=resultsDb)
group1b.runAll()
mask1b = ddf_nvisit1b.metricValues.mask
data1b = ddf_nvisit1b.metricValues.data[~mask1b]<jupyter_output><empty_output><jupyter_text>Do the same as above, now modifying the above contraint to keep only the **DDF** observations ("not DDF").<jupyter_code># metric = the "maximum" of the "airmass" for each group of visits in the slicer
#metric1b = LogTgapsMetric(allGaps=True)

# slicer = a grouping or subdivision of visits for the simulated survey
# based on their position on the sky (using a Healpix grid)
#slicer1b = slicers.UserPointsSlicer(ra=[150.11], dec=[2.14])

# constraint = the sql query (or 'select') that selects all visits in r band
constraint1c = 'filter = "r"'
constraint1c += ' and note like "DD%"'
# MetricBundle = combination of the metric, slicer, and sqlconstraint
ddf_nvisit1c = metricBundles.MetricBundle(metric1, slicer1, constraint1c)
bundleDict1c = {'logtgaps': ddf_nvisit1c}
group1c = metricBundles.MetricBundleGroup(bundleDict1c, opsdb1, outDir=outDir, resultsDb=resultsDb)
group1c.runAll()
mask1c = ddf_nvisit1c.metricValues.mask
data1c = ddf_nvisit1c.metricValues.data[~mask1c]<jupyter_output><empty_output><jupyter_text>Now compare the result for all the observations, the DDF observations only, and the WFD observations.<jupyter_code>fig = plt.figure(figsize=(10,8))
plt.hist(np.log10(data1[0]), bins=50, log=True, histtype='step', label='All XMM-LSS')
plt.hist(np.log10(data1b[0]), bins=50, log=True, histtype='step', label='WFD')
plt.hist(np.log10(data1c[0]), bins=50, log=True, histtype='step', label='DDF')

plt.xlim(-3, 3.5)
plt.ylim(bottom=1)
plt.xlabel('log10(days)')
plt.ylabel('Number of pairs')
plt.legend(loc='upper left')
#plt.savefig('XMMLSS_r_WFD_DDF_baseline.png')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10877403846153846, 'fraction_numerical': 0.026642628205128204, 'mean_word_length': 4.269880365939479, 'pattern_counts': {'":': 25, '<': 55, '<?xml version=': 0, '>': 55, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '649948', 'license': 'MIT', 'n_tokens_mistral': 5242, 'n_tokens_neox': 4907, 'n_words': 1513, 'path': 'contrib/04_Log_Delta_T_MetricWFDvsDDF.ipynb', 'repo_name': 'RichardsGroup/LSST_OpSim', 'size': '14995'}