#! /usr/bin/perl -w
# XcompL [-r] [-s] [<left> <right> <count>]
#
#  Compares and plots the measured and calculated behavior of
#   a stateless component-based system.
#
#  The XcompL version uses a SystemCode file that has a built-in
#   loop, driven by an input series from here
#
#  Assumes that the SYN script has been run in the current directory without
#   the -X option, so that "SystemCode" is available to execute and "system.scf"
#   has the name of the calculated file in its first line.
#   These are the functions compared.
#
#  If any of the non-switch parameters are present, all must be.  Otherwise the default
#    is the interval defined by the calculated file and a count of about 
#    three times the number of subdomains it has.
#    The -r option causes the sampling to be random
#    and default 3 times as frequent as the count described above.
#
#  The normal sampling is equispaced across the interval; with the -s option
#    the subdomains within the interval are instead sampled the number of
#    times in the system.scf file (which means that -s overrides <count>).
#
# -G causes the graph to be displayed on the screen.

use sampling;
use component;
use util;

$Goption = 1;
$useRand = 0;
$useSubd = 0;
$maybeParm = 0;
A:
for ($n=0; defined($ARGV[$n]); $n++) {
  if ($ARGV[$n] eq "-r") {
    $useRand = 1;
    next A;
  }
  if ($ARGV[$n] eq "-s") {
    $useSubd = 1;
    next A;
  }
  if ($ARGV[$n] eq "-G") { #shouldn't be there...
    $Goption = 1;
    next A;
  }
#can't do this, conflict with $left negative
  #if ($ARGV[$n] =~ /^-/) { # strange option, ignore
    #next A;
  #}
  #not one of options
  $maybeParm = 1;
  last A;
}
if ($maybeParm) { #possibly some parameters
  unless (defined($ARGV[$n+2])) { die "Parameters: <left> <right> <count> must all be given";}
  $left = $ARGV[0+$n];
  $right = $ARGV[1+$n];
  $count = $ARGV[2+$n];
}
#debug
#print "PARAM:  $left,$right,$count\n";

# Get the input range from the theory.ccf file
# and allow override from command line
open(ST,"system.scf") or die "must have a system.scf file";
chop($lastheory = <ST>);
if ($useSubd) { #get samples/subd from first line of system.scf (fix later?)
  $firstline = <ST>;
  @firstsub = split(" ",$firstline);
  $sampersub = $firstsub[2];
}
close(ST); #leave open if to be read later; see $subhits below
open(TCCF,$lastheory) or die "A calculated ", $lastheory, " file must be present";
chop($line = <TCCF>);
@codes = split(" ",$line);
if ($line !~ /^theory/) {die "doesn't seem to be a theory file";}
#debug
#print "$line \n";
$Rely = 0;
if (defined($codes[3]) && $codes[2] eq "r") { #reliability file
  $Rely = $codes[3]; #confidence %
}

open(B,">subd.plt"); #for marking boundaries
open(C,">pred.plt"); #calculation values
open(CR,">runpred.plt"); #calculated run times
$k = 0;
while ($line = <TCCF>) { #read and store boundaries, store & plot theory graph.
  @subd = split(" ",$line);
  $LB[$k] = $subd[0];
  $UB[$k] = $subd[1];
#print "$k   "; #debug
    print B "$subd[0] 0\n";
    $slope[$k] = $subd[2];
    $intcpt[$k] = $subd[3];
    $v1 = $slope[$k]*$LB[$k] + $intcpt[$k];
    $v2 = $slope[$k]*$UB[$k] + $intcpt[$k];
    $runslope[$k] = $subd[4];
    $runintcpt[$k] = $subd[5];
    $v3 = $runslope[$k]*$LB[$k] + $runintcpt[$k];
    $v4 = $runslope[$k]*$UB[$k] + $runintcpt[$k];
    print C "$LB[$k] $v1\n";
    print C "$UB[$k] $v2\n";
    print CR "$LB[$k] $v3\n";
    print CR "$UB[$k] $v4\n";
    $verr[$k] = 0;
    $rerr[$k] = 0;
    $diverge[$k] = "";
  $subhits[$k] = 0;
  $k++;
}
$numSubd = $k;
close(C);
close(CR);
print B "$UB[$numSubd-1] 0\n";
close(B);
unless (defined($right)) { #no parameters; past here cannot test $right to see if params
  $left = $LB[0];
  $right = $UB[$numSubd-1];
  $count = 3*$numSubd;
  if ($useRand) {  #increase default count
    $count *= 3;
  }
}
$rs = "";
if ($useRand) {
  $rs = " (random)";
}
$princt = $count;
if ($useSubd) {
  $rs = $rs." (by subdomain)";
  $princt = $numSubd*$sampersub;
}
printf "Comparing on [%.2f, %.2f) for %d points %s\n",$left,$right,$princt,$rs;

unless (-e "SystemCode") {die "No code `SystemCode' to measure";}
unlink "trace_path"; unlink "tmp_hash"; #for trace option
open(M,">meas.plt");
open(MR,">runmeas.plt");
if ($useRand) {
  $input = $left + rand($right-$left);
}
else {
  &sampling::init($left,$right,$count);
  $input = &sampling::next;
}

$vsqrs = 0;
$rsqrs = 0;
SUB:
for ($h=0; $h<$numSubd; $h++) { #for all subdomains 
  #must jump out immediately if not -s
  #set up samples/subdomain
  if ($useSubd) { #override the global sampling
    #print "sub: $h; $LB[$h] $UB[$h] $sampersub\n"; if ($LB[$h] >= $UB[$h]) { next SUB; } #debug
    &sampling::init($LB[$h],$UB[$h],$sampersub);
    $input = &sampling::next;
    $verr[$h] = 0;
    $rerr[$h] = 0;
    $count = $sampersub*$numSubd;
  }
$K = 0;  #subdomain index left off
$firstime = 1;
while ($input ne ".") {
  $thdiverge = 0;
  $codiverge = 0;  #assume a functional value for both
  #find the subdomain it's in 
  if ($useSubd) { #easy! 
    $k = $h;
  } else { #(should be binary search for random case!)
    if ($input < $LB[0] || $input > $UB[$numSubd-1]) { #out of range
      $k = $numSubd;
    } else {
      ITER: 
      for ($k=$K; $k < $numSubd; $k++) {
        if ($input >= $LB[$k] && $input < $UB[$k]) {
          unless ($useRand) {
            $K = $k;  #remember where found
          }
          last ITER;
        }
      }
    }
  }
  if ($k >= $numSubd || $runslope[$k] == 0 && $runintcpt[$k] < 0) { #theory undefined
    $thdiverge = 1;
    $diverge[$k] = "theory^ ";
  }
  if ($firstime) {
    ($meas, $rmeas, $pid, $cin, $cout, $cerr) = component::excomp($input, "perl SystemCode", "begin");
#print "out: $meas; pid: $pid; CIN: $cin   (from excomp)\n";  #debug
    $firstime = 0;
  } else {
    ($meas, $rmeas, $pid, $cin, $cout, $cerr) = component::excomp($input, "perl SystemCode", "cont",$pid, $cin, $cout, $cerr);
  }
#print "out: $meas; pid: $pid; CIN: $cin   (from excomp)\n";  #debug
  if ($meas =~ /LOOPING/) {
    die "System may be looping on input ", $input, "; stopped";
  }
  unless ($rmeas =~ /\d/ ) {#no digits, usually means error for code
#print "Code stderr:  $rmeas\n"; #debug
    $codiverge = 1;
    $meas = -1;
    $rmeas = -1; #phony values
    if ($k < $numSubd) { #a subdomain was found
      $diverge[$k] .= "code^ ";
    }
  } 
  print M "$input $meas \n";
  print MR "$input $rmeas \n";
  #compare with theory
  if ($k < $numSubd) {
    #debug
    #print "$k	";
    $v = $slope[$k]*$input + $intcpt[$k];
    $r = $runslope[$k]*$input + $runintcpt[$k];
    $vsqrs += $meas*$meas;
    $rsqrs += $rmeas*$rmeas;
    $subhits[$k]++;
    $verr[$k] += ($meas - $v)**2;
    #debug
    #if ($verr[$k] > .1) {print "$k [$LB[$k], $UB[$k]) error $input @ $slope[$k]-$intcpt[$k]\n";}
    if ($Rely) {
      $rerr[$k] = $r;
    } else {
      $rerr[$k] += ($rmeas - $r)**2;
    }
  }
  if ($useRand) {
    $input = $left + rand($right-$left);
  }
  else {
    $input = &sampling::next;
  }
} #end sampling loop
unless ($useSubd) {
  last SUB; #must exit after one "subdomain" iteration if regular, not -s
}
} #end loop over subdomains
($meas, $rmeas, $pid, $cin, $cout, $cerr) = component::excomp(".", "perl SystemCode", "end",$pid, $cin, $cout, $cerr); #close out looping code
close(M);
close(MR);
$vrms = sqrt($vsqrs/$count);
$rrms = sqrt($rsqrs/$count);

if ($Rely) { #reliability
  print "normalized rms error    $Rely-% confidence\n";
  print "{count}	interval     func            reliability\n";
} else { #run time
  print "normalized rms errors\n";
  print "{count}	interval	func	run\n";
}

$wverr = 0;
$wrerr = 0;
$totalen = 0;
if ($Rely) {
  $fform = "%6.4f";
} else {
  $fform = "%6.2f";
}
$maxverr = 0;
$maxrerr = 0;
for ($k=0; $k < $numSubd; $k++) {
  if ($diverge[$k] ne "") {
    printf "{%d}	[%.2f, %.2f)    %s \n", $subhits[$k], $LB[$k], $UB[$k], $diverge[$k];
  } else {
    if ($subhits[$k] > 0) {
      $verr[$k] = 100*sqrt($verr[$k]/$subhits[$k])/(abs($vrms)+.0000001);
      if ($maxverr < $verr[$k]) {
        $maxverr = $verr[$k];
      }
      unless ($Rely) { #don't fool with reliability value
        $rerr[$k] = 100*sqrt($rerr[$k]/$subhits[$k])/(abs($rrms)+.0000001);
        if ($maxrerr < $rerr[$k]) {
          $maxrerr = $rerr[$k];
        }
      }
      $thislen = $UB[$k] - $LB[$k];
      $wverr += $thislen*$verr[$k];
      $wrerr += $thislen*$rerr[$k];
      $totalen += $thislen;
      printf "{%d}	[%.2f, %.2f)	%6.2f%%	$fform%% \n", $subhits[$k],$LB[$k], $UB[$k], $verr[$k], $rerr[$k];
    }
  }
} #end loop over subdomains
printf "Weighted ave:		%6.2f%%	$fform%% \n", $wverr/$totalen, $wrerr/$totalen;
printf "  Maximum:		%6.2f%% $fform%% \n", $maxverr, $maxrerr;

if ($useRand) { #measured files out of order
  #gnuplot doesn't care?
  #rename "runmeas.plt", "/tmp/t";
  #`mv runmeas.plt /tmp/t`;
  #`sort -g /tmp/t >runmeas.plt`;
  #rename "meas.plt", "/tmp/t";
  #`mv meas.plt /tmp/t`;
  #`sort -g /tmp/t >meas.plt`;
}

open(PF,">complot");
  print PF "set terminal postscript color\n";
  print PF qq(set output "plot.ps"\n);
  print PF qq(plot "meas.plt", "runmeas.plt", "subd.plt" with points pt 2 pointsize 3, "pred.plt" with lines, "runpred.plt" with lines\n);
  close(PF);
system("gnuplot complot"); #plot measured run time even if reliability case

if ($Goption) {
  open(PF,">complot");
  #print PF "set terminal x11\n";
  #print PF "set nokey\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);
  print PF qq(set xrange [$left:$right]\n);
  if ($Rely) {
    print PF qq(set y2label "$Rely%-confidence reliability bound" font "Helvetica,16"\n);
    print PF qq(plot "runpred.plt" axes x1y2 with lines lt 3,"subd.plt" axes x1y1 with points pointsize 3, "pred.plt" axes x1y1 with lines lt -1, "meas.plt" axes x1y1 with lines lt 4\n);
  } else {
    print PF qq(set y2label "Run time values" font "Helvitica,16"\n);
    print PF qq(plot "runpred.plt" axes x1y2 with lines lt 3,"subd.plt" axes x1y1 with points pointsize 3, "pred.plt" axes x1y1 with lines lt -1, "meas.plt" axes x1y1 with lines lt 4, "runmeas.plt" axes x1y2 with lines lt 5\n);
  }
  #print PF "pause -1\n";
  close(PF);
  util::GNUplotit("complot");
}
