#!/usr/bin/perl
#
# Given a dataset of measurements of a star's position, 
#   do the following many times:
#
#     a) generate an "alternate" version of the measurements, using
#           random gaussian noise added to the actual values
#     b) fit the "alternate" data with a model to determine values for
#                    parallax, mu_ra, mu_dec
#  
# Use all those fits to the "alternate" datasets to estimate the 
#    uncertainty in the fitted parameters
#
# MWR 6/27/2021

$debug = 1;

srand();


# name of file in which we write (and overwrite) alternate input data
$alternate_file = "./run_alternate_input.dat";
# name of Perl script to run on the alternate input
$script = "./matrix_gxand.pl";
#$script = "./fit_gxand.pl";



#####################################################################
#    Parse args
#
if ($#ARGV != 1) {
  printf "usage: run_alternate.pl   num_iter  positions.dat \n";
  exit(1);
}
my $num_iter = $ARGV[0];
my $pos_data_file = $ARGV[1];
if ($debug > 0) {
  printf "#  num_iter        %d \n", $num_iter;
  printf "#  pos_data_file   %s \n", $pos_data_file;
}



##########################################################
#    Read data from the given datafile
#
$num_pos = 0;
open(POSITION, $pos_data_file) || die("can't open file $pos_data_file");
while (<POSITION>) {
  $line = $_;
  if ($line =~ /^#/) {
    next;
  }
  @words = split(/\s+/, $line);
  if ($#words != 4) {
    printf " invalid input line has %d words, not expected 5 \n", $#words;
    printf "%s", $line;
    exit(1);
  }
  $jd_array[$num_pos] = $words[0];
  $ra_pos_array[$num_pos] = $words[1];
  $ra_uncert_array[$num_pos] = $words[2];
  $dec_pos_array[$num_pos] = $words[3];
  $dec_uncert_array[$num_pos] = $words[4];
  $num_pos++;
}
close(POSITION);
if ($debug > 0) {
  printf " read %d entries in file %s \n", $num_pos, $pos_data_file;
}


##################################################################
#  debugging only
if (0 == 1) {
  for ($i = 0; $i < 10; $i++) {
    $x = gaussian_rand();
    printf " rand value %3d is %6.3f \n", $i, $x;
  }
}


#########################################################################
#  First, we run the script on the actual measurements.
#      This will give us the "best" values for each parameter
#
$cmd = "perl $script $pos_data_file";
$ret = exec_cmd($cmd);
# parse the output to extract the values of the fitted parameters
($best_parallax, $best_pmra, $best_pmdec, $best_K, $best_L) = extract_params($ret);
if ($debug > 0) {
  printf "#  best_parallax = %10.2f \n", $best_parallax;
  printf "#  best_pmra     = %10.2f \n", $best_pmra;
  printf "#  best_pmdec    = %10.2f \n", $best_pmdec;
  printf "#  best_K        = %10.2f \n", $best_K;
  printf "#  best_L        = %10.2f \n", $best_L;
}



##########################################################################
#     main loop: create alternate data, fit model, record fitted values
#

for ($i_iter = 0; $i_iter < $num_iter; $i_iter++) {
  if ($debug > 0) {
    printf "  top of loop, iter %d \n", $i_iter;
  }

  # create a file with alternate measurements
  open(ALTERNATE_FILE, ">$alternate_file") || die("can't open file $alternate_file");
  for ($i_pos = 0; $i_pos < $num_pos; $i_pos++) {
    
    $alt_jd = $jd_array[$i_pos];
    $alt_ra_pos = $ra_pos_array[$i_pos] + 
                       $ra_uncert_array[$i_pos]*gaussian_rand();
    $alt_ra_uncert = $ra_uncert_array[$i_pos];
    $alt_dec_pos = $dec_pos_array[$i_pos] + 
                       $dec_uncert_array[$i_pos]*gaussian_rand();
    $alt_dec_uncert = $dec_uncert_array[$i_pos];

    printf ALTERNATE_FILE "%12.5f %10.1f %10.1f %10.1f %10.1f \n",
           $alt_jd, $alt_ra_pos, $alt_ra_uncert, 
           $alt_dec_pos, $alt_dec_uncert;

  }
  close(ALTERNATE_FILE);

  # and now run the script which fits a model to the alternate 
  #    measurements
  $cmd = "perl $script $alternate_file";
  $ret = exec_cmd($cmd);

  # parse the output to extract the values of the fitted parameters
  ($this_parallax, $this_pmra, $this_pmdec, $this_K, $this_L) = extract_params($ret);
  if ($debug > 0) {
    printf " iter %3d yields parallax %8.2f \n", $i_iter, $this_parallax;
    printf " iter %3d yields pmra     %8.2f \n", $i_iter, $this_pmra;
    printf " iter %3d yields pmdec    %8.2f \n", $i_iter, $this_pmdec;
    printf " iter %3d yields K        %8.2f \n", $i_iter, $this_K;
    printf " iter %3d yields L        %8.2f \n", $i_iter, $this_L;
  }

  $parallax_array[$i_iter] = $this_parallax;
  $pmra_array[$i_iter] = $this_pmra;
  $pmdec_array[$i_iter] = $this_pmdec;
  $K_array[$i_iter] = $this_K;
  $L_array[$i_iter] = $this_L;

}


######################################################################
#     Now determine the mean and stdev of all those values

$parallax_mean = calc_mean($num_iter, @parallax_array);
$parallax_stdev = calc_sigma($num_iter, @parallax_array);
$pmra_mean = calc_mean($num_iter, @pmra_array);
$pmra_stdev = calc_sigma($num_iter, @pmra_array);
$pmdec_mean = calc_mean($num_iter, @pmdec_array);
$pmdec_stdev = calc_sigma($num_iter, @pmdec_array);
$K_mean = calc_mean($num_iter, @K_array);
$K_stdev = calc_sigma($num_iter, @K_array);
$L_mean = calc_mean($num_iter, @L_array);
$L_stdev = calc_sigma($num_iter, @L_array);

if ($debug > 0) {
  printf "#  parallax     mean %10.2f  stdev %10.2f \n", 
                $parallax_mean, $parallax_stdev;
  printf "#  pm_ra        mean %10.2f  stdev %10.2f \n", 
                $pmra_mean, $pmra_stdev;
  printf "#  pm_dec       mean %10.2f  stdev %10.2f \n", 
                $pmdec_mean, $pmdec_stdev;
  printf "#  K            mean %10.2f  stdev %10.2f \n", 
                $K_mean, $K_stdev;
  printf "#  pm_dec       mean %10.2f  stdev %10.2f \n", 
                $L_mean, $L_stdev;
}


######################################################################
#     Print out final results

printf " parallax   %10.2f  +/-  %10.2f \n", 
               $best_parallax, $parallax_stdev;
printf " pm_ra      %10.2f  +/-  %10.2f \n", 
               $best_pmra, $pmra_stdev;
printf " pm_dec     %10.2f  +/-  %10.2f \n", 
               $best_pmdec, $pmdec_stdev;
printf " K          %10.2f  +/-  %10.2f \n", 
               $best_K, $K_stdev;
printf " L          %10.2f  +/-  %10.2f \n", 
               $best_L, $L_stdev;








exit 0;









######################################################################
# PROCEDURE: gaussian_rand
#
# DESCRIPTION: Returns a random number drawn from a gaussian (normal)
#              distribution with mean 0.0 and standard deviation 1.0.
#
#              Stolen from 
#
#                 http://www.unix.com.ua/orelly/perl/cookbook/ch02_11.htm
#
#
sub gaussian_rand {
    my ($u1, $u2);  # uniformly distributed random numbers
    my $w;          # variance, then a weight
    my ($g1, $g2);  # gaussian-distributed numbers

    do {
        $u1 = 2 * rand() - 1;
        $u2 = 2 * rand() - 1;
        $w = $u1*$u1 + $u2*$u2;
    } while ( $w >= 1 );

    $w = sqrt( (-2 * log($w))  / $w );
    $g2 = $u1 * $w;
    $g1 = $u2 * $w;
    # could return both if wanted
    #   return wantarray ? ($g1, $g2) : $g1;
    # but I need just one
    return($g1);
}







############################################################################
# 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 > 0) {
    printf "cmd is ..$cmd.. \n";
  }
  $ret = `$cmd`;
  if ($debug > 0) {
    printf "ret is ..$ret.. \n";
  }

  return($ret);
}



############################################################################
# PROCEDURE: extract_params
#
# DESCRIPTION: Given a set of lines which is the output of a script,
#              which looks like this:
#
#                  # solution   parallax     283.5  
#                  # solution   pm_ra       2842.9  
#                  # solution   pm_dec       413.2 
#                  # solution   K             20.1 
#                  # solution   L            -53.8 
#
#              read the values of the parallax, pm_ra, pm_dec, K, L.
#              Return them in that order.
#
# RETURNS:
#              (parallax, pm_ra, pm_dec, K, L)
#
#
sub extract_params {

  my($result_str);
  my($parallax, $pm_ra, $pm_dec, $K, $L);
  my(@lines, @words);
  my($bad_value);

  $result_str = $_[0];

  # initialize the three parameters to values that are impossible
  $bad_value = -99999.9;
  $parallax = $bad_value;
  $pm_ra = $bad_value;
  $pm_dec = $bad_value;
  $K = $bad_value;
  $L = $bad_value;


  @lines = split(/\n/, $result_str);
  foreach $line (@lines) {
    @words = split(/\s+/, $line);
    if ($words[2] eq "parallax") {
      $parallax = $words[3];
    }
    if ($words[2] eq "pm_ra") {
      $pm_ra = $words[3];
    }
    if ($words[2] eq "pm_dec") {
      $pm_dec = $words[3];
    }
    if ($words[2] eq "K") {
      $K = $words[3];
    }
    if ($words[2] eq "L") {
      $L = $words[3];
    }
  }

  # check to see if we read values for each param
  if (($parallax == $bad_value) || ($pm_ra == $bad_value) ||
      ($pm_dec == $bad_value) || ($K == $bad_value) || ($L == $bad_value))  {
    printf "extract_params: unable to read one or more params ... quitting \n";
    exit(1);
  }

  return($parallax, $pm_ra, $pm_dec, $K, $L);

}






########################################################################
# 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);
}

