<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Examples of Stacking BOSS Spectra using Speclite
Examples of using the speclite package to perform basic operations on spectral data accessed with the bossdata package.  To keep the examples small, we use data from a single BOSS plate (6641 observed on MJD 56383) and show how to work with both the individual spec-lite files and the combined spPlate file (see here for details on the different data products).
Package Initialization
Step1:  Stacked Sky
Get a list of sky spectra on plate 6641
Step2:  Plot a stacked spectrum
Step3:  Stack individual Spec-lite files
Loop over all sky spectra on the plate.  The necessary spec-lite files will be automatically downloaded, if necessary, which will take several minutes.
Step4:  Stack Spectra from one Plate file
Accumulate the sky spectra from a Plate file, which will be automatically downloaded if necessary.
Step5:  Stacked Quasars
Get a list of sky spectra on plate 6641, observed on MJD 56383
Step6:  Plot the redshift distribution of the selected quasars
Step7:  Stack spectra from individual Spec-lite files
Loop over all quasar spectra on the plate. The necessary spec-lite files will be automatically downloaded, if necessary, which will take several minutes.
Step8:  Stack spectra from one Plate file
Step9:  Transform each spectrum to its quasar rest frame.  We perform this operation in place (re-using the memory of the input array) and in parallel on all spectra.
Step10:  Resample each spectrum to a uniform rest wavelength grid and stack them together to calculate the mean rest-frame quasar spectrum.  The resample() and accumulate() operations re-use the same memory for each input spectrum, so this loop has fixed (small) memory requirements, independent of the number of spectra being stacked.
<ASSISTANT_TASK:>
Python Code:
%pylab inline
import speclite
print(speclite.version.version)
import bossdata
print(bossdata.__version__)
finder = bossdata.path.Finder()
mirror = bossdata.remote.Manager()
Explanation: Examples of Stacking BOSS Spectra using Speclite
Examples of using the speclite package to perform basic operations on spectral data accessed with the bossdata package.  To keep the examples small, we use data from a single BOSS plate (6641 observed on MJD 56383) and show how to work with both the individual spec-lite files and the combined spPlate file (see here for details on the different data products).
Package Initialization
End of explanation
spAll = bossdata.meta.Database(lite=True)
sky_table = spAll.select_all(where='PLATE=6641 and OBJTYPE="SKY"')
print('Found {0} sky fibers for plate 6641.'.format(len(sky_table)))
Explanation: Stacked Sky
Get a list of sky spectra on plate 6641:
End of explanation
def plot_stack(data, truncate_percentile):
    valid = data['ivar'] > 0
    wlen = data['wavelength'][valid]
    flux = data['flux'][valid]
    dflux = data['ivar'][valid]**(-0.5)
    plt.figure(figsize=(12,5))
    plt.fill_between(wlen, flux, lw=0, color='red')
    plt.errorbar(wlen, flux, dflux, color='black', alpha=0.5, ls='None', capthick=0)
    plt.xlim(np.min(wlen), np.max(wlen))
    plt.ylim(0, np.percentile(flux + dflux, truncate_percentile))
    plt.xlabel('Wavelength ($\AA$)')
    plt.ylabel('Flux $10^{-17}$ erg/(s cm$^2 \AA$)')
    plt.tight_layout();
Explanation: Plot a stacked spectrum:
End of explanation
spec_sky = None
for row in sky_table:
    filename = finder.get_spec_path(plate=row['PLATE'], mjd=row['MJD'], fiber=row['FIBER'], lite=True)
    spectrum = bossdata.spec.SpecFile(mirror.get(filename))
    data = spectrum.get_valid_data(include_sky=True, use_ivar=True, fiducial_grid=True)
    spec_sky = speclite.accumulate(spec_sky, data, data_out=spec_sky, join='wavelength',
                                   add=('flux', 'sky'), weight='ivar')
spec_sky['flux'] += spec_sky['sky']
plot_stack(spec_sky, truncate_percentile=97.5)
Explanation: Stack individual Spec-lite files
Loop over all sky spectra on the plate.  The necessary spec-lite files will be automatically downloaded, if necessary, which will take several minutes.
End of explanation
plate_sky = None
filename = finder.get_plate_spec_path(plate=6641, mjd=56383)
plate = bossdata.plate.PlateFile(mirror.get(filename))
plate_data = plate.get_valid_data(sky_table['FIBER'], include_sky=True, use_ivar=True, fiducial_grid=True)
for data in plate_data:
    plate_sky = speclite.accumulate(plate_sky, data, data_out=plate_sky, join='wavelength',
                                    add=('flux', 'sky'), weight='ivar')
plate_sky['flux'] += plate_sky['sky']
plot_stack(plate_sky, truncate_percentile=97.5)
Explanation: Stack Spectra from one Plate file
Accumulate the sky spectra from a Plate file, which will be automatically downloaded if necessary.
End of explanation
DR12Q = bossdata.meta.Database(finder, mirror, quasar_catalog=True)
qso_table = DR12Q.select_all(where='PLATE=6641 and ZWARNING=0', what='PLATE,MJD,FIBER,Z_VI')
print('Found {0} QSO targets for plate 6641.'.format(len(qso_table)))
Explanation: Stacked Quasars
Get a list of sky spectra on plate 6641, observed on MJD 56383:
End of explanation
plt.hist(qso_table['Z_VI'], bins=25);
plt.xlabel('Redshift z')
plt.ylabel('Quasars')
plt.tight_layout();
Explanation: Plot the redshift distribution of the selected quasars:
End of explanation
fiducial_grid = np.arange(1000.,3000.)
rest_frame, resampled, spec_qso = None, None, None
for row in qso_table:
    filename = finder.get_spec_path(plate=row['PLATE'], mjd=row['MJD'], fiber=row['FIBER'], lite=True)
    spectrum = bossdata.spec.SpecFile(mirror.get(filename))
    data = spectrum.get_valid_data(use_ivar=True, fiducial_grid=True)
    rest_frame = speclite.redshift(z_in=row['Z_VI'], z_out=0, data_in=data, data_out=rest_frame, rules=[
        dict(name='wavelength', exponent=+1),
        dict(name='flux', exponent=-1),
        dict(name='ivar', exponent=+2)])
    resampled = speclite.resample(rest_frame, x_in='wavelength', x_out=fiducial_grid, y=('flux', 'ivar'),
                                  data_out=resampled)
    spec_qso = speclite.accumulate(spec_qso, resampled, data_out=spec_qso, join='wavelength',
                                   add='flux', weight='ivar')
plot_stack(spec_qso, truncate_percentile=99.5)
Explanation: Stack spectra from individual Spec-lite files
Loop over all quasar spectra on the plate. The necessary spec-lite files will be automatically downloaded, if necessary, which will take several minutes.
End of explanation
filename = finder.get_plate_spec_path(plate=6641, mjd=56383)
plate = bossdata.plate.PlateFile(mirror.get(filename))
plate_data = plate.get_valid_data(qso_table['FIBER'], use_ivar=True, fiducial_grid=True)
zorder = np.argsort(qso_table['Z_VI'])
Explanation: Stack spectra from one Plate file
End of explanation
z_in = qso_table['Z_VI'][:,np.newaxis]
plate_data = speclite.redshift(z_in=z_in, z_out=0, data_in=plate_data, data_out=plate_data, rules=[
        dict(name='wavelength', exponent=+1),
        dict(name='flux', exponent=-1),
        dict(name='ivar', exponent=+2)
    ])
Explanation: Transform each spectrum to its quasar rest frame.  We perform this operation in place (re-using the memory of the input array) and in parallel on all spectra.
End of explanation
resampled, plate_qso = None, None
for data in plate_data:
    resampled = speclite.resample(data, x_in='wavelength', x_out=fiducial_grid, y=('flux', 'ivar'), data_out=resampled)
    plate_qso = speclite.accumulate(spec_qso, resampled, data_out=plate_qso, join='wavelength',
                                   add='flux', weight='ivar')
plot_stack(plate_qso, truncate_percentile=99.5)
Explanation: Resample each spectrum to a uniform rest wavelength grid and stack them together to calculate the mean rest-frame quasar spectrum.  The resample() and accumulate() operations re-use the same memory for each input spectrum, so this loop has fixed (small) memory requirements, independent of the number of spectra being stacked.
End of explanation
<END_TASK>