<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Simulating DESI Spectra
The goal of this notebook is to demonstrate how to generate some simple DESI spectra using the quickgen utility.  For simplicity we will only generate 1D spectra and skip the more computationally intensive (yet still instructive!) step of extracting 1D spectra from simulated 2D spectra (i.e., so-called "pixel-level simulations").
For additional (albeit somewhat outdated) information and documentation about quickgen see
https
Step1:  Check and set your environment
Next, make sure we have all the right environment variables set (assuming the bash shell).  If any of these environment variables are missing please set them in your .bashrc file (and then restart this notebook) or create them for just this notebook session using the %set_env magic command, as we demonstrate below.
Step2:  Let's reassign the $SPECPROD environment to something other than dailytest so that we don't conflict with the outputs of the standard DESI integration test.  In addition, we need to make raw data input $DESI_SPECTO_DATA match $DESI_SPECTRO_SIM/$PIXPROD where the simulated data will be written.
Step3:  Specify the parameters of the simulation.
Next, let's specify the number and spectral type distribution of spectra we want to simulate, and the random seed.  Setting the seed here (which can be any number at all!) ensures that your simulations are reproducible.  Let's also explicitly set the night of the "observations" (the default is to use the current date) and the expid or exposure ID number (which would allow you to simulate more than one DESI exposure).
The flavor option is used to choose the correct sky-brightness model and it also determines the distribution of targets for a given flavor.  For example, flavor='dark' returns the right relative sampling density of ELGs, LRGs, and QSOs.  The other available (science target) options for flavor are 'dark', 'gray', 'grey', 'bright', 'bgs', 'mws', 'lrg', 'elg', 'qso', and 'std'.  (You can also set flavor to either 'arc' or 'flat' but that would be boring!)
Step4:  Generating noiseless spectra.
The first step is to generate the fibermap and simspec files needed by quickgen.  The fibermap table contains (simulated) information about the position of each target in the DESI focal plane, while the simspec table holds the "truth" spectra and the intrinsic properties of each object (redshift, noiseless photometry, [OII] flux, etc.).
In detail, the simspec and fibermap data models are described at
* http
Step5:  Reading the fibermap and spectral metadata
First, let's see what got written to the raw data directory as a result of that last command.
Step6:  Let's go a step further and read the fibermap and simspec files from on-disk.
Note that in general code should not generate filepaths by hand, but rather call desispec.io.findfile to find the file it needs.  If you need to override the standard environment variable locations, you can use the outdir option, while still letting it construct the canonical filename for each type of file.
Step7:  Make a simple plot
Here's a fun simple plot of the redshift histogram distributions.  Now you try!
Step8:  Simulating spectra using quickgen.
We're now ready to simulate DESI spectra using quickgen!  Since we're calling quickgen from within this notebook (rather than from the command line in a terminal) we have to parse the simspec and fibermap filename inputs first.
quickgen generates four types of files and writes them to the \$DESI_SPECTRO_REDUX/\$SPECPROD/exposures directory
Step9:  Inspect the output cframe files
Let's briefly look at one of the cframe files for the blue camera using desispec.io.frame.read_frame, which returns a Frame class with all the attributes you might want.
Step10:  Let's make a quick plot of the zeroth spectrum.
Step11:  Regrouping the spectra
As you can imagine, working with cframe files is pretty tedious, especially across three cameras, 10 spectrographs, and more than 35 million targets!  Therefore, let's combine and reorganize the individual cframe files into spectra files grouped on the sky. Spectra are organized into healpix pixels (here chosen to have nside=64).  If you're interested, you can read more about the healpix directory structure here
Step12:  Inspect the output (combined and regrouped) spectra
So what did we end up with in the redux output directory?
* exposures/{night}/{expid}/
Step13:  As a quick example, let's plot up the zeroth spectrum in healpix pixel 19435.
<ASSISTANT_TASK:>
Python Code:
import os
import numpy as np
import matplotlib.pyplot as plt
from astropy.io import fits
from astropy.table import Table
import desispec.io
import desisim.io
from desisim.obs import new_exposure
from desisim.scripts import quickgen
from desispec.scripts import group_spectra
%pylab inline
Explanation: Simulating DESI Spectra
The goal of this notebook is to demonstrate how to generate some simple DESI spectra using the quickgen utility.  For simplicity we will only generate 1D spectra and skip the more computationally intensive (yet still instructive!) step of extracting 1D spectra from simulated 2D spectra (i.e., so-called "pixel-level simulations").
For additional (albeit somewhat outdated) information and documentation about quickgen see
https://desi.lbl.gov/DocDB/cgi-bin/private/ShowDocument?docid=1429
The heart of quickgen is the SpecSim package, which you can read out here:
http://specsim.readthedocs.io/en/stable
If you identify any errors or have requests for additional functionality please create a new issue on
https://github.com/desihub/desisim/issues
or send a note to &#100;&#101;&#115;&#105;&#45;&#100;&#97;&#116;&#97;&#64;&#100;&#101;&#115;&#105;&#46;&#108;&#98;&#108;&#46;&#103;&#111;&#118;.
Getting started.
First, import all the package dependencies.
End of explanation
def check_env():
    for env in ('DESIMODEL', 'DESI_ROOT', 'DESI_SPECTRO_SIM', 'DESI_SPECTRO_DATA', 
            'DESI_SPECTRO_REDUX', 'SPECPROD', 'PIXPROD'):
        if env in os.environ:
            print('{} environment set to {}'.format(env, os.getenv(env)))
        else:
            print('Required environment variable {} not set!'.format(env))
check_env()
Explanation: Check and set your environment
Next, make sure we have all the right environment variables set (assuming the bash shell).  If any of these environment variables are missing please set them in your .bashrc file (and then restart this notebook) or create them for just this notebook session using the %set_env magic command, as we demonstrate below.
End of explanation
%set_env SPECPROD=example
%set_env PIXPROD=example
rawdata_dir = desisim.io.simdir()
%set_env DESI_SPECTRO_DATA=$rawdata_dir
print('Simulated raw data will be written to {}'.format(desisim.io.simdir()))
print('Pipeline will read raw data from {}'.format(desispec.io.rawdata_root()))
print('    (without knowing that it was simulated)')
print('Pipeline will write processed data to {}'.format(desispec.io.specprod_root()))
Explanation: Let's reassign the $SPECPROD environment to something other than dailytest so that we don't conflict with the outputs of the standard DESI integration test.  In addition, we need to make raw data input $DESI_SPECTO_DATA match $DESI_SPECTRO_SIM/$PIXPROD where the simulated data will be written.
End of explanation
nspec = 100
seed = 555
flavor = 'dark'
night = '20170615'
expid = 0
Explanation: Specify the parameters of the simulation.
Next, let's specify the number and spectral type distribution of spectra we want to simulate, and the random seed.  Setting the seed here (which can be any number at all!) ensures that your simulations are reproducible.  Let's also explicitly set the night of the "observations" (the default is to use the current date) and the expid or exposure ID number (which would allow you to simulate more than one DESI exposure).
The flavor option is used to choose the correct sky-brightness model and it also determines the distribution of targets for a given flavor.  For example, flavor='dark' returns the right relative sampling density of ELGs, LRGs, and QSOs.  The other available (science target) options for flavor are 'dark', 'gray', 'grey', 'bright', 'bgs', 'mws', 'lrg', 'elg', 'qso', and 'std'.  (You can also set flavor to either 'arc' or 'flat' but that would be boring!)
End of explanation
fibermap, truth = new_exposure(flavor=flavor, nspec=nspec, seed=seed, night=night, 
                               expid=expid, tileid=None, exptime=None)
Explanation: Generating noiseless spectra.
The first step is to generate the fibermap and simspec files needed by quickgen.  The fibermap table contains (simulated) information about the position of each target in the DESI focal plane, while the simspec table holds the "truth" spectra and the intrinsic properties of each object (redshift, noiseless photometry, [OII] flux, etc.).
In detail, the simspec and fibermap data models are described at
* http://desidatamodel.readthedocs.io/en/latest/DESI_SPECTRO_SIM/PIXPROD/NIGHT/simspec-EXPID.html
* http://desidatamodel.readthedocs.io/en/latest/DESI_SPECTRO_DATA/NIGHT/fibermap-EXPID.html
To generate these files we'll use new_exposure, a convenience function for generating random typical exposures of various types for testing.  However, note that new_exposure isn't intended for every possible analysis; if you want to use your own mix of objects, you just need to write your own fibermap and simspec files following that format instead of calling new_exposure.
Note that in our call to new_exposure, the tileid and exptime (exposure time) optional inputs are shown for demonstration purposes but do not need to be explicitly set.  In particular, the default exposure time is based on the value specified in the $DESIMODEL/data/desi.yaml parameter file.
Note: The simspec file format may change in the near future, so structure your code to separate generating spectra from the code for writing these particular formats.
End of explanation
rawdata_dir = desispec.io.rawdata_root()
!find $rawdata_dir | sort
Explanation: Reading the fibermap and spectral metadata
First, let's see what got written to the raw data directory as a result of that last command.
End of explanation
fiberfile = desispec.io.findfile('fibermap', night=night, expid=expid)
simspecfile = desisim.io.findfile('simspec', night=night, expid=expid)
print('Reading fibermap file {}'.format(fiberfile))
hdu = fits.open(fiberfile)
hdu.info()
fibermap = Table(hdu['FIBERMAP'].data)
hdu.close()
fibermap[:3]
print('Reading simspec file {}.'.format(simspecfile))
hdu = fits.open(simspecfile)
hdu.info()
meta = Table(hdu['METADATA'].data)
hdu.close()
meta[:3]
Explanation: Let's go a step further and read the fibermap and simspec files from on-disk.
Note that in general code should not generate filepaths by hand, but rather call desispec.io.findfile to find the file it needs.  If you need to override the standard environment variable locations, you can use the outdir option, while still letting it construct the canonical filename for each type of file.
End of explanation
allobjtype = meta['OBJTYPE']
redlim = (-0.2, 1.1*meta['REDSHIFT'].max())
fig, ax = plt.subplots()
for objtype in sorted(set(allobjtype)):
    indx = objtype == allobjtype
    hh = ax.hist(meta['REDSHIFT'][indx], bins=nspec//3, 
                 label=objtype, alpha=0.5, range=redlim)
ax.set_xlabel('Redshift')
ax.set_ylabel('Number of Simulated Spectra')
ax.legend(loc='upper right', ncol=3)
ax.margins(0.2)
ax.set_xlim(redlim)
Explanation: Make a simple plot
Here's a fun simple plot of the redshift histogram distributions.  Now you try!
End of explanation
args = quickgen.parse([
    '--simspec', simspecfile,
    '--fibermap', fiberfile
])
quickgen.main(args)
Explanation: Simulating spectra using quickgen.
We're now ready to simulate DESI spectra using quickgen!  Since we're calling quickgen from within this notebook (rather than from the command line in a terminal) we have to parse the simspec and fibermap filename inputs first.
quickgen generates four types of files and writes them to the \$DESI_SPECTRO_REDUX/\$SPECPROD/exposures directory: calib, sky, cframe, and frame files.  We will use the cframe, or calibrated frame files, which contain the flux-calibrated and sky-subtracted DESI spectra (one file per brz camera and spectrograph).  
The data model and the other files and their contents are documented here:
http://desidatamodel.readthedocs.io/en/latest/DESI_SPECTRO_REDUX/PRODNAME/exposures/NIGHT/EXPID/index.html
The code in the following cell calls the equivalent of the command line:
quickgen --simspec {simspecfile} --fibermap {fiberfile}
End of explanation
cframefile = desispec.io.findfile('cframe', night=night, expid=expid, camera='b0')
print('Reading {}'.format(cframefile))
cframe = desispec.io.frame.read_frame(cframefile)
dir(cframe)
Explanation: Inspect the output cframe files
Let's briefly look at one of the cframe files for the blue camera using desispec.io.frame.read_frame, which returns a Frame class with all the attributes you might want.
End of explanation
print(cframe.wave.shape, cframe.flux.shape)
fig, ax = plt.subplots()
ax.errorbar(cframe.wave, cframe.flux[0, :], 1/np.sqrt(cframe.ivar[0, :]))
ax.set_xlabel('Wavelength (A)')
ax.set_ylabel('Flux ($10^{-17}$ erg/s/cm$^2$)')
Explanation: Let's make a quick plot of the zeroth spectrum.
End of explanation
nside = 64
args = group_spectra.parse(['--hpxnside', '{}'.format(nside)])
group_spectra.main(args)
Explanation: Regrouping the spectra
As you can imagine, working with cframe files is pretty tedious, especially across three cameras, 10 spectrographs, and more than 35 million targets!  Therefore, let's combine and reorganize the individual cframe files into spectra files grouped on the sky. Spectra are organized into healpix pixels (here chosen to have nside=64).  If you're interested, you can read more about the healpix directory structure here:
https://github.com/desihub/desispec/blob/master/doc/nb/Intro_to_DESI_spectra.ipynb
Regrouping is especially important for real observations with overlapping tiles where the same object could be reobserved on different exposures separated by short or large amounts of time.
TODO: Add the spectra files to desidatamodel.
To regroup the spectra, we will run the (notebook) equivalent of the command:
desi_group_spectra --hpxnside 64
End of explanation
reduxdir = desispec.io.specprod_root()
!find $reduxdir | sort
Explanation: Inspect the output (combined and regrouped) spectra
So what did we end up with in the redux output directory?
* exposures/{night}/{expid}/: individual spectrograph camera spectra ("frames")
  grouped by night/expid
* spectra-64/: spectra grouped by healpix location on the sky
End of explanation
specfilename = desispec.io.findfile('spectra', groupname=19435, nside=nside)
print('Reading {}'.format(specfilename))
specobj = desispec.io.read_spectra(specfilename)
dir(specobj)
specobj.wave.keys(), specobj.flux.keys()
thisone = 0
fig, ax = plt.subplots()
for camera, color in zip( ('b', 'r', 'z'), ('blue', 'red', 'magenta') ):
    ax.plot(specobj.wave[camera], specobj.flux[camera][thisone], color=color)
ax.set_xlabel('Wavelength (A)')
ax.set_ylabel('Flux ($10^{-17}$ erg/s/cm$^2$)')
Explanation: As a quick example, let's plot up the zeroth spectrum in healpix pixel 19435.
End of explanation
<END_TASK>