#!/usr/bin/perl
use Getopt::Std;
require 'ctime.pl';
require './fixnum.pl';
use strict;

my $version = 3.0;

print STDERR "\n\n\n#####\n#####\n";
my $bondrad = 0.65;
my ($type1,$type2,$atomdef,$residef,$stdvol,$outfile,
    $min,$max,$incr,$debug,$spec,$listfile,$type2) = getINPUT();
writeTime($outfile);

if($debug) { print STDOUT "Start $min\nStop $max\nIncr $incr\nSpec $spec\nlist $listfile\n"; }
#foreach my $phos (1.54,1.64,1.74) {
#foreach my $phos2 (1.54,1.64,1.74) {
for(my $phos=1.40;$phos<=1.60;$phos+=0.01) {
for(my $phos2=1.52;$phos2<=1.72;$phos2+=0.01) {
   my $start = getTime();
   if($debug) { print STDERR "Creating Atom Defs file\n"; }
   my $types = editAtomDefs($atomdef,$phos,$phos2);

   if($debug) { print STDERR "Getting PDB IDs\n"; }
   my @pdbs = getPdbIDs();

   my $outf = "PHOS-".fixnum12($phos)."-temp-$spec\.out";
   if($debug) { print STDERR "Getting Volumes into $outf\n"; }
   foreach my $pdbid (@pdbs) {
     my $temp = "temp-$spec\.out";
     if($debug) { print STDERR "  get volume for ",lc($pdbid),"\n"; }
     if($debug) { print STDERR "perl getVolumes.pl -a $types -r $pdbid -o $temp\n"; }
     system("perl getVolumes.pl -a $types -r $pdbid -o $temp -u $spec");
     if(!stat("$temp")) { die "no $temp file\n"; }
     system("cat $temp >> $outf");
     system("rm -f $temp");
   }

   if($debug) { print STDERR "Starting analysis\n"; }
   my ($totdev,$totperdev,$rnaatoms,$numatoms,$count,$totatoms) 
     = AnalyzeVol($outf);

   if($debug) { print STDERR "Outputing Info\n"; }
   
   printOutInfo($outfile,$phos,$start,$totdev,$totperdev,
                $rnaatoms,$numatoms,$count,$totatoms,$phos2);
}}


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

sub getINPUT {
  my %args;
  getopt("spotnxidulv",\%args);
  if($args{v}) { print STDERR "\n\nOptimize, version $version\n\n"; die; }
  my $min = $args{n} || 1.2;
  my $max = $args{x} || 2.4;
  my $incr = $args{i} || 0.05;
  my $debug = $args{d} || 1;
  my $spec = $args{u} || int(rand(900000)+100000);
  my $listfile = $args{l} || "pdbs.txt";
  #my $type = $args{t} || "P4H0u";
  my $type1 = "O2H0s";
  my $type2 = "O2H0b";
  my $path = getPath();
  my $outfile = $args{o} || $path."type-vs-stdev.fst.".$spec;
  my $atomdef = $path."atom-defs.dat";
  if(!stat($atomdef)) { die; }
  my $residef = $path."resi-defs.dat";
  if(!stat($residef)) { die; }
  my $stdvol = $path."stdvol173.dat";
  if(!stat($stdvol)) { die; }
  return ($type1,$type2,$atomdef,$residef,$stdvol,$outfile,$min,$max,
	$incr,$debug,$spec,$listfile,$type2);
};


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

sub isRNA {
  my ($resi) = @_;
  $resi = uc($resi);
  $resi =~ s/^\s+//; $resi =~ s/\s+$//;
  my %hash = ("ADE",1,"GUA",1,"URI",1,"CYT",1);
  if($resi && $hash{$resi}) { return 1; }
  return 0;
};


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

sub getPath {
  use Cwd;
  my $dir = cwd();
  $dir =~ s/$/\//;
  my $path = "/home/vossman/phos/";
  if(!stat($path)) { die; }
  if($path eq $dir) { $path = ""; }
  return $path;
};


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

sub getTime {
  my @times = split(' ', &ctime(time));
  my $time = $times[4].' '.$times[1].' '.$times[2]." ".$times[3];
  return $time;
};


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

sub subtTime {
  my ($time1,$time2) = @_;
#time2 > time1
  $time1 =~ s/:/ /g; $time2 =~ s/:/ /g;
  my @times1 = split(' ',$time1);
  my @times2 = split(' ',$time2);
  $times1[1] = mon2num($times1[1]);
  $times2[1] = mon2num($times2[1]);
  my ($year,$month,$day,$hour,$min,$sec);
  #sec
    $sec = $times2[5] - $times1[5];
    if($sec < 0) { $times2[4]--; $sec+=60; }
  #min
    $min = $times2[4] - $times1[4];
    if($min < 0) { $times2[3]--; $min+=60; }
  #hour
    $hour = $times2[3] - $times1[3];
    if($hour < 0) { $times2[2]--; $hour+=24; }
  #days - grossly inaccurate
    $day = $times2[2] - $times1[2];
    if($day < 0) { $times2[1]--; $day+=30; }
  #months
    $month = $times2[1] - $times1[1];
    if($month < 0) { $times2[0]--; $month+=12; }
  #year
    $month = $times2[0] - $times1[0];
  if($debug > 1) { print STDERR "$year,$month,$day,$hour,$min,$sec"; }
  if($year) {
    return "$year year(s), $month month(s)";
  } elsif($month) {
    return "$month month(s), $day day(s)";
  } elsif($day) {
    return "$day day(s), $hour hour";
  } elsif($hour) {
    return "$hour hour, $min min";
  } elsif($min) {
    return "$min min, $sec sec";
  } elsif($sec) {
    return "$sec sec";
  }
  return 0;
};


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

sub mon2num {
  my ($mon) = @_;
  my @mons = ("Jan","Feb","Mar","Apr","May","Jun",
              "Jul","Aug","Sep","Oct","Nov","Dec",);
  for(my $i=0; $i<=$#mons; $i++) {
    if($mon eq $mons[$i]) {
      return $i+1;
    }
  }
  return 0;
};


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

sub getPdbIDs { 
  my @pdbs;
  if($debug) { print STDERR "Reading $listfile for PDBs\n"; }
  open(PDB, "< $listfile") or die;
  while(defined(my $line = <PDB>)) {
    $line =~ s/^\s+//; $line =~ s/\s+$//;
    if($line && length($line) == 4) {
      push(@pdbs,$line);
    }
  }
  if($#pdbs <= 0) { die "not enough pdbid in pdbs.txt"; }
  return @pdbs;
};


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

sub AnalyzeVol {
   my ($outf) = @_;
#PASS ONE -- get central limits 5% and 95%
   my %min; my %max; my $totatoms =0;
  {
   open(OUTF, "< $outf") or die "Can't open $outf for reading:\n";
   my @keys; my %isKEY;
   my %list;
   while(defined(my $line = <OUTF>) ) {
     my @chunks = split('\t',$line);
     foreach my $test (@chunks) { $test =~ s/\s+$//; $test =~ s/^\s+//; }
     my $ATOM = $chunks[4];
     my $RESI = $chunks[3];
     if(isRNA($RESI) && $chunks[13] eq "\[OK\]") {
        $totatoms++;
        my $key = getKey($RESI,$ATOM);
        my $vol = $chunks[10];
        push(@{$list{$key}},$vol);
        if(!$isKEY{$key}) {
          push(@keys,$key);
          $isKEY{$key} = 1;
        }
     }
   }
   close(OUTF);
   if($debug) { print STDOUT  "RES-ATM count   min     max     mincut  maxcut   median\n"; }
   foreach my $key (@keys) {
     @{$list{$key}} = sort { $a <=> $b } @{$list{$key}};
     $min{$key} = ${$list{$key}}[int(0.05*($#{$list{$key}})+0.5)];
     $max{$key} = ${$list{$key}}[int(0.95*($#{$list{$key}}))];
     if($debug) {
       #print STDOUT "$key\t",fixnum50($#{$list{$key}}+1),"  ",fixnum22(${$list{$key}}[0]),
       #           "\t",fixnum22(${$list{$key}}[$#{$list{$key}}]),"\t",fixnum22($min{$key}),
       #           "\t",fixnum22($max{$key}),"\t",fixnum22(${$list{$key}}[int($#{$list{$key}}/2)]),
       #           "\n";
     }
   }
  } #close PASS ONE

#PASS TWO
   open(OUTF, "< $outf") or die "Can't open $outf for reading:\n";
   my @keys; my %ovol; my %okay; my %ovolsq; my %isKEY;
   my $numatoms =0; my $rnaatoms = 0;
   while(defined(my $line = <OUTF>) ) {
     my @chunks = split('\t',$line);
     foreach my $test (@chunks) { $test =~ s/\s+$//; $test =~ s/^\s+//; }
     my $ATOM = $chunks[4];
     my $RESI = $chunks[3];
     if(isRNA($RESI)) { $rnaatoms++; }
     if(isRNA($RESI) && $chunks[13] eq "\[OK\]") {
       my $key = getKey($RESI,$ATOM);
       my $vol = $chunks[10];
       if($vol > $min{$key} && $vol < $max{$key}) {
         $numatoms++;
         if(!$isKEY{$key}) {
           push(@keys,$key);
           $isKEY{$key} = 1;
           $okay{$key} = 1;
           $ovol{$key} = $vol;
           $ovolsq{$key} = ($vol**2);
         } else {
           $okay{$key}++;
           $ovol{$key} += $vol;
           $ovolsq{$key} += ($vol**2);
         }
       } 
     }
   }
   close(OUTF);

    system("rm -f $outf");

   my %stdev; my %mean; my $totdev=0; my $totperdev=0; my $count=0;
   foreach my $key (@keys) {
     if($okay{$key} > 10) {
       $count++;
       $stdev{$key} = sqrt(($ovolsq{$key}-($ovol{$key}**2)/$okay{$key})/($okay{$key}-1));
       $mean{$key} = $ovol{$key} / $okay{$key};
       $totdev += $stdev{$key};
       $totperdev += $stdev{$key} / $mean{$key};
       #print STDERR "$key\t",fixnum32($mean{$key}),"\t",fixnum22($stdev{$key}),"\t",
       #      fixnum32($totdev),"\n"; 
     } 
     print STDERR ".";
   }
  return ($totdev,$totperdev,$rnaatoms,$numatoms,$count,$totatoms);
};


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

sub writeTime {
  ($outfile) = @_;
  open(OUT, ">> $outfile")
    or die "Couldn't open $outfile for writing: $!\n";
  print OUT getTime(),"\n";
  close(OUT);
};


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

sub editAtomDefs {
   my ($atomdef,$phos,$phos2) = @_;
   open(ATOM, "< $atomdef")
      or die "Couldn't open $atomdef for reading: $!\n";
   my $types = $atomdef;
   substr($types,-3) = "temp.$spec";
   open(TYPES, "> $types")
      or die "Couldn't open $types for writing: $!\n";
   while(defined(my $line = <ATOM>)) {
      my @chunks = split(' ', $line);
      if($chunks[0] eq $type1) {
         print TYPES "$type1   ",fixnum12($phos),"    $bondrad\n"; 
      } elsif($chunks[0] eq $type2) {
         print TYPES "$type2   ",fixnum12($phos2),"    $bondrad\n";
      } else {
         print TYPES "$line";
      }
   }
   close(ATOM);
   close(TYPES);
   return $types;
};


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

sub printOutInfo {
   my ($outfile,$phos,$start,$totdev,$totperdev,
       $rnaatoms,$numatoms,$count,$totatoms,$phos2) = @_;
   open(OUT, ">> $outfile")
      or die "Couldn't open $outfile for writing: $!\n";
   print OUT fixnum12($phos)," ",fixnum12($phos2)," \t",fixnum37(100*$numatoms/$rnaatoms),"\tperok\n";
   print OUT fixnum12($phos)," ",fixnum12($phos2),"  \t",fixnum37($totperdev),"\ttotper\n";
   print OUT "det  ",fixnum12($phos),"  ",fixnum30($count);
   print OUT " totdev ",fixnum35($totdev),"  ",fixnum35(100*$totdev/$count)," tot\%dev ";
   print OUT fixnum25($totperdev),"\%  ",fixnum25(100*$totperdev/$count),"\%\n";
   print OUT "Time:  ",subtTime($start,getTime()),"\tatoms= $numatoms\/$totatoms\/$rnaatoms";
   print OUT "\t",fixnum32(100*$numatoms/$rnaatoms),"\% xx $type1 xx $type2\n";
   close(OUT);
};


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

sub getKey {
  my ($res,$atm) = @_;
  if($atm eq "P" || $atm eq "O1P" || $atm eq "O2P") {
    return "SUG-".$atm;
  } elsif(length($atm) == 3 && substr($atm,2,1) eq "\'") {
    return "SUG-".$atm;
  }
  return $res."-".$atm;
};