#! /usr/bin/perl -w
# Xcompt  [-T] [-S] [<.ccfc file> [<Otherfile>]]
#
#  Compares the result in a .ccfc calculated file (with state) with the <Otherfile>
#   which has the ".datt" structure of lines with:  
#     input state1 state2 ... stateN [output | runtime | N states]
#
#  Default (or -O) is output behavior; -T for run time; -S for state.  
#   The default Otherfile is "output.datt"; for -T | -S, runtime.datt | state.datt
#   Otherfiles are created by calling Xsys, unless
#   the existing Otherfile is newer than SystemCode.
#
#  The values from Otherfile are placed in the subdomains given by the .ccfc file
#   and averaged, then compared with the single .ccfc-file value
#
#  The comparison is graphed to show subdomains as floating rectangles
#    if a single state; else mag vector used for state.
#
#

use util;

$plotOpt = "O"; #default is to plot output
$highwater = -1; #last "-" arg
$other = "output.datt";
for ($opt=0; defined($ARGV[$opt]); $opt++) {
  if ($ARGV[$opt] eq "-T") {
    $plotOpt = "T";
    $other = "runtime.datt";
    $highwater = $opt;
  } elsif ($ARGV[$opt] eq "-S") {
    $plotOpt = "S";
    $other = "state.datt";
    $highwater = $opt;
  } elsif ($ARGV[$opt] eq "-O") {
    $plotOpt = "O";
    $other = "output.datt";
    $highwater = $opt;
  }
}
if (defined($ARGV[$highwater+1])) {
  $ccf = $ARGV[$highwater+1];
  unless (-e $ccf) {
    die "given file doesn't exist\n";
  } 
} else { #look in system.scf
  open (SCF, "<system.scf") || die "name of final theory file must be supplied\n";
  chop($ccf = <SCF>);
  close SCF;
}
if (defined($ARGV[$highwater+2])) { 
  $other = $ARGV[$highwater+2];
}
if (defined($ARGV[$highwater+2])) { #other file given
  die "2nd file to compare doesn't exist" unless (-e $other);
}
unless (-e $other && (-M $other < -M "SystemCode")) {
  warn "Running SystemCode...";
  system "perl Xsys -V";
}
open(CCFC, "<$ccf");
$line = <CCFC>;
@lane = split(" ",$line); # theory number-of-states number-input-subdomains numbers-state-subdomains state-initials
die "bad .ccfc file format in $ccf" unless ($lane[0] eq "theory");
$NumStates = $lane[1];
$NumInSubd = $lane[2];
for ($i=0; $i<$NumStates; $i++) {
  $NumStSubd[$i] = $lane[3+$i];
  #$initial[$i] = $lane[3+$NumStates+$i]; Never used?  OK, I think
}
for ($i=0; $i<$NumInSubd; $i++) {
  $line = <CCFC>;
  @lane = split(" ",$line);
  $LBin[$i] = $lane[0];
  $UBin[$i] = $lane[1];
}
for ($i=0; $i<$NumStates; $i++) {
  for ($j=0; $j<$NumStSubd[$i]; $j++) {
    $line = <CCFC>;
    @lane = split(" ",$line);
    $LBst[$i][$j] = $lane[0];
    $UBst[$i][$j] = $lane[1];
  } 
}
open(C,">pred.plt");
open(CP,">plat.plt") if ($NumStates == 1);
open(D,">pred.datt");
$inS = 0;
for ($i=0; $i<$NumStates; $i++) {
  $stS[$i] = 0;
}
%result = ();
$feas = 0;
while ($line = <CCFC>) {
  @lane = split(" ",$line); # count out run s1 s2 ... sN
  if ($lane[0]) { #data count > 0
    $feas++;
    $input = ($LBin[$inS]+$UBin[$inS])/2;
    $hash = "$inS";  # hash is a list of the subdomain subscripts, used with "$result"
    print C "$input";
    print D "$input";
    $vec = 0;
    for ($i=0; $i<$NumStates; $i++) {
      $j = $stS[$i];
#debug
#print "[$inS,$j): $line\n";
      $hash .= " $j";
      $stave = ($LBst[$i][$j] + $UBst[$i][$j])/2;
      print D " $stave";
      $vec += $stave**2;
    }
    if ($NumStates == 1) {
      $vec = $stave;
    } else {
      $vec = sqrt($vec);
    }
    print C " $vec";
    if ($plotOpt eq "O") {
      $val = $lane[1];
      print D " $lane[1]\n"; #output
      $result{$hash} = $lane[1];
    } elsif ($plotOpt eq "T") {
      $val = $lane[2];
      print D " $lane[2]\n"; #run time
      $result{$hash} = $lane[2];
    } else { #state
      $vec = 0;
      $statelist = "";
      for ($i=0; $i<$NumStates; $i++) {
        print D " $lane[3+$i]";
        $statelist .= "$lane[3+$i] ";
        $vec += ($lane[3+$i])**2;
      }
      if ($NumStates == 1) {
        $val = $lane[3];
      } else {
        $val = sqrt($vec);
      }
      print D "\n"; 
      chop($result{$hash} = $statelist);  #remove final blank
    }
    print C " $val\n";
    if ($NumStates == 1) { #plot plateaux
      print CP "$LBin[$inS] $LBst[0][$stS[0]] $val\n";
      print CP "$LBin[$inS] $UBst[0][$stS[0]] $val\n";
      print CP "\n";
      print CP "$UBin[$inS] $LBst[0][$stS[0]] $val\n";
      print CP "$UBin[$inS] $UBst[0][$stS[0]] $val\n";
      print CP "\n\n";
    }
  }
  if (++$inS >= $NumInSubd) {
    $inS = 0;
    BS:
    for ($s=0; $s<$NumStates; $s++ ) {
      if (++$stS[$s] >= $NumStSubd[$s]) {
        $stS[$s] = 0;
      } else { #leave rest alone
        last BS;
      }
    }
  }
}
close(C);
if ($NumStates == 1) {
  close(CP);
  $plts[scalar(@plts)] = "\"$other\" with points";
  $plts[scalar(@plts)] = "\"plat.plt\" with lines";
} else {
  $plts[scalar(@plts)] = "\"meas.plt\" with points";
  $plts[scalar(@plts)] = "\"pred.plt\" with points";
}
close(D);

#for ($h,$v) (%result) { print "Subs $h; Cal: $v\n"; }  #debug

################################################
sub find { #locate and return subdomain index for value
  my ($val, $code, $i);
  $val = $_[0];
  $code = $_[1]; # -1 means input
                 # 0..NumStates-1 means that state
  if ($code == -1) { #find input subdomain
    for ($i=0; $i<$NumInSubd; $i++) {
      if ($val >= $LBin[$i] && $val < $UBin[$i]) {
        return $i;
      }
    }
    warn "Input $val not found in subdomains";
  } else { #find state
    for ($i=0; $i<$NumStSubd[$code]; $i++) {
      if ($val >= $LBst[$code][$i] && $val < $UBst[$code][$i]) {
        return $i;
      }
    }
    warn "State $val not found in state-$code list of subdomains";
  }
}
##############################################

%accum = ();
%howmany = ();
%noneN = ();
for ($s=0; $s<$NumStates; $s++) {
  $mstot[$s] = 0;
}
$msN = 0;
open (CF, "<$other"); #checked for existence above
if ($NumStates > 1) { #need a 'vector' file for measurements
  open (DF, ">meas.plt") or die "can't open 'meas.plt'";
}
while ($line = <CF>) { #
  @lane = split(" ",$line);
  #locate subdomain
  $I = find($lane[0],-1); 
  $hash = "$I";
  $insvec = 0;  #state vector (more than 1 state)
  $outsvec = 0;
  for ($s=0; $s<$NumStates; $s++) {
    $insvec += $lane[1+$s]**2;
    if ($plotOpt eq "S") {
      $outsvec += $lane[$NumStates+1+$s]**2;
    } else {
      $outsvec = $lane[$NumStates+1]**2;
    }
    $S = find($lane[1+$s],$s);
    $hash .= " $S";
    if ($plotOpt eq "S") {
    $mstot[$s] += $lane[$NumStates+1+$s]**2;
    } else {
      $mstot[0] += $lane[$NumStates+1]**2;
    }
  }
  $insvec = sqrt($insvec);
  $outsvec = sqrt($outsvec);
  $msN++;
  if ($NumStates > 1) {
    print DF "$lane[0] $insvec $outsvec\n";
  }
  if (defined($result{$hash})) {
    if ($plotOpt eq "S") { #result is a list if more than one state
      @R = split(" ",$result{$hash});
      for ($s=0; $s<$NumStates; $s++) {
        $accum{$hash}[$s] += ($lane[$NumStates+1+$s] - $R[$s])**2;
      }
    } else {
      $accum{$hash}[0] += ($lane[$NumStates+1] - $result{$hash})**2;
    }
    $howmany{$hash} += 1;
  } else {
    $R = "(none)";
    $noneN{$hash} = 1;
  }
}
close (CF);
if ($NumStates > 1) {
  close DF;
}
if ($plotOpt eq "S") {
  $upb = $NumStates;
} else {
  $upb = 1;
}
for ($s=0; $s<$upb; $s++) {
  $rms[$s] = sqrt($mstot[$s]/$msN);
}

###################################################
sub inproc {
@spa = split(" ",$a);
@spb = split(" ",$b);
for ($n=scalar(@spa)-1; $n>=0; $n--) {
  if ($spa[$n] < $spb[$n]) {return -1;}
  if ($spa[$n] > $spb[$n]) {return 1;}
}
return 0;
}
######################################################

$which = "Output" if ($plotOpt eq "O");
$which = "Run time" if ($plotOpt eq "T");
$which = "Result state" if ($plotOpt eq "S");
print "Relative % r-m-s errors in $which \n";

$totalsubds = 0;
if ($NumStates == 1) { #print errs in table format
print "RMS to normalize: $rms[0]\n";  #debug
  $tabval = ();
  print "	";
  foreach $b (@LBin) {
    printf "[%2.1f	", $b; #abbreviate the [ , ) to [
  }
  print "\n";
  #retrieve all the hash values
  for $h (sort inproc keys %result) {
    if ($howmany{$h}) { 
      $av = sqrt($accum{$h}[0]/$howmany{$h});
      @hs = split(" ",$h);
#print "in: $hs[0]; st: $hs[1];  $av/$howmany{$h}\n";  #debug
      $tabval[$hs[0]][$hs[1]] = $av/$rms[0]; # %err; undefined if none
    }
  }
  $totArea = 0;
  $totErr = 0;
  $maxErr = 0;
  for ($j=0; $j<$NumStSubd[0]; $j++) { #loop over states
    printf "[%2.1f	", $LBst[0][$j]; #abbreviate the [ , ) to [
    for ($i=0; $i<$NumInSubd; $i++) { #loop over inputs
      $thisArea = ($UBin[$i]-$LBin[$i])*($UBst[0][$j]-$LBst[0][$j]);
      if (defined($tabval[$i][$j])) {
        $totArea += $thisArea;
        $thisErr = $tabval[$i][$j];
        if ($thisErr > $maxErr) {
          $maxErr = $thisErr;
        }
        printf "%3.1f	", 100.0*$thisErr;
        $totErr += $thisErr*$thisArea;
        $totalsubds++;
      } else {
        print " **	";
      }
    }
  print "\n";
  }
  if ($totArea > 0) {
    $totErr *= 100.0/$totArea;
    printf "Weighted r-m-s (max) error: %4.2f%% (%4.2f%%)\n",$totErr,100*$maxErr;
  }
} else { #explicit x-prod states
  for ($s=0; $s<$NumStates; $s++) {
    $averr[$s] = 0;
    $maxerr[$s] = 0;
  }
  for $h (sort inproc keys %howmany) {
    if ($plotOpt eq "S") {
      $upb = $NumStates;
    } else {
      $upb = 1;
    }
    for ($s=0; $s<$upb; $s++) {
      $ave[$s] = sqrt($accum{$h}[$s]/$howmany{$h});
    }
    @hs = split(" ",$h); # input  state(s)
    printf("[%4.2f,%4.2f)",$LBin[$hs[0]],$UBin[$hs[0]]);
    for ($i=0; $i<$NumStates; $i++) {
      $s = $hs[1+$i];
      if ($s eq "u") { #should never happen...
        printf(" X (undefined)");
      } else {
        printf(" X [%4.2f,%4.2f)",$LBst[$i][$s],$UBst[$i][$s]);
      }
    }
    printf(": Err ");
    if ($plotOpt eq "S") {
      $upb = $NumStates;
    } else {
      $upb = 1;
    }
    for ($s=0; $s<$upb; $s++) {
      $err = 100*$ave[$s]/$rms[$s];
      printf(" %3.2f%%", $err);
      $averr[$s] += $err;
      if ($err > $maxerr[$s]) {
        $maxerr[$s] = $err;
      }
    }
    printf("\n");
    $totalsubds++;
  }
  if ($plotOpt eq "S") {
    $upb = $NumStates;
    $vecer = 0;
    $mvecer = 0;
  } else {
    $upb = 1;
  }
  printf("Average (Max) error ");
  for ($s=0; $s<$upb; $s++) {
    $averr[$s] /= $totalsubds;
    printf(" %3.2f%% (%3.2f%%)",$averr[$s],$maxerr[$s]);
    $vecer += $averr[$s]**2;
    $mvecer += $maxerr[$s]**2;
  }
  if ($plotOpt eq "S" && $NumStates>1) {
    $vecer = sqrt($vecer);
    $mvecer = sqrt($mvecer);
    printf(" vec: %3.2f%% (%3.2f%%)",$vecer,$mvecer);
  }
  printf("\n");
} #end of multiple-state case
$tnone = 0;
for $h (keys %noneN) {
  if ($NumStates == 1) { #break out only if just one state
    if ($tnone == 0) {
      print "Subdomains executed, but no calculated value:\n";
    }
    @two = split(' ',$h);
    $bisub = $two[0];
    $bssub = $two[1];
    $iL = $LBin[$bisub];
    $iR = $UBin[$bisub];
    $sL = $LBst[0][$bssub];
    $sR = $UBst[0][$bssub];
    printf(" [%4.2f,%4.2f)x[%4.2f,%4.2f)\n",$iL,$iR,$sL,$sR);
  } else {
    print " (indicies) $h\n";
  }
  $tnone++;
#debug
#print "key: $h\n";
}
if ($tnone) { #some ccfc subd values were missing
  print "$tnone subdomains with actual but no calculated values\n";
}
$totalsubs = $NumInSubd;
for ($i=0; $i<$NumStates; $i++) {
  $totalsubs *= $NumStSubd[$i];
}
$infeas = $totalsubs - $feas;
print "Execution fell in $totalsubds subdomains out of $feas\n";
if ($infeas) {
  print "$infeas subdomains were calculated to be infeasible out of $totalsubs\n";
}

$to_plt = join(",",@plts);
open(PF,">complot");
#print PF "set terminal X11\n";
print PF qq(set xlabel "Input Space"\n);
print PF qq(set ylabel "State Space"\n);
print PF qq(set zlabel "Outputs"\n) if ($plotOpt eq "O");
print PF qq(set zlabel "Run Time"\n) if ($plotOpt eq "T");
print PF qq(set zlabel "Result State"\n) if ($plotOpt eq "S");
#print PF qq(set nokey\n);
if ($NumStates == 1) {
  print PF qq(set hidden3d\n);
}
print PF qq(set xtics nomirror\n);
print PF qq(set ytics nomirror\n);
print PF qq(set ticslevel .05\n);
print PF qq(splot $to_plt \n);
#print PF qq(pause -1 \n);
close(PF);  

util::GNUplotit("complot");
#system("gnuplot complot");
