#/usr/bin/perl
use strict;
require "/home/vossman/lib/fixnum.pl";

srand (12);

my $min=10000;
my $rna=1;
my %isRNA;
my $max=100;
my %isType;
my %fix;
my @types;
my @rules;

open(IN, "< atom-defs.fix") or die;
while(defined(my $line = <IN>)) {
  if(substr($line,0,4) eq "!pro") { $rna = 0; }
  if(substr($line,0,4) eq "BOND") {
    my @chunks = split(' ',$line);
    #print STDERR "BOND  ", fill($chunks[1],7), fill($chunks[2],7), 
	#fixnum43($chunks[3]),"  ", fixnum43($chunks[4]),"\n";
    my $type = substr($chunks[1],0,4); my $type2 = substr($chunks[2],0,4);
    my $rule = $type."_".$type2."_".$chunks[3]."_".$chunks[4];
    push(@rules,$rule);
    if($rna) { $isRNA{$type} = 1; }
    if(!$isType{$type}) {
      push(@types,$type);
      $isType{$type}=rand(1);
    }
    if($rna) { $isRNA{$type2} = 1; }
    if(!$isType{$type2}) {
      push(@types,$type2);
      $isType{$type2}=rand(1);
    }
  }
}
close(IN);

my $olderr=1; my $toterr=0;
for(my $i=0;$i<=$max;$i++) {
  $olderr = $toterr;
  $toterr=0;
  print STDERR "\niter $i\n";
  @types = shuffle(@types);
  if($i==$max) { @types = sort {$a cmp $b } @types; }
  foreach my $type (@types) {
    if($isRNA{$type}) { print STDERR "!"; }
    my $sum=0; my $count=0; my $sumsq =0; my $num=0;
    print STDERR "$type\t";
    foreach my $rule (@rules) {
      my @bits = split("_",$rule);
      if($bits[0] eq $type || $bits[1] eq $type) {
        my $oth;
        if($bits[0] eq $type) { $oth = $bits[1]; }
        else { $oth = $bits[0]; }
        $sum += ($bits[3]-$isType{$oth})*$bits[2];
        $sumsq += ($isType{$type}+$isType{$oth}-$bits[3])**2 *$bits[2];
        $count+=$bits[2]; #no energy
        $num++;
      }
    }
    $isType{$type} = $fix{$type} || ($sum/$count+$isType{$type})/2;
    my $err = sqrt($sumsq/$count)*100;
    print STDERR fixnum25($isType{$type}),"\t",
  	fixnum25($err),"\tN=$num\n";
    $toterr += $err;
  }
  if($min > $toterr) { $min = fixnum35($toterr); }
  $toterr = fixnum35($toterr);
  print STDERR "total error ---------- ",$toterr,"  best $min\n";
}

print STDERR "\#\#\#\# FINAL RESULT \#\#\#\#\n\n";

@types = sort {$a cmp $b } @types;
foreach my $type (@types) {
  if($isRNA{$type}) {
    my $sum=0; my $count=0; my $sumsq =0; my $num=0;
    print STDERR "$type\t";
    foreach my $rule (@rules) {
      my @bits = split("_",$rule);
      if($bits[0] eq $type || $bits[1] eq $type) {
        my $oth;
        if($bits[0] eq $type) { $oth = $bits[1]; }
        else { $oth = $bits[0]; }
        $sum += ($bits[3]-$isType{$oth})*$bits[2];
        $sumsq += ($isType{$type}+$isType{$oth}-$bits[3])**2 *$bits[2];
        $count+=$bits[2];
        $num++;
      }
    }
    my $err = sqrt($sumsq/$count)*100;
    print STDERR fixnum25($isType{$type}),"\t",
        fixnum25($err),"\tN=$num\n";
    $toterr += $err;
  }
  if($min > $toterr) { $min = fixnum35($toterr); }
  $toterr = fixnum35($toterr);
}
print STDERR "total error ---------- ",$toterr,"  best $min\n";

sub shuffle {
  my (@list) = @_;
  my @newlist;
  while($#list >= 0) {
     my $num = int(rand($#list+1));
     push(@newlist,$list[$num]);
     $list[$num]=$list[$#list];
     $#list--;
  }
  return @newlist;
};

sub fill {
  my ($str,$num) = @_;
  while(length($str) < $num) {
    $str .= " ";
  }
  return $str;
};