479 origVariance =
False, nptsFull = 1e6, keepPlots =
True, titleFs=14):
480 """Plot diffim residuals for LOCAL and SPATIAL models"""
485 for cell
in kernelCellSet.getCellList():
486 for cand
in cell.begin(
True):
488 if not (cand.getStatus() == afwMath.SpatialCellCandidate.GOOD):
491 cand = diffimLib.cast_KernelCandidateF(cand)
492 diffim = cand.getDifferenceImage(diffimLib.KernelCandidateF.ORIG)
493 orig = cand.getScienceMaskedImage()
495 ski = afwImage.ImageD(kernel.getDimensions())
496 kernel.computeImage(ski,
False, int(cand.getXCenter()), int(cand.getYCenter()))
498 sbg =
background(int(cand.getXCenter()), int(cand.getYCenter()))
499 sdiffim = cand.getDifferenceImage(sk, sbg)
502 bbox = kernel.shrinkBBox(diffim.getBBox())
503 tdiffim = diffim.Factory(diffim, bbox)
504 torig = orig.Factory(orig, bbox)
505 tsdiffim = sdiffim.Factory(sdiffim, bbox)
508 candidateResids.append(np.ravel(tdiffim.getImage().getArray() /
509 np.sqrt(torig.getVariance().getArray())))
510 spatialResids.append(np.ravel(tsdiffim.getImage().getArray() /
511 np.sqrt(torig.getVariance().getArray())))
513 candidateResids.append(np.ravel(tdiffim.getImage().getArray() /
514 np.sqrt(tdiffim.getVariance().getArray())))
515 spatialResids.append(np.ravel(tsdiffim.getImage().getArray() /
516 np.sqrt(tsdiffim.getVariance().getArray())))
518 fullIm = diffExposure.getMaskedImage().getImage().getArray()
519 fullMask = diffExposure.getMaskedImage().getMask().getArray()
521 fullVar = exposure.getMaskedImage().getVariance().getArray()
523 fullVar = diffExposure.getMaskedImage().getVariance().getArray()
526 bitmaskBad |= afwImage.MaskU.getPlaneBitMask(
'NO_DATA')
527 bitmaskBad |= afwImage.MaskU.getPlaneBitMask(
'SAT')
528 idx = np.where((fullMask & bitmaskBad) == 0)
529 stride = int(len(idx[0]) // nptsFull)
530 sidx = idx[0][::stride], idx[1][::stride]
531 allResids = fullIm[sidx] / np.sqrt(fullVar[sidx])
533 testFootprints = diffimTools.sourceToFootprintList(testSources, warpedTemplateExposure,
534 exposure, config, pexLog.getDefaultLog())
535 for fp
in testFootprints:
536 subexp = diffExposure.Factory(diffExposure, fp[
"footprint"].getBBox())
537 subim = subexp.getMaskedImage().getImage()
539 subvar = afwImage.ExposureF(exposure, fp[
"footprint"].getBBox()).getMaskedImage().getVariance()
541 subvar = subexp.getMaskedImage().getVariance()
542 nonfitResids.append(np.ravel(subim.getArray() / np.sqrt(subvar.getArray())))
544 candidateResids = np.ravel(np.array(candidateResids))
545 spatialResids = np.ravel(np.array(spatialResids))
546 nonfitResids = np.ravel(np.array(nonfitResids))
550 from matplotlib.font_manager
import FontProperties
551 except ImportError, e:
552 print "Unable to import pylab: %s" % e
558 fig.canvas._tkcanvas._root().lift()
562 fig.suptitle(
"Diffim residuals: Normalized by sqrt(input variance)", fontsize=titleFs)
564 fig.suptitle(
"Diffim residuals: Normalized by sqrt(diffim variance)", fontsize=titleFs)
566 sp1 = pylab.subplot(221)
567 sp2 = pylab.subplot(222, sharex=sp1, sharey=sp1)
568 sp3 = pylab.subplot(223, sharex=sp1, sharey=sp1)
569 sp4 = pylab.subplot(224, sharex=sp1, sharey=sp1)
570 xs = np.arange(-5, 5.05, 0.1)
571 ys = 1. / np.sqrt(2 * np.pi) * np.exp( -0.5 * xs**2 )
573 sp1.hist(candidateResids, bins=xs, normed=
True, alpha=0.5, label=
"N(%.2f, %.2f)"
574 % (np.mean(candidateResids), np.var(candidateResids)))
575 sp1.plot(xs, ys,
"r-", lw=2, label=
"N(0,1)")
576 sp1.set_title(
"Candidates: basis fit", fontsize=titleFs-2)
577 sp1.legend(loc=1, fancybox=
True, shadow=
True, prop = FontProperties(size=titleFs-6))
579 sp2.hist(spatialResids, bins=xs, normed=
True, alpha=0.5, label=
"N(%.2f, %.2f)"
580 % (np.mean(spatialResids), np.var(spatialResids)))
581 sp2.plot(xs, ys,
"r-", lw=2, label=
"N(0,1)")
582 sp2.set_title(
"Candidates: spatial fit", fontsize=titleFs-2)
583 sp2.legend(loc=1, fancybox=
True, shadow=
True, prop = FontProperties(size=titleFs-6))
585 sp3.hist(nonfitResids, bins=xs, normed=
True, alpha=0.5, label=
"N(%.2f, %.2f)"
586 % (np.mean(nonfitResids), np.var(nonfitResids)))
587 sp3.plot(xs, ys,
"r-", lw=2, label=
"N(0,1)")
588 sp3.set_title(
"Control sample: spatial fit", fontsize=titleFs-2)
589 sp3.legend(loc=1, fancybox=
True, shadow=
True, prop = FontProperties(size=titleFs-6))
591 sp4.hist(allResids, bins=xs, normed=
True, alpha=0.5, label=
"N(%.2f, %.2f)"
592 % (np.mean(allResids), np.var(allResids)))
593 sp4.plot(xs, ys,
"r-", lw=2, label=
"N(0,1)")
594 sp4.set_title(
"Full image (subsampled)", fontsize=titleFs-2)
595 sp4.legend(loc=1, fancybox=
True, shadow=
True, prop = FontProperties(size=titleFs-6))
597 pylab.setp(sp1.get_xticklabels()+sp1.get_yticklabels(), fontsize=titleFs-4)
598 pylab.setp(sp2.get_xticklabels()+sp2.get_yticklabels(), fontsize=titleFs-4)
599 pylab.setp(sp3.get_xticklabels()+sp3.get_yticklabels(), fontsize=titleFs-4)
600 pylab.setp(sp4.get_xticklabels()+sp4.get_yticklabels(), fontsize=titleFs-4)
607 if keepPlots
and not keptPlots:
610 print "%s: Please close plots when done." % __name__
615 print "Plots closed, exiting..."
617 atexit.register(show)