<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Simple Error Model for Twinkles
This notebook will calculate a simple error model for the Twinkles data as a function of the seeing and depth.
Requirements
You will need the DESC Monitor and its dependencies.
You will also need the truth database and OpSim database for the desired Twinkles run. This example uses the run 1.1 settings and thus we use kraken_1042 as our OpSim database, but further investigation will be done when 1.3 results are ready and we will move to the minion_1016 database. LSST OpSim databases can be downloaded here.
Motivation for an Error Model
The Error Model will be a way to estimate the uncertainties in the Twinkles data from the properties of a visit and an object. There are three major motivations behind this
Step1:  Load necessary database connections.
The Monitor has the ability to load the true fluxes and the observed fluxes (stored in MySQL databases on NERSC) for us. If you are not running this on jupyter-dev you will have to open an ssh-tunnel as described in the Monitor package setup.
Let's start by specifying the location of our truth database and setting up a truth database connection.
Since we are using the same field and chip in run 1.1 and 1.3 we can use the run 1.1 stars database without issues when 1.3 is ready.
Step2:  Then we'll establish a database connection to the NERSC MySQL database for the observed data from Twinkles.
Step3:  And finally we'll establish the connection the Opsim database.
Step4:  Load the visit info from DM processed PhoSim images and from Opsim
Our simplest error model wants to model flux errors as a function of visit seeing and 5-sigma depth. Therefore, we will need to pull this info from Opsim and from the DM databases to compare. The Monitor is set up to do this too.
Step5:  Use Monitor to assemble depth and seeing curves
We use the CcdVisit table to populate our 5-sigma depth curves and the seeing curves from the PhoSim images for each visit. We name our instance of the Monitor a worker since it does all the work of gathering data for us.
Step6:  Comparing 5-sigma depth
Step7:  And here we plot by filter.
Step8:  It looks like there are discrepancies between the values measured using the PhoSim images and the Opsim values. We need to look into this. It looks like the DM calculated values in the PhoSim images are consistently deeper than the Opsim depth. Also, the effect looks the worst in the u and y filters. This indicates that there are differences that need to be explained between the PhoSim sky and the OpSim sky.
Seeing Comparison
Step9:  There are also discrepencies here as well. It seems like the PhoSim+DM values are consistently under the Opsim values. 
The large outliers were a result of cosmic ray repair not being turned on in the processing and was discussed in this issue. We therefore do not expect to see these particular outliers in the 1.3 data.
Comparing measured fluxes to the simulation "truth"
In the rest of this notebook we compare the flux values from our PhoSim images to the "true" values that we know from the input catalogs to the simulation. We start with a position matching method that finds matches between the objects of the observed catalogs and the truth catalogs based upon angular distance.
Then the worker calculates the flux residuals (observed - truth) and finds the mean and variance of this value among the objects in each visit.
Note
Step10:  Plotting Bias and Sigma
Now that we have the data we need in the form of mean flux residuals and mean squared flux residuals for each visit we can combine this with our depth and seeing information to construct plots that show the bias and sigma as functions of these values.
<ASSISTANT_TASK:>
Python Code:
import os
import desc.monitor
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from lsst.sims.photUtils import calcNeff
%matplotlib inline
%load_ext autoreload
%autoreload 2
Explanation: Simple Error Model for Twinkles
This notebook will calculate a simple error model for the Twinkles data as a function of the seeing and depth.
Requirements
You will need the DESC Monitor and its dependencies.
You will also need the truth database and OpSim database for the desired Twinkles run. This example uses the run 1.1 settings and thus we use kraken_1042 as our OpSim database, but further investigation will be done when 1.3 results are ready and we will move to the minion_1016 database. LSST OpSim databases can be downloaded here.
Motivation for an Error Model
The Error Model will be a way to estimate the uncertainties in the Twinkles data from the properties of a visit and an object. There are three major motivations behind this:
Develop a method to rapidly simulate large amounts of catalog level datasets for time-variable astronomical objects without the need for image simulation.
Understand the uncertainties that will be present in observed data to improve the likelihoods that we will use for Bayesian estimates of other properties such as the time delays for strongly lensed AGN or the detections of Supernovae.
Look for evidence of systematic errors related to image processing. Also look for major differences between OpSim expected values and those calculated from PhoSim images. Comparisons between OpSim and PhoSim will be useful for other simulation tasks throughout DESC.
So far, we have worked on modeling the relationships between flux uncertainties and the observed seeing and 5 sigma depth of a visit.
NOTE: As mentioned above the current analysis in this notebook is done using Twinkles Run 1.1. This notebook will be updated with our latest findings when Run 1.3 is ready.
Import necessary modules
Don't worry about the warnings below.
End of explanation
star_db_name = '../../twinkles_run1.1.db'
truth_dbConn = desc.monitor.TruthDBInterface(database=star_db_name, driver='sqlite')
Explanation: Load necessary database connections.
The Monitor has the ability to load the true fluxes and the observed fluxes (stored in MySQL databases on NERSC) for us. If you are not running this on jupyter-dev you will have to open an ssh-tunnel as described in the Monitor package setup.
Let's start by specifying the location of our truth database and setting up a truth database connection.
Since we are using the same field and chip in run 1.1 and 1.3 we can use the run 1.1 stars database without issues when 1.3 is ready.
End of explanation
twinkles_dbConn = desc.monitor.DBInterface(database='DESC_Twinkles_Level_2',
                                           #host='127.0.0.1', port='3307',  ##if not running jupyter-dev
                                           host='scidb1.nersc.gov', port=3306,
                                           driver='mysql', project='Twinkles Run1.1')
Explanation: Then we'll establish a database connection to the NERSC MySQL database for the observed data from Twinkles.
End of explanation
opsim_dbConn = desc.monitor.OpsimDBInterface('../../kraken_1042_sqlite.db')  ##Run 1.1 OpSim database
#opsim_dbConn = desc.monitor.OpsimDBInterface('../../minion_1016_sqlite.db')
Explanation: And finally we'll establish the connection the Opsim database.
End of explanation
dm_visit_info = twinkles_dbConn.get_all_visit_info()
opsim_info = opsim_dbConn.get_summary_depth_info_for_visits(1427) #1427 is the Twinkles field ID
#Find entries with the proper obsHistIds in the Opsim data to match the Twinkles run
obs_list = []
for visit_id in dm_visit_info['visit_id']:
    obs_list.append(np.where(opsim_info['obsHistID'] == visit_id)[0][0])
opsim_info = opsim_info[obs_list]
Explanation: Load the visit info from DM processed PhoSim images and from Opsim
Our simplest error model wants to model flux errors as a function of visit seeing and 5-sigma depth. Therefore, we will need to pull this info from Opsim and from the DM databases to compare. The Monitor is set up to do this too.
End of explanation
worker = desc.monitor.Monitor(twinkles_dbConn, truth_dbConn=truth_dbConn)
depth_curve = worker.measure_depth_curve()
seeing_curve = worker.measure_seeing_curve()
Explanation: Use Monitor to assemble depth and seeing curves
We use the CcdVisit table to populate our 5-sigma depth curves and the seeing curves from the PhoSim images for each visit. We name our instance of the Monitor a worker since it does all the work of gathering data for us.
End of explanation
fig = plt.figure()
bins = 15
n,bins,p = plt.hist(depth_curve.lightcurve['mag'], histtype='step', lw=4, bins=15, label='From CcdVisit', range=(21.5,26.5))
plt.hist(opsim_info['fiveSigmaDepth'], histtype='step', bins=bins, lw=4, label='Opsim Values')
plt.legend()
plt.xlabel('5 sigma depth (mags)')
plt.ylabel('Number of Visits')
plt.title('Twinkles 1.1 5-sigma Depth')
#plt.ylim(0, 6500)
Explanation: Comparing 5-sigma depth:
Data Management Calculated Depth from PhoSim images
We use values from CcdVisit table in our Twinkles pserv database along with the following adapted from LSE-40, http://www.astro.washington.edu/users/ivezic/Teaching/Astr511/LSST_SNRdoc.pdf. 
seeing is a value direct from the CcdVisit table.
$Noise$ is equivalent to the sky_noise values in the CcdVisit table.
$zero_point$ below is the zero_point value in the CcdVisit table.
We have that for a single Gaussian profile and in the $\textbf{background dominated limit}$:
$n_{eff} = 2.266*(\frac{seeing}{pixel_scale})^{2}$
$SNR = \frac{Signal}{Noise}$
$SNR = \frac{Source_Counts}{\sigma_{tot}*\sqrt{n_{eff}}}$
So that, SNR = 5 gives
$Source_Counts_{5\sigma} = 5 * \sigma_{tot} * \sqrt{n_{eff}}$
and to convert to flux in maggies we scale by the zeropoint and then multiply by 1e9 to get nanomaggies.
$F_{5} = \frac{10^{9}* Source_Counts_{5\sigma}}{zero_point}$
But, the version we use in the monitor is modified a bit to account for the Poisson noise of the source and thus derive a general form outside of the background dominated limit.
$SNR = \frac{Signal}{\sqrt{Signal + Noise^2}}$
$SNR = 5 = \frac{Source_Counts_{5\sigma}}{\sqrt{Source_Counts_{5\sigma} + (\sigma_{tot}*\sqrt{n_{eff}})^2}}$
We solve this as a quadratic equation and take the positive root.
$Source_Counts_{5\sigma} = \frac{25 + 25\sqrt{1 + \frac{4}{25}(\sigma_{tot}*\sqrt{n_{eff}})^2}}{2}$
and once again convert to flux in nanomaggies using the zeropoint and 1e9 conversion factor.
These results are the $\color{red}{red\ line}$ in the plot below.
Opsim Depth
We compare the results from PhoSim images with the FWHMeff values from the Summary table in the Opsim database.
Now we will plot and compare the two values for 5-sigma depth we get.
End of explanation
fig = plt.figure(figsize=(18, 12))
fig_num = 1
for filter_val in ['u', 'g', 'r', 'i', 'z', 'y']:
    fig.add_subplot(2,3,fig_num)
    n,bins,p = plt.hist(depth_curve.lightcurve['mag'][depth_curve.lightcurve['bandpass'] == str('lsst'+filter_val)], 
                        histtype='step', lw=4, bins=15, label='CcdVisit', range=(21.5,26.5))
    plt.hist(opsim_info['fiveSigmaDepth'][opsim_info['filter'] == filter_val], histtype='step', bins=bins, lw=4, label='Opsim Values')
    if fig_num == 1:
        plt.legend()
    plt.xlabel('5 sigma depth (mags)')
    plt.ylabel('Number of Visits')
    #if fig_num == 1:
    #    plt.ylim(0, 2800)
    plt.title(filter_val)
    fig_num += 1
plt.suptitle('Twinkles 1.1 5-sigma Depth by filter')
plt.subplots_adjust(top=0.93)
Explanation: And here we plot by filter.
End of explanation
plt.hist(opsim_info['FWHMeff'] - seeing_curve.seeing_curve['seeing'], range=(-0.2, 0.4), bins=20)
plt.title(r'Opsim $FWHM_{\bf{eff}}$ - DM seeing')
plt.xlabel(r'Opsim $FWHM_{\bf{eff}}$ - DM seeing (arcsec)')
plt.ylabel('# of visits')
fig = plt.figure(figsize=(8,6))
plt.scatter(opsim_info['FWHMeff'], opsim_info['FWHMeff'] - seeing_curve.seeing_curve['seeing'])
l1, = plt.plot(np.arange(0, 1.8, 0.01), np.zeros(len(np.arange(0, 1.8, 0.01))), c='r', label='DM seeing = Opsim seeing')
plt.xlim(0, 1.8)
#plt.ylim(0, 1.8)
plt.xlabel(r'Opsim $FWHM_{\bf{eff}}$ (arcsec)')
plt.ylabel(r'Opsim $FWHM_{\bf{eff}}$ - DM seeing (arcsec)')
plt.legend([l1], ['PhoSim+DM seeing = Opsim seeing'], loc=2)
plt.title('Twinkles 1.1 Seeing Comparison')
Explanation: It looks like there are discrepancies between the values measured using the PhoSim images and the Opsim values. We need to look into this. It looks like the DM calculated values in the PhoSim images are consistently deeper than the Opsim depth. Also, the effect looks the worst in the u and y filters. This indicates that there are differences that need to be explained between the PhoSim sky and the OpSim sky.
Seeing Comparison
End of explanation
distances = worker.match_catalogs(return_distance=True, within_radius=1./3600.)
worker.calc_flux_residuals(depth_curve, seeing_curve)
Explanation: There are also discrepencies here as well. It seems like the PhoSim+DM values are consistently under the Opsim values. 
The large outliers were a result of cosmic ray repair not being turned on in the processing and was discussed in this issue. We therefore do not expect to see these particular outliers in the 1.3 data.
Comparing measured fluxes to the simulation "truth"
In the rest of this notebook we compare the flux values from our PhoSim images to the "true" values that we know from the input catalogs to the simulation. We start with a position matching method that finds matches between the objects of the observed catalogs and the truth catalogs based upon angular distance.
Then the worker calculates the flux residuals (observed - truth) and finds the mean and variance of this value among the objects in each visit.
Note: The following is currently Run 1.1 data since we have not finished the run 3 processing. As soon as it is available we will update these plots with 1.3 data.
End of explanation
fig = plt.figure(figsize=(18,12))
i=1
for band in ['u', 'g', 'r', 'i', 'z', 'y']:
    fig.add_subplot(2,3,i)
    worker.plot_bias_map(with_bins=12, in_band=band, use_existing_fig=fig)
    i+=1
plt.tight_layout()
fig = plt.figure(figsize=(18,12))
i=1
for band in ['u', 'g', 'r', 'i', 'z', 'y']:
    fig.add_subplot(2,3,i)
    worker.plot_bias_map(with_bins=12, in_band=band, use_existing_fig=fig,
                         normalize=True)
    i+=1
plt.tight_layout()
fig = plt.figure(figsize=(18,12))
i=1
for band in ['u', 'g', 'r', 'i', 'z', 'y']:
    fig.add_subplot(2,3,i)
    worker.plot_bias_scatter(in_band=band, use_existing_fig=fig)
    i+=1
    plt.gca().set_axis_bgcolor('k')
plt.tight_layout()
fig = plt.figure(figsize=(18,12))
i=1
for band in ['u', 'g', 'r', 'i', 'z', 'y']:
    fig.add_subplot(2,3,i)
    worker.plot_bias_scatter(in_band=band, use_existing_fig=fig, normalize=True)
    i+=1
    plt.gca().set_axis_bgcolor('k')
plt.tight_layout()
fig = plt.figure(figsize=(18,12))
i=1
for band in ['u', 'g', 'r', 'i', 'z', 'y']:
    fig.add_subplot(2,3,i)
    worker.plot_variance_map(with_bins=12, in_band=band, use_existing_fig=fig)
    i+=1
plt.tight_layout()
fig = plt.figure(figsize=(18,12))
i=1
for band in ['u', 'g', 'r', 'i', 'z', 'y']:
    fig.add_subplot(2,3,i)
    worker.plot_variance_scatter(in_band=band, use_existing_fig=fig)
    i+=1
plt.tight_layout()
Explanation: Plotting Bias and Sigma
Now that we have the data we need in the form of mean flux residuals and mean squared flux residuals for each visit we can combine this with our depth and seeing information to construct plots that show the bias and sigma as functions of these values.
End of explanation
<END_TASK>