#!/usr/bin/perl -w

#use strict;


#多基因的遗传算法

use AI::Genetic;


#        benz_furan_MMTORSION
open FH,"benz_furan_mmenergyGA.tab" or die "can't find the file";
my @qmenergies=<FH>;

chomp(@qmenergies);


map($_=~s/^\d+\s+//,@qmenergies);
#
#map($_*=627.51,@qmenergies);

open FH,"benz_furan_mmenergyGA.tab";
my @mmenergies=<FH>;
chomp(@mmenergies);
map($_=~s/^\d+\s+//s,@mmenergies);


my @angles=(0 .. 72);

@angles=map($_*5,@angles);


my $ga = new AI::Genetic(
       -fitness    => \&fitnessFunc,
       -type       => 'rangevector',
       -population => 1000,
       -crossover  => 0.85,   #一个基因的话 没有交叉的必要0.9
       -mutation   => 0.1,
       -terminate  => \&terminateFunc,
      );



$ga->init([
[1,300],
[1,8],
[0,36],
 
[0,300],
[1,8],
[0,36]
 
 
 
]);#基因的数目



$ga->evolve('rouletteSinglePoint', 1000);     # 选择的方法                   一共进化的代数








my $individual=$ga->getFittest(1);

my $resultx=$individual->genes();
my @resugenes=@$resultx;
 
print  join "**",@resugenes;
print "\n";
 
 my $resulty=$individual->score();
my $rmsdd=1/$resulty;

print "   fittness: $resulty   rmsd:  $rmsdd\n";

 
 




sub fitnessFunc
{
    my $genes = shift;
   # print join "     ", @$genes;
  #  die "111";
    my @valuesmm=();#to save new mm energies
   
    my $gene_pk1=@$genes[0]/10;
    my $gene_pk2=@$genes[3]/10;
     my $gene_pn1=@$genes[1];
     my $gene_pn2=@$genes[4];
     my $gene_phase1=@$genes[2]*10;
     my $gene_phase2=@$genes[5]*10;
     for(0..72)
     {
         
#          print "$mmenergies[$_]";
     #    print "$gene_pk1*(1+cos(($gene_pn1*$angles[$_]-$gene_phase1)/180*3.14159))\n";
     #    die "dd";
            my $energydihe1=$gene_pk1*(1+cos(($gene_pn1*$angles[$_]-$gene_phase1)/180*3.14159));
            my $energydihe2=$gene_pk2*(1+cos(($gene_pn2*$angles[$_]-$gene_phase2)/180*3.14159));
            push @valuesmm,$energydihe2+$energydihe1+$mmenergies[$_];
     }
     my @valuesmm_sort=sort{$a <=> $b} @valuesmm;
       
         @valuesmm=map($_-$valuesmm_sort[0],@valuesmm);
         $rmsd=&calcrmsd(\@valuesmm,\@qmenergies);
         
         #能用一个余弦函数表示就不用两个余弦函数表示；进行罚分处理
         if($gene_pk1==$gene_pk2 and $gene_pn1 == $gene_pn1 and  $gene_phase1 ==  $gene_phase2)
         {
               $rmsd*=2;
         }
     
     
     
     
   

    my $fitness;
   
    $fitness=1/$rmsd;
        print "$gene_pk1 $gene_phase1  $gene_pn1      $gene_pk2   $gene_pn1    $gene_phase2
            fitness: $fitness      $rmsd\n";
    return $fitness;

 
}

sub calcrmsd
{
   my @sander=@{$_[0]};
   my @gauss=@{$_[1]};
   my $rmsd=0;
   if(@sander-0==@gauss-0)
   {
       for(0..@sander-1)
       {
           $rmsd+=($sander[$_]-$gauss[$_])**2;
       }
         return $rmsd/(@sander-0);
       
       
   }
   else
   {
       
       print "error ,the sander and qm have different energy\n";
     exit 0;
   }
}



 sub terminateFunc
 {
    my $ga = shift;
     
    # terminate if reached some threshold.
    my $THRESHOLD=10000;
    return 1 if $ga->getFittest->score > $THRESHOLD;
    return 0;
 }