<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 XID+ Example Run Script
(This is based on a Jupyter notebook, available in the XID+ package and can be interactively run and edited)
XID+ is a probababilistic deblender for confusion dominated maps. It is designed to
Step1:  Work out what small tiles are in the test large tile file for PACS
Step2:  You can fit with the numpyro backend.
<ASSISTANT_TASK:>
Python Code:
import numpy as np
from astropy.io import fits
from astropy import wcs
import pickle
import dill
import sys
import os
import xidplus
import copy
from xidplus import moc_routines, catalogue
from xidplus import posterior_maps as postmaps
from builtins import input
Explanation: XID+ Example Run Script
(This is based on a Jupyter notebook, available in the XID+ package and can be interactively run and edited)
XID+ is a probababilistic deblender for confusion dominated maps. It is designed to:
Use a MCMC based approach to get FULL posterior probability distribution on flux
Provide a natural framework to introduce additional prior information
Allows more representative estimation of source flux density uncertainties
Provides a platform for doing science with the maps (e.g XID+ Hierarchical stacking, Luminosity function from the map etc)
Cross-identification tends to be done with catalogues, then science with the matched catalogues.
XID+ takes a different philosophy. Catalogues are a form of data compression. OK in some cases, not so much in others, i.e. confused images: catalogue compression loses correlation information. Ideally, science should be done without compression.
XID+ provides a framework to cross identify galaxies we know about in different maps, with the idea that it can be extended to do science with the maps!!
Philosophy: 
build a probabilistic generative model for the SPIRE maps
Infer model on SPIRE maps
Bayes Theorem
$p(\mathbf{f}|\mathbf{d}) \propto p(\mathbf{d}|\mathbf{f}) \times p(\mathbf{f})$
In order to carry out Bayesian inference, we need a model to carry out inference on.
For the SPIRE maps, our model is quite simple, with likelihood defined as:
    $L = p(\mathbf{d}|\mathbf{f}) \propto |\mathbf{N_d}|^{-1/2} \exp\big{ -\frac{1}{2}(\mathbf{d}-\mathbf{Af})^T\mathbf{N_d}^{-1}(\mathbf{d}-\mathbf{Af})\big}$
where:
    $\mathbf{N_{d,ii}} =\sigma_{inst.,ii}^2+\sigma_{conf.}^2$
Simplest model for XID+ assumes following:
All sources are known and have positive flux (fi)
A global background (B) contributes to all pixels 
PRF is fixed and known
Confusion noise is constant and not correlated across pixels
Because we are getting the joint probability distribution, our model is generative i.e. given parameters, we generate data and vica-versa
Compared to discriminative model (i.e. neural network), which only obtains conditional probability distribution i.e. Neural network, give inputs, get output. Can't go other way'
Generative model is full probabilistic model. Allows more complex relationships between observed and target variables
End of explanation
from healpy import pixelfunc
order_large=6
order_small=10
tile_large=21875
output_folder='../../../test_files/'
outfile=output_folder+'Tile_'+str(tile_large)+'_'+str(order_large)+'.pkl'
with open(outfile, 'rb') as f:
    obj=pickle.load(f)
priors=obj['priors']
theta, phi =pixelfunc.pix2ang(2**order_large, tile_large, nest=True)
tile_small = pixelfunc.ang2pix(2**order_small, theta, phi, nest=True)
priors[0].moc.write(output_folder+'Tile_'+str(tile_large)+'_'+str(order_large)+'_moc_.fits')
from astropy.table import Table, join, vstack,hstack
Table([priors[0].sra,priors[0].sdec]).write(output_folder+'Tile_'+str(tile_large)+'_'+str(order_large)+'_table.fits')
moc=moc_routines.get_fitting_region(order_small,tile_small)
for p in priors:
    p.moc=moc
    p.cut_down_prior()
    p.prior_bkg(0.0,1)
    p.get_pointing_matrix()
moc.write(output_folder+'Tile_'+str(tile_small)+'_'+str(order_small)+'_moc_.fits')
print('fitting '+ str(priors[0].nsrc)+' sources \n')
print('there are '+ str(priors[0].snpix)+' pixels')
%%time
from xidplus.stan_fit import PACS
#priors[0].upper_lim_map()
#priors[0].prior_flux_upper=(priors[0].prior_flux_upper-10.0+0.02)/np.max(priors[0].prf)
fit=PACS.all_bands(priors[0],priors[1],iter=1000)
Took  13205.7 seconds (3.6 hours)
outfile=output_folder+'Tile_'+str(tile_small)+'_'+str(order_small)
posterior=xidplus.posterior_stan(fit,priors)
xidplus.save(priors,posterior,outfile)
post_rep_map = postmaps.replicated_maps(priors, posterior, nrep=2000)
band = ['PACS_100', 'PACS_160']
for i, p in enumerate(priors):
    Bayesian_Pval = postmaps.make_Bayesian_pval_maps(priors[i], post_rep_map[i])
    wcs_temp = wcs.WCS(priors[i].imhdu)
    ra, dec = wcs_temp.wcs_pix2world(priors[i].sx_pix, priors[i].sy_pix, 0)
    kept_pixels = np.array(moc_routines.sources_in_tile([tile_small], order_small, ra, dec))
    Bayesian_Pval[np.invert(kept_pixels)] = np.nan
    Bayes_map = postmaps.make_fits_image(priors[i], Bayesian_Pval)
    Bayes_map.writeto(outfile + '_' + band[i] + '_Bayes_Pval.fits', overwrite=True)
cat = catalogue.create_PACS_cat(posterior, priors[0], priors[1])
kept_sources = moc_routines.sources_in_tile([tile_small], order_small, priors[0].sra, priors[0].sdec)
kept_sources = np.array(kept_sources)
cat[1].data = cat[1].data[kept_sources]
cat.writeto(outfile + '_PACS_cat.fits', overwrite=True)
Explanation: Work out what small tiles are in the test large tile file for PACS
End of explanation
%%time
from xidplus.numpyro_fit import PACS
fit_numpyro=PACS.all_bands(priors)
outfile=output_folder+'Tile_'+str(tile_small)+'_'+str(order_small)+'_numpyro'
posterior_numpyro=xidplus.posterior_numpyro(fit_numpyro,priors)
xidplus.save(priors,posterior_numpyro,outfile)
post_rep_map = postmaps.replicated_maps(priors, posterior_numpyro, nrep=2000)
band = ['PACS_100', 'PACS_160']
for i, p in enumerate(priors):
    Bayesian_Pval = postmaps.make_Bayesian_pval_maps(priors[i], post_rep_map[i])
    wcs_temp = wcs.WCS(priors[i].imhdu)
    ra, dec = wcs_temp.wcs_pix2world(priors[i].sx_pix, priors[i].sy_pix, 0)
    kept_pixels = np.array(moc_routines.sources_in_tile([tile_small], order_small, ra, dec))
    Bayesian_Pval[np.invert(kept_pixels)] = np.nan
    Bayes_map = postmaps.make_fits_image(priors[i], Bayesian_Pval)
    Bayes_map.writeto(outfile + '_' + band[i] + '_Bayes_Pval_numpyro.fits', overwrite=True)
cat = catalogue.create_PACS_cat(posterior_numpyro, priors[0], priors[1])
kept_sources = moc_routines.sources_in_tile([tile_small], order_small, priors[0].sra, priors[0].sdec)
kept_sources = np.array(kept_sources)
cat[1].data = cat[1].data[kept_sources]
cat.writeto(outfile + '_PACS_cat_numpyro.fits', overwrite=True)
moc.area_sq_deg
100.0*(20.0*np.pi*(1.0/3600.0)**2)/moc.area_sq_deg
Explanation: You can fit with the numpyro backend.
End of explanation
<END_TASK>