#!/usr/bin/perl
#
# Read a datafile.
#    Figure out the noise in the data.
#
# Iterate many times
#    Generate a known signal
#    Add random noise to the signal
#    Attempt to fit the noisy signal with a model
#    Note the result
#
# Compute stdev in all the results


$debug = 0;

use POSIX;

# read data from this file
#
#  No, we just estimate the stdev of noise here
$noise_stdev = 0.3;

# iterate this many times
$N_iter = 300;

$amp = 1.0;
$fwhm = 1.0;
$sigma = $fwhm / 2.35;


# Generate a signal with a known time
$signal_time = 12.89;
$start_t = 0;
$end_t = 20;
$dt = 0.1;

$ndata = 0;
for ($t = $start_t; $t < $end_t; $t += $dt) {
     
  $signal_t_array[$ndata] = $t;
  $arg = 0.0 - ($t - $signal_time)*($t - $signal_time)/(2.0*$sigma*$sigma);
  $signal_v_array[$ndata] = $amp*exp($arg);
  $ndata++;
}
if ($debug > 0) {
  for ($i = 0; $i < $ndata; $i++) {
    printf "  i %5d  t %8.3f  signal %8.3f \n", 
         $i, $signal_t_array[$i], $signal_v_array[$i];
  }
}


# Now, we iterate N_iter times
for ($i_iter = 0; $i_iter < $N_iter; $i_iter++) {

  # create an instance of the signal, with random gaussian noise added
  for ($i = 0; $i < $ndata; $i++) {
    $instance_t_array[$i] = $signal_t_array[$i];
    $instance_v_array[$i] = $signal_v_array[$i] +
                                $noise_stdev*gaussian_rand();
  }

  # now, we try to fit the instance
  $best_time = -1;
  $min_k = 1e6;

    # we loop over possible guessed times
    $start_guess = 1;
    $end_guess = 19;
    $d_guess = 0.1;
    for ($guess_t = $start_guess; $guess_t < $end_guess; $guess_t += $d_guess) {
    
         if ($debug > 0) {
           printf " guess_t is now %8.3f \n", $guess_t;
         }
         $k_sum = 0;
    
         # we now generate a model with this guessed time
         #    and compare to the data, compute k (like chi)
         $start_t = 0;
         $end_t = 100;
         $dt = 1.0;
    
         for ($i = 0; $i < $ndata; $i++) {
         
           $t = $instance_t_array[$i];
           $arg = 0.0 - ($t - $guess_t)*($t - $guess_t)/(2.0*$sigma*$sigma);
           $model_v = $amp*exp($arg);
         
           $diff = $model_v - $instance_v_array[$i];
           $diff_sq = $diff*$diff;
           $k_sum += $diff_sq;
    
           if ($debug > 0) {
             printf " %8.3f %8.3f  %8.3f    k now %8.3f \n", 
                $t, $model_v, $data_v_array[$i], $k_sum;
           }
         
         }
    
         # now we can spit out the value of guess and the k_sum
         if ($debug > 0) {
           printf "  guess_t %8.3f  k_sum %8.3f \n", $guess_t, $k_sum;
         }

         if ($k_sum < $min_k) {
           $min_k = $k_sum;
           $best_time = $guess_t;
           if ($debug > 0) {
             printf " new best guess %8.3f  yields k %8.3f \n",
                    $best_time, $k_sum;
           }
         }
    
    }

    # what was the best guess?
    printf " iter %3d   best t = %8.3f    k = %8.3f \n", 
             $i_iter, $best_time, $min_k;


}



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

