#!/usr/bin/perl -w

#calculation script CalcF for one construct
# floating-point intervals version
# Uses "xxx.ccft" files for component specs

use dec15;

$realeps = &dec15::geteps();
$eps = $realeps;  #larger than actual float granularity 

sub warnmess{ #print warning based on V option
  unless ($Voption && defined($_[0])) {return;}
  $mess = "";
  MLOOP: for ($wm=0;;$wm++) {
    unless (defined($_[$wm])) {last MLOOP;}
    $mess = $mess.$_[$wm];
  }
warn $mess, "\n";
}

#get options passed by SYNF
$Voption = 0;  #default is silent
$allOptions = "";
OPT: for ($i = 0;;$i++) { #over all options
  unless (defined($ARGV[$i])) {last OPT;}
  $allOptions = $allOptions." ".$ARGV[$i];
  if ($ARGV[$i] eq "-V") {
    $Voption = 1;
  }
}
$thecurrentflag = "L"; #to mark the output theory files
$Rely = 0; # global may be set in the process_ccf subroutine

#########################################
#   read configuration file system .sscf
#########################################
$i = 0;
open(SSCF, "system.sscf");
chop($type = <SSCF>);
#debug
#print "|$type|\n";
if($type eq "conditional" || $type eq "loop") {
  chop ($ccf0 = <SSCF>);
}
chop($ccf1 = <SSCF>);  #The TRUE branch for conditional
unless ($type eq "loop") {
  chop($ccf2 = <SSCF>);
}
close SSCF;

########################################################
#Subroutine definition:
#return the vector of subdomain values for parameter file
#value is an array of strings, the string at subscript K
#being the line containing the subdomain intervals.

sub process_ccf { #one parameter: file name to read
  my ($i, $filePre, $line, @r);
#debug
#print "ccf file: $_[0] \n";
  @filePre = split(/\./,$_[0]);
  $ProgFile = $filePre[0];  # use .ccf file name w/o .ccf
  #debug
  #print "Program file: $ProgFile \n";
  open(CCF,$_[0]) or die "no component control file ",$_[0];
  chop($line = <CCF>);
  if ($line !~ /^theory/) { #must get the file made by COMPF
    close(CCF);  #use the same name
    #the following isn't the best way to get the ccft name...
    open(CCF,$filePre[0].".ccft") or die "component ",$filePre[0]," has no measured .ccft file";
    $line = <CCF>;  #get the code line to parallel .ccft case above
  }
  @code = split(" ",$line);
  if (defined($code[3]) && $code[2] eq "r") { #reliability file
    $Rely = $code[3]; #set by any file being "r", value is the confidence % >0
  }
  $i = 0;
  while ($line = <CCF>) {
    $r[$i] = $line;
    $i++;
  }
  close(CCF);
  return @r; #and globals $ProgFile, $Rely
}
#end of subroutine
##########################################################

#the "rangex" arrays contain line strings for the subdomain intervals
#and must be split out to be used

#Test component is range0
if ($type eq "conditional" || $type eq "loop") {
  @range0 = process_ccf($ccf0);
  $range0cnt = $#range0+1;
  $p0 = $ProgFile;
}

#range1 is 1st in series, TRUE for conditional, body for loop
@range1 = process_ccf($ccf1);
$range1cnt = $#range1+1;
$p1 = $ProgFile;

#range2 is 2nd in series, FALSE for conditional
unless ($type eq "loop") {
  @range2 = process_ccf($ccf2);
  $range2cnt = $#range2+1;
  $p2 = $ProgFile;
}

if ($Rely) { #some file was reliability, do that
  $thecurrentflag = "LC r $Rely";
}

##########################################
# component 0 is the test, if a conditional or loop
# component 1 first in a series, "true" part for conditional, body for a loop
# then component 2 if a series or "false" part of conditional
##########################################

#Read and store values from component measurements
#TBD do something about this ugly code

if ($type eq "conditional" || $type eq "loop") {
  for ($i=0; $i < $range0cnt; $i++) {
    @line = split(" ",$range0[$i]);
    #debug
    #print "@line\n";
    $subDleft0[$i] = $line[0];
    $subDright0[$i] = $line[1];
    $aveval = $line[2]*($line[0] + $line[1])/2.0 + $line[3];
    if ($aveval > -0.1 && $aveval < 0.2) {
      $compromise = 0;
    } elsif ($aveval >= 0.2 && $aveval < 1.2) {
      $compromise = 1;
    } else { #way off, quit
      die "for component ",$p0," test value ",$aveval, "on subdomain [",$line[0],", ",$line[1],") out of line";
    }
    unless (abs($aveval - $compromise) < 0.1) {  #had to fudge pretty badly
      warn "test component ",$p0," does not have clear binary value on subdomain [",$line[0],", ",$line[1],") -- $aveval replaced by $compromise";
    }
  #$valSlope0[$i] = 0;  not used...
  $valIntcpt0[$i] = $compromise;
  $runSlope0[$i] = $line[4];
  $runIntcpt0[$i] = $line[5];
  }
}
for ($i=0; $i < $range1cnt; $i++) { 
  @line = split(" ",$range1[$i]);
  $subDleft1[$i] = $line[0];
  $subDright1[$i] = $line[1];
  $valSlope1[$i] = $line[2];
  $valIntcpt1[$i] = $line[3];
  $runSlope1[$i] = $line[4];
  $runIntcpt1[$i] = $line[5];
}

unless ($type eq "loop") {
  for ($i=0; $i < $range2cnt; $i++) { 
    @line = split(" ",$range2[$i]);
    $subDleft2[$i] = $line[0];
    $subDright2[$i] = $line[1];
       $valSlope2[$i] = $line[2];
       $valIntcpt2[$i] = $line[3];
       $runSlope2[$i] = $line[4];
       $runIntcpt2[$i] = $line[5];
  }
}

#################################################
# Subroutine definitions:
#
# take two intervals as arrays (left 0, right 1, left 2, right 3) 
#as input, and return # the intersection interval as a similar array.
# If the intersection is empty, return [1,0)
#

sub intersect { 
my(@int1,@int2,$t,@nul);
@nul = (1, 0);
@int1 = ($_[0], $_[1]);
@int2 = ($_[2], $_[3]);
#debug
#print "@int1\n";
#print "@int2\n";
# Make @int1 lie to the left
if ($int1[0] > $int2[0]) { 
  @int2 = ($_[0], $_[1]);
  @int1 = ($_[2], $_[3]);
}
if ($int1[0] < $int2[0]) { #strictly to the left
  if ($int1[1] <= $int2[0]) { #no intersection
    return @nul;
  }
  else {
    if ($int1[1] <= $int2[1]) { 
      return ($int2[0],$int1[1]);
    }
    else { #must be that 1 completely covers 2
      return ($int2[0],$int2[1]);
    }
  }
}
else { #left ends coincide
  if ($int1[1] <= $int2[1]) {$t = $int1[1];} else {$t = $int2[1];}
  return ($int1[0],$t);
}
}

##################################################
# test for empty interval
#   parameters are endpoints
sub notnull {
my($left, $right);
 $left = shift;
 $right = shift;
 #debug
 #print "empty? [$left $right)\n";
 if ($left >= $right) { return 0;  #empty
 }
 else { return 1;
 }
}
#####################################################

$itf = 6;
$theorytest = "once.ccf";
$theoryseries = "another.ccf";
$again = "again.ccf";
$theory = "theory.ccf";

if ($type eq "conditional") {
  $k = 0; #counter for new subdomains created in @range3 structure
  $jstartT = 0;
  $jstartF = 0; #left-off indices for 2nd component subdomains
  for ($i=0; $i<$range0cnt; $i++) { #intersect subdomains
    $testTrue = $valIntcpt0[$i];
    $foundany = 0;
    if ($testTrue) {  # test is true on this subdomain
#TBD clean up repetitive code...  (careful of ...1 and ...2)
      ITL: for ($j=$jstartT; $j<$range1cnt; $j++) { # find any intersections with the proper component
        @int = intersect(($subDleft0[$i],$subDright0[$i]),($subDleft1[$j],$subDright1[$j]));
#debug
#print "@int \n";
        if (notnull(@int)) { #add to list of new subdomains
#debug
#if ($int[1]-$int[0] < .00001) {
#print "small intersection [$int[0], $int[1])\n";
#}
          $foundany = 1;
          $range3[$k] = dec15::dump15s($int[0]).dec15::dump15s($int[1]);
          if ($Rely) { #reliability; just the intercept
            $thtimeSlope[$k] = 0;
            $thtimeIntcpt[$k] = $runIntcpt1[$j] * $runIntcpt0[$i];
          } else {
            $thtimeSlope[$k] = $runSlope1[$j] + $runSlope0[$i];
            $thtimeIntcpt[$k] = $runIntcpt1[$j] + $runIntcpt0[$i];
          }
          $thvalSlope[$k] = $valSlope1[$j];
          $thvalIntcpt[$k] = $valIntcpt1[$j];
          $range3[$k] = $range3[$k] .dec15::dump15s($thvalSlope[$k]).dec15::dump15s($thvalIntcpt[$k]).dec15::dump15s($thtimeSlope[$k]).dec15::dump15s($thtimeIntcpt[$k])."TRUE";
          $k++;
        }
        else { #null intersection
          if ($foundany) {
            $jstartT = $j - 1;  #next time start where one was found
            last ITL;  #can't be any more this time
          }
        }
      }
    }
    else { # test false on this subdomain
      IFL: for ($j=$jstartF; $j<$range2cnt; $j++) { # find any intersections with the proper component
        @int = intersect(($subDleft0[$i],$subDright0[$i]),($subDleft2[$j],$subDright2[$j]));
#debug
#print "@int \n";
        if (notnull(@int)) { #add to list of new subdomains
          $foundany = 1;
          $range3[$k] = dec15::dump15s($int[0]).dec15::dump15s($int[1]);
          if ($Rely) { #reliability; just the intercept
            $thtimeSlope[$k] = 0;
            $thtimeIntcpt[$k] = $runIntcpt2[$j] * $runIntcpt0[$i];
          } else {
            $thtimeSlope[$k] = $runSlope2[$j] + $runSlope0[$i];
            $thtimeIntcpt[$k] = $runIntcpt2[$j] + $runIntcpt0[$i];
          }
          $thvalSlope[$k] = $valSlope2[$j];
          $thvalIntcpt[$k] = $valIntcpt2[$j];
          $range3[$k] = $range3[$k] .dec15::dump15s($thvalSlope[$k]).dec15::dump15s($thvalIntcpt[$k]).dec15::dump15s($thtimeSlope[$k]).dec15::dump15s($thtimeIntcpt[$k])."FALSE";
          $k++;
        }
        else { #null intersection
          if ($foundany) {
            $jstartF = $j - 1;  #next time start where one was found
            last IFL;  #can't be any more this time
          }
        }  
      }
    }
  }
  $range3cnt = $k;
}
elsif ($type eq "series") {

#sub locateOutput { #compare serial & binary search
#$Serial = locateOutputS($_[0]);
#$Binary = locateOutputB($_[0]);
#if ($Serial != $Binary) {
  #warn "Searches disagreed: ", $Serial, " != ",$Binary;
#}
#return $Serial;
#}

#sub locateOutput{ #parameter is the point to locate in 2nd component domains
sub locateOutputS{ #parameter is the point to locate in 2nd component domains
 #return is the index in $range2 structure, or -1 if not found 
  my($find,$j);
  $find = $_[0];
  for ($j = 0; $j < $range2cnt; $j++) { #serial search for now
    @rec = split(" ", $range2[$j]);
    if ($find >= $rec[0] && $find < $rec[1]) { #found here
      return $j;
    }
  }
  return -1;
}

#sub locateOutputB{ #binary search version
sub locateOutput{ #binary search version
  my($needed,$L,$R,$j);
  $needed = $_[0];
  $L = 0;
  $R = $range2cnt-1;  #extremes of subscripts of intervals
  INNER: while (1) {  #must jump out
    #subdomains are in order and disjoint, but may not cover the domain
    $j = int(($L + $R)/2);  #careful, must look at two intervals at the end
                            # if two, then $j is $L to start
  #debug
  #$str = "binary search for ".$needed." in [".$L.",". $R ."], index: ".$j;
  #print "$str \n";
  #end debug
    if ($R - $L + 1 <= 2) { #must be in these subdomains if here
      #try first one
      #debug
      #printf "needed: %.18e in [%.18e,%.18e)\n", $needed, $subDleft2[$j],$subDright2[$j];
      if ($needed >= $subDleft2[$j] && $needed < $subDright2[$j]) { #found
        last INNER
      }
      if ($L != $R) { # there were two, try the other
        $j = $R;
      #debug
      #printf "needed: %.18e in [%.18e,%.18e)\n", $needed, $subDleft2[$j],$subDright2[$j];
        if ($needed >= $subDleft2[$j] && $needed < $subDright2[$j]) { #found
          last INNER
        }
      }
      #not in either one, so can't be found
      return -1;
    }
    else { # [L,R] more than two subdomains; which way to shift?
      if ($needed < $subDleft2[$j]) { #it's left of here
        $R = $j;
      }
      else { #right
        $L = $j;
      }
    }
  } #end while(1)
  return $j;
} #end binary search version

    @rec0 = split(" ",$range2[0]);
    @reck = split(" ",$range2[$range2cnt-1]);
    $k = 0;  #accumulate new subdomains in $range3[$k]
    SUBD: for($i=0; $i<$range1cnt; $i++) { #over all first subdomains
      $begin1 = $subDleft1[$i];
      $end1 = $subDright1[$i];
      #1st component subdomain is [$begin1, $end1]
      if ($end1 < $begin1) { # should not happen
        die "interval [,",$begin1,", ",$end1,") reversed?";
      }
      #debug
      #print "Subd: [$begin1, $end1)\n";
      $slope = $valSlope1[$i];
      $intcpt = $valIntcpt1[$i];
      if (abs($slope) <= $eps) { #special case for faster step-function approx
        $thvalSlope[$k] = 0.0;
        $j = locateOutput($intcpt);
        if ($j < 0 ) { #not found
          warn "Functional value on subdomain [",$subDleft1[$i],", ",$subDright1[$i],") undefined because output $needed from ",$p1, " falls in no ",$p2," subdomain";
          $thvalIntcpt[$k] = 0; #phony value
          $thtimeSlope[$k] = 0;
          $thtimeIntcpt[$k] = -1; #"undefined!"
          $TF = "SERIES";
        } else {
          $thvalIntcpt[$k] = $valSlope2[$j]*$intcpt + $valIntcpt2[$j];
          if ($Rely) { #reliability; just the intercept
            $thtimeSlope[$k] = 0;
            $thtimeIntcpt[$k] = $runIntcpt2[$j] * $runIntcpt1[$i];
          } else {
            $thtimeSlope[$k] = $runSlope1[$i] + $slope*$runSlope2[$j];
            $thtimeIntcpt[$k] = $runSlope2[$j]*$intcpt + $runIntcpt2[$j] + $runIntcpt1[$i];
          }
          @check = split(" ",$range2[$j]);
          if (defined($check[6])) {
            $TF = $check[6];
          }
          else {
            $TF = "SERIES";
          }
        }
        #debug
        #print "$i/$j:$thvalIntcpt[$k]($TF) ";
        $range3[$k] = dec15::dump15s($subDleft1[$i]).dec15::dump15s($subDright1[$i]). dec15::dump15s($thvalSlope[$k]).dec15::dump15s($thvalIntcpt[$k]).dec15::dump15s($thtimeSlope[$k]).dec15::dump15s($thtimeIntcpt[$k]).$TF;

        $k++; #bump count of new intervals
        next SUBD; #done with 1st-comp subdomain
      } # end special case of zero slope
      #get here only if $slope != 0
      $L = 0;  # count of intervals into which this subdomain is broken
      $reltiny = &dec15::epsplace($end1 - $begin1); #float error relative to interval
      $left = $begin1;
      $lastj = -99; #the 2nd comp subdomain found on the previous time
      $sev = 1; #number of epsilons to ignore for float fuzz
      #$flag = 10; #see debug below
      SUBSUB: while ($left < $end1) { #jump out when all intervals obtained
        $diverge = 0; #assume the current interval will have a functional value
        #Assumes contiguous subdomains in 2nd component
        $needed = $slope*$left + $intcpt; 
        #debug
        #print "in: $left, out: $needed\n";
        $j = locateOutput($needed); 
        $inj = $j;  # save for adjusting below
        if ($left > $begin1 && $j == $lastj) {#found same 2nd-comp interval, must move on
          if ($slope > 0) {
            $j++;
          } else { # should really fix locateOutput to know about negative slopes...
            $j--; 
            $inj = $j; #so don't report it...
          }
        }
        @check = split(" ",$range2[$j]);
        if ($slope < 0 && abs($needed - $check[0]) <= $eps) {
          $j--;
        }
        #uncomment to eliminate some tiny subdomains
        #if ($slope > 0 && abs($needed - $check[1]) <= $eps) {
        #  $j++;
        #}
        if (abs($j - $inj) > 0) { #did adjustment, OK?
          #debug and declare "flag" above
          #if ($flag-- > 0) { warn "adjust start from $inj to $j"; @check = split(" ",$range2[$inj]); warn "in (1)[",$begin1,", ",$end1,"), (2)[",$check[0],", ",$check[1],"): subdomain [",$left,", ? ); slope: ",$slope,"; value: ",$needed,"; count: ",$L; }
        }
          if ($j < 0 || $j >= $range2cnt) { #should also check if the proper endpoint in here
            $j = -1; #fake not found
          }
        if ($j < 0) { #not found
          if ($L) { #but some found already, so rest uncovered
            $right = $end1;
            warn "Functional value on subdomain [",$left,", ",$end1,") undefined because output $needed from ",$p1, " falls in no ",$p2," subdomain";
            $diverge = 1;
          } else {
            # not-found case if none found so far
            if ($slope > 0) { #reflect back the first interval start
              $right = ($rec0[0] - $intcpt)/$slope;
            } else { #slope negative, reflect back the last interval end
              $right = ($reck[1] - $intcpt)/$slope
            }
            if ($slope == 0 || $left > $end1 || $right <= $left) {
              $right = $end1;
            }
            #try again with next bit (if any), but this interval uncovered
            warn "Functional value on subdomain [",$begin1,", ",$right,") undefined because output $needed from ",$p1, " falls in no ",$p2," subdomain";
            $diverge = 1;
          }
            #continue below with "undefined" interval
        } else {# $left's output found at $j
          $jj = locateOutput($slope*$end1 +$intcpt);
          if ($jj == $j) { #same second-comp interval (covers slope 0 or almost 0 case)
            #whole interval retained
            $right = $end1;
          }
          #reflect back 2nd-component interval boundary to 1st domain
          elsif ($slope > 0) {
            $right = ($subDright2[$j] - $intcpt)/$slope; 
          }
          else { #slope < 0
            $right = ($subDleft2[$j] - $intcpt)/$slope;  #worry about open end?
          } #end slope cases  
          if ($right >= $end1) {
            $right = $end1;
          }
        } # end of the found case
        @check = split(" ",$range2[$j]); #get 2nd component values
        if ($right - $left <= $sev*$reltiny|| $right - $left <= $eps) {# tiny interval
          warnmess "in (1)[",$begin1,", ",$end1,"), (2)[",$check[0],", ",$check[1],"): tiny subdomain created [",$left,", ",$right,"); slope: ",$slope,"; value: ",$needed,"; count: ",$L;
            #$left = $right + $sev*&dec15::epsplace($right); #shift and omit this interval 
            #next SUBSUB; #fix when the theory file is written
        }
        if ($diverge) { #functional value undefined on this subdomain
          #mark by negative runtime 
          $thvalSlope[$k] = 0;
          $thvalIntcpt[$k] = 0; #phony values
          $thtimeSlope[$k] = 0;
          $thtimeIntcpt[$k] = -1;
        } else { #regular case
          if (defined($check[6])) {
            $TF = $check[6];
          }
          else {
            $TF = "SERIES";
          }
          $thvalSlope[$k] = $slope*$valSlope2[$j];
          $thvalIntcpt[$k] = $valSlope2[$j]*$intcpt + $valIntcpt2[$j];
          #if (abs($thvalSlope[$k]) > 1000 && $right - $left < .01*($reck[1] - $rec0[0])) { #change to constant!
          if ($end1-$begin1 < .001 && abs($thvalSlope[$k]) > 1000 && abs($thvalSlope[$k])*($right - $left)  > .3*($reck[1] - $rec0[0])) { #change to constant!
            warnmess "Made step-function interval [",$left,", ", $right,"), slope was ",$thvalSlope[$k];
            $thvalIntcpt[$k] = $thvalSlope[$k]*($right + $left)/2 + $thvalIntcpt[$k];
            $thvalSlope[$k] = 0;
          }
          if ($Rely) { #reliability; just the intercept
            $thtimeSlope[$k] = 0;
            $thtimeIntcpt[$k] = $runIntcpt2[$j] * $runIntcpt1[$i];
          } else {
            # run times: mx+b + m'x'+b', x' = output of C1 on input x
            $thtimeSlope[$k] = $runSlope1[$i] + $slope*$runSlope2[$j];
            $thtimeIntcpt[$k] = $runSlope2[$j]*$intcpt + $runIntcpt2[$j] + $runIntcpt1[$i];
          }
        } #end regular case
        $range3[$k] = dec15::dump15s($left).dec15::dump15s($right). dec15::dump15s($thvalSlope[$k]).dec15::dump15s($thvalIntcpt[$k]).dec15::dump15s($thtimeSlope[$k]).dec15::dump15s($thtimeIntcpt[$k]).$TF;

        #debug
        #print "$range3[$k] \n";
        $L++;
        $k++;
        if ($right >= $end1) { #test for done with covering this subdomain
          last SUBSUB;
        }
        $left = $right;
        $lastj = $j;
      } #end of SUBSUB loop
    }  # end for $i (first comp subdomain loop)
    $range3cnt = $k;
}
elsif ($type eq "loop") {
  #create subdirectory, do single calculation there, and bring back to this
  # directory the computed "theory" file 
  $thidx = 0;
  $subdir = ".SubCalcLoop";  #name of the subdirectory
  `rm -r -f $subdir`; #destroy directory so it can be used repeatedly
  mkdir $subdir;
  `cp * $subdir`;  #copy needed files
  chdir $subdir; #shift to new directory
#debug -- where are we?
#print `pwd`;
  #build a proper .sscf file
  open (SSCF, ">system.sscf") || die "can't create .sscf file for calculation";
  if ($p1 =~ /^theory/) {
    @tmp = split(/\./,$ccf1);
    $p1 = $tmp[0];  #no "program" name, use file name w/o .ccf
  }
  $p2 = "(ident)";
  $Cfalse = "ident.ccft";
  $Ctrue = $ccf1;
  $Ctest = $ccf0;
  warnmess "  Conditional:  IF ",$p0," THEN ",$p1," ELSE ",$p2," FI -> once  -> again";
  print SSCF "conditional\n";
  print SSCF "$Ctest\n";
  print SSCF "$Ctrue\n";
  print SSCF "$Cfalse\n";
  close SSCF;
#debug
#print `cat system.sscf`;
  if (-1 == system "./CalcF ".$allOptions) {die "`CalcF' script failed"} ;
  $cmd = "mv -f ".$theory." ../".$theorytest;
  `$cmd`;
  chdir ".."; #back to parent directory
  #debug 
  #system "cp $theorytest /tmp/looptest.ccf";

  # copy FALSE subdomains into "theory.ccf"
  #  and the TRUE subdomains into "again.ccf"
  open (THEO, $theorytest) || die "Can't open ", $theorytest;
  $tmp = <THEO>;  #skip first line (file name)
  open (THE, ">".$theory) || die "Can't open ",$theory;
  print THE "theory $thecurrentflag\n";
  open (AGAIN, ">".$again) || die "Can't open ",$again;
  print AGAIN "theory $thecurrentflag\n";
  $anytrue = 0;
  $totalsubs = 0;
  while ($line = <THEO>) {
    @isit = split(" ", $line);
    $totalsubs++;
      if ($isit[$itf] eq "TRUE")  {
        print AGAIN $line;
        $anytrue += 1;
      }
      else { #came out FALSE
        if ($isit[$itf] ne "FALSE") {
        `cat $theorytest`;
        die "Non-Aristotlian!";
        }
        print THE $line;
      }
  }
  close(THEO);
  close(THE);
  close(AGAIN);
  $anytrue1 = $anytrue; #to compare number of subdomains still true after one unwinding
  if ($anytrue) {
    warnmess "  (`again' is stripped of false subdomains -- ",$anytrue,"/",$totalsubs," active)";
    $againame = $again."0";
    system "cat $again > /tmp/$againame";
  }

FL:
  while ($anytrue) { #not all the subdomains came out FALSE
    #check for unending theory loop
    open(ONE,">again1");
    open(AGAIN,$again);
    $tmp = <AGAIN>; #discard name
    while ($line = <AGAIN>) { #copy line
      print ONE $line;
    }
    close(ONE);
    close(AGAIN);
  
    #create subdirectory, do calculation there, and bring back the result
    $subdir = ".SubCalcLoop";  #name of the subdirectory
    `rm -r -f $subdir`; #destroy directory so it can be used repeatedly
    mkdir $subdir;
    `cp * $subdir`;  #copy needed files
    chdir $subdir; #shift to new directory
  #debug -- where are we?
  #print `pwd`;
    #build a proper .sscf file
    open (SSCF, ">system.sscf") || die "can't create .sscf file for calculation";
    #$p1 = "again";
    #$p2 = "theorytest";
    print SSCF "series\n";
    print SSCF "$again \n";
    print SSCF "$theorytest \n";
    close SSCF;
  #debug
  #print `cat system.sscf`;
  #system "cat $again"; system "cat $theorytest"; die "stop.";
    if (-1 == system "./CalcF ".$allOptions) {die "`CalcF' script failed"} ;
    $cmd = "mv -f ".$theory." ../".$theoryseries;
    `$cmd`;
    chdir ".."; #back to parent directory
  
    #copy FALSE subdomains into "theory.ccf", TRUE into "again.ccf"
    #check if all subdomains have taken FALSE branch
  
    open (THEO, $theoryseries) || die "can't open $theoryseries";
    $tmp = <THEO>;  #skip first line (file name)
    open (THE, ">>".$theory) || die "can't append to $theory";
    open (AGAIN, ">".$again) || die "can't open $again";
    print AGAIN "theory $thecurrentflag\n";
    $anytrue = 0;
    while ($line = <THEO>) {
      @isit = split(" ", $line);
        if ($isit[$itf] eq "TRUE")  {
          print AGAIN $line;
          $anytrue += 1;
        }
        else { #came out FALSE
          if ($isit[$itf] ne "FALSE") {
            `cat $theoryseries>/tmp/foo`;
            die "Non-Aristotlian!";
          }
          print THE $line;
        }
    }
    close(THEO);
    close(THE);
    close(AGAIN);
    #check for non-termination
    if ($anytrue == $anytrue1) { #no subdomain eliminated
      # not sure if this algorithm is correct
        open(AGAIN,$again);
        $aga = <AGAIN>; #skip name
        open (ONE,"again1"); 
        CMPE: while (1) {
          $one = <ONE>;
          $aga = <AGAIN>;
          unless ($one) {last CMPE;}
          @part1 = split(" ",$one);
          @part2 = split(" ",$aga);
          if ($part1[0] != $part2[0] || $part1[1] != $part2[1]) { #subdomain differs
            #false alarm, keep going
            #debug
            #print "$one\n"; print "$aga \n";
            last CMPE;
          }
        }
        unless ($one) { #exited loop at end, same subdomains
          die "  Loop calculation will never terminate";
        }
        close(AGAIN);
        close(ONE);
    }
    $thidx++;
    #debug 
    #$theoryname = $theory.$thidx; system "cp $theory /tmp/$theoryname";

    #sort the theory.ccf file in place  -- surely there must be a better way?!
    open(CCF,$theory);
    open (T1,">t1");
    $tmp = <CCF>;
    print T1 $tmp;
    close(T1);
    open (T2,">t2");
    while ($tmp = <CCF>) {
      print T2 $tmp;
    }
    close(T2);
    close(CCF);
    `rm -f $theory`;
    `cp t1 $theory`;
    `sort -g t2 >> $theory`;
    if ($anytrue) {
      warnmess "  Series: again;once  -> again ",$anytrue," subdomains still active";
      $anytrue1 = $anytrue; #ready for next iteration
      system "rm -f again1";
      system "cp $again again1";
      #debug 
      #$againdex++; $againame = $again.$againdex; `cp $again /tmp/$againame`;
    }
    else {
      warnmess "  Series: again;once -> again, termination";
    }
  } #end loop over $anytrue
}
else { 
  die "Construct type ", $type, " not implemented";
}

####################################################
#  create "theory" .ccf file of results 
####################################################

unless ($type eq "loop") {  #loop file already written above
# write the "theory" equivalent file
  open(TCCF, ">".$theory) || die "Could not open `theory' file";
  print TCCF "theory $thecurrentflag\n";
  if ($type eq "conditional") {
    @t1 = @t2 = ();
    #must sort conditional subdomains
    @range3 = sort {(@t1=split(" ",$a))[0] <=> (@t2=split(" ",$b))[0]} @range3;
  } #else is "series," already sorted
  $R = -99; #phony value to start
  K3:
  for($k=0; $k<$range3cnt; $k++) { #check for contiguous as we go
    $lastR = $R; 
    @rline = split(" ",$range3[$k]);
    $L = $rline[0];
    $R = $rline[1]; #subdomain [L,R)
    if ($k > 0 && $lastR != $L && $p1.".ccf" ne $again) { #not contiguous, fix
      # but not if we are doing a loop, OK there
      warnmess "Gap ",$lastR,") [",$L," in ",$type," subdomains\n";
      $rline[0] = $lastR;
      $fixed = "";
      $kk = 0;
      while (defined($rline[$kk])) {
        $fixed = $fixed.$rline[$kk++]." "; 
      }
      $range3[$k] = $fixed;
    }
    @rline = split(" ",$range3[$k]);
    if ($rline[1] - $rline[0] <= $eps) {warn "Empty interval [",$rline[0],", ".$rline[1],") dropped\n"; next K3; }
    print TCCF "$range3[$k] \n";
  }
  close(TCCF);
}
