Index: /branches/eam_branches/ipp-20211108/ippTests/ippCheck/smfcheck.pro
===================================================================
--- /branches/eam_branches/ipp-20211108/ippTests/ippCheck/smfcheck.pro	(revision 41988)
+++ /branches/eam_branches/ipp-20211108/ippTests/ippCheck/smfcheck.pro	(revision 41989)
@@ -1,3 +1,7 @@
 # compare smf-to-smf
+if (not($?SAVEPLOT)) set SAVEPLOT = 0
+if (not($?PLOTROOT))
+  $PLOTROOT = smfcompare
+end  
 
 input tap.sh
@@ -31,109 +35,286 @@
       echo ippCheckKeys $1 $Next $2 $Next
       exec ippCheckKeys $1 $Next $2 $Next
-      break
     end
-  end
-end
-
-macro check.keywords
+
+    # if extname is XYnn.psf, compare the headers AND compare the table values
+    strstr -q $w0:0 .psf -start s1
+    if ($s1 == 4)
+      $Next = $i - 1
+      echo $w0:0 == $Next
+      echo ippCheckKeys $1 $Next $2 $Next
+      exec ippCheckKeys $1 $Next $2 $Next
+    end
+  end
+end
+
+
+# for a single chip, I want outputs with names like:
+# expname.chipID.smfcmp.field.png
+# expname.chipID.smfcmp.txt
+
+### example of plotting all outputs for a single chip:
+macro go.sample
+  $SAVEPLOT = 1
+  $KAPA = kapa -noX
+  smf.compare.onechip o9451g0319o/o9451g0319o.XY01 o9451g0319o.1833521.cm.2467375.smf o9451g0319o.1833521.cm.2467377.smf 2
+end
+
+# generate plot and comparison txt file for one chip (extension)
+macro smf.compare.onechip
   if ($0 != 5)
-    echo "USAGE: check.keywords (smf1) (ext1) (smf2) (ext2)"
-    break
-  end
-
-  exec fhead -x $2 $1 > t1.hdr
-  exec fhead -x $4 $3 > t2.hdr
-  
-  
-  # this comparison compares FITS Table fields
-  exec grep "^TTYPE" t1.hdr | prcol 3 > f1.hdr
-  exec grep "^TTYPE" t2.hdr | prcol 3 > f2.hdr
-
-  break -auto off
-  echo "fields added to $3"
-  exec diff f1.hdr f2.hdr | grep '>'
-
-  echo "fields removed from $1"
-  exec diff f1.hdr f2.hdr | grep '<'
-  break -auto on
-end
-
-macro match.keywords
+    echo "USAGE: smf.compare.onechip (output) (smf1) (smf2) (extension)"
+    break
+  end
+
+  local output smf1 smf2 extnum
+  $output = $1
+  $smf1 = $2
+  $smf2 = $3
+  $extnum = $4
+
+  dirname $output -var mydir
+  mkdir $mydir
+
+  load.smf.fields $smf1 $extnum t0
+  load.smf.fields $smf2 $extnum t1
+
+  smf.pos.compare t0 t1
+  smf.compare.fields $output t0 t1
+  output stdout
+end
+
+# load in the data for a pair of smf fields for a single chip 
+macro load.smf.fields
+  if ($0 != 4)
+    echo "USAGE: load.smf.fields (smffile) (extension) (version)"
+    echo "  loads all fields and saves them with version appended to the names"
+    break
+  end
+
+  local version extnum filename
+  $filename = $1
+  $extnum = $2
+  $version = $3
+
+  data $filename
+  read -fits $extnum -extnum -list-fields -q; list tfields -join allfields
+  read -fits $extnum -extnum $allfields
+  
+  foreach f $allfields
+    set $f\_$version = $f
+    delete $f
+  end
+end
+
+$MATCH_RADIUS = 1; # match radius in pixels
+macro smf.pos.compare
+  if ($0 != 4)
+    echo "USAGE: smf.pos.compare (output) (v1) (v2)"
+    echo "generates matched position lists & stats"
+    break
+  end
+
+  local output v1 v2 i
+  $output = $1
+  $v1     = $2
+  $v2     = $3
+
+  match2d X_PSF_$v1 Y_PSF_$v1 X_PSF_$v1 Y_PSF_$v1 $MATCH_RADIUS -closest
+
+  # initial questions to answer:
+  # * how many sources from #1 are not in #2?
+  # * how many sources from #2 are not in #1?
+  reindex Xm_$v1 = X_PSF_$v1 using index2 -keep-unmatched
+  reindex Ym_$v1 = Y_PSF_$v1 using index2 -keep-unmatched
+
+  reindex Xm_$v1 = X_PSF_$v1 using index1 -keep-unmatched
+  reindex Ym_$v1 = Y_PSF_$v1 using index1 -keep-unmatched
+ 
+  vstat -q Xm_$v1
+  $N1_keep = $NPTS
+  $N1_miss = Xm_$v1[] - $N1_keep
+  
+  vstat -q Xm_$v1
+  $N2_keep = $NPTS
+  $N2_miss = Xm_$v1[] - $N2_keep
+  
+ 
+  echo "$N1_keep & $N1_miss vs $N2_keep & $N2_miss"
+
+  # positional consistency:
+  set dX = X_PSF_$v1 - Xm_$v1
+  set dY = Y_PSF_$v1 - Ym_$v1
+
+  vstat dX
+  vstat dY
+end
+
+macro smf.compare.fields
+  if ($0 != 4)
+    echo "USAGE: smf.pos.compare (output) (v1) (v2)"
+    echo "generates matched position lists & stats"
+    break
+  end
+
+  local output v1 v2 i
+  $output = $1
+  $v1     = $2
+  $v2     = $3
+
+  smf.fields.magrange.init
+  exec rm -f $output.txt
+
+  # XXX this will fail if the two smfs do not have the same fields
+  list fields -split $allfields
+  for i 0 $fields:n
+    $f = $fields:$i
+    smf.fields.magrange.plot $output $f $v1 $v2
+  end  
+end
+
+macro smf.fields.magrange.init
+  clear -s
+  antialias 0.4
+  label -fn helvetica 24
+  resize 2000 800
+  section a 0.0 0.0 0.7 1.0
+  section b 0.7 0.0 0.3 1.0
+  $BSTY = -tickpad 0.1 -labelpadx 3.0 -xpad 3.5 +xpad 1.0
+end
+
+macro smf.fields.magrange.plot
   if ($0 != 5)
-    echo "USAGE test.keywords (smf1) (ext1) (smf2) (ext2)"
-    break
-  end
-
-  # smf1/ext1 holds the referece list
-  list keys -x "smfkeys.sh $1 $2"
-
-  rd h1 $1 -x $2
-  rd h2 $3 -x $4
-  
-  for i 0 $keys:n
-    if ("$keys:$i" == "XTENSION") continue; # opihi transforms XTENSION to SIMPLE
-
-    $mykey = $keys:$i
-    keyword h1 $mykey val1
-    keyword h2 $mykey val2
-
-    # XXX I do not have a good way to test the
-    # type of the header value. I can just check
-
-    # empty string needs to be specially handled
-    if (("$val1" == "") || ("$val2" == ""))
-      if ("$val1" == "$val2") 
-        # echo "matched values $keys:$i : $val1 vs $val2"
-      else
-        echo "mismatched values $keys:$i : $val1 vs $val2"
-      end
-      continue
-    end
-
-    # multi-word strings need to be specially handled
-    getchr "$val1" " " s1
-    getchr "$val2" " " s2
-    if (($s1 > -1) || ($s2 > -1))
-      # echo "multiword: $val1 : $s1"
-      # echo "multiword: $val2 : $s2"
-      if ("$val1" == "$val2") 
-        # echo "matched multivalues $keys:$i : $val1 vs $val2"
-      else
-        echo "mismatched multivalues $keys:$i : $val1 vs $val2"
-      end
-      continue
-    end
-
-    # strings with +,-,*,/,^ need to be specially handled
-    $hasmath1 = 0
-    $hasmath2 = 0
-    foreach mychar + - * / ^ : ,
-      getchr "$val1" "$mychar" s1
-      if ($s1 > -1)
-        $hasmath1 = 1
-      end
-      getchr "$val2" "$mychar" s2
-      if ($s2 > -1)
-        $hasmath2 = 1
-      end
-    end
-
-    if ($hasmath1 || $hasmath2)
-      #echo "math: $val1 : $s1"
-      #echo "math: $val2 : $s2"
-      if ("$val1" == "$val2") 
-        # echo "matched multivalues $keys:$i : $val1 vs $val2"
-      else
-        echo "mismatched multivalues $keys:$i : $val1 vs $val2"
-      end
-      continue
-    end
-
-    if (isword($val1)) 
-      if ("$val1" == "$val2") continue
-    else
-      if ($val1 == $val2) continue
-    end
-    echo "mismatched values $keys:$i : $val1 vs $val2"
-  end
-end
+    echo "USAGE: smf.fields.magrange.plot (output) {field) (v1) (v2)"
+    echo " NOTE: run smf.pos.compare (or equiv) first to generate index1, index2"
+    break
+  end
+  
+  local output field v1 v2
+  $output = $1
+  $field  = $2
+  $v1     = $3
+  $v2     = $4
+
+  # plot the new values vs the old values (all old values)
+  reindex tmp_$v2 = $field\_$v2 using index1 -keep-unmatched
+  set dtmp = tmp_$v2 - $field\_$v1
+  hist.range dtmp Ntmp dx
+  vstat -q dtmp
+
+  # save the summary line to be written to a file by the calling function
+  output $output.txt
+  fprintf "%24s : %12.6g %12.6g : %12.6g : %12.6g - %12.6g" $field $MEAN $SIGMA $hist:median $hist:minVal $hist:maxVal
+  output stdout
+
+  if (isnan($hist:median))
+    continue
+  end
+
+  clear
+  section a
+  lim PSF_INST_MAG_$v1 $hist:minVal $hist:maxVal; box $BSTY; plot PSF_INST_MAG_$v1 dtmp -pt cir -op 0.3
+  section b
+  if (dx[] == 1)
+    lim dx Ntmp; box $BSTY; plot -pt cir dx Ntmp; label -x $field
+  else
+    lim dx Ntmp; box $BSTY; plot -x hist dx Ntmp; label -x $field
+  end
+ 
+  if ($SAVEPLOT)
+    png -name $output.$field.png
+  end
+end
+
+macro hist.range
+  if ($0 != 4)
+    echo "USAGE: hist.range (value) (Nvalue) (dvalue)"
+    break
+  end
+
+  # assume we have a small enough vector to sort and select
+  # the min and max percentile explicitly
+
+  local Nm Np No value10 value50 value90 range
+
+  set _tmp_sort = $1
+  sort _tmp_sort
+  $Nm = 0.1*_tmp_sort[]
+  $Np = 0.9*_tmp_sort[]
+  $No = 0.5*_tmp_sort[]
+
+  # 10 and 90 percentile values
+  $value10 = _tmp_sort[$Nm]
+  $value90 = _tmp_sort[$Np]
+  $value50 = _tmp_sort[$No]
+
+  # set the range to be 2x the 10 - 90 percentile range
+  $range = 3*($value90 - $value10)
+  $hist:minVal = $value50 - 0.5*$range
+  $hist:maxVal = $value50 + 0.5*$range
+  $hist:delVal = $range / 250
+  $hist:median = $value50
+
+  if ($range == 0.0)
+    # if all values are the same, generate a vector with a single value of all points
+    create $3 $hist:median {$hist:median + 1.5}
+    set $2 = zero($3) + 0.0
+    $2[0] = _tmp_sort[]
+    return    
+  end
+
+  if (isnan($range))
+    # if all values are the same, generate a vector with a single value of all points
+    create $3 0 1
+    set $2 = zero($3) + 0.0
+    $2[0] = _tmp_sort[]
+    return    
+  end
+
+  # XXX add histogram option to ignore out-of-range values
+  histogram $1 $2 $hist:minVal $hist:maxVal $hist:delVal -range $3
+  $2[0] = 0; $2[-1] = 0
+end
+
+# NOTE: in this method, match2d + reindex gives two pairs of vectors: the matched sources
+macro test.m1
+
+  $v1 = t0
+  $v2 = t1
+
+  match2d X_PSF_$v1 Y_PSF_$v1 X_PSF_$v2 Y_PSF_$v2 $MATCH_RADIUS
+  reindex Xm_$v1 = X_PSF_$v1 using index1
+  reindex Ym_$v1 = Y_PSF_$v1 using index1
+  
+  reindex Xm_$v2 = X_PSF_$v2 using index2
+  reindex Ym_$v2 = Y_PSF_$v2 using index2
+
+  set dX = Xm_$v1 - Xm_$v2
+  set dY = Ym_$v1 - Ym_$v2
+
+  vstat dX
+  vstat dY
+end
+
+# NOTE: in this method, match2d -closest + reindex -keep-unmatched gives two pairs of vectors which match the original
+# source lists (unmatched are NAN)
+macro test.m2
+
+  $v1 = t0
+  $v2 = t1
+
+  match2d X_PSF_$v1 Y_PSF_$v1 X_PSF_$v2 Y_PSF_$v2 $MATCH_RADIUS -closest
+  reindex Xm_$v1 = X_PSF_$v1 using index2 -keep-unmatched
+  reindex Ym_$v1 = Y_PSF_$v1 using index2 -keep-unmatched
+  set dX1 = Xm_$v1 - X_PSF_$v2
+  set dY1 = Ym_$v1 - Y_PSF_$v2
+  vstat dX1
+  vstat dY1
+  
+  reindex Xm_$v2 = X_PSF_$v2 using index1 -keep-unmatched
+  reindex Ym_$v2 = Y_PSF_$v2 using index1 -keep-unmatched
+  set dX2 = X_PSF_$v1 - Xm_$v2
+  set dY2 = Y_PSF_$v1 - Ym_$v2
+  vstat dX2
+  vstat dY2
+end
+