2013-04-24

The following images review the statistics for the SAS g.00000 data in four processings:

• SAS_v8. This is the current public release version of SAS, and was generated on 2012-07-19 22:32:40 with software_ver 34160.
• SAS123.ipp. This is a reprocessing of SAS using the current IPP tag (2013-04-19 23:59:00 with software_ver 35266). This is taken as the default against which the others are measured.
• standard. This is the initial proposed set of stacking changes, including SAFE=F, circularization fix, and bimodal input rejection. (2013-04-24 03:41:11/35307M)
• unrestricted. This is identical to the standard reduction above, but sets the minimum input threshold to 10.0pxl (only inputs above that are rejected with the bimodal code). (2013-04-24 03:41:11/35307M)

First, fraction of input images accepted into the stack (positive good). The SAS_v8 reduction had no input clipping, and so all excluded inputs are rejected at the PSF matching phase. The proposed standard reduction is somewhat more aggressive than the current tag, likely as a result of the "year 3" data having slightly different seeing than the "year 1+2" data. This induces a bimodality that the rejection code picks up on. The unrestricted reduction is more inclusive of things with slightly larger sizes, and due to the fixes in the chi^2 matching, excludes fewer images than the SAS_v8.

Fill factor (positive good) is not very surprising. Turning safety off results in much better fill factors compared with the previous reductions.

The .num files (positive good) contain the number of input images that contributed to a given pixel. This is largely a reflection of the acceptance plot above, showing that as we exclude more input images, the average contributions to a given pixel drops. The unrestricted reduction offers a median improvement of about 0.4 inputs/pixel, with the other reductions having little change.

The above data suggests that switching to the unrestricted option (which will only clip the worst outliers) is probably a good idea. The concern for this is that we will increase the size of the target PSF. Correcting the circularization bug seems to help constrain the targets reasonably well. The SAS_v8 reduction did not have this correction, and did not restrict any inputs, and has the largest PSF size. The standard reduction clips strongly to keep the target PSF size down, and as such, does not increase it relative to the current tag. The unrestricted PSFs are larger than the current tag, but looking at the histogram of FWHM values, it does not seem to push to extreme values, with nearly all targets smaller than 8pxl.

2013-04-23

The circularization bug appears to be a result of an issue in the pmPSF_AxesToModel and pmPSF_ModelToAxes functions that resulted in a factor of ~2 error in the sizes. With this bug resolved, the Target PSF for the example below drops to 4.815245 pxl. Removing any restriction on the input FWHM sizes (-Df PSF.INPUT.THRESH 10.0 so all inputs with FWHM < 10 are accepted) produces a target PSF of 6.571348. As even this is a factor of 2 smaller than the target PSF created before the bug fix, I'm beginning to suspect that we may have been aggressively clipping the inputs as a result of this bug.

The two metrics that are most likely to clarify whether a given stacking recipe is better than another are the average value of the .num image produced (the average number of inputs for a given pixel in the stack) and the S/N distribution of detected objects. For the initial test, I've run the rejection I've been developing (the standard recipe), along with the rejection with PSF.INPUT.THRESH set to 10.0 (the unrestricted recipe). The following image shows the number image on the bottom row, along with a zoom in of a section of the unconvolved image (with finding chart displayed). In the region of this zoom, the additional inputs that are not rejected at input (and have therefore survived the PSF matching phase as well) seem to help fill in gaps and correct some defects (dark spot, bottom middle of image; linear feature, middle right of image). Below this is a plot of the ratio of S/N values for all objects detected by psphot on these images as a function of y-position. As the change in depth is largely a function of y-position, this shows the difference. Note that at y < 1500pxl, the S/N ratio bifurcates, with the samples from areas of added inputs yielding a higher S/N (largely consistent with an improvement of 2 inputs: sqrt(8)/sqrt(10) ~ 0.9). The scatter in the S/N ratio increases at high y, correlating with the other region of increased depth.

This test suggests that restricting inputs based on FWHM may have been necessary due to a bug in the circularization process, and may no longer be needed. The follow-up test will be to measure these two metrics for the sample of:

• Stacks created with current IPP tag.
• Stacks created with proposed "standard" recipe, that does exclude inputs based on FWHM.
• Stacks created with proposed "unrestricted" recipe, that sets the input threshold at 10.0.

2013-04-22

Circularization

When ppStack generates a target PSF, it does so by taking the envelope of a "circularized" realization of the input PSFs. This appears to work fine for PSFs that are already circular, but for highly elliptical PSFs, this breaks down for unknown reasons. As an example, the following figure shows the PSF grid for input 14 for a certain SAS stack. The target grid shown in the middle is dominated by the large size of input 14. However, disabling circularization creates the field shown in the right panel. The circularized form seems much larger than expected from this elongated PSF. Looking at the profiles along the major and minor axis and comparing this to the circularized form shows that there is an error somewhere in this calculation.

The effect of this circularization error is that the input PSF cuts work on a PSF radial profile corresponding to the elliptical form. This produces a FWHM on the order of 4. This is below the threshold, and is therefore accepted. However, upon circularization, the FWHM balloons to ~11, resulting in a target PSF of ~11. This likely explains some of the unexpected jumps in target FWHM from a "well behaved" set of inputs.

Armin's questions/points

• Q1: "One good test to see if everything is propagated correctly is to check the noise"
• A1: I chose six horizontal stripes across the image, and calculated the image noise and variance median for each 100-pixel tall sample (partially because I have that code done and easy to use). To remove the effects of stars (which appear as horizontal excursions in the following plot), I used the median absolute deviation as the image noise statistic, converted in the plot to a Gaussian sigma. The variance image was measured using a simple median. No covariance term was applied, although a quick by-eye estimate suggests a covariance term of ~1.1. Each horizontal stripe sample is offset from the previous by 200 counts on both axes to reduce the confusion of overlaps.

• Q2: "If you look at the data, it is clear that for low flux, the X2 levels out at a constant value, and for large flux value it increases. This means that c1=0.0"
• A2: Agreed. Some consideration was made to switch to fitting a chi^2 floor with an exponential rise with flux, but that should yield effectively the same result as the log-log fit used.
• Q3: "What kind of errors do you assume when you fit stamp_chisqr versus stamp_flux?"
• A3: We use no errors on this fit, so they are implicitly weighted equivalently. This is certainly why the fit is biased when low-flux/high chi^2 points are included (such as in fit 0 here stacking_coverage.20130307). The quick answer would be to use Poissonian error estimates, which would reduce the weight of these outliers. I think the log-log fits reduce the weights of these outliers more (log10(~200)/log10(~30) < sqrt(~200)/sqrt(~30)).

2013-04-19

Comparison of SAS processing with current tag and with proposed updates. The first two columns show the apples-to-apples comparison of the current IPP tag and the proposed changes to the trunk. This shows that the improvement in the input rejection is roughly equal in the two versions of the code, with the SAS123 showing a slight favor to the proposed changes (30% improved input fraction compared with 20% that did not). The increase in the fill factor is largely a result of switching the safety flag to false, and not due to the increase in input number. The final column shows the comparison between the proposed changes and the original SAS_v8 processing from 2012-07-19. This processing did not have the stack input filtering enabled to the degree it currently is, and as a result, only failures in the PSF-matching phase contribute to the input rejection. Without this filtering, the fraction of inputs used is obviously higher. This raises the question of whether input filtering is really helpful for 3pi stacks.

	SAS12 tag vs proposed	SAS123 tag vs proposed	SAS123 SAS_v8 vs proposed
Fill factor
Fraction of inputs accepted
Difference in input fraction

2013-04-16 A

I was able to track down the quadratic fit problem as a result of fitting data with outliers that span many orders of magnitude. I put together a quick hack by taking the logarithm of the chi^2 and flux values, and fitting the quadratic to that. This greatly reduces the scatter found, and provides better fits to the actual data. The following comparison image shows the original stack combination, the reduction with the improved input FWHM cuts (which introduced the chi^2 cut issue), the reduction with this logarithm hack, and finally a result with all of the previous improvements but now with SAFE = F.

2013-04-16

While testing the input FWHM cut, I discovered that an image that appeared valid by eye was being flagged during the PSF matching stage as having a overly deviant chi^2 value. A study of the code revealed that the chi^2 value used is calculated in the following way:

1. For all valid stamps in the PSF match, the flux and squared residual are fit with a quadratic function (x0 + x1 * flux + x2 * flux^2).
2. This quadratic function is subtracted from the squared residual, providing a flux-neutral measurement of the squared residuals.
3. The RMS scatter of these remaining residuals is saved as the chi^2 value for the PSF match (there's a sum over all used kernels, but the example case has only a single kernel).

The following table shows the set of PSF match values, along with the parameters of the quadratic function. This table shows that we seem to be needlessly matching PSFs for images already rejected at the input stage. In addition, there appear to be 2 iterations of this function (pmSubtractionRejectStamps) being called for images 7+.

Nstamp	mean	rms	x0	x1	x2	Image	Rejected?
27	-1e-06	60.7719	19.8792	0.000116288	-9.31452e-13	0
119	1e-06	64.4908	57.5848	-8.45089e-05	3.59486e-11	1
26	-1e-06	43.3212	27.068	5.53666e-05	1.7642e-13	2
107	-2e-05	2756.68	1273.49	-0.00181565	1.34383e-10	3	At Input Stage
106	-9e-06	1105.5	614.941	-0.000711782	1.01758e-10	4	Yes
102	-1e-06	105.827	18.7579	0.000123878	-3.35511e-12	5	At Input Stage
101	-0	51.924	27.2631	7.72215e-05	1.07841e-12	6
100	-6e-06	738.86	222.832	-0.000305871	2.12201e-10
99	8e-06	273.861	95.194	-0.000156691	2.03263e-10	7	Yes
124	4e-06	226.235	18.6949	0.000113144	4.15584e-12
122	2e-06	92.0136	18.5325	0.000110472	3.66559e-13	8
123	1e-06	135.964	-1.48911	0.000144198	-2.92106e-12
121	1e-06	77.5818	7.60706	0.000104179	7.44259e-14	9
103	3e-06	330.224	-70.9388	0.000411208	-5.42577e-12
101	-3e-06	53.8796	-7.55995	0.000204406	5.83814e-12	10
116	0	55.2009	15.2849	4.40914e-05	6.08214e-12
116	-0	23.573	32.1768	1.87319e-06	8.4501e-12	11

We can now plot up the data and quadratic fits to validate their quality. These plots clearly point that the large chi^2 values being used to reject input exposures at the PSF matching stage are caused by poor quadratic fits to relatively well constrained data. This supports the "by-eye" assessment that the images being rejected do not look like they should be rejected.

Comparing polynomial fits from IPP with external code, it appears that the assumption that chisq(flux) is a well-behaved quadratic function is not necessarily true. Fitting a quadratic to log10(chisq(log10(flux))) looks to provide a better fit to the data trend.

Image	Plot	Logfit Plot
0
1
2
3
4
5
6
7a
7b
8a
8b
9a
9b
10a
10b
11a
11b

2013-04-08

The new proposed algorithm for the initial ppStack FWHM cuts uses a Gaussian mixture model analysis to attempt to determine the best separation between a possible two component distribution.

• If Punimodal > 0.5 (this distribution is best described by a single Gaussian mode); limit = MAX(FWHM_i)
• Otherwise, assume two Gaussians best describe it
  • If pi1 > pi2 (the smaller FWHM group has the majority of the samples); limit = m1 + k * sigma1
  • If pi2 > pi1 (the larger FWHM group has the majority of the samples); limit = m2 + k * sigma2
  • If abs(pi1 - pi2) < separation threshold (taken as ~0.2. assume this is nearly unimodal); limit = m2 + k * sigma2
• If N <= 4; limit = MAX(FWHM_i)
• If limit > high_threshold; limit = high_threshold
• If limit < low_threshold; limit = low_threshold

The results of this algorithm change are shows below, displaying the prior results along with this proposed change. The first figure shows the fraction of input images accepted into the stack for each filter (g,r,i,z,y). The second shows the distribution of the FWHM limits.

2013-04-04

I have the first set of ppStack changes in place and merged to the trunk. This includes the fix for the PSF target that was resulting in excess detections. I have run a test on the same SAS stack that Ken and Heather were using for IRAF tests. This test has the SAFE flag set to FALSE, allowing single inputs to contribute.

2013-03-07

For the Arp 220 example, of the six input warps that contain the galaxies, two are excluded due to large FWHM, with another warp having a FWHM close to the hard limit of 7.5. This input was excluded in the original processing, but not in my reprocessing of the warps (FWHM for the reprocessing was 7.48). Because of this, only three (four) warps are available. We currently enforce the "safe combination" rules, which do not allow a pixel to be populated from only one input exposure. The following images show the input warps that contain Arp 220, with the top showing the reprocessed set, and the bottom showing the original warps. The red X show which inputs were excluded due to their FWHM being large. Below that is a comparison of the stacks obtained for these warps with the default combination (left panels) and with the FWHM limit set to 100, which allows all warps to contribute (right panels).