
list tests
 test1
 test2
 test3
 test4
 test5
 test6
 test7
 test8
end

## NOTE: this is not a TAP test
macro mean.test
  if ($0 != 2)
    echo "USAGE: mean.test (errorbar)"
    break
  end

  create x 0 100

  delete C0fit dC0fit Cnfit Mfit dMfit
  for i 0 500
    gaussdev y x[] 0 1.0   ; # mean of 0, sigma of 1
    
    set dy = 1.0 + zero(x) ; # set error to be correct
    fit1d_irls -q x y 0 -dy dy -wt wt -mask mask -use-median -Nboot 1000

    concat $C0    C0fit
    concat $dC0  dC0fit

    vstat -q y
    concat $MEAN Mfit
    concat {$SIGMA / sqrt($NPTS)} dMfit
  end

  vstat -q C0fit;  set MeanC  = $MEAN; set StdevC = $SIGMA
  vstat -q dC0fit; set MeanDC = $MEAN
  vstat -q Mfit;   set MeanM  = $MEAN; set StdevM = $SIGMA
  vstat -q dMfit;  set MeanDM = $MEAN

  fprintf "%7.3f %7.3f" $MeanC $MeanM
  fprintf "%7.3f %7.3f" $StdevC $StdevM
end

macro cliptest.plot
  if ($0 != 3)
    echo "USAGE: cliptest.plot (errorbar) (output)"
    break
  end

  set dy = $1 + zero(x) ; # set error to be correct
  fit1d_irls x y 0 -dy dy -wt wt -mask mask -use-median

  subset ygood = y where not(mask)

  vstat y
  vstat ygood

  histogram wt Nwt 0 1 0.005 -range dw
  vstat wt

  subset wts = wt where (wt < 0.3*$MEDIAN)
  echo "0.3 of MEDIAN : wts[]"

  subset wts = wt where (wt < 0.3*$MEAN)
  echo "0.3 of MEAN : wts[]"

  dev -n 0; resize 1200 600
  label -fn courier 24; 
  section a 0.4 0.0 0.6 1.0; lim x y; clear; box -ypad 0.1 -xpad 3.2 -labelpadx 2.5 -labels 1001; label -x sequence +y value +x "red: clipped, blue: unclipped"; 
  plot x y; plot -c blue -pt 7 x y where not(mask); plot -c red -pt 7 -lw 2 x y where mask

  section b 0.0 0.0 0.4 1.0; lim dw Nwt; box +ypad 0.1 -xpad 3.2 -labelpadx 2.5; label -x "weight" -y "Npoints"; plot -x 1 dw Nwt
  png -name $2
end

## NOTE: this is not a TAP test
macro cliptest

  create x 0 100
  gaussdev y x[] 0 1.0   ; # mean of 0, sigma of 1

  if (1)
    y[10] = 5.0
    y[20] = 7.0
    y[30] = 4.0
  end

  cliptest.plot 1.0 irls.sample.v0.png

  cliptest.plot 0.1 irls.sample.v1.png

  cliptest.plot 0.01 irls.sample.v2.png
end

## NOTE: this is not a TAP test
macro outlier.test
  if ($0 != 2)
    echo "USAGE: outlier.test (errorbar)"
    break
  end

  create x 0 100

  delete C0fit dC0fit Cnfit Mfit dMfit
  for i 0 300
    gaussdev y x[] 0 1.0   ; # mean of 0, sigma of 1
    
    if (0)
      y[10] = 5.0
      y[20] = 7.0
      y[30] = 4.0
    end
    
    set dy = $1 + zero(x) ; # set error to be correct
    fit1d_irls -q x y 0 -dy dy -wt wt -mask mask -use-median -Nboot 1000
    concat $C0 C0fit
    concat $dC0 dC0fit
    concat $Cnfit Cnfit

    vstat -q y
    concat $MEAN Mfit
    concat {$SIGMA / sqrt($NPTS)} dMfit
  end

  vstat C0fit
  vstat dC0fit
  vstat Cnfit
  vstat Mfit
  vstat dMfit
end

# fit a line without errors
macro test1
 $PASS = 1
 break -auto off

 delete -q x y

 create x 0 100
 set y = 3 + 5*x
 fit -q x y 1

 if ($Cn != 1)
   $PASS = 0
 end
 if (abs($C0 - 3) > 1e-5)
   $PASS = 0
 end
 if (abs($C1 - 5) > 1e-5)
   $PASS = 0
 end
end

# fit a line with errors
macro test2
 $PASS = 1
 break -auto off

 delete -q x y dy

 create x 0 100
 set dy = 0.1*rnd(x) - 0.05
 set y = 3 + 5*x + dy
 fit -q x y 1

 if ($Cn != 1)
   $PASS = 0
   echo "wrong number of elements : $Cn"
 end
 if (abs($C0 - 3) > 0.06)
   $PASS = 0
   echo "wrong value for C0 : $C0"
 end
 if (abs($C1 - 5) > 0.06)
   $PASS = 0
   echo "wrong value for C1 : $C1"
 end
end

# fit a line with errors and weights
macro test3
 $PASS = 1
 break -auto off

 delete -q x y dy

 create x 0 100
 set dy = 0.1*rnd(x) - 0.05
 set y = 3 + 5*x + dy
 set dy = 0.1 + zero(x)
 fit -q x y 1 -dy dy

 if ($Cn != 1)
   $PASS = 0
 end
 if (abs($C0 - 3) > 0.06)
   $PASS = 0
 end
 if (abs($C1 - 5) > 0.06)
   $PASS = 0
 end
end

# fit a line with errors, weights, and outliers 
macro test4
 $PASS = 1
 break -auto off

 delete -q x y dy

 create x 0 100
 set dy = 0.1*rnd(x) - 0.05
 set y = 3 + 5*x + dy
 set dy = 0.1 + zero(x)
 y[5] = 23
 y[20] = -10
 y[50] = 0.0
 fit -q x y 1 -dy dy -clip 3 3

 if ($Cn != 1)
   $PASS = 0
 end
 if ($Cnv != 97)
   $PASS = 0
 end
 if (abs($C0 - 3) > 0.06)
   $PASS = 0
 end
 if (abs($C1 - 5) > 0.06)
   $PASS = 0
 end
end

# fit a quadratic without errors
macro test5
 $PASS = 1
 break -auto off

 delete -q x y

 create x 0 100
 set y = 3 + 5*x - 4*x^2
 fit -q x y 2

 if ($Cn != 2)
   $PASS = 0
 end
 if (abs($C0 - 3) > 1e-5)
   $PASS = 0
 end
 if (abs($C1 - 5) > 1e-5)
   $PASS = 0
 end
 if (abs($C2 + 4) > 1e-5)
   $PASS = 0
 end
end

# fit a quadratic with errors
macro test6
 $PASS = 1
 break -auto off

 delete -q x y dy

 create x 0 100
 set dy = 0.1*rnd(x) - 0.05
 set y = 3 + 5*x - 4*x^2 + dy
 fit -q x y 2

 if ($Cn != 2)
   $PASS = 0
 end
 if (abs($C0 - 3) > 0.06)
   $PASS = 0
 end
 if (abs($C1 - 5) > 0.06)
   $PASS = 0
 end
 if (abs($C2 + 4) > 0.06)
   $PASS = 0
 end
end

# fit a quadratic with errors and weights
macro test7
 $PASS = 1
 break -auto off

 delete -q x y dy

 create x 0 100
 set dy = 0.1*rnd(x) - 0.05
 set y = 3 + 5*x - 4*x^2 + dy
 set dy = 0.1 + zero(x)
 fit -q x y 2 -dy dy

 if ($Cn != 2)
   $PASS = 0
 end
 if (abs($C0 - 3) > 0.06)
   $PASS = 0
 end
 if (abs($C1 - 5) > 0.06)
   $PASS = 0
 end
 if (abs($C2 + 4) > 0.06)
   $PASS = 0
 end
end

# fit a quadratic with errors, weights, and outliers 
macro test8
 $PASS = 1
 break -auto off

 delete -q x y dy

 create x 0 100
 set dy = 0.1*rnd(x) - 0.05
 set y = 3 + 5*x - 4*x^2 + dy
 set dy = 0.1 + zero(x)
 y[5] = 23
 y[20] = -10
 y[50] = 0.0

 # it takes 4 iterations to successfully reject the outliers above...
 fit -q x y 2 -dy dy -clip 3 4

 if ($Cn != 2)
   $PASS = 0
 end
 if ($Cnv != 97)
   $PASS = 0
 end
 if (abs($C0 - 3) > 0.06)
   $PASS = 0
 end
 if (abs($C1 - 5) > 0.06)
   $PASS = 0
 end
 if (abs($C2 + 4) > 0.06)
   $PASS = 0
 end
end
