#!/usr/bin/perl -w 
use Getopt::Std;
require 'ctime.pl';
require "/home/vossman/lib/fixnum.pl";

print STDERR "\n";

$gap=0.01;
$grange=50;
 
getopt("ior", \%args);

$isresi{'CYT'}=1;
$isresi{'URI'}=1;
$isresi{'GUA'}=1;
$isresi{'ADE'}=1;
$isresi{'SUG'}=1;

@times = split(' ', &ctime(time));
$time = $times[2].$times[1].fixnum20b($times[4]%100);

#Set input variables
$in = $args{i} || "in.out";
$resi = $args{r} || "../../spec-defs.dat";   #designed for resi-defs.dat
if($args{o}) {
   $out = $args{o};
} else {
   $out = "./stdvol258-ext_$time.dat";
}
if(stat($out)) { 
  print STDERR "deleting $out...";
  sleep(5); 
  system("rm -f $out");
  print STDERR "done\n\n";
}


#print OUT "<P><FONT SIZE=-2>Page generated on ",$time," </FONT>\n";
$header = " atom\ttype\tcount\t\% OK\t mode\t mean\t sd\t% SD\t  min\t  max\n";
$header2 = "------\t-----\t-----\t----\t-----\t------\t----\t-----\t------\t------\n";

#Open Files






open(RESI, "< $resi")
   or die "Couldn't open $resi for reading: $!\n";
my %isSug;
print STDERR "reading $resi file\n";
my $curresi; my $i=0;
while(defined(my $line = <RESI>)) {
   if($line && $line =~ /^[A-Za-z]+/ && !$i) {
      $line =~ s/\s+$//; $line =~ s/^\s+//;
      @chunks = split(' ', $line);
      foreach $test (@chunks) { $test =~ s/\s+$//; $test =~ s/^\s+//; }
      if($chunks[1]) {
         push(@resi,$chunks[0]);
         print STDERR $chunks[0];
         $curresi = $chunks[0];
         $i=$chunks[1];
      }
   } elsif($curresi && $line && $line =~ /^  [A-Za-z]+/) {
      $line =~ s/\s+$//; $line =~ s/^\s+//;
      @chunks = split(' ', $line);
      foreach $test (@chunks) { $test =~ s/\s+$//; $test =~ s/^\s+//; }
      push(@{$curresi},$chunks[0]);
      my $label = $curresi."_".$chunks[0];
      $TYPE{$label}=$chunks[1];
      print STDERR ".";
      $i--;
      if($curresi eq "SUG") { $isSug{$chunks[0]}=1; }
   }
}
close(RESI);
print STDERR "\n\n";








print STDERR "parsing infile...";
my %lines;
open(IN, "< $in")
  or die "Couldn't open $in for reading: $!\n";
my $count=0;
while(defined(my $line = <IN>)) {
  $count++;
  if($count%5000 == 0) { print STDERR "."; }
  my @chunks = split('\t',$line);
  my $RES = $chunks[3];
  my $ATM = $chunks[4];
  $RES =~ s/\s+//g;
  $ATM =~ s/\s+//g;
  if($isSug{$ATM}) {
    push(@{"SUG__".$ATM},$line);
  } else {
    push(@{$RES."__".$ATM},$line);
  }
  #print STDERR $RES."__".$ATM;
}
close(IN);
print STDERR "done\n\n";








open(OUT, ">> $out")
      or die "Couldn't open $out for writing: $!\n";
print OUT 
"# NucProt type set of 258 Volumes\n#\n# Generated: ",&ctime(time)
,"# Input file $in\n# Resi file $resi\n#\n\n";
close(OUT);

print STDERR "Writing to $out...\n";

foreach $RESI (@resi) {
   open(OUT, ">> $out")
      or die "Couldn't open $out for writing: $!\n";
   $tvol = 0;
   $tvolsq = 0;
   $tokay = 0; $tposs = 0; $tbad = 0;
   $tmodevol=0; $tminvol=0; $tmaxvol=0;
   if($out ne "stdout") { print STDERR "\n$RESI\t[ "; }
   $k = -1;
   print OUT "$RESI\n",$header,$header2;
   foreach $ATOM (@{$RESI}) {
      $k++;
      if($out ne 'stdout') { print STDERR "$ATOM "; }
      $okay = 0; $poss=0; $bad=0;
      $ovol = 0;
      $ovolsq = 0;
      $maxvol = 0;
      $minvol = 1500;
      #open(IN, "< $in")
      #   or die "Couldn't open $in for reading: $!\n";
      #while( defined($line = <IN>) ) {
      foreach my $line (@{$RESI."__".$ATOM}) {
        @chunks = split('\t',$line);
        foreach $test (@chunks) { $test =~ s/\s+$//; $test =~ s/^\s+//; }
        if($ATOM eq $chunks[4] && ( ($RESI eq "SUG" && $isresi{$chunks[3]})
            || $RESI eq $chunks[3])) {
           #if(1) {
           #if($chunks[13] ne "bad") {
           if($chunks[13] eq "\[OK\]") {
              $okay++;
              $vol = $chunks[10];
              $ovol += $vol;
              $volt = int(int($vol/$gap)*$gap*10000+$gap/10)/10000;
              $mode{$volt}++;
              $ovolsq = $ovolsq + $vol*$vol;
              if($vol > $maxvol) { $maxvol = $vol; }
              if($vol < $minvol) { $minvol = $vol; }
           } 
           elsif($chunks[13] eq "possible") { $poss++; }
           elsif($chunks[13] eq "bad") { $bad++; }
        } 
     }
     #close(IN);
     $tokay += $okay;
     $tbad += $bad;
     $tposs += $poss;

     if($okay > 1) {
        $ovolsq = sqrt( ($ovolsq - $ovol*$ovol/$okay) / ($okay - 1) );
        $ovol = $ovol / $okay;
        $tvol += $ovol;
        $tvolsq += $ovolsq*$ovolsq;
     } elsif ($okay == 1) {
        $ovolsq = 0;
     }

     if($okay > 1) {
#FIND MODE
       $modemax=0;
       for($l=$gap*($grange+1);$l<45;$l+=$gap) {
          $avgmode = 0;
          $l=int($l*10000+$gap/100)/10000;
          for($k=-$grange;$k<=$grange;$k++) {
             $l2 = int(($l - $k*$gap)*1000 + $gap/100)/1000;
             if($mode{$l2}) {
                if(abs($l2 - $l) > 0.3) {
                   #$factor = int(2*$gap/abs($l-$l2)*1000+$gap/100)/1000
                   $factor = int(exp(-5*(abs($l-$l2)-0.1)**2)*1000+$gap/100)/1000
                } else { $factor = 1; }
                $avgmode += $mode{$l2}*$factor;
             }
          }
          if($avgmode > $modemax) {
             $modemax=$avgmode;
             $modevol = $l;
          }
          $l2=$l-$gap*($grange+1);
          $mode{$l2}=0;
        }
        $tmodevol+=$modevol;
        $tmaxvol+=$maxvol;
        $tminvol+=$minvol;

        my $label = $RESI."_".$ATOM;
        $okay2=$okay*100.0/($okay+$poss+$bad);
        print OUT  "  $ATOM\t",$TYPE{$label},"\t",fixnum50($okay),"\t",fixnum21($okay2),"\t",
                     fixnum22($modevol),"\t",fixnum23($ovol),"\t",fixnum12($ovolsq),"\t",
		     fixnum22($ovolsq/$ovol*100),"\t",
                     fixnum23($minvol),"\t",fixnum23($maxvol),"\n";
     } else {
        print STDERR "\be ";
     }
  }
  if($out ne 'stdout') { print STDERR "]"; }
  if($tokay > 1) {
     $tvolsq = sqrt($tvolsq);
     $tokay2=$tokay*100.0/($tokay+$tposs+$tbad);
     print OUT $header2;
     print OUT "TOTAL\t\t",fixnum50($tokay),"\t",fixnum21($tokay2),"\t",
                fixnum31($tmodevol),"\t",fixnum32($tvol),"\t",fixnum12($tvolsq),"\t\t",
                fixnum32($tminvol),"\t",fixnum32($tmaxvol),"\n\n";
  }
  close(OUT);
}

print STDERR "\n\nIf an \"e\" follows the character; no \[OK\] atoms were found.\n\n";