<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Fitting PRFs in K2 TPFs from Campaign 9.1
In this simple tutorial we will show how to perform PRF photometry in a K2 target pixel file using PyKE and oktopus.
This notebook was created with the following versions of PyKE and oktopus
Step1:  1. Importing the necessary packages
As usual, let's start off by importing a few packages from the standard scientific Python stack
Step2:  Since we will perform PRF photometry in a Target Pixel File (TPF), let's import KeplerTargetPixelFile and KeplerQualityFlags classes from PyKE
Step3:  It's always wise to take a prior look at the data, therefore, let's import photutils so that we can perform aperture photometry
Step4:  Additionally, we will need a model for the Pixel Response Function and for the scene. We can import those from PyKE, as well.
Step5:  Finally, we will use oktopus to handle our statistical assumptions
Step6:  2. Actual PRF Photometry
Let's start by instantiating a KeplerTargetPixelFile object (you may either give a url or a local path to the file)
Step7:  Note that we set quality_bitmask=KeplerQualityFlags.CONSERVATIVE_BITMASK, which means that cadences that have specific QUALITY flags such as Attitude tweak, Safe mode, etc, will be ignored.
Let's take a look at the pixel data using the plot method from KeplerTargetPixelFile
Step8:  Now, let's create circular apertures around the targets using photutils
Step9:  We can also use photutils to create aperture photometry light curves from the drawn apertures
Step10:  Let's visualize the light curves
Step11:  Looking at the data before performing PRF photometry is important because it will give us insights on how to construct our priors on the parameters we want to estimate.
Another important part of PRF photometry is the background flux levels. We can either estimate it beforehand and subtract it from the target fluxes or we can let the background be a free parameter during the estimation process. We will choose the latter here, because we will assume that our data comes from Poisson distributions. Therefore, we want the pixel values to be positive everywhere. And, more precisely, remember that the subtraction of two Poisson random variables is not a Poisson rv. Therefore, subtracting any value would break our statistical assumption.
In any case, let's take a look at the background levels using the estimate_background_per_pixel method from KeplerTargetPixelFile class.
Step12:  Ooops! Looks like there is something funny happening on that frame because the background levels are way bellow zero.
Let's plot the pixel data to see what's going on
Step13:  Ok, we see that there is something unusual here. Let's just ignore this cadence for now and move on to our PRF photometry.
Let's redraw the background light curve using more meaningful bounds
Step14:  Now, let's create our PRF model using the SimpleKeplerPRF class
Step15:  This is a simple model based on the PRF calibration data. It depends on the channel and the dimensions of the target pixel file that we want to model. This model is parametrized by stellar positions and flux.
To combine one or more PRF models and a background model, we can use the SceneModel class
Step16:  Note that this class takes a list of PRF objects as required inputs. Additionally, a parameter named bkg_model can be used to model the background variations. The default is a constant value for every frame.
Now that we have taken a look at the data and created our model, let's put our assumptions on the table by defining a prior distribution for the parameters.
Let's choose a uniform prior for the whole parameter space. We can do that using the UniformPrior class
Step17:  This class takes two parameters
Step18:  Now, we can feed both the SceneModel and the UniformPrior objects to the PRFPhotometry class
Step19:  Finally, we use the fit method in which we have to pass the pixel data tpf.flux.
Step20:  Note that our Poisson likelihood assumption is embedded in the PRFPhotometry class. That can be changed while creating PRFPhotometry through the loss_function parameter.
Now, let's retrieve the fitted parameters which are store in the opt_params attribute
Step21:  Let's visualize the parameters as a function of time
Step22:  Oops! That outlier is probably the funny cadence we identified before
Step23:  We can retrieve the residuals using the get_residuals method
Step24:  We can also get the pixels time series for every single model as shown below
Step25:  Let's then visualize the single models
<ASSISTANT_TASK:>
Python Code:
import pyke
pyke.__version__
import oktopus
oktopus.__version__
Explanation: Fitting PRFs in K2 TPFs from Campaign 9.1
In this simple tutorial we will show how to perform PRF photometry in a K2 target pixel file using PyKE and oktopus.
This notebook was created with the following versions of PyKE and oktopus:
End of explanation
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rc
rc('text', usetex=True)
font = {'family' : 'serif',
        'size'   : 22,
        'serif'  : 'New Century Schoolbook'}
rc('font', **font)
Explanation: 1. Importing the necessary packages
As usual, let's start off by importing a few packages from the standard scientific Python stack:
End of explanation
from pyke import KeplerTargetPixelFile, KeplerQualityFlags
Explanation: Since we will perform PRF photometry in a Target Pixel File (TPF), let's import KeplerTargetPixelFile and KeplerQualityFlags classes from PyKE:
End of explanation
import photutils.aperture as apr
Explanation: It's always wise to take a prior look at the data, therefore, let's import photutils so that we can perform aperture photometry:
End of explanation
from pyke.prf import SimpleKeplerPRF, SceneModel, PRFPhotometry
Explanation: Additionally, we will need a model for the Pixel Response Function and for the scene. We can import those from PyKE, as well.
End of explanation
from oktopus import UniformPrior
Explanation: Finally, we will use oktopus to handle our statistical assumptions:
End of explanation
tpf = KeplerTargetPixelFile('https://archive.stsci.edu/missions/k2/target_pixel_files/c91/'
                            '224300000/64000/ktwo224364733-c91_lpd-targ.fits.gz',
                            quality_bitmask=KeplerQualityFlags.CONSERVATIVE_BITMASK)
Explanation: 2. Actual PRF Photometry
Let's start by instantiating a KeplerTargetPixelFile object (you may either give a url or a local path to the file):
End of explanation
tpf.plot(frame=100)
Explanation: Note that we set quality_bitmask=KeplerQualityFlags.CONSERVATIVE_BITMASK, which means that cadences that have specific QUALITY flags such as Attitude tweak, Safe mode, etc, will be ignored.
Let's take a look at the pixel data using the plot method from KeplerTargetPixelFile:
End of explanation
tpf.plot(frame=100)
apr.CircularAperture((941.5, 878.5), r=2).plot(color='cyan')
apr.CircularAperture((944.5, 875.5), r=2).plot(color='cyan')
Explanation: Now, let's create circular apertures around the targets using photutils:
End of explanation
lc1, lc2 = np.zeros(len(tpf.time)), np.zeros(len(tpf.time))
for i in range(len(tpf.time)):
    lc1[i] = apr.CircularAperture((941.5 - tpf.column, 878.5 - tpf.row), r=2).do_photometry(tpf.flux[i])[0]
    lc2[i] = apr.CircularAperture((944.5 - tpf.column, 875.5 - tpf.row), r=2).do_photometry(tpf.flux[i])[0]
Explanation: We can also use photutils to create aperture photometry light curves from the drawn apertures:
End of explanation
plt.figure(figsize=[17, 4])
plt.plot(tpf.time, lc1, 'o', markersize=2)
plt.xlabel("Time")
plt.ylabel("Flux")
plt.figure(figsize=[17, 4])
plt.plot(tpf.time, lc2, 'ro', markersize=2)
plt.xlabel("Time")
plt.ylabel("Flux")
Explanation: Let's visualize the light curves:
End of explanation
bkg = tpf.estimate_bkg_per_pixel(method='mode')
plt.figure(figsize=[17, 4])
plt.plot(tpf.time, bkg, 'ko', markersize=2)
plt.xlabel("Time")
plt.ylabel("Flux")
Explanation: Looking at the data before performing PRF photometry is important because it will give us insights on how to construct our priors on the parameters we want to estimate.
Another important part of PRF photometry is the background flux levels. We can either estimate it beforehand and subtract it from the target fluxes or we can let the background be a free parameter during the estimation process. We will choose the latter here, because we will assume that our data comes from Poisson distributions. Therefore, we want the pixel values to be positive everywhere. And, more precisely, remember that the subtraction of two Poisson random variables is not a Poisson rv. Therefore, subtracting any value would break our statistical assumption.
In any case, let's take a look at the background levels using the estimate_background_per_pixel method from KeplerTargetPixelFile class.
End of explanation
funny_cadence = np.argwhere(bkg < 0)[0][0]
tpf.plot(frame=funny_cadence)
Explanation: Ooops! Looks like there is something funny happening on that frame because the background levels are way bellow zero.
Let's plot the pixel data to see what's going on:
End of explanation
plt.figure(figsize=[17, 4])
plt.plot(tpf.time, bkg, 'ko', markersize=2)
plt.xlabel("Time")
plt.ylabel("Flux")
plt.ylim(2150, 2400)
Explanation: Ok, we see that there is something unusual here. Let's just ignore this cadence for now and move on to our PRF photometry.
Let's redraw the background light curve using more meaningful bounds:
End of explanation
sprf = SimpleKeplerPRF(tpf.channel, tpf.shape[1:], tpf.column, tpf.row)
Explanation: Now, let's create our PRF model using the SimpleKeplerPRF class:
End of explanation
scene = SceneModel(prfs=[sprf] * 2)
Explanation: This is a simple model based on the PRF calibration data. It depends on the channel and the dimensions of the target pixel file that we want to model. This model is parametrized by stellar positions and flux.
To combine one or more PRF models and a background model, we can use the SceneModel class:
End of explanation
prior = UniformPrior(lb=[10e3, 940., 877., 10e3, 943., 874., 1e3],
                     ub=[60e3, 944., 880., 60e3, 947., 877., 4e3])
Explanation: Note that this class takes a list of PRF objects as required inputs. Additionally, a parameter named bkg_model can be used to model the background variations. The default is a constant value for every frame.
Now that we have taken a look at the data and created our model, let's put our assumptions on the table by defining a prior distribution for the parameters.
Let's choose a uniform prior for the whole parameter space. We can do that using the UniformPrior class:
End of explanation
scene.plot(*prior.mean)
Explanation: This class takes two parameters: lb, ub. lb stands for lower bound and ub for upper bound.
The order of those values should correspond to the order of the parameters in our model. For example, an object from
SimpleKeplerPRF takes flux, center_col, and center_row. Therefore, we need to define the prior values on that same order. And since we have two targets, that results in six parameters that have to be defined. The last parameter is the background constant.
Let's visualize our model evaluated at the mean value given by our prior probability:
End of explanation
phot = PRFPhotometry(scene_model=scene, prior=prior)
Explanation: Now, we can feed both the SceneModel and the UniformPrior objects to the PRFPhotometry class:
End of explanation
opt_params = phot.fit(tpf.flux)
Explanation: Finally, we use the fit method in which we have to pass the pixel data tpf.flux.
End of explanation
flux_1 = opt_params[:, 0]
xcenter_1 = opt_params[:, 1]
ycenter_1 = opt_params[:, 2]
flux_2 = opt_params[:, 3]
xcenter_2 = opt_params[:, 4]
ycenter_2 = opt_params[:, 5]
bkg_hat = opt_params[:, 6]
Explanation: Note that our Poisson likelihood assumption is embedded in the PRFPhotometry class. That can be changed while creating PRFPhotometry through the loss_function parameter.
Now, let's retrieve the fitted parameters which are store in the opt_params attribute:
End of explanation
plt.figure(figsize=[18, 4])
plt.plot(tpf.time, flux_1, 'o', markersize=2)
plt.ylabel("Flux")
plt.xlabel("Time")
plt.ylim()
Explanation: Let's visualize the parameters as a function of time:
End of explanation
outlier = np.argwhere(flux_1 > 25000)[0][0]
outlier == funny_cadence
plt.figure(figsize=[18, 4])
plt.plot(tpf.time, flux_1, 'o', markersize=2)
plt.ylabel("Flux")
plt.xlabel("Time")
plt.ylim(20500, 22000)
plt.figure(figsize=[18, 4])
plt.plot(tpf.time, xcenter_1, 'o', markersize=2)
plt.ylabel("Column position")
plt.xlabel("Time")
plt.figure(figsize=[18, 4])
plt.plot(tpf.time, ycenter_1, 'o', markersize=2)
plt.ylabel("Row position")
plt.xlabel("Time")
plt.figure(figsize=[18, 4])
plt.plot(tpf.time, flux_2, 'ro', markersize=2)
plt.ylim(17000, 19000)
plt.ylabel("Flux")
plt.xlabel("Time")
plt.figure(figsize=[18, 4])
plt.plot(tpf.time, xcenter_2, 'ro', markersize=2)
plt.ylabel("Column position")
plt.xlabel("Time")
plt.figure(figsize=[18, 4])
plt.plot(tpf.time, ycenter_2, 'ro', markersize=2)
plt.ylabel("Row position")
plt.xlabel("Time")
plt.figure(figsize=[18, 4])
plt.plot(tpf.time, bkg_hat, 'ko', markersize=2)
plt.ylabel("Background Flux")
plt.xlabel("Time")
plt.ylim(2350, 2650)
Explanation: Oops! That outlier is probably the funny cadence we identified before:
End of explanation
residuals = phot.get_residuals()
plt.imshow(residuals[100], origin='lower')
plt.colorbar()
Explanation: We can retrieve the residuals using the get_residuals method:
End of explanation
prf_1 = np.array([scene.prfs[0](*phot.opt_params[i, 0:3])
                  for i in range(len(tpf.time))])
prf_2 = np.array([scene.prfs[1](*phot.opt_params[i, 3:6])
                  for i in range(len(tpf.time))])
Explanation: We can also get the pixels time series for every single model as shown below:
End of explanation
plt.imshow(prf_1[100], origin='lower', extent=(940, 949, 872, 880))
plt.colorbar()
plt.imshow(prf_2[100], origin='lower', extent=(940, 949, 872, 880))
plt.colorbar()
Explanation: Let's then visualize the single models:
End of explanation
<END_TASK>