#!/usr/bin/perl -w
# CmeasF script (measure component outputs and run times)
#  Use .pscf file to process all components and create .ccft files
#  that contain their specs
#
use samplingf;
use dec15;
$eps = &dec15::geteps();

#get options
$Coption = 0;
$Loption = 1;  #assume linear
$Voption = 0;
OPT: for ($i = 0;;$i++) { #over all options
  unless (defined($ARGV[$i])) {last OPT;}
  if ($ARGV[$i] eq "-V") {
    $Voption = 1;
  }
  if ($ARGV[$i] eq "-C") {
    $Coption = 1; $Loption = 0;
  }
  if ($ARGV[$i] eq "-L") {
    $Coption = 0; $Loption = 1;
  }
}

#process system configuration file
$sys_desc = "system.pscf" ;
open(SYSTEM, $sys_desc ) || die "could not open the file", $sys_desc ;
$comp_name = <SYSTEM>; #discard polish line

#store component .ccf names
@components = ();
$i = 0;
while ($comp_name = <SYSTEM>) { #read to end
  chop($comp_name);
  $components[$i] = $comp_name;
  $i++
}
if ($Coption) { #must also do ident file
  $components[$i++] = "ident.ccf";
}
close SYSTEM ;

#debug
#foreach $comp (@components) {print $comp."\n"; }

########################################subroutine#######################
sub fitline {
# best-fit line in the sense of least squares of differences
# parameters:  $_[0] - X-vector; $_[1] - Y-vector; $_[2] - X-ave; $_[3] - Y-ave; $_[4] - Count
my($xvec, $yvec, $xave, $yave, $N, $Slope, $Inter, $dev, $cross, $var, $i);
$xvec = $_[0];
$yvec = $_[1];
$xave = $_[2];
$yave = $_[3];
$N = $_[4];
$cross = 0;
$var = 0;
for ($i=0; $i<$N; $i++) {
  $cross += @$xvec[$i] * @$yvec[$i];
  $var += @$xvec[$i] * @$xvec[$i];
}
if ($N == 1) {
  $Slope = 0;
}
else {
  $Slope = ($cross/$N - $yave*$xave)/($var/$N - $xave*$xave);
}
if (abs($cross/$N - $yave*$xave) <= $eps) {
  #debug
  #$a = $cross/$N;
  #$b = $yave*$xave;
  #$c = $var/$N;
  #$d = $xave*$xave;
  #print "$Slope $a - $b / $c - $d\n";
  #end debug
  $Slope = 0;
}
$Inter = $yave - $Slope*$xave;
$dev = 0;
for ($i=0; $i<$N; $i++) {
  # It should be possible to do this in the first pass as in std dev -- fix TBD
  $dev += ($Slope*@$xvec[$i] + $Inter -@$yvec[$i])**2;
}
$dev = sqrt($dev/$N)/($yave+$eps); #rms deviation from the line
return $Slope, $Inter, $dev;
#  return (slope,intercept,deviation) of the fitted line
}
#########################################################################

sub printf15 { # print line from a ccf file (9 values)
  #first parameter is a file name to write to
  #needed to avoid loss of precision when a float is near an integer in Perl
  my($file, $t, $N);
  $file = shift; #file handle or 0 if none
  for ((1..2)) {
    $t = shift;
    printf $file "%s",&dec15::dump15s($t);
  }
  $t = shift;
  if ($t eq "L") {
    $N = 6;
  }
  else {
    $N = 4;
  }
  printf $file "%s ",$t;
  for ((1..$N)) {
    $t = shift;
    printf $file "%s",&dec15::dump15s($t);
  }
  printf $file "\n";
}
#########################################################################

# get ready to print results
#for each component process its file and execute test
CCFL: foreach $comp (@components) {
  open(COMP, $comp) || die "could not open the file ", $comp;
  # get the name of the executable
  chop($comp_name = <COMP>) ;
  if ($comp_name eq "theory") { #not a real component, skip
    next CCFL;
  }
  #real component, measure it
  #see if the executable is there
  die "component ", $comp_name, " can't be executed." unless (-x $comp_name);
  @wholename = split(/\./,$comp);
  #debug
  #print "@wholename\n";
  $comptheory = $wholename[0].".ccft";
  open(COMPTH, ">".$comptheory)  || die "could not create theory file ", $comptheory;
  print COMPTH "theory\n";
  #read subdomain information
  while($linearr = <COMP>) {
    #read line from .ccf file
    @line = split(' ',$linearr);
    $sub_start = $line[0] ;
    $sub_end = $line[1] ;
    $test_cnt = $line[2] ;
    &samplingf::init($sub_start, $sub_end, $test_cnt);
    $cnt = 0 ;
    $input = &samplingf::next ;
    #for each subdomain get equi-spaced input and call components with this
    @x_vals = ();  #inputs
    $x_bar = 0;  #average
    @y_vals = ();  #functional outputs
    $y_bar = 0;
    $y_sqrs = 0;
    @t_vals = ();  #run times
    $t_bar = 0;
    $t_sqrs = 0;
    while($input ne ".") {
      $x_vals[$cnt] = $input;
      $x_bar += $input;
      $sed = "sed 's/^/out:  /'";  # command to reproduce stdin to stdout with a leading "out:  "
      $cmd = "(echo ".$input."|".$comp_name."|".$sed.") 2>&1|" ;
      open(PI,$cmd);
      for($k=1; $k<=2; $k++) {
        @outstrm = split(' ',<PI>);
        #The output stream is: 1) RT of component; 2) "out:  ".output of component
        # just to be safe, check for the "out" and don't use the order
        #debug
        #print "$outstrm[0] \n";
        #print "$outstrm[1] \n";
        if ($outstrm[0] =~ /out:/) { #functional output, record
          $y = $outstrm[1];
          $y_vals[$cnt] = $y;
          $y_bar += $y;
          $y_sqrs += $y*$y;
        }
        else { #run time; record 
          $t = $outstrm[0];
          $t_vals[$cnt] = $t;
          $t_bar += $t;
          $t_sqrs += $t*$t;
        }
      }
      close(PI);
      $input = &samplingf::next ;
      $cnt++ ;
    } #end loop over inputs
    $x_bar /= $cnt;
    $y_bar /= $cnt;
    $t_bar /= $cnt;
    #debug
    #print "$cnt $y_bar\n";
    if ($Loption || $cnt == 1) { #don't compute stdev
      $y_sdev = 0;
      $t_sdev = 0;
    }
    else {
      $y_sdev = sqrt(abs(($y_sqrs - $cnt*$y_bar*$y_bar)/($cnt -1)))/($y_bar+$eps);
      $t_sdev = sqrt(abs(($t_sqrs - $cnt*$t_bar*$t_bar)/($cnt -1)))/($t_bar+$eps);
    }
#debug
#print "N: $cnt, xave: $x_bar, yave: $t_bar, stdev: $y_sdev\n";
    if ($Loption) {
      ($R_slope, $R_inter, $R_dev) = fitline(\@x_vals, \@t_vals, $x_bar, $t_bar, $cnt);
      #debug
      #print "line:  $R_slope x + $R_inter\n";
      ($V_slope, $V_inter, $V_dev) = fitline(\@x_vals, \@y_vals, $x_bar, $y_bar, $cnt);
      #debug
      #print "line:  $V_slope x + $V_inter\n";
      printf15(COMPTH, $sub_start, $sub_end,"L",$V_slope,$V_inter,$R_slope,$R_inter,$V_dev,$R_dev)
;
    }
    if ($Coption) {
      printf15(COMPTH, $sub_start, $sub_end, "C", $y_bar, $t_bar, $y_sdev, $t_sdev);
    }
    #debug
    #print "sigma: $y_sdev, rsigma: $t_sdev ";
    #if ($Loption) {
      #print "rms: $V_dev, Rrms: $R_dev";
    #}
    #print "\n";
    #end debug
  } #end loop over subdomains, one component
  close COMP ;
} #end loop over components
