Index: trunk/psphot/test/tap_psphot_varmodel.pro
===================================================================
--- trunk/psphot/test/tap_psphot_varmodel.pro	(revision 33963)
+++ trunk/psphot/test/tap_psphot_varmodel.pro	(revision 34086)
@@ -2,85 +2,111 @@
 # -*-sh-*-
 
-# $KAPA = kapa -noX
-
-# PSF.CONVOLVE : if true, we insert delta functions (and optionally
-#                galaxies) and smooth the image with the psf model
-#                (uses a GAUSS regardless of the model). Note that
-#                PSF.CONVOLVE = T is faster than F, but (a) only
-#                allows Gauss models and (b) only yields quantized
-#                locations
-
-# config for ppImage to generate chip, mask, weight
-$ppImageConfig = -recipe PPIMAGE PPIMAGE_N
-$ppImageConfig = $ppImageConfig -Db BACKGROUND T
-$ppImageConfig = $ppImageConfig -Db CHIP.FITS T
-$ppImageConfig = $ppImageConfig -Db CHIP.MASK.FITS T
-$ppImageConfig = $ppImageConfig -Db CHIP.VARIANCE.FITS T
-$ppImageConfig = $ppImageConfig -Db BASE.FITS F
-$ppImageConfig = $ppImageConfig -Db VARIANCE.BUILD T
-$ppImageConfig = $ppImageConfig -Db PHOTOM F
-
-# basic options for the these images (filter, location, obstype)
-$BaseOptions = -type OBJECT -filter r -skymags 20.86 -ra 270.70 -dec -23.70 -pa 0.0
-$BaseOptions = $BaseOptions -Df PSASTRO:DVO.GETSTAR.MAX.RHO 50000.0
-
-# create an image with fake sources and insert the resulting cmf file into a dvodb
-$RefConfig = -camera SIMTEST 
-$RefConfig = $RefConfig -recipe PPSIM STACKTEST.MAKE 
-$RefConfig = $RefConfig -D PSASTRO:PSASTRO.CATDIR catdir.ref 
-$RefConfig = $RefConfig -Db PSF.CONVOLVE F
-
-# options for the reference image
-$RefOptions = $BaseOptions -exptime 100.0 
-$RefOptions = $RefOptions -seeing 1.0 
-$RefOptions = $RefOptions -D PSF.MODEL PS_MODEL_GAUSS 
-$RefOptions = $RefOptions -Df STARS.DENSITY 10.0 
-$RefOptions = $RefOptions -Df STARS.SIGMA.LIM 0.5
-
-# basic config for ppSim with randomly distributed stars and NO galaxies
-$RealConfig = -camera SIMTEST 
-$RealConfig = $RealConfig -recipe PPSIM STACKTEST.RUN 
-$RealConfig = $RealConfig -D PSASTRO:PSASTRO.CATDIR catdir.ref
-$RealConfig = $RealConfig -Db STARS.FAKE F
-$RealConfig = $RealConfig -Db STARS.REAL T 
-$RealConfig = $RealConfig -Db MATCH.DENSITY F 
-$RealConfig = $RealConfig -Db PSF.CONVOLVE F
-$RealConfig = $RealConfig -Df STARS.DENSITY 10.0
-$RealConfig = $RealConfig -Df STARS.SIGMA.LIM 2.5
-$RealConfig = $RealConfig -Db GALAXY.FAKE F 
-$RealConfig = $RealConfig -Db GALAXY.GRID F 
-
-# options for the repeated images
-$RealOptions = $BaseOptions -exptime 30.0
-  
-$ExtraOptions = -D PSF.MODEL PS_MODEL_GAUSS
-
-# sample alternate options:
-# $ppSimOptions = $FakeOptions -D PSF.MODEL PS_MODEL_PS1_V1
-# $ppSimOptions = $FakeOptions -Df PSF.ARATIO 1.2
-# $ppSimOptions = $FakeOptions -Df PSF.THETA +30.0
-# $ppSimOptions = $FakeOptions -D PSF.MODEL PS_MODEL_GAUSS
-
-list fwhm 
- 1.0 
- 1.1 
- 1.2 
- 1.5
+# This script includes a set of tests to demonstrate the dependence of the faint-end bias on the weighting scheme
+# We have 3 weighting schemes :
+# CONSTANT -- per-pixel weight is fixed (disadvantage: lower S/N, especially for higher sky?)
+# IMAGE_VAR - per-pixel weight is Poisson from image (disadvantage: faint-end bias)
+# MODEL_VAR - per-pixel weight is Poisson from model (disadvantage: 2 linear-fit passes)
+
+# Functions:
+# 
+# init  : initialize variables
+# mkref : generate a fake reference catalog in a DVO database
+# mkexp : generate a fake exposure from fake catalog and detrend it (saves basename.in.cmf & basename.fits)
+# runphot : run psphot on an exposure (saves basename.cmf)
+# ckchip.mags : generate a set of summary plots for the psphot analysis
+
+# I would like to produce a grid of tests:
+# sky (bright, middle, dark)
+# fwhm (1.0, 1.3, 1.6)
+# PSF_MODEL input (GAUSS PS1_V1)
+# PSF_MODEL apply (GAUSS PS1_V1)
+# variance mode: CONSTANT, IMAGE_VAR, MODEL_VAR
+
+macro ppsimtest
+  data $1
+  # sf1 : measured flux inserted in image
+  # sm2 : theoretical magnitude for star
+  read x 1 y 2 sf1 3 sm2 5 Io 12
+  set sm1 = -2.5*log(sf1)
+  set dsm = sm1 - sm2
+  lim -n 2 sm1 -0.2 0.2; clear; box; plot sm1 dsm
+end
+
+macro go.one
+  mkdir test
+
+  local sky fwhm psf_in psf_out
+
+  # mkref creates refimage.* and catdir.ref
+  file catdir.ref found
+  if (not($found))
+    mkref
+  end
+
+  $KAPA = kapa -noX
+  if (not($?RefConfig)) init
+
+  foreach sky 20.0
+    foreach fwhm 1.0
+      foreach psf_in PS1_V1
+
+        sprintf name "test/test.%02d.%02d.%s" $sky {10*$fwhm} $psf_in
+        mkexp $name $sky $fwhm $psf_in
+	data $name.dat
+	read fluxIn 3 MagPredIn 5
+	set MagRealIn = -2.5*log(fluxIn)
+	set dMagIn = MagRealIn - MagPredIn
+        
+        foreach psf_out GAUSS PS1_V1
+          foreach mode CONSTANT IMAGE_VAR MODEL_VAR
+            runphot $name $name.$psf_out.$mode "-D LINEAR_FIT_VARIANCE_MODE $mode -D PSF_MODEL PS_MODEL_$psf_out"
+            ckchip.mags $name.in.cmf $name.$psf_out.$mode.cmf $name.$psf_out.$mode 0.0
+          end
+        end
+      end
+    end
+  end
 end
 
 macro go
-
   mkdir test
 
-  $ExtraOptions = -D PSF.MODEL PS_MODEL_GAUSS
-  mkexp test/image.00 1.0
-  $ExtraOptions = -D PSF.MODEL PS_MODEL_PGAUSS
-  mkexp test/image.01 1.0
-  $ExtraOptions = -D PSF.MODEL PS_MODEL_PS1_V1
-  mkexp test/image.02 1.0
+  local sky fwhm psf_in psf_out
+
+  # mkref creates refimage.* and catdir.ref
+  file catdir.ref found
+  if (not($found))
+    mkref
+  end
+
+  $KAPA = kapa -noX
+  if (not($?RefConfig)) init
+
+  foreach sky 19.0 20.0 21.0
+    foreach fwhm 1.0 1.3 1.6
+      foreach psf_in GAUSS PS1_V1
+
+        sprintf name "test/test.%02d.%02d.%s" $sky {10*$fwhm} $psf_in
+        # mkexp $name $sky $fwhm $psf_in
+	data $name.dat
+	read fluxIn 3 MagPredIn 5
+	set MagRealIn = -2.5*log(fluxIn)
+	set dMagIn = MagRealIn - MagPredIn
+        
+        foreach psf_out GAUSS PS1_V1
+          foreach mode CONSTANT IMAGE_VAR MODEL_VAR
+            # runphot $name $name.$psf_out.$mode "-D LINEAR_FIT_VARIANCE_MODE $mode -D PSF_MODEL PS_MODEL_$psf_out"
+            ckchip.mags $name.in.cmf $name.$psf_out.$mode.cmf $name.$psf_out.$mode 0.0
+          end
+        end
+      end
+    end
+  end
 end
 
 # create a reference database of fake stars to be used by ppSim below
 macro mkref
+  if (not($?RefConfig)) init
+
   exec rm -rf catdir.ref
   exec rm -f refimage.fits
@@ -100,4 +126,29 @@
 # create a realistic distribution of fake stars, GAUSS PSF
 macro mkexp
+  if ($0 != 5)
+    echo "USAGE: mkexp basename sky fwhm psf_model"
+    break
+  end
+
+  local fwhm basename psf_model
+  $basename = $1
+  $sky = $2
+  $fwhm = $3
+  $psf_model = $4
+
+  $ExtraOptions = -D PSF.MODEL PS_MODEL_$psf_model
+
+  # create the raw image
+  echo ppSim -seeing $fwhm -skymags $sky -nx 3000 -ny 3000 $RealConfig $RealOptions $ExtraOptions $basename
+  exec ppSim -seeing $fwhm -skymags $sky -nx 3000 -ny 3000 $RealConfig $RealOptions $ExtraOptions $basename
+  exec /bin/mv -f $basename.cmf $basename.in.cmf
+
+  # create the chip output
+  echo ppImage $ppImageConfig -file $basename.fits $basename
+  exec ppImage $ppImageConfig -file $basename.fits $basename
+end
+
+# create a realistic distribution of fake stars, GAUSS PSF
+macro mkexp.deep
   if ($0 != 3)
     echo "USAGE: mkexp basename fwhm"
@@ -109,7 +160,35 @@
   $fwhm = $2
 
+  $RealOptionsDeep = $BaseOptions -exptime 100
+
   # create the raw image
-  echo ppSim -seeing $fwhm -nx 3000 -ny 3000 $RealConfig $RealOptions $ExtraOptions $basename
-  exec ppSim -seeing $fwhm -nx 3000 -ny 3000 $RealConfig $RealOptions $ExtraOptions $basename
+  echo ppSim -seeing $fwhm -nx 3000 -ny 3000 $RealConfig $RealOptionsDeep $ExtraOptions $basename
+  exec ppSim -seeing $fwhm -nx 3000 -ny 3000 $RealConfig $RealOptionsDeep $ExtraOptions $basename
+  exec /bin/mv -f $basename.cmf $basename.in.cmf
+
+  # create the chip output
+  echo ppImage $ppImageConfig -file $basename.fits $basename
+  exec ppImage $ppImageConfig -file $basename.fits $basename
+end
+
+# create a realistic distribution of fake stars, GAUSS PSF
+macro mkexp.bright
+  if ($0 != 3)
+    echo "USAGE: mkexp basename fwhm"
+    break
+  end
+
+  local fwhm basename
+  $basename = $1
+  $fwhm = $2
+
+  # basic options for the these images (filter, location, obstype)
+  $BaseOptions = -type OBJECT -filter r -skymags 20.0 -ra 270.70 -dec -23.70 -pa 0.0
+  $BaseOptions = $BaseOptions -Df PSASTRO:DVO.GETSTAR.MAX.RHO 50000.0
+  $RealOptionsDeep = $BaseOptions -exptime 100
+
+  # create the raw image
+  echo ppSim -seeing $fwhm -nx 3000 -ny 3000 $RealConfig $RealOptionsDeep $ExtraOptions $basename
+  exec ppSim -seeing $fwhm -nx 3000 -ny 3000 $RealConfig $RealOptionsDeep $ExtraOptions $basename
   exec /bin/mv -f $basename.cmf $basename.in.cmf
 
@@ -133,4 +212,109 @@
   echo psphot -threads 4 -file $basename.ch.fits -mask $basename.ch.mk.fits -variance $basename.ch.wt.fits $outname $options
   exec psphot -threads 4 -file $basename.ch.fits -mask $basename.ch.mk.fits -variance $basename.ch.wt.fits $outname $options
+end
+
+# compare two cmf files with extname Chip.psf 
+# things to compare:
+# * completeness (which sources in (1) are not detected in (2)
+# * positions (X_PSF, Y_PSF) 
+# * instrumental psf mags
+# * position errors (no input errors; use a model?)
+# * measured FWHM?
+# * kron mags (fluxes)
+# * etc, etc
+macro ckchip.mags
+  if ($0 != 5)
+    echo "USAGE: ckchip.mags (raw) (out) (output) (zpt_off)"
+    break
+  end
+
+list pairs
+  PSF_INST_MAG_out      M_raw                 PSF_INST_MAG_raw 2 -0.21 0.21 V
+  AP_MAG_out            M_raw                 PSF_INST_MAG_raw 2 -0.21 0.21 V
+end
+
+  load.cmf $1 Chip.psf raw
+  load.cmf $2 Chip.psf out
+
+  # images generated with convolution will not have the right output positions
+  set X_raw = int(X_PSF_raw) + 0.5
+  set Y_raw = int(Y_PSF_raw) + 0.5
+  set M_raw = PSF_INST_MAG_raw + $4
+  set K_out = -2.5*log(KRON_FLUX_out)
+  match2d X_PSF_out Y_PSF_out X_PSF_raw Y_PSF_raw 1.5 -index1 index1 -index2 index2
+
+  local i NX NY nx ny N
+
+  resize 1000 1000
+
+  clear
+  section a0 0.0 0.0 1.0 0.5
+  label -fn courier 14; 
+  style -pt 0 -sz 0.4
+  show.pair 0
+  reindex ap = AP_MAG_out using index1
+  set dap = ap - v2
+  plot -c red -pt 2 -sz 0.5 rv dap
+
+  delete -q imag_V dmag_V smag_V Dmag_V dmag_Vraw
+  for imag -11 -6 0.2
+    subset dmsub = delta if (rv > $imag) && (rv < $imag + 0.2)
+    vstat -q dmsub
+    concat $imag imag_V
+    concat $MEAN dmag_V
+    concat $MEDIAN Dmag_V
+    concat $SIGMA smag_V
+  end
+
+  vstat -q dmag_V
+  $offset = $MEDIAN
+  set dmag_V = dmag_V - $offset
+  set Dmag_V = Dmag_V - $offset
+
+  section a1 0.0 0.50 1.0 0.25
+  lim rv -0.025 0.075; label -fn courier 14; box; 
+  # plot -c black -pt 0 -sz 0.3 MagPredIn dMagIn
+  plot -c red   -pt 7 -sz 2.0 imag_V dmag_V
+  plot -c darkgreen -pt 3 -sz 2.0 imag_V Dmag_V
+  # plot -c gold -pt 4 -sz 2.0 imag_V dmag_Vraw
+  plot -c blue  -pt 2 -sz 2.0 imag_V smag_V
+  # label -y "mean (red), median (green), unclipped mean (gold), sigma (blue) of delta"
+  label -y "mean (red), median (green), sigma (blue) of delta"
+  sprintf line "OFFSET: %6.3f" $offset
+  textline -frac 0.1 0.8 -fn courier 24 "$line"
+
+  subset dmsub = delta if (rv < -13) && (rv > -16)
+  vstat -q dmsub
+  sprintf line "STDEV: %6.3f" $SIGMA
+  textline -frac 0.1 0.7 -fn courier 24 "$line"
+
+
+  set lChiNorm = log(PSF_CHISQ_out / PSF_NDOF_out)
+  reindex chi = lChiNorm using index1
+  section a2 0.0 0.75 1.0 0.25
+  label -fn courier 14; 
+  lim rv chi; box; 
+  plot -c red -pt 0 -sz 0.5 rv chi
+
+  label -y "chisq"
+  png -name $3.png
+
+  # section a1 0.0 0.5 1.0 0.5
+  # style -pt 0 -sz 0.4
+  # show.pair 1
+end
+
+macro go.phot
+ runphot test.00 test.00.varmode.C "-Db SAVE.RESID T -D LINEAR_FIT_VARIANCE_MODE CONSTANT"
+ runphot test.00 test.00.varmode.I "-Db SAVE.RESID T -D LINEAR_FIT_VARIANCE_MODE IMAGE_VAR"
+ runphot test.00 test.00.varmode.M "-Db SAVE.RESID T -D LINEAR_FIT_VARIANCE_MODE MODEL_VAR"
+ runphot test.00 test.00.varmode.S "-Db SAVE.RESID T -D LINEAR_FIT_VARIANCE_MODE MODEL_SKY"
+end
+
+macro go.vars
+ dev -n varC; ckchip.mags test.00.in.cmf test.00.varmode.C.cmf test.00.varmode.C 0.0
+ dev -n varI; ckchip.mags test.00.in.cmf test.00.varmode.I.cmf test.00.varmode.I 0.0
+ dev -n varM; ckchip.mags test.00.in.cmf test.00.varmode.M.cmf test.00.varmode.M 0.0
+ dev -n varS; ckchip.mags test.00.in.cmf test.00.varmode.S.cmf test.00.varmode.S 0.0
 end
 
@@ -347,5 +531,6 @@
     lim rv $word:4 $word:5; box; plot rv delta
   end
-  label -y '$word:0' -x '$word:2'
+  $line = '$word:0' - '$word:1'
+  label -y "$line" -x '$word:2'
 end
 
@@ -637,4 +822,71 @@
 end
 
+macro init
+  # config for ppImage to generate chip, mask, weight
+  $ppImageConfig = -recipe PPIMAGE PPIMAGE_N
+  $ppImageConfig = $ppImageConfig -Db BACKGROUND T
+  $ppImageConfig = $ppImageConfig -Db CHIP.FITS T
+  $ppImageConfig = $ppImageConfig -Db CHIP.MASK.FITS T
+  $ppImageConfig = $ppImageConfig -Db CHIP.VARIANCE.FITS T
+  $ppImageConfig = $ppImageConfig -Db BASE.FITS F
+  $ppImageConfig = $ppImageConfig -Db VARIANCE.BUILD T
+  $ppImageConfig = $ppImageConfig -Db PHOTOM F
+  
+  # basic options for the these images (filter, location, obstype)
+  $BaseOptions = -type OBJECT -filter r -ra 270.70 -dec -23.70 -pa 0.0
+  $BaseOptions = $BaseOptions -Df PSASTRO:DVO.GETSTAR.MAX.RHO 50000.0
+  
+  # PSF.CONVOLVE : if true, we insert delta functions (and optionally
+  #                galaxies) and smooth the image with the psf model
+  #                (uses a GAUSS regardless of the model). Note that
+  #                PSF.CONVOLVE = T is faster than F, but (a) only
+  #                allows Gauss models and (b) only yields quantized
+  #                locations
+
+  # create an image with fake sources and insert the resulting cmf file into a dvodb
+  $RefConfig = -camera SIMTEST 
+  $RefConfig = $RefConfig -recipe PPSIM STACKTEST.MAKE 
+  $RefConfig = $RefConfig -D PSASTRO:PSASTRO.CATDIR catdir.ref 
+  $RefConfig = $RefConfig -Db PSF.CONVOLVE F
+  
+  # options for the reference image
+  $RefOptions = $BaseOptions
+  $RefOptions = $RefOptions -exptime 100.0 
+  $RefOptions = $RefOptions -seeing 1.0 
+  $RefOptions = $RefOptions -skymags 21.0  
+  $RefOptions = $RefOptions -D PSF.MODEL PS_MODEL_GAUSS 
+  $RefOptions = $RefOptions -Df STARS.DENSITY 10.0 
+  $RefOptions = $RefOptions -Df STARS.SIGMA.LIM 0.5
+
+  # basic config for ppSim with randomly distributed stars and NO galaxies
+  $RealConfig = -camera SIMTEST 
+  $RealConfig = $RealConfig -recipe PPSIM STACKTEST.RUN 
+  $RealConfig = $RealConfig -D PSASTRO:PSASTRO.CATDIR catdir.ref
+  $RealConfig = $RealConfig -Db STARS.FAKE F
+  $RealConfig = $RealConfig -Db STARS.REAL T 
+  $RealConfig = $RealConfig -Db MATCH.DENSITY F 
+  $RealConfig = $RealConfig -Db PSF.CONVOLVE F
+  $RealConfig = $RealConfig -Df STARS.DENSITY 10.0
+  $RealConfig = $RealConfig -Df STARS.SIGMA.LIM 1.0
+  $RealConfig = $RealConfig -Db GALAXY.FAKE F 
+  $RealConfig = $RealConfig -Db GALAXY.GRID F 
+  
+  # options for the repeated images
+  $RealOptions = $BaseOptions -exptime 30.0
+    
+  # sample alternate options:
+  # $ppSimOptions = $FakeOptions -D PSF.MODEL PS_MODEL_PS1_V1
+  # $ppSimOptions = $FakeOptions -Df PSF.ARATIO 1.2
+  # $ppSimOptions = $FakeOptions -Df PSF.THETA +30.0
+  # $ppSimOptions = $FakeOptions -D PSF.MODEL PS_MODEL_GAUSS
+  
+  list fwhm 
+   1.0 
+   1.1 
+   1.2 
+   1.5
+  end
+end
+
 if ($SCRIPT)
   fulltest 4