Studies to help understand what spectroperfectionism provides as input for coaddition, tests of the coadd algorithm, and plots for the documentation.
In [1]:
%pylab inline
In [2]:
import os
import os.path
import sys
import copy
In [3]:
import scipy.linalg
import scipy.interpolate
In [4]:
import galsim
In [5]:
assert 'DESIMODEL' in os.environ, '$DESIMODEL is not defined.'
sys.path.append(os.path.join(os.environ['DESIMODEL'],'py'))
import desimodel.simulate
In [6]:
assert 'DESI_SPECTRO_DATA' in os.environ, '$DESI_SPECTRO_DATA is not defined.'
assert 'DESI_SPECTRO_REDUX' in os.environ, '$DESI_SPECTRO_REDUX is not defined.'
assert 'PRODNAME' in os.environ, '$PRODNAME is not defined.'
sys.path.append('../py')
import desispec.io
import desispec.coaddition
import desispec.resolution
Use 2 Angstrom binning for the true SED.
In [7]:
wlen = np.linspace(5000.,6000.,501,endpoint=True)
Create a unit-area Gaussian line:
In [8]:
def line(wlen,wlen0,sigma):
# Allow for non-uniform wavelength binning.
dwlen = np.empty_like(wlen)
dwlen[:-1] = np.diff(wlen)
dwlen[-1] = dwlen[-2]
return dwlen*np.exp(-0.5*((wlen-wlen0)/sigma)**2)/(np.sqrt(2*np.pi)*sigma)
In [9]:
assert 1.==np.sum(line(wlen,5500.,25.)),'line() is not correctly normalized.'
Combine a slowly varying continuum with two emission lines. The second line is chosen to be below the detector resolution.
In [10]:
sed = (5.-0.5*((wlen-wlen[0])/(wlen[-1]-wlen[0]))**2 +
300.*line(wlen,5325.,50.) + 30.*line(wlen,5700.,5.))
In [11]:
plt.fill_between(wlen,sed,facecolor=(0,0,1,0.25),color='k');
plt.xlabel('Wavelength (Angstrom)')
plt.ylabel('Flux (arb.units / Angstrom / second)');
Build simple models of fiber optics and a focal plane CCD for toy studies.
In [12]:
class Fiber(object):
"""
A simple model for the optics of a single fiber illuminating the focal plane.
The PSF is modeled as an elliptical Gaussian core and a circular Moffat (beta=3) tail.
Each constructor arg can be provided a either a constant value or else a tuple (begin,end)
of values that are linearly interpolated along the trace in the direction of increasing
wavelength. Linear dimensions for all parameters are in dimensionless units where [0,1]
gives the full width (in the nominal wavelength direction) of the focal plane.
Args:
x: focal-plane position of the trace in the nominal wavelength direction, in
dimensionless units.
y: focal-plane position of the trace in the nominal fiber-number direction, in
dimensionless units.
sigma: Sigma for the Gaussian PSF in dimensionless units.
ellipticity: ratio of the PSF ellipse minor to major axes lengths.
angle: clockwise rotation of the PSF ellipse major axes from the +x direction
in degrees.
tail_fraction: fraction of flux in the PSF tail component.
tail_size: scale parameter for the PSF tail component in dimensionless units.
"""
def __init__(self,x,y,sigma,ellipticity,angle,tail_fraction=0.1,tail_size=0.2):
self.x = self._param(x)
self.y = self._param(y)
self.sigma = self._param(sigma)
self.ellipticity = self._param(ellipticity)
self.angle = self._param(angle)
self.tail_fraction = self._param(tail_fraction)
self.tail_size = self._param(tail_size)
def get_psf(self,t):
"""
Calculates the normalized PSF at one point along the trace. The PSF centroid
and size are in dimensionless units where [0,1] gives the full width of the
focal plane.
"""
frac = self.tail_fraction(t)
return galsim.Add(
galsim.Gaussian(flux = 1-frac,sigma = self.sigma(t))
.shear(q = self.ellipticity(t),beta = self.angle(t)*galsim.degrees),
galsim.Moffat(flux = frac,beta = 3,scale_radius = self.tail_size(t))
).shift(dx = self.x(t),dy = self.y(t))
def offset(self,dy):
fiber = copy.copy(self)
fiber.y = lambda t: self.y(t)+dy
return fiber
@staticmethod
def _param(p):
"""
Returns a function of the dimensionless trace position 0 <= t <= 1 that is
either constant or linearly interpolated.
"""
if isinstance(p,tuple):
p1,p2 = p
dp = p2-p1
return lambda t: p1+dp*t
else:
return lambda t: p
In [13]:
class Detector(object):
"""
A simple model for the 2D response of a multi-fiber spectrograph.
Args:
fibers(iterable): fibers illuminating the focal plane.
min_wlen(float): wavelength corresponding to start of each trace in Angstroms.
max_wlen(float): wavelength corresponding to end of each trace in Angstroms.
num_wlen_bins(int): number of wavelength bins along each trace used to calculate
the detector response matrix.
flux_calib(float): calibration factor that converts from flux units to total
detected electrons per unit time on the CCD.
threshold(float): used to determine the mask of relevant pixels for each fiber.
ccd_width(int): width of simulated CCD in pixels.
ccd_height(int): height of simulated CCD in pixels.
mean_dark_current(float): Mean dark current per pixel in units of electrons
per unit time.
readout_noise_rms(float): RMS Gaussian noise of the readout in electron units.
"""
def __init__(self,fibers,min_wlen=5000.,max_wlen=6000.,num_wlen_bins=81,
flux_calib=200.,threshold=5e-3,ccd_width=100,ccd_height=40,
mean_dark_current=0.05,readout_noise_rms=0.1):
self.fibers = fibers
self.flux_calib = flux_calib
dwlen = (max_wlen - min_wlen)/num_wlen_bins
self.min_wlen = min_wlen
self.max_wlen = max_wlen
self.wlen = 0.5*dwlen + np.linspace(min_wlen,max_wlen,num_wlen_bins,endpoint=False)
self.ccd = galsim.Image(ccd_width,ccd_height,scale = 1./ccd_width)
self.mean_dark_current = mean_dark_current
self.readout_noise_rms = readout_noise_rms
# Calculate the response matrix for each fiber.
responses = np.empty((len(fibers),num_wlen_bins,ccd_height*ccd_width))
dx,dy = -0.5,-0.5*ccd_height/ccd_width
for i,fiber in enumerate(self.fibers):
for j in range(num_wlen_bins):
t = (j+0.5)/num_wlen_bins
# Factor of dwlen/wlen converts from flux/Angstrom/sec to photons/sec.
psf = self.flux_calib*dwlen/self.wlen[j]*fiber.get_psf(t).shift(dx,dy)
psf.drawImage(image = self.ccd)
responses[i,j] = self.ccd.array.flat
# Assign pixels to each fiber.
relevance = np.sum(responses,axis=1)
above_threshold = (relevance > threshold*np.max(relevance,axis=1)[:,np.newaxis])
priority = np.argsort(relevance,axis=0)
highest_priority = (priority[-1] == np.arange(len(fibers))[:,np.newaxis])
self.masks = above_threshold & highest_priority
# Transpose and trim each fiber's response matrix to its pixel mask.
self.responses = [ ]
for i,mask in enumerate(self.masks):
self.responses.append(responses[i].T[mask])
# Reshape masks to the 2D CCD geometry.
self.masks = self.masks.reshape((len(fibers),ccd_height,ccd_width))
def simulate(self,wlen,SEDs,exposure_time,oversampling=5):
"""
Simulate the mean detector response to the SEDs illuminating each fiber.
Args:
wlen(numpy.ndarray): Array of wavelengths in Angstroms where the SEDs are tabulated.
Should cover (at least) the full detector wavelength coverage.
SEDs(iterable): A list of SEDs for each fiber. Each SED is represented by
a numpy array of fluxes per angstrom per second in arbitrary flux units.
exposure_time(float): Simulated exposure time in seconds.
oversampling(int): SED response is simulated as a sum of oversampling*num_wlen_bins
equally spaced PSFs along the fiber trace.
Returns:
numpy.ndarray: Simulated mean CCD pixel values in units of detected electrons,
including the mean dark current.
"""
assert len(SEDs) == len(self.fibers),'Expected %d SEDs' % len(self.fibers)
# Initialize a CCD image with the mean dark current.
image = self.ccd.copy()
image.array[:] = exposure_time*self.mean_dark_current
h,w = image.array.shape
# Calculate oversampled wavelength grid.
num_samples = oversampling*len(self.wlen)
dw = (self.max_wlen - self.min_wlen)/num_samples
wvec = 0.5*dw + np.linspace(self.min_wlen,self.max_wlen,num_samples,endpoint = False)
tvec = (wvec - self.min_wlen)/(self.max_wlen - self.min_wlen)
SEDarray = np.empty((len(self.fibers),len(self.wlen)))
for i,fiber in enumerate(self.fibers):
# Build a linear interpolation of this fiber's SED.
SED_func = scipy.interpolate.interp1d(wlen,SEDs[i],kind='linear')
# Factor of dw/w converts from flux/Angstrom/sec to photons/sec.
flux_norm = exposure_time*self.flux_calib*dw/wvec*SED_func(wvec)
for j,t in enumerate(tvec):
psf = flux_norm[j]*fiber.get_psf(t).shift(-0.5,-0.5*h/w)
psf.drawImage(image = image,add_to_image = True)
# Resample this fiber's SED to our response wavelength grid.
SEDarray[i] = SED_func(self.wlen)
# Return a simulated exposure dictionary.
return { 'exptime': exposure_time, 'image': image.array,
'wlen_in': wlen, 'SEDs_in': SEDs, 'wlen_out': self.wlen, 'SEDs_out': SEDarray }
def add_noise(self,mean_image,num_realizations,seed=123):
"""
Add random noise realizations to a simulated mean detector image.
Includes both shot noise and readout noise.
Args:
mean_image(numpy.ndarray): Mean CCD image in detected electrons, including
dark current, with shape (height,width).
num_realizations(int): Number of random noise realizations to generate.
seed(int): Random seed to use.
"""
assert mean_image.shape == self.ccd.array.shape,'Invalid mean_image shape'
h,w = self.ccd.array.shape
generator = np.random.RandomState(seed)
# Generate Poisson shot noise realizations.
noisy = generator.poisson(lam=mean_image,size=(num_realizations,h,w)).astype(float)
# Add readout noise.
noisy += generator.normal(scale=self.readout_noise_rms,size=(num_realizations,h,w))
return noisy
Implement the reduction algorithm described in Bolton & Schlegel 2009 (BS) http://arxiv.org/abs/0911.2689. The decorrelation algorithm that BS is based on was probably first used in cosmology by Hamilton & Tegmark 2000 http://arxiv.org/abs/astro-ph/9905192.
In [14]:
class Pipeline(object):
def __init__(self,detector):
self.detector = detector
@staticmethod
def decorrelate(Cinv):
"""
Implements the BS / HT decorrelation algorithm, which uses the matrix square root of Cinv
to form a diagonal basis. This is generally a better choice than the eigenvector or
Cholesky bases since it leads to more localized basis vectors.
Args:
Cinv(numpy.ndarray): Square array of inverse covariance matrix elements. The matrix
is assumed to be positive definite but we do not check this.
Returns:
tuple: Tuple ivar,R of uncorrelated flux inverse variances and the corresponding
resolution matrix. These have shapes (nflux,) and (nflux,nflux) respectively.
The rows of R give the resolution-convolved responses to unit flux for each
wavelength bin.
"""
# Calculate the matrix square root of Cinv to diagonalize the flux errors.
L,X = scipy.linalg.eigh((Cinv+Cinv.T)/2.)
# Check that all eigenvalues are positive.
assert np.all(L > 0), 'Found some negative Cinv eigenvalues.'
# Check that the eigenvectors are orthonormal so that Xt.X = 1
assert np.allclose(np.identity(len(L)),X.T.dot(X))
Q = X.dot(np.diag(np.sqrt(L)).dot(X.T))
# Check BS eqn.10
assert np.allclose(Cinv,Q.dot(Q))
# Calculate the corresponding resolution matrix and diagonal flux errors.
s = np.sum(Q,axis=1)
R = Q/s[:,np.newaxis]
ivar = s**2
# Check BS eqn.14
assert np.allclose(Cinv,R.T.dot(np.diag(ivar).dot(R)))
return ivar,R
@staticmethod
def deconvolve(pixel_values,pixel_ivar,response):
"""
Calculate the weighted linear least-squares flux solution for an observed trace.
Args:
pixel_values(numpy.ndarray): 1D array of dark-current subtracted observed image pixel
values in units of detected electrons. Only pixels within the relevance mask for
this fiber should be included here.
pixel_ivar(numpy.ndarray): Uncorrelated pixel inverse-variances to use for weighting.
response(numpy.ndarray): Linear response matrix that transforms a flux vector in
a vector of mean predicted pixel values.
Returns:
tuple: Tuple deconvolved,Cinv,best_fit of the deconvolved best-fit flux (which
will generally include high-frequency oscillations), the correlated
inverse-variance matrix for the decorrelated flux bins, and the best fit pixel
values (dark-current subtracted).
"""
# Calculate the flux inverse covariance for the deconvolved least-squares solution.
ATN = response.T.dot(np.diag(pixel_ivar))
Cinv = ATN.dot(response)
# Solve the linear least-squares problem.
deconvolved,res,rank,sing = scipy.linalg.lstsq(Cinv,ATN.dot(pixel_values))
if rank < len(deconvolved):
print 'WARNING: deconvolved inverse-covariance is not positive definite.'
# Calculate the best-fit pixel values.
best_fit = response.dot(deconvolved)
return deconvolved,Cinv,best_fit
def reduce(self,exposure,refit=True):
"""
Calculate the BS spectro-perfect reduction of the observed image for each fiber.
Use two passes to reduce the bias in the ivar weighting due to statistical fluctuations.
Args:
exposure(dict): Dictionary of simulated exposure info.
refit(bool): Perform a second pass using best fit from first pass to improve ivar.
Returns:
list: A list of reduction objects for each fiber.
"""
image = exposure['image']
dark = exposure['exptime']*self.detector.mean_dark_current
noise_var = self.detector.readout_noise_rms**2
reductions = [ ]
for i,(response,mask) in enumerate(zip(self.detector.responses,self.detector.masks)):
pixel_values = image[mask] - dark
# Use the noisy pixel values for an initial estimate of the pixel inverse-variance.
pixel_ivar = (pixel_values + dark + noise_var)**-1
deconvolved,Cinv,best_fit = self.deconvolve(pixel_values,pixel_ivar,response)
if refit:
# Use the best-fit image for an improved (less biased) estimate of the pixel ivar.
pixel_ivar = (best_fit + dark + noise_var)**-1
deconvolved,Cinv,best_fit = self.deconvolve(pixel_values,pixel_ivar,response)
# Calculate the decorrelated errors and resolution matrix.
ivar,resolution = self.decorrelate(Cinv)
# Normalize to the exposure time to get fluxes.
ivar *= exposure['exptime']**2
deconvolved /= exposure['exptime']
# Convolve the reduced and true fluxes with the resolution.
flux = resolution.dot(deconvolved)
truth = resolution.dot(exposure['SEDs_out'][i])
# Save this reduction.
reductions.append(Reduction(self.detector.wlen,flux,ivar,resolution,truth))
return reductions
The co-addition algorithms to combine exposures from a single camera and to combine camera coadds resampled onto a global wavelength grid are both implemented as methods of the Reduction class:
In [15]:
class Reduction(object):
"""
Represents the results of a pipeline reduction for a single fiber.
"""
def __init__(self,wlen,flux,ivar,resolution,truth):
self.wlen = wlen
self.flux = flux
self.ivar = ivar
self.resolution = resolution
self.truth = truth
self.initialize()
def initialize(self):
# Initialize the quantities we will accumulate during co-addition.
self.Cinv = self.resolution.T.dot(np.diag(self.ivar).dot(self.resolution))
#self.Cinv_f = self.resolution.T.dot(np.diag(self.ivar).dot(self.flux))
#self.Cinv_t = self.resolution.T.dot(np.diag(self.ivar).dot(self.truth))
self.Cinv_f = self.resolution.T.dot(self.ivar*self.flux)
self.Cinv_t = self.resolution.T.dot(self.ivar*self.truth)
def mask_out(self,mask):
self.flux[mask] = 0.
self.ivar[mask] = 0.
self.resolution[mask] = 0.
self.initialize()
def __iadd__(self,other):
"""
Coadd this reduction with another reduction that uses the same wavelength grid.
"""
assert np.array_equal(self.wlen,other.wlen),'Cannot coadd different wavelength grids.'
# Accumulate weighted deconvolved fluxes.
self.Cinv += other.Cinv
self.Cinv_f += other.Cinv_f
self.Cinv_t += other.Cinv_t
# Recalculate the deconvolved solution and resolution.
self.ivar,self.resolution = Pipeline.decorrelate(self.Cinv)
R_it = scipy.linalg.inv(self.resolution.T)
self.flux = R_it.dot(self.Cinv_f)/self.ivar
self.truth = R_it.dot(self.Cinv_t)/self.ivar
return self
@staticmethod
def coadd(global_grid,*reductions):
"""
Combine a list of reductions using a global wavelength grid.
"""
num_global = len(global_grid)
# Step 1: accumulated weighted sums of deconvolved flux estimates.
Cinv = np.zeros((num_global,num_global))
Cinv_f = np.zeros_like(global_grid)
Cinv_t = np.zeros_like(global_grid)
for reduction in reductions:
resampler = Reduction.get_resampling_matrix(global_grid,reduction.wlen)
Cinv += resampler.T.dot(reduction.Cinv.dot(resampler))
Cinv_f += resampler.T.dot(reduction.Cinv_f)
Cinv_t += resampler.T.dot(reduction.Cinv_t)
# Step 2: check for global wavelength bins with no information available.
mask = (np.diag(Cinv) > 0)
keep = np.arange(num_global)[mask]
keep_t = keep[:,np.newaxis]
# Step 3: find decorrelated basis.
ivar = np.zeros_like(global_grid)
resolution = np.zeros_like(Cinv)
ivar[mask],resolution[keep_t,keep] = Pipeline.decorrelate(Cinv[keep_t,keep])
# Step 4: calculate decorrelated flux vectors.
flux = np.zeros_like(global_grid)
truth = np.zeros_like(global_grid)
R_it = scipy.linalg.inv(resolution[keep_t,keep].T)
flux[mask] = R_it.dot(Cinv_f[mask])/ivar[mask]
truth[mask] = R_it.dot(Cinv_t[mask])/ivar[mask]
return Reduction(global_grid,flux,ivar,resolution,truth)
@staticmethod
def get_resampling_matrix(global_grid,local_grid):
"""
Build the rectangular matrix that linearly resamples from the global grid to a local grid.
The local grid range must be contained within the global grid range.
Args:
global_grid(numpy.ndarray): Sorted array of global grid wavelengths.
local_grid(numpy.ndarray): Sorted array of local grid wavelengths.
"""
assert np.all(np.diff(global_grid) > 0),'Global grid is not strictly increasing.'
assert np.all(np.diff(local_grid) > 0),'Local grid is not strictly increasing.'
# Locate each local wavelength in the global grid.
global_index = np.searchsorted(global_grid,local_grid)
assert local_grid[0] >= global_grid[0],'Local grid extends below global grid.'
assert local_grid[-1] <= global_grid[-1],'Local grid extends above global grid.'
# Lookup the global-grid bracketing interval (xlo,xhi) for each local grid point.
# Note that this gives xlo = global_grid[-1] if local_grid[0] == global_grid[0]
# but this is fine since the coefficient of xlo will be zero.
global_xhi = global_grid[global_index]
global_xlo = global_grid[global_index-1]
# Create the rectangular interpolation matrix to return.
alpha = (local_grid - global_xlo)/(global_xhi - global_xlo)
local_index = np.arange(len(local_grid),dtype=int)
matrix = np.zeros((len(local_grid),len(global_grid)))
matrix[local_index,global_index] = alpha
matrix[local_index,global_index-1] = 1 - alpha
return matrix
In [16]:
def plot_response(detector,num_spots,save=None):
# Calculate the spot spacing in wavelength bins.
spot_spacing = (len(detector.wlen)-1)//(num_spots-1)
spot_image = np.zeros_like(detector.ccd.array)
mask_image = np.zeros_like(detector.ccd.array)
for i,(response,mask) in enumerate(zip(detector.responses,detector.masks)):
for spot in response.T[::spot_spacing]:
spot_image[mask] += spot
mask_image += (i+1)*mask
plt.figure(figsize=(7,8))
plt.subplot(2,1,1)
h,w = detector.ccd.array.shape
for fiber in detector.fibers:
plt.plot((w*fiber.x(0),w*fiber.x(1)),(w*fiber.y(0),w*fiber.y(1)),'k-')
assert len(detector.fibers) == 3,'Mask colormap is harded coded for 3 fibers.'
cmap = matplotlib.colors.ListedColormap(['lightgray','r','g','b'])
plt.imshow(mask_image,interpolation='none',origin='lower',cmap=cmap,vmin=-0.5,vmax=3.5)
plt.colorbar(orientation='horizontal',ticks=np.arange(4),
pad=0.1).set_label('Fiber number')
plt.axis('off')
plt.subplot(2,1,2)
plt.imshow(spot_image,interpolation='none',origin='lower')
plt.colorbar(orientation='horizontal',pad=0.1).set_label('Unit flux PSF spots')
plt.axis('off')
plt.tight_layout()
if save: plt.savefig(save)
In [17]:
def pipeline_test(detector,simulated,num_exposures,save=None):
global redux1,redux2
wlen_out = simulated['wlen_out']
pipeline = Pipeline(detector)
noisy = detector.add_noise(simulated['image'],num_exposures);
exposure = copy.copy(simulated)
coadd = [ ]
#
plt.figure(figsize=(7,8))
plt.subplot(2,1,1)
plt.imshow(noisy[0],interpolation='none',origin='lower')
plt.colorbar(orientation='horizontal',pad=0.1).set_label('Detected electrons')
plt.axis('off')
#
plt.subplot(2,1,2)
for i in range(num_exposures):
exposure['image'] = noisy[i]
redux = pipeline.reduce(exposure)
for j in range(len(detector.fibers)):
if i == 0:
# Show errors from first reduction to compare with scatter.
truth = redux[j].truth
dflux = redux[j].ivar**(-0.5)
plt.fill_between(redux[j].wlen,truth-dflux,truth+dflux,alpha=0.2,lw=0,color='r')
# Show the true input SED
plt.plot(exposure['wlen_in'],exposure['SEDs_in'][j],'g-')
# Plot the reduced fluxes for odd-numbered fibers.
if j%2 == 1:
plt.scatter(redux[j].wlen,redux[j].flux,color='b',s=2.,alpha=0.5)
# Coadd all realizations onto the first exposure.
if i == 0:
coadd.append(redux[j])
else:
coadd[j] += redux[j]
# Plot coadd of all exposures for even-numbered fibers.
if j%2 == 0 and i+1 == num_exposures:
plt.errorbar(coadd[j].wlen,coadd[j].flux,coadd[j].ivar**-0.5,ls='None',color='k')
plt.xlim(detector.min_wlen,detector.max_wlen)
plt.ylim(0,1.05*np.max(simulated['SEDs_in']))
plt.xlabel('Wavelength (Angstrom)')
plt.ylabel('Flux (arb.units / Angstrom / second)')
plt.tight_layout()
if save: plt.savefig(save)
In [18]:
def coadd_test(sim1,sim2,detector,save=None):
exp1 = copy.copy(sim1)
exp2 = copy.copy(sim2)
exp1['image'] = detector.add_noise(sim1['image'],1)[0]
exp2['image'] = detector.add_noise(sim2['image'],1)[0]
pipeline = Pipeline(detector)
redux1 = pipeline.reduce(exp1)
redux2 = pipeline.reduce(exp2)
fiber1 = redux2[0]
fiber2 = redux1[2]
fiber12 = copy.copy(fiber1)
fiber12 += fiber2
#
plt.figure(figsize=(7,8))
plt.subplot(2,1,1)
plt.errorbar(fiber12.wlen,fiber12.flux,fiber12.ivar**-0.5,color='k',fmt='.',ls='None')
#plt.plot(fiber1.wlen,fiber1.truth,'r-')
#plt.plot(fiber2.wlen,fiber2.truth,'b-')
plt.plot(fiber12.wlen,fiber12.truth,'k-')
err1 = fiber1.ivar**-0.5
plt.fill_between(fiber1.wlen,fiber1.truth-err1,fiber1.truth+err1,color='r',lw=0,alpha=0.2)
err2 = fiber2.ivar**-0.5
plt.fill_between(fiber2.wlen,fiber2.truth-err2,fiber2.truth+err2,color='b',lw=0,alpha=0.2)
plt.xlabel('Wavelength (Angstrom)')
plt.ylabel('Flux (arb.units / Angstrom / second)')
#
plt.subplot(2,1,2)
plt.errorbar(fiber1.wlen,fiber1.flux,fiber1.ivar**-0.5,color='r',ls='None',alpha=0.5)
plt.errorbar(fiber2.wlen,fiber2.flux,fiber2.ivar**-0.5,color='b',ls='None',alpha=0.5)
plt.plot(fiber12.wlen,fiber12.truth,'k-')
err12 = fiber12.ivar**-0.5
plt.fill_between(fiber12.wlen,fiber12.truth-err12,fiber12.truth+err12,color='k',lw=0,alpha=0.2)
plt.xlabel('Wavelength (Angstrom)')
plt.ylabel('Flux (arb.units / Angstrom / second)')
plt.tight_layout()
if save: plt.savefig(save)
Build a 3-fiber spectrograph with different PSF variations (unrealistically bad) along each trace:
In [19]:
detector = Detector([
Fiber(x=(0.1,0.9),y=(0.15,0.05),sigma=0.020,ellipticity=0.75,angle=(-90.,90.)),
Fiber(x=(0.1,0.9),y=(0.25,0.15),sigma=0.015,ellipticity=0.50,angle=(-45.,135.)),
Fiber(x=(0.1,0.9),y=(0.35,0.25),sigma=0.010,ellipticity=0.25,angle=(0.,180.)),
],ccd_width=100,ccd_height=40)
In [20]:
plot_response(detector,num_spots=9)#,save='../tex/fig/detector.pdf')
Simulate illuminating the fibers by the same SED with normalizations scaled by factors of 2,1,0.5:
In [21]:
exposure = detector.simulate(wlen,[2.0*sed,sed,0.5*sed],exposure_time = 200.);
Compare the results of reducing independent noise realizations of this simulation:
In [22]:
pipeline_test(detector,exposure,4)#,save='../tex/fig/repeats.pdf')
Simulate short and long (4x) exposures with all fibers illuminated by the same SED:
In [23]:
short_exposure = detector.simulate(wlen,[sed,sed,sed],exposure_time = 200.);
In [24]:
long_exposure = detector.simulate(wlen,[sed,sed,sed],exposure_time = 800.);
In [25]:
coadd_test(short_exposure,long_exposure,detector)#,save='../tex/fig/stacked.pdf')
Simulate three cameras with overlapping wavelength coverage:
In [26]:
b_camera = Detector([
Fiber(x=(0.1,0.9),y=(0.15,0.05),sigma=0.015,ellipticity=0.75,angle=(-90.,90.)),
],ccd_width=50,ccd_height=10,min_wlen=5000.,max_wlen=5400.,num_wlen_bins=20)
r_camera = Detector([
Fiber(x=(0.1,0.9),y=(0.15,0.05),sigma=0.010,ellipticity=0.50,angle=(-45.,135.)),
],ccd_width=50,ccd_height=10,min_wlen=5200.,max_wlen=5800.,num_wlen_bins=60)
z_camera = Detector([
Fiber(x=(0.1,0.9),y=(0.15,0.05),sigma=0.015,ellipticity=0.25,angle=(0.,180.)),
],ccd_width=50,ccd_height=10,min_wlen=5600.,max_wlen=6000.,num_wlen_bins=25)
In [27]:
b_exposure = b_camera.simulate(wlen,[sed],exposure_time=200.)
r_exposure = r_camera.simulate(wlen,[sed],exposure_time=200.)
z_exposure = z_camera.simulate(wlen,[sed],exposure_time=200.)
b_exposure['image'] = b_camera.add_noise(b_exposure['image'],1,seed=1)[0]
r_exposure['image'] = r_camera.add_noise(r_exposure['image'],1,seed=2)[0]
z_exposure['image'] = z_camera.add_noise(z_exposure['image'],1,seed=3)[0]
In [28]:
def global_coadd_test(save=None):
b_frame = Pipeline(b_camera).reduce(b_exposure)[0]
r_frame = Pipeline(r_camera).reduce(r_exposure)[0]
z_frame = Pipeline(z_camera).reduce(z_exposure)[0]
# Define our global wavelength grid. We deliberately align with the blue camera
# grid so we can check that coadd = blue where expected.
wlen = np.linspace(4990.,6010.,52)
# Mask out some pixels in the red camera to test how this performs.
r_frame.mask_out(np.arange(6,12))
r_frame.mask_out(np.arange(28,36))
r_cuts = r_frame.wlen[(6,12,28,36),]
# Perform the coaddition across cameras.
coadd = Reduction.coadd(wlen,b_frame,r_frame,z_frame)
#
plt.figure(figsize=(7,8))
plt.subplot(2,1,1)
for frame,color in zip((b_frame,r_frame,z_frame),('b','r','g')):
mask = (frame.ivar > 0)
truth = frame.truth[mask]
err = frame.ivar[mask]**-0.5
plt.fill_between(frame.wlen[mask],truth-err,truth+err,color=color,lw=0,alpha=0.2)
mask = (coadd.ivar > 0)
plt.plot(coadd.wlen[mask],coadd.truth[mask],'k-')
err = coadd.ivar[mask]**-0.5
plt.errorbar(coadd.wlen[mask],coadd.flux[mask],err,color='k',fmt='.',ls='None')
plt.xlim(wlen[0]-10,wlen[-1]+10)
plt.ylim(0.85*np.min(sed),1.1*np.max(sed))
plt.xlabel('Wavelength (Angstrom)')
plt.ylabel('Flux (arb.units / Angstrom / second)')
for r_cut in r_cuts:
plt.axvline(r_cut,color='r',ls=':')
#
plt.subplot(2,1,2)
for frame,color in zip((b_frame,r_frame,z_frame),('b','r','g')):
mask = (frame.ivar > 0)
err = frame.ivar[mask]**-0.5
plt.errorbar(frame.wlen[mask],frame.flux[mask],err,color=color,ls='None')
mask = (coadd.ivar > 0)
truth = coadd.truth[mask]
err = coadd.ivar[mask]**-0.5
plt.fill_between(coadd.wlen[mask],truth-err,truth+err,color='gray',lw=0,alpha=0.2)
plt.plot(coadd.wlen[mask],coadd.truth[mask],'k-')
plt.xlim(wlen[0]-10,wlen[-1]+10)
plt.ylim(0.85*np.min(sed),1.1*np.max(sed))
plt.xlabel('Wavelength (Angstrom)')
plt.ylabel('Flux (arb.units / Angstrom / second)')
for r_cut in r_cuts:
plt.axvline(r_cut,color='r',ls=':')
#
plt.tight_layout()
if save: plt.savefig(save)
In [29]:
global_coadd_test()#save='../tex/fig/global.pdf')
In [30]:
desi_instrument = desimodel.simulate.Instrument(basePath=os.environ['DESIMODEL'])
In [31]:
def plot_grids(s=200,first=0.6,save=None):
# Define the linear wavelength grids used in each camera.
bgrid = np.linspace(3579.0,5938.8,3934)
rgrid = np.linspace(5635.0,7730.8,3494)
zgrid = np.linspace(7445.0,9824.0,3966)
# Build a linear grid with 1A spacing.
linear_grid = np.arange(3579.0,9826.0,1.0)
# Build a log-lambda grid covering the full range with the specified first bin size.
ratio = zgrid[-1]/bgrid[0]
ngrid = 1 + np.floor(np.log(ratio)/np.log((bgrid[0]+first)/bgrid[0])).astype(int)
log_grid = bgrid[0]*np.power(ratio,np.arange(ngrid)/(ngrid-1.))
assert bgrid[0]==log_grid[0] and zgrid[-1] == log_grid[-1],'Bad logspace parameters'
# Calculate bin sizes.
bsize = np.diff(bgrid)
rsize = np.diff(rgrid)
zsize = np.diff(zgrid)
linear_size = np.diff(linear_grid)
log_size = np.diff(log_grid)
#
plt.figure(figsize=(7,4))
plt.plot(bgrid[::s],bsize[::s],'b+')
plt.plot(rgrid[::s],rsize[::s],'rx')
plt.plot(zgrid[::s],zsize[::s],'g+')
plt.plot(linear_grid[::s],linear_size[::s],'k+')
plt.plot(log_grid[::s],log_size[::s],'kx')
# Superimpose the PSF FWHM values from $DESIMODEL.
for band,color in zip('brz',('blue','red','green')):
wlen = desi_instrument.psfFWHMWavelength[band].wavelength
fwhm = desi_instrument.psfFWHMWavelength[band].values
plt.plot(wlen,fwhm,color=color,ls='-')
# Superimpose the BOSS global wavelength grid.
boss_n1 = np.floor(1e4*np.log10(log_grid[0])).astype(int)
boss_n2 = np.ceil(1e4*np.log10(log_grid[-1])).astype(int)
boss_grid = np.power(10.,1e-4*np.arange(boss_n1,boss_n2+1))
boss_size = np.diff(boss_grid)
plt.plot(boss_grid[::s],boss_size[::s],'k--')
#
plt.xlim(log_grid[0],log_grid[-1])
plt.xlabel('Wavelength $\lambda$ (Angstrom)')
plt.ylabel('$\Delta\lambda$ (Angstrom)')
plt.tight_layout()
if save: plt.savefig(save)
In [32]:
plot_grids()#save='../tex/fig/grids.pdf')
In [33]:
os.uname()
Out[33]:
In [34]:
brick_path = desispec.io.meta.findfile('brick',brickid = '3582p000',band = 'r')
In [35]:
brick_file = desispec.io.brick.Brick(brick_path)
In [36]:
def brick_load(brick_file,index=0):
flux = brick_file.hdu_list[0].data[index]
ivar = brick_file.hdu_list[1].data[index]
wlen = brick_file.hdu_list[2].data
R = desispec.resolution.Resolution(brick_file.hdu_list[3].data[index])
return desispec.coaddition.Spectrum(wlen,flux,ivar,R)
In [37]:
spectrum = brick_load(brick_file,0)
In [38]:
def benchmark(spectrum,stage=0):
if stage == 1:
print 'scipy.linalg.eigh'
scipy.linalg.eigh(spectrum.Cinv.todense())
elif stage == 2:
print 'scipy.linalg.inv'
scipy.linalg.inv(spectrum.resolution.todense().T)
elif stage >= 3:
print 'coaddition.decorrelate'
desispec.coaddition.decorrelate(spectrum.Cinv)
elif stage >= 4:
print 'coaddition.Spectrum.finalize'
spectrum.finalize()
In [39]:
%timeit -n10 -r10 benchmark(spectrum,stage=0)
In [40]:
%timeit -n3 -r1 benchmark(spectrum,stage=1)
In [41]:
%timeit -n3 -r1 benchmark(spectrum,stage=2)
In [42]:
%timeit -n3 -r1 benchmark(spectrum,stage=3)
In [43]:
%timeit -n3 -r1 benchmark(spectrum,stage=4)
Look at how much quantities change going through a round trip of decorrelation algegra:
In [44]:
flux1 = np.copy(spectrum.flux)
ivar1 = np.copy(spectrum.ivar)
R1 = spectrum.resolution.todense()
spectrum.finalize()
flux2 = spectrum.flux
ivar2 = spectrum.ivar
R2 = spectrum.resolution.todense()
In [45]:
def round_trip_plot(ndiag=5,save=None):
plt.figure(figsize=(8,8))
#
plt.subplot(3,1,1)
plt.scatter(flux1,flux2-flux1,alpha=0.5,s=1.,color='r')
plt.xlim(-0.5,2.)
plt.ylim(-0.06,0.06)
plt.xlabel('Flux')
plt.ylabel('Change in Flux')
plt.plot((-0.5,2.),(-0.005,+0.02),'b--')
plt.plot((-0.5,2.),(+0.005,-0.02),'b--')
plt.grid()
#
plt.subplot(3,1,2)
plt.scatter(ivar1,ivar2-ivar1,alpha=0.5,s=1.,color='r')
plt.xlim(0.,6.)
plt.ylim(-0.25,0.25)
plt.xlabel('Inverse Variance')
plt.ylabel('Change in Inverse Variance')
plt.plot((0,6),(0,+0.06),'b--')
plt.plot((0,6),(0,-0.06),'b--')
plt.grid()
#
plt.subplot(3,1,3)
colors = iter(matplotlib.cm.rainbow(np.linspace(0, 1, 2*ndiag+1)))
for k in range(-ndiag,ndiag+1):
d1 = np.diag(R1,k=k)
d2 = np.diag(R2,k=k)
plt.scatter(d1,d2-d1,s=2.,alpha=0.5,color=next(colors))
plt.xlim(-0.05,0.5)
plt.ylim(-0.025,0.025)
plt.xlabel('Resolution Value')
plt.ylabel('Change in Resolution Value')
plt.plot((-0.05,0.5),(-0.0005,+0.005),'b--')
plt.plot((-0.05,0.5),(-0.0005,-0.005),'b--')
plt.grid()
#
plt.tight_layout()
if save:
plt.savefig(save)
In [46]:
round_trip_plot()#save='../tex/fig/roundtrip.png')
In [46]: