#! /usr/bin/perl -w
# XcompN  [-T] [-S] <.ccfc file> [<Other file>]
#  compares the result in a .ccfc calculated file with the <Other file>
#   which is presumed to have the same structure
#
#  Default is output behavior; -T for run time; -S for state.  
#   The default Other is output.datt; for -T and -S runtime.datt | state.datt
#     The line format is:
#     input N states [output | runtime | N states]
#
#  The values from Other are placed in the subdomains given by the .ccfc file
#   and averaged, then compared with the single .ccfc-file value

use samplings;

$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") {
    warn "Comparison not implemented for multiple states";
    $plotOpt = "S";
    $other = "state.datt";
    $highwater = $opt;
  }
}
die "A .ccf[tc] file must be given" unless (defined($ARGV[$highwater+1]));
$ccf = $ARGV[$highwater+1];
die "given file ",$ccf," not found" unless (-e $ccf);
if (defined($ARGV[$highwater+2])) { #other file given
  $other = $ARGV[$highwater+2];
  die "2nd file to compare doesn't exist" unless (-e $other);
} else { #default to "output.datt" etc., assume that SystemCode must create
  system "./compR -O -I SystemCode";
}
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];
}
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(D,">pred.datt");
$inS = 0;
for ($i=0; $i<$NumStates; $i++) {
  $stS[$i] = 0;
}
%result = ();
while ($line = <CCFC>) {
  @lane = split(" ",$line); # count out run s1 s2 ... sN
  if ($lane[0]) { #data count > 0
    $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;
    }
    $vec = sqrt($vec);
    print C " $vec";
    if ($plotOpt eq "O") {
      print C " $lane[1]\n"; #output
      print D " $lane[1]\n"; #output
      $result{$hash} = $lane[1];
    } elsif ($plotOpt eq "T") {
      print C " $lane[2]\n"; #run time
      print D " $lane[2]\n"; #run time
      $result{$hash} = $lane[2];
    } else {
      $vec = 0;
      $statelist = "";
      for ($i=0; $i<$NumStates; $i++) {
        print D " $lane[2+$i]";
        $statelist .= "$lane[2+$i] ";
        $vec += ($lane[2+$i])**2;
      }
      print C " $vec\n";  #state
      print D "\n"; 
      chop($result{$hash} = $statelist);  #remove final blank
    }
  }
  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);
close(D);
#debug
#foreach ($h,$r) (%result) {print "Subs: $h; Calc: $r\n"; }

################################################
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;
      }
    }
    die "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;
      }
    }
    die "State $val not found in state-$code list of subdomains";
  }
}
##############################################

%accum = ();
%howmany = ();
$mstot = 0;
$msN = 0;
$noneN = 0;
open (CF, "<$other"); #checked for existence above
while ($line = <CF>) { #
  @lane = split(" ",$line);
  #locate subdomain
  $I = find($lane[0],-1); 
  $hash = "$I";
  for ($s=0; $s<$NumStates; $s++) {
    if ($lane[1+$s] eq "u") { #state not initialized yet
      $S = find($initial[$s],$s);
    } else {
      $S = find($lane[1+$s],$s);
    }
    $hash .= " $S";
  }
  $mstot += $lane[$NumStates+1]**2;
  $msN++;
  $accum{$hash} += $lane[$NumStates+1];
  $howmany{$hash} += 1;
  #TBD:  this isn't right for multiple states; fix
}
close (CF);
$rms = sqrt($mstot/$msN);

$totalsubds = 0;
$averr = 0;
while ($h = each(%howmany)) {
  $ave = $accum{$h}/$howmany{$h};
  if (defined($result{$h})) {
    $R = $result{$h};
  } else {
    $R = "(none)";
    $noneN++;
  }
#debug
#print "Subs: $h; Calc: $R; Other: $ave\n";
  @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]);
    }
  }
  if ($R ne "(none)") {
    $err = 100*abs($R - $ave)/$rms;
    printf(": Err %3.2f%%\n", $err);
    $averr += $err;
    $totalsubds++;
  } else {
    printf(": Not Calculated\n");
  }
}
if ($noneN) { #some ccfc subd values were missing
  print "There were $noneN subdomains with actual but no calculated values\n";
}
$totalsubs = $NumInSubd;
for ($i=0; $i<$NumStates; $i++) {
  $totalsubs *= $NumStSubd[$i];
}
print "Execution fell in $totalsubds subdomains out of $totalsubs\n";
$averr /= $totalsubds;
printf("Average error %3.2f%%\n",$averr);
