#!/pkg/bin/perl -w
# hmisstest.pl
#July 31, 2005 modify to speed up the case where a row, col. pair has a 2 and a constant. 
#
# July 11, 2005 modify for large-scale execution.
# June21, 2005 modifying to solve the haplotyping problem so as to minimize the number of
# conflicts.
# June 4, 2005 changed use of stringlen, right and left. This should have no affect on
# the output when all columns are selected, but fixed a bug when only an interval is
#selected, and also the logic is clearer now.
#
#December 5, 2004
# use diagnostics
# Generates the ILP to solve the problem of
# setting the missing values in order to minimize the number of
# resulting incompatible pairs of sites.
#
open IN, "$ARGV[0]";
open OUT, ">constraints";
open OUT1, ">objective";
open OUT2, ">trace"; #for debugging
$line1 = <IN>;
$line2 = <IN>;
chomp $line1;
chomp $line2;
$lpoutfile = "hlp.lp";
@lines = <IN>;
$line = $lines[1];
#print "$line";
chomp $line;
$line =~ tr/ //d;
$stringlen = length ($line);
print "$stringlen ***\n";

$li = 0;
foreach $line (@lines) {
      $hli = 2*$li;
      $hlines[$hli] = $line;
print "$hli, $hlines[$hli]\n";
      $hli++;
      $hlines[$hli] = $line;
print "$hli, $hlines[$hli]\n";
$li++;
}


%count = ();
$left = 0;
$right = $stringlen-1; # right and left now define the range of the indices into the data. 

#print "do the whole range or an interval - input d for whole range or two numbers for an interval\n";
#$line = <STDIN>;
#if ($line =~ /d/) {
#$left = 0;
#$right = $stringlen-1; # right and left now define the range of the indices into the data. 
#}
#elsif ($line =~ /(\d+) (\d+)/) {
#  $left = $1-1;
#  $right = $2-1;
# print "got those numbers\n";
#}

   foreach $i ($left...$right -1) {
     foreach $j ($i+1...$right) {
     $ii = $i+1;
     $jj = $j+1;
     $key = $ii . ',' . $jj;
 #   print OUT "$key\n";

     %{$hash{$key}} = ();     
     $count{$key} = 0;
     }
   }


print OUT1 "minimize \n";
# print OUTM "maximize \n";
     
$linecount = 1;
foreach $line (@hlines) {
 $line =~ tr/ //d;
 chomp $line;

   @chars = split(//, $line);
   foreach $i ($left...$right -1) {    
     foreach $j ($i+1...$right) {
     $ii = $i+1;
     $jj = $j+1;
     $key = $ii . ',' . $jj;
     $value  = $chars[$i] . ',' . $chars[$j];
     if (($chars[$i] eq '2') or ($chars[$j] eq '2')) {
    #    print OUT "$i, $j \n";
       $variable{$key} = 1;
       $vhash{$linecount}{$key} = $value;
   #   print OUT "Another $key; $value, $count{$key}\n";
       if ($chars[$i] ne '2') {              # July 31, 2005 this is a kluge and a bass-ackwards way to do this,
                                             # but given the existing code, it is a quick hack and seems to
					     # work - anyway, the speed of this code is not critical. Rather
					     # it is the solution speed of the ILP produced that matters.
          $value0 = "$chars[$i],0";
          $value1 = "$chars[$i],1";
             if (! defined $hash{$key}{$value0}) {
             $hash{$key}{$value0} = 1;
             $count{$key}++;
             }

             if (! defined $hash{$key}{$value1}) {
             $hash{$key}{$value1} = 1;
             $count{$key}++;
             }
       }
       if ($chars[$j] ne '2') {
          $value0 = "0,$chars[$j]";
          $value1 = "1,$chars[$j]";
             if (! defined $hash{$key}{$value0}) {
             $hash{$key}{$value0} = 1;
             $count{$key}++;
             }

             if (! defined $hash{$key}{$value1}) {
             $hash{$key}{$value1} = 1;
             $count{$key}++;
             }
       }
     }
     elsif (! defined $hash{$key}{$value}) {
       $hash{$key}{$value} = 1;
       $count{$key}++;
   #   print OUT "Another $key; $value, $count{$key}\n";
     }
  } 
 }
  $linecount++;
}


 $count = 0;
 foreach $ekey (sort keys %count) {
  if ($count{$ekey} == 4) {
   print OUT2 "XX $ekey  $count{$ekey} \n";
   print OUT "LB$ekey = 1\n";           # this puts in an LB variable set to 1 for
                                        # every pair that is already in conflict. Remove this to
                                        # just compute the extra conflicts rather than the total conflicts.
       print OUT1 "+ LB$ekey ";
       if (0 == $count % 8) {
       print OUT1 "\n";
       }
  $count++;
  }
 }
    print OUT1 "\n";


$linecount = 0;
foreach $line (@hlines) {
$linecount++;

 $line =~ tr/ //d;
 chomp $line;
 @chars = split(//, $line);
 foreach $i ($left...$right) {
   $ii = $i+1;
     if ($chars[$i] eq '2') {
#        print OUT "Y$linecount,$ii <= 1\n";
        $bin{"Y$linecount,$ii"} = 1;
	if (0 == $linecount %2 ) {     # if this is the second line of a variable pair
	                               # then add in the constraint that exactly one
				       # of the variables must be set to 1 and the other to 0.
	$lcm1 = $linecount - 1;
	print OUT "Y$linecount,$ii + Y$lcm1,$ii = 1\n";
	}
     }
     elsif ($chars[$i] eq '1') {
         print OUT "Y$linecount,$ii = 1\n";
     }
     else {
         print OUT "Y$linecount,$ii = 0\n";
     }
 }
}


 $k = 0;
 $useful = 0;
 $paircount = $subcount = $sub22count = 0;
   foreach $i ($left...$right -1) {
     foreach $j ($i+1...$right) {
     $ii = $i+1;
     $jj = $j+1;
       $key = $ii . ',' . $jj;
      if ((defined $variable{$key}) and ($count{$key} < 4)) {           # we only further examine a pair of columns if they have some 
                                                                         # variable element(s) are not already incompatible
        $string11 =  $string10 = $string01 = $string00 = "";
        $c11 = $c10 = $c01 = $c00 = 0;
        $I11 = $I10 = $I01 = $I00 = 0;
        $stringforpair = "";

        foreach $lin (1..$linecount){
         if (defined $vhash{$lin}{$key}) {
         $paircount++; 
 #        print "$paircount\n";

             $value = $vhash{$lin}{$key};
             if ($vhash{$lin}{$key} eq '2,2') {   

#             print OUT "X$lin,$key,1,1 + X$lin,$key,1,0 +  X$lin,$key,0,1 + X$lin,$key,0,0 = 1\n";
#
#            print OUT "X$lin,$key,1,1 - Y$lin,$ii <= 0\n";  #if X is set to 1,1 then both Y variables
#            print OUT "X$lin,$key,1,1 - Y$lin,$jj <= 0\n";  # must be set to 1
#
#            print OUT "X$lin,$key,1,0 - Y$lin,$ii <= 0\n";  #if X is set to 1,0, then the first Y variable
#            print OUT "X$lin,$key,1,0 + Y$lin,$jj <= 1\n";  #is set to 1 and the second to 0
#
#            print OUT "X$lin,$key,0,1 + Y$lin,$ii <= 1\n";
#            print OUT "X$lin,$key,0,1 - Y$lin,$jj <= 0\n";
#
#            print OUT "X$lin,$key,0,0 + Y$lin,$ii <= 1\n";
#            print OUT "X$lin,$key,0,0 + Y$lin,$jj <= 1\n";
#
             $stringforpair .= "X$lin,$key,1,1 + X$lin,$key,1,0 +  X$lin,$key,0,1 + X$lin,$key,0,0 = 1\n";

             $stringforpair .= "X$lin,$key,1,1 - Y$lin,$ii <= 0\n";  #if X is set to 1,1 then both Y variables
             $stringforpair .= "X$lin,$key,1,1 - Y$lin,$jj <= 0\n";  # must be set to 1

             $stringforpair .= "X$lin,$key,1,0 - Y$lin,$ii <= 0\n";  #if X is set to 1,0, then the first Y variable
             $stringforpair .= "X$lin,$key,1,0 + Y$lin,$jj <= 1\n";  #is set to 1 and the second to 0

             $stringforpair .= "X$lin,$key,0,1 + Y$lin,$ii <= 1\n";  #this is the symmetric case when X is set to 0,1
             $stringforpair .= "X$lin,$key,0,1 - Y$lin,$jj <= 0\n";

             $stringforpair .= "X$lin,$key,0,0 + Y$lin,$ii <= 1\n";  # this is the case when  X is set to 0,0
             $stringforpair .= "X$lin,$key,0,0 + Y$lin,$jj <= 1\n";

             $bin{"X$lin,$key,1,1"} = $bin{"X$lin,$key,1,0"} =  $bin{"X$lin,$key,0,1"} = $bin{"X$lin,$key,0,0"} = 1;

#            print "2,2\n";

             if (!defined $hash{$key}{'1,1'}) {
              $I11 = 1;
              $c11++;                       # this records that 1,1 is missing in the pair of columns but can be 
                                              # created in this row.
              $string11 .= "+ X$lin,$key,1,1 ";    
                                                 # what we want to do here is to count the number of missing cases in
                                                 # this column pair, and then check that we even can create
                                                 # an incompatibility. If we can, then we need to create the
                                                 # inequalities to check if the appropriate X's were set to
                                                 # a new incompatibility. The string1,1 is accumulating the
                                                 # various X variables that would contribute to the 1,1 case,
                                                 # if they will later be used in the inequalities that check
                                                 # if the appropriate variables have been set in order to
                                                 # create a new incompatibility.
             }

              if (!defined $hash{$key}{'1,0'}) {
              $I10 = 1;
              $c10++;
              $string10 .= "+ X$lin,$key,1,0 ";

             }
              if (!defined $hash{$key}{'0,1'}) {
              $I01 = 1;
              $c01++;
              $string01 .= "+ X$lin,$key,0,1 ";
             }
              if (!defined $hash{$key}{'0,0'}) {
              $I00 = 1;
              $c00++;
              $string00 .= "+ X$lin,$key,0,0 ";
             }
             $sub22count++;
             $subcount++;
             }





   #     print OUT "$lin, $key $vhash{$lin}{$key}\n";  
   #     print OUT "Another $key; $value, $count{$key}\n";
           }
         }

   #     print OUT "$key, $count{$key}\n";
         if ($I11 + $I10 + $I01 + $I00 + $count{$key}  == 4) {  # here we check whether an incompatibility is possible for that column pair.
         $useful++;

     #   print OUT "Column pair $ii, $jj is useful\n";
         print OUT "$stringforpair\n";
         $rstring = "";

         if ($c11 > 0) {                          # need to create inequalities that set D11 = 1 iff one of the X11 variables are
                                                  # set to 1. That is, the 11 deficiency is rectified.

                        print OUT "$string11 - $c11 D$key,1,1 <= 0\n";  # This says that if any of the X11 variables are set to 1, then
                                                              # D11 will be set to 1.
                        print OUT "$string11 - D$key,1,1 => 0\n";  # This says that D11 is set to 1 only if at least one of the 
                                                          # X11 variables are set to 1.
                        $rstring .= "+ D$key,1,1 ";
                        $bin{"D$key,1,1"} = 1;

         }               
         if ($c10 > 0) {
                        print OUT "$string10 - $c10 D$key,1,0 <= 0\n";  
                        print OUT "$string10 - D$key,1,0 => 0\n";  
                        $rstring .= "+ D$key,1,0 ";
                        $bin{"D$key,1,0"} = 1;
         }               
         if ($c01 > 0) {
                        print OUT "$string01 - $c01 D$key,0,1 <= 0\n";  
                        print OUT "$string01 - D$key,0,1 => 0\n";  
                        $rstring .= "+ D$key,0,1 ";
                        $bin{"D$key,0,1"} = 1;
         }               
         if ($c00 > 0) {
                        print OUT "$string00 - $c00 D$key,0,0 <= 0\n";  
                        print OUT "$string00 - D$key,0,0 => 0\n";  
                        $rstring .= "+ D$key,0,0 ";
                        $bin{"D$key,0,0"} = 1;
         }               

                        $rhs = 3 - $count{$key};
                        print OUT "$rstring - 4 LB$key => -$count{$key}\n";  # these inequalities set LB to 1 if and only if
                        print OUT "$rstring - LB$key <= $rhs\n";    # the Y values are set so that the column pair is
                                                                                # incompatible.
                        $bin{"LB$key"} = 1;
                        print OUT1 "+ LB$key ";
                        $k++;
                        if (0 == $k % 8) {
                        print OUT1 "\n";
                        } 
                        
         }
#        else {
#            print OUT "Column pair $ii, $jj is not useful\n";
#        } 
      }
     }
  }

 print OUT1 "\n st \n";
# print OUTM "\n st \n";

print OUT "binaries\n";
 foreach $index (sort keys %bin) {
  print OUT "$index\n";
}

 print OUT "end\n";

 print "The number of non-perfect pairs is $count\n";
 print "The potential number of subpairs is $paircount\n";
 print "The subpotential number of subpairs is $subcount\n";
 print "The number of useful pairs is $useful\n";

close (IN);
close (OUT);
close (OUT1);

system ("cat objective constraints> $lpoutfile");
print "The output is in file $lpoutfile\n";