<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Before beginning this exercise you must download some data files, which can be retrieved from here
Step1:  We have provided three images containing stars, taken with 3 different CCDs, in "stars_X.npz" where X = 0, 1, 2.
The true magnitudes of all the stars are also known, taken from a reliable reference catalogue.
OK, I'll be honest; these are simulations.  The stars' centers are all at centres of pixels
make your lives easier; there is no spatial structure in the PSF; and there is also no noise added of any kind.  You may buy me a drink later.
To read the image data and the calibs say something like
Step2:  Just to make sure that we're all on the same page, here's code to display the image and stars (imshow is a simple utility imported from rhlUtils -- feel free to use plt.imshow if you'd rather)
Step3:  Time for you to do some work.  Write some code to estimate a PSF model by simply averaging all the objects, giving each equal weight.  You should add the option to use only some subset of the stars (e.g. the faintest 25%).
Your model should simply be an image of the PSF, normalised to have a maximum value of 1.0
You could use the calibs object to find all the stars, but in the real world there are stars in the data that are
not in the catalogue so you'll have to write a simple object finder (i.e. don't use the calibs!).  It's sufficient to find pixels that are larger than any neighbours to the left, right, top, or bottom; no, this isn't quite how you'd do it in the real world but it's not far off.
I told you that there were three data sets, taken in different places on the sky and with CCDs with different properties.  Choose one set (e.g. "stars_0.npz") to carry out the following activities, then go back and look at the other two.
Step4:  OK, now use your PSF model to create an image of the residuals created by subtracting the scaled PSF from the stars.
How does it look?  Do you see what's going on?  Remember, I said that the PSF isn't a function of position, but CCDs aren't perfect.
<ASSISTANT_TASK:>
Python Code:
import os
import numpy as np
import matplotlib.pyplot as plt
from rhlUtils import BBox, CCD, Image, imshow
%matplotlib notebook 
%config InlineBackend.figure_format = 'retina'
#%matplotlib qt
#%gui qt
dataDir = # complete
mag0 = 33   # Magnitude of an object with 1 detected photon
Explanation: Before beginning this exercise you must download some data files, which can be retrieved from here: 
https://northwestern.box.com/s/rsb5wjb7dztg2128flzgex4fmq19havu
Be sure to move the corresponding files to the directory where you are running this notebook. Furthermore, you will need to provide the path to these data as the variable dataDir, below.
End of explanation
data = np.load(os.path.join(dataDir, "stars_0.npz"))
image, calibs = data["image"], data["calibs"]
image0 = image.copy()    # Keep a copy
Explanation: We have provided three images containing stars, taken with 3 different CCDs, in "stars_X.npz" where X = 0, 1, 2.
The true magnitudes of all the stars are also known, taken from a reliable reference catalogue.
OK, I'll be honest; these are simulations.  The stars' centers are all at centres of pixels
make your lives easier; there is no spatial structure in the PSF; and there is also no noise added of any kind.  You may buy me a drink later.
To read the image data and the calibs say something like:
End of explanation
plt.figure(1)
plt.clf()
imshow(image, vmin=0, vmax=1000)
plt.title("Data")
plt.plot(calibs[:, 0], calibs[:, 1], '+')   # calibs[:, 2] contains the object's magnitude (not flux)
plt.show()
Explanation: Just to make sure that we're all on the same page, here's code to display the image and stars (imshow is a simple utility imported from rhlUtils -- feel free to use plt.imshow if you'd rather)
End of explanation
image = image0.copy()
#...
plt.figure(2)
plt.clf()
imshow(psfIm, vmin=0, vmax=1.1)
plt.title("PSF model")
plt.show();
Explanation: Time for you to do some work.  Write some code to estimate a PSF model by simply averaging all the objects, giving each equal weight.  You should add the option to use only some subset of the stars (e.g. the faintest 25%).
Your model should simply be an image of the PSF, normalised to have a maximum value of 1.0
You could use the calibs object to find all the stars, but in the real world there are stars in the data that are
not in the catalogue so you'll have to write a simple object finder (i.e. don't use the calibs!).  It's sufficient to find pixels that are larger than any neighbours to the left, right, top, or bottom; no, this isn't quite how you'd do it in the real world but it's not far off.
I told you that there were three data sets, taken in different places on the sky and with CCDs with different properties.  Choose one set (e.g. "stars_0.npz") to carry out the following activities, then go back and look at the other two.
End of explanation
image = image0.copy()
# ...
plt.figure(3)
plt.clf()
imshow(image, vmin=image.min(), vmax=image.max()) # , vmin=0, vmax=100)
plt.title("Residuals")
plt.show();
Explanation: OK, now use your PSF model to create an image of the residuals created by subtracting the scaled PSF from the stars.
How does it look?  Do you see what's going on?  Remember, I said that the PSF isn't a function of position, but CCDs aren't perfect.
End of explanation
<END_TASK>