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

########################################

my %args;
getopt("ior", \%args);

my %isRNA = ('CYT',1,'URI',1,'GUA',1,'ADE',1);
my %isAtom;
my %isSUG;
my %atomVol;

my @times = split(' ', &ctime(time));
my $time = $times[1].'-'.$times[2].'-'.$times[4];

########################################

#Set input variables
my @in = ($args{i}) || <*.out>;
my $resi = $args{r} || "../../spec-defs.dat";   #designed for resi-defs.dat
my $out = $args{o} || "base-sizes-$time.dat";
open(OUT, "> $out")
  or die "Couldn't open $out for writing: $!\n";
print "Writing to $out...\n";

print OUT "\# NucProt volume set of base sizes\n\#\n\# Generated: ",
  &ctime(time),"\# Resi file $resi\n\#\n\n";

########################################

#Open Files
my @res = parseResiFile($resi);

print OUT "file\tADE\tCYT\tURI\tGUA\tSUG\tADE2\tCYT2\tURI2\tGUA2\n";
foreach my $inf (@in) {
  cleanAtomVol(@res);
  parseOutFile($inf);
  print STDERR "Calculating Statistics...\n";
  my %mean;
  foreach my $key (@res) {
    my ($r,$a) = split('_',$key);
    my ($m,$stdev) = meanstdev(@{$atomVol{"$key"}});
    $mean{$r}+=$m;
    $mean{$r.2}+=$m;
    if($r eq "SUG") {
      foreach my $r2 ("ADE","CYT","GUA","URI") {
        my $key2 = $r2."_".$a; 
        ($m,$stdev) = meanstdev(@{$atomVol{"$key2"}});
        $mean{$r2.2}+=$m;
      }
    }
  }
  my $inf2 = $inf;
  $inf2 =~ s/-edit.out$//; 
  $inf2 =~ s/-\d\d\d\d\d\d//; 
  $inf2 =~ s/-voronoi//;
  print OUT "$inf2\t",fixnum32($mean{"ADE"}),"\t",
    fixnum32($mean{"CYT"}),"\t",
    fixnum32($mean{"URI"}),"\t",
    fixnum32($mean{"GUA"}),"\t",
    fixnum32($mean{"SUG2"}),"\t",
    fixnum32($mean{"ADE2"}),"\t",
    fixnum32($mean{"CYT2"}),"\t",
    fixnum32($mean{"URI2"}),"\t",
    fixnum32($mean{"GUA2"}),"\n";
}
#end MAIN

########################################
########################################

sub parseResiFile {
#opens resi file
#returns list of residues and atoms in form of resi-atom
#allow fills global hash %isAtom
  my ($resifile) = @_;
  my @res;
  open(RESI, "< $resifile")
    or die "Couldn't open $resi for reading: $!\n";
  my $i=0; my $curresi;
  while(defined(my $line = <RESI>)) {
    if($line && substr($line,0,1) && !$i) {
      $line =~ s/\s+$//; $line =~ s/^\s+//;
      my @chunks = split(' ', $line);
      foreach my $test (@chunks) { $test =~ s/\s+$//; $test =~ s/^\s+//; }
      if($chunks[1]) {
        $curresi = $chunks[0];
        $i=$chunks[1];
      }
    } elsif($curresi && $line && substr($line,2,1)) {
      $line =~ s/^\s+//;
      my @chunks = split(' ', $line);
      my $test = $chunks[0];
      $test =~ s/\s+$//; 
      $test =~ s/^\s+//;
      push(@res,$curresi."_".$test);
      if($curresi eq "SUG") {
        $isAtom{"ADE_".$test} = 1; 
        $isAtom{"CYT_".$test} = 1; 
        $isAtom{"GUA_".$test} = 1; 
        $isAtom{"URI_".$test} = 1; 
        $isSUG{$test}=1; 
      }
      $isAtom{$curresi."_".$test} = 1;
      $i--;
    }
  }
  close(RESI);
  return @res;
};

########################################

sub parseOutFile {
  my ($in) = @_;
  my $total = getNumLines($in);
  print STDERR "\nreading the file ($total) $in\n";
  my $count; my $carrot=0;
  printBar(76);
  open(IN, "< $in") or die;
  while(defined(my $line = <IN>)) {
    $count++;
    my @chunks = split('\t',$line);
    foreach my $test (@chunks) { 
      $test =~ s/\s+$//; 
      $test =~ s/^\s+//; 
    }
    my $atom = $chunks[4]; my $resi = $chunks[3];
    my $status = $chunks[13];
    my $key = $chunks[3]."_".$chunks[4];
    if($status eq "\[OK\]" && $isAtom{$key}) {
      my $vol = $chunks[10];
      if($isSUG{$atom} && $isRNA{$resi}) {
        my $key2 = "SUG_".$chunks[4];
        push(@{$atomVol{$key2}},$vol);
      }  
      push(@{$atomVol{$key}},$vol);
    }
    if($count > $carrot) {
      print STDERR "^";
      $carrot = $carrot + $total / 76;
    }
  }
  print STDERR "\n";
  close(IN);
  return;
};

########################################

sub meanstdev {
  my (@nums) = @_;
  if($#nums < 1) { return 0; } 
  my ($sum,$sumsq) = (0,0);
  foreach my $num (@nums) {
    $sum += $num;
    $sumsq += ($num**2);
  }
  my $mean = ($sum/($#nums+1));
  my $stdev = sqrt(($sumsq-$sum**2/($#nums+1))/$#nums);
  return ($mean,$stdev);
};

########################################

sub cleanAtomVol {
#clean out %atomVol;
  my (@res) = @_;
  foreach my $key (@res) {
    my ($r,$a) = split('_',$key);
    @{$atomVol{$key}} = 0;
    if($r eq "SUG") {
      foreach my $r2 ("ADE","CYT","GUA","URI") {
        my $key2 = $r2."_".$a; 
         @{$atomVol{$key2}} = 0;
      }
    }
  }
};

########################################

sub getNumLines {
  my ($in) = @_;
  system("wc -l $in > temp.txt");
  open(IN, "< temp.txt") or die;
  my $num = <IN>;
  close(IN);
  my @chunks = split(' ',$num);
  $num = $chunks[0];
  $num =~ s/\s+$//; $num =~ s/^\s+//;
  return $num;
};

########################################

sub printBar {
  my ($num) = @_;
  my $mult=1;
  while($num/$mult > 80) {
    $mult++;
  }
  $num = int($num/$mult+0.5);
  print STDERR "|";
  for(my $i=1;$i<=$num-2;$i++) {
    if($i==int(($num-2)/2+0.501)) { print STDERR "\!"; }
    elsif($i==int(($num-2)/4+0.501)) { print STDERR "\#"; }
    elsif($i==int(3*($num-2)/4+0.501)) { print STDERR "\#"; }
    elsif($i%5==0) { print STDERR "+"; }
    else { print STDERR "-"; }
    if($i%100==0) { print STDERR "\n"; }
  }
  print STDERR "|\n";
  return $mult;
};

########################################

sub clean {
  my ($str) = @_;
  $str =~ s/\s+$//; $str =~ s/^\s+//;
  return $str;
};

########################################

sub list2vec {
  my (@list) = @_;
  my $vec;
  foreach my $item (@list) {
    $vec .= $item.",";
  }
  return $vec;
};

########################################

sub vec2list {
  my ($vec) = @_;
  return split(',',$vec);
};

########################################

sub vec2pdb {
#HETATM91176  O   HOH   496
  my ($vec,$count,$num,$resnum) = @_;
  my ($x,$y,$z,$r) = vec2list($vec);
  my $line = "HETATM";
  $line .= fixnum50($count);
  $line .= "  O   HOH  ";
  $line .= fixnum40($resnum);
  $line .= "    ";
  $line .= fixnum43($x+0.012);
  $line .= fixnum43($y+0.034);
  $line .= fixnum43($z+0.056);
  $line .= "  1.00";
  if($num) {
    $line .= fixnum32($num);
  } else {
    $line .= fixnum32(0.0);
  }
  $line .= "\n";
  return $line;
};