#!/usr/bin/perl
#
# Run the "fit_gxand.pl" script repeatedly, while making modifications
#   to its input file, in order to perform a jackknife analysis
#   of its output.
#
# MWR 6/16/2021

$debug = 0;

# name of script which finds best-fit parameters
$script = "./fit_gxand.pl";
# file with all the measurements of stellar position, one per line
$star_file = "./gxand_pos.dat";
# file with one measurement removed
$temp_file = "./run_jackknife.tmp1";
# we send output of the fitting script into this temp output file
$fit_output = "./run_jackknife.tmp2";


###################################################################
#     Count the number of lines in input datafile
open(STAR_FILE, "$star_file") || die("can't open file $star_file");
$num_measures = 0;
while (<STAR_FILE>) {
  $line = $_;
  if ($line =~ /^#/) {
    next;
  }
  $line_array[$num_measures] = $line;
  $num_measures++;
}
close(STAR_FILE);
if ($debug > 1) {
  printf " input file $star_file has %d lines \n", $num_measures;
}


#################################################################
#     Loop over the input file, removing one line at a time
#
for ($i_measure = 0; $i_measure < $num_measures; $i_measure++) {

  open(TEMP_FILE, ">$temp_file") || die("can't create file $temp_file");
  for ($j = 0; $j < $num_measures; $j++) {
    if ($j != $i_measure) {
      printf TEMP_FILE "%s", $line_array[$j];
    }
  }
  close(TEMP_FILE);

  # run the fitting program on the modified input
  $cmd = "perl $script $temp_file > $fit_output";
  $ret = exec_cmd($cmd);

  # read the fitted parameters and save them for later
  open(FIT_OUTPUT, "$fit_output") || die("can't open file $fit_output");
  while (<FIT_OUTPUT>) {
    $line = $_;
    @words = split(/\s+/, $line);
    if (($words[1] eq "param") && ($words[3] eq "val")) {
      $param_index = $words[2];
      $param_val = $words[4];
      $param_array[$i_measure][$param_index] = $param_val;
    }
  }
  close(FIT_OUTPUT);
  if ($debug > 1) {
    for ($i = 0; $i < 5; $i++) {
      printf "  iter %2d  param %2d has val %7.1f \n", 
            $i_measure, $i, $param_array[$i_measure][$i];
    }
  }

} #  end of loop over measurements


#################################################################
#     Analyze the results of the individual trials
#

for ($i = 0; $i < 5; $i++) {

  for ($i_measure = 0; $i_measure < $num_measures; $i_measure++) {
    $sigma_array[$i_measure] = $param_array[$i_measure][$i];
    if ($debug > 1) {
      printf "   param %2d  measure %2d   val %7.1f \n", 
            $i, $i_measure, $sigma_array[$i_measure];
    }
  }
  $sigma = calc_sigma($num_measures, @sigma_array);
  printf "  param %2d has sigma %7.1f \n", $i, $sigma;

}





exit 0;






##############################################################################
# PROCEDURE: exec_cmd
#
# DESCRIPTION: Execute the given shell command line.  If the $debug flag
#              is set, we print to stdout the command line, and
#              also print out the result string.
#
# RETURNS:
#              the result of the command
#
#
sub exec_cmd {

  my($cmd, $ret);

  $cmd = $_[0];

  if ($debug > 1) {
    printf "cmd is ..$cmd.. \n";
  }
  $ret = `$cmd`;
  if ($debug > 1) {
    printf "ret is ..$ret.. \n";
  }

  return($ret);
}




########################################################################
# PROCEDURE: calc_mean
#
# DESCRIPTION: given the number of elements in an array, and the
#              array itself, calculate the mean value of the array.
#
# RETURNS:
#              the mean value
#
sub calc_mean {

  my($num, @array);
  my($i, $sum);
  my($mean);

  $num = $_[0];
  for ($i = 0; $i < $num; $i++) {
    $array[$i] = $_[$i+1];
  }

  if ($num == 0) {
    $mean = 0.0;
    return($mean);
  }
 
  # printf "%2d ", $num;
  for ($i = 0; $i < $num; $i++) {
  #   printf "%6.3f ", $array[$i];
  }

  $sum = 0.0;
  for ($i = 0; $i < $num; $i++) {
    $sum += $array[$i];
  }
  $mean = $sum/$num;
  # printf(" %7.4f \n", $mean);

  return($mean);
}


########################################################################
# PROCEDURE: calc_sigma
#
# DESCRIPTION: given the number of elements in an array, and the
#              array itself, calculate the standard deviation from the
#              mean of an array
#
# RETURNS:
#              the stdev
#
sub calc_sigma {

  my($num, @array);
  my($i, $sum);
  my($sumsq);
  my($mean, $stdev);

  $num = $_[0];
  for ($i = 0; $i < $num; $i++) {
    $array[$i] = $_[$i+1];
  }

  if ($num < 2) {
    $stdev = 0.0;
    return($stdev);
  }
 
  # printf "%2d ", $num;
  for ($i = 0; $i < $num; $i++) {
  #   printf "%6.3f ", $array[$i];
  }

  $sum = 0.0;
  $sumsq = 0.0;
  for ($i = 0; $i < $num; $i++) {
    $sum += $array[$i];
    $sumsq += $array[$i]*$array[$i];
  }
  $mean = $sum/$num;
  $stdev = sqrt(($sumsq - $num*$mean*$mean)/($num - 1));
  # printf(" %7.4f \n", $stdev);

  return($stdev);
}

