#! /usr/bin/perl -w
# Xcute  <.ccf[t] file> [<start> <stop> <count>]  runs and plots a component behavior.
#
#  If no <file> extension is given, ".ccf" is assumed.
#  The file may point to any of:
#    -By-hand descriptive .ccf file (for which the .ccft measured file must exist)
#    -'theoryN.ccf' calculated file
#    -Measured .ccft file
#  In the latter two cases, only the approximate files are plotted.
#
#  In the first case, measurements are made and plotted as well as the
#    approximation, and rms errors displayed (for this set of
#    executions, not the ones that made the .ccft file).  
#    Measurements not repeated if files current.
#
#  Default for the optional count is about 5 times the number of subdomains.
#    The three optional values must all be given if any is.
# 
#  Measurements recorded in a file "<ccf name>.dat" in lines:  input output runtime
#    This file is used to date "current" measurements
#
use sampling;
use component;
use util;

unless (defined($ARGV[0])) {
  die "A .ccf[t] file must be given";
}
$ccf = $ARGV[0];
@spfile = split('\.',$ccf);
if (not defined($spfile[1])) { #bare file, assume .ccf
  $ccf .= ".ccf";
}
$base = $spfile[0];  #for naming measurement data file
$ccftFile = $base.".ccft";
open(CCF, $ccf) || die "given file ",$ccf," not found";
$nameline = <CCF>;
@nameand = split(' ',$nameline);
$name = $nameand[0];
$parmsupplied = 0;
if (defined($ARGV[1])) {
  die "All <start> <stop> <count> must be given" unless (defined($ARGV[3]));
  $start = $ARGV[1];
  $stop = $ARGV[2];
  $count = $ARGV[3];
  $parmsupplied = 1;
}
#record the subdomains and get the domain extent
#All file formats have this data in the same place
$k = 0;
while ($line = <CCF>) {
  @vals = split(" ",$line);
  $LB[$k] = $vals[0];
  $UB[$k] = $vals[1];
  $k++;
}
$left = $LB[0];
$right = $UB[$k-1];
$numSubd = $k;
if ($parmsupplied) { #check reasonableness
  die "given range [$start,$stop) for $count points doesn't match .ccf values"
    unless ($start >= $left && $stop <= $right);
} else { #no params given; start/stop adjusted to all subdomains
  $start = $left;
  $stop = $right;
  $count = 5*$numSubd;
}

if ($name ne "theory") {  #.ccf for a real component, measure
  unless (-e $name) {die "No code file $name to measure";}
  # BUT: the corresponding xxx.ccft file must also exist 
  unless (-e $ccftFile) {
    die "an approximation file (.ccft) must exist for the component";
  }
  $theoryonly = 0;
  $testfile = $base.".dat";  #holds measured data
  unless ($parmsupplied) {
    $testfile = $base.".alldat";
  }
  if (!$parmsupplied && -e $testfile && (-M $testfile < -M $ccftFile) && (-M $testfile < -M $name) && (-M $testfile < -M $ccf)) { #current, no measure unless explicit range given
    $goodata = 0;
    goto PLOT;  #cheap shot
  }
#print "[$left, $right)  $k $count\n";  #debug
  
  $goodata = 1;
  printf "Sampling [%.2f, %.2f) for %d points\n",$start,$stop,$count;
  open(M,">$testfile"); #"$base.dat" set above
  $y_sqrs = 0;
  $y_bar = 0;
  $t_sqrs = 0;
  $t_bar = 0;
  &sampling::init($start,$stop,$count);
  $mcount = 0;
  while (($input = &sampling::next) ne ".") {
    ($meas, $rmeas) = component::excomp($input, $name);
    $inputs[$mcount] = $input;
    $outputs[$mcount] = $meas;
    $runtimes[$mcount] = $rmeas;
    $y_bar += $meas;
    $y_sqrs += $meas*$meas;
    $t_sqrs += $rmeas*$rmeas;
    $t_bar += $rmeas;
    print M "$input $meas $rmeas\n";
    $mcount++;
  } #end loop over inputs
  $y_bar /= $count;
  $yrms = sqrt($y_sqrs/$count);
  $t_bar /= $count;
  $trms = sqrt($t_sqrs/$count);
  close(M);
} else { #not real component file, must be theory only (.ccf or .ccft)
  $ccftFile = $ccf;  # put name back to original
  $theoryonly = 1;
}

PLOT:
# May or may not have measurements; $theoryonly is the flag
unless ($theoryonly) { #measurement stuff
  unless ($goodata) { #read in previous measurement data
    open (M,"<$testfile"); #exists from above
    $km = 0;
    $yrms = 0;
    $trms = 0;
    while ($line = <M>) {
      @sline = split(' ',$line);  #input output runtime
      $inputs[$km] = $sline[0];
      $outputs[$km] = $sline[1];
      $yrms += $sline[1]**2;
      $runtimes[$km] = $sline[2];
      $trms += $sline[2]**2;
      $km++
    }
    close M;
    $mcount = $km;  #number of measured points
    $yrms = sqrt($yrms/$mcount);
    $trms = sqrt($trms/$mcount);
  } #end reading measurements from file
  # data in arrays is now valid unless $theoryonly
#print "RMS: $yrms\n";  #debug
  #write measured plot files
  open (OUT, ">output.dat");
  open (RUN, ">runtime.dat");
  for ($k=0; $k<$mcount; $k++) {
    print OUT "$inputs[$k] $outputs[$k]\n";
    print RUN "$inputs[$k] $runtimes[$k]\n";
  }
  close OUT;
  close RUN;
} #end measurement stuff

open(CCF,"<".$ccftFile) or die "No theory file $ccftFile"; 
$nameline = <CCF>;
unless ($nameline =~ /^theory/) {die "The `theory' file doesn't appear to be one";}
@codes = split(" ",$nameline);
$rely = 0; #  assume it is run time in the theory file
if (defined($codes[3]) && $codes[2] eq "r") { #reliability instead of run time
  $rely = 1;
  $conf = $codes[3];
}
unless ($theoryonly) {
  if ($rely) { #reliability
    print "normalized rms func error	$conf-% confidence\n";
    print "interval	func	reliability\n";
  } else { #run time
    print "normalized rms errors\n";
    print "interval	func	run\n";
  }
}
open(B,">subd.dat"); #for marking boundaries
open(C,">pred.dat"); #calculation output values
open(CR,">runpred.dat"); #calc runtimes
$WtdVerr = 0;
$WtdRerr = 0;
$km = 0;  #measurements counter, advances along with subdomains
$k = 0; #counter in theory file
#print "theory: $theoryonly; range: $start--$stop\n";  #debug
MLOOP:
while ($line = <CCF>) { # plot theory, print error table if measured.
  @subd = split(" ",$line);
  if ($parmsupplied && ($UB[$k] <= $start || $LB[$k] >= $stop)) { #outside range given
    $k++;
    next MLOOP;
  }
  print B "$subd[0] 0\n";
  $slope = $subd[2];
  $intcpt = $subd[3];
  $rslope = $subd[4];
  $rintcpt = $subd[5];
  if ($parmsupplied && $LB[$k] < $start) {
    $Lp = $start;
  } else {
    $Lp = $LB[$k];
  }
  if ($parmsupplied && $UB[$k] > $stop) {
    $Rp = $stop;
  } else {
    $Rp = $UB[$k];
  }
  $v1 = $slope*$Lp + $intcpt;
  $r1 = $rslope*$Lp + $rintcpt;
  $v2 = $slope*$Rp + $intcpt;
  $r2 = $rslope*$Rp + $rintcpt;
  print C "$Lp $v1\n";
  print C "$Rp $v2\n";
  print CR "$Lp $r1\n";
  print CR "$Rp $r2\n";
  unless ($theoryonly) {
    #$verr = 100*$subd[6]/$yrms;
    #$rerr = 100*$subd[7]/$trms;
    $verr = 0;
    $rerr = 0;
    $startkm = $km;
    while ($km < $mcount && $inputs[$km] < $UB[$k]) { #add errors for subdomain $k
      $x = $inputs[$km];
      $v = $slope*$x + $intcpt;
      $r = $rslope*$x + $rintcpt;
      $verr += ($v - $outputs[$km])**2;
      $rerr += ($r - $runtimes[$km])**2;
      $km++;
    }
    if ($km == $startkm) { #no data in subdomain
      $verr = 100; #show 100% error
      $rerr = 100;
    } else {
      $nk = $km - $startkm;
#print "subd: $k; $nk points; $v approx; $verr errs2\n";  #debug
      if ($yrms && $trms) {
        $verr = 100*sqrt($verr/$nk)/$yrms;
        $rerr = 100*sqrt($rerr/$nk)/$trms;
      }
    }
    $fform = "%5.2f";
    if ($rely) { #relability, correct "run time" error to reliability value
      $rerr = $rintcpt;
      $fform = "	%6.4f";
    }
  } #end measured
    
  unless ($theoryonly) {
    $WtdVerr += $verr*($UB[$k] - $LB[$k]);
    $WtdRerr += $rerr*($UB[$k] - $LB[$k]);
    printf "[%.2f, %.2f)	%5.2f%%	$fform%%\n",$LB[$k], $UB[$k], $verr, $rerr;
  }
  $k++;
} #end MLOOP
unless ($theoryonly) {
  $WtdVerr /= $stop - $start;
  $WtdRerr /= $stop - $start;
  printf "Weighted ave:	%5.2f%%	$fform%%\n", $WtdVerr, $WtdRerr;
}
close CCF;
close(C);
close(CR);
print B "$UB[$k-1] 0\n";
close(B);

  $meastring = "\"output.dat\" axes x1y1 with lines lt 4,";
  $runmeastring = "\"runtime.dat\" axes x1y2 with lines lt 5,";
  if ($rely) {
    $runmeastring = "";
  }
  if ($theoryonly) {
    $meastring = "";
    $runmeastring = "";
  }
  open(PF,">complot");  #same "complot" used by everyone!
  #print PF "set nokey\n";
  #print PF "set size 1.3,1.3\n";
  if ($parmsupplied) {
    print PF qq(set xrange [$start:$stop]\n);
  }
  print PF qq(set nox2tics\n);
  print PF qq(set border 11\n);
  print PF qq(set xlabel "Input showing subdomains" font "Helvetica,16"\n);
  print PF qq(set ylabel "Functional output values" font "Helvetica,16"\n);
  print PF qq(set xtics nomirror font "Helvetica,14"\n);
  print PF qq(set ytics nomirror font "Helvetica,14"\n);
  print PF qq(set y2tics font "Helvetica,14"\n);
  if ($rely) {
    print PF qq(set y2label "$conf%-confidence reliability bound" font "Helvetica,16"\n);
  } else {
    print PF qq(set y2label "Run time values" font "Helvitica,16"\n);
  }
  print PF "plot \"runpred.dat\" axes x1y2 with lines lt 3,\"subd.dat\" axes x1y1 with points pointsize 3,$meastring $runmeastring \"pred.dat\" axes x1y1 with lines lt -1\n";
  #print PF qq(plot "runpred.plt" axes x1y2 with lines lt 3,"subd.plt" axes x1y1 with points pointsize 3,$meastring $runmeastring "pred.plt" axes x1y1 with lines lt -1\n);
  #print PF "pause -1\n";
  close(PF);
  util::GNUplotit("complot");
  #system("gnuplot complot");
