#!/usr/local/perl
use Storable qw(dclone);
use List::Util qw(shuffle);
use IO::File;
$| = 1;

$InputParameter = $ARGV[0];
open (SQIU, $InputParameter)or die 'Could not open InputParameter.txt for writing';
while (<SQIU>) {
        if ($. == 2){
                print $_;
                chomp;
                @Parameters = split(/\t/,$_);
        }
}
#print $Parameters[1];
$shu1 = $Parameters[1];
$result = $Parameters[2];
mkdir "$result", 0777 unless -d "$result";
open (OUT, ">$result/result_OUT");
OUT -> autoflush;
$Ni = $Parameters[3];
$Nch = $Parameters[4];
$Nmut = $Parameters[5];
$Nrec = $Parameters[6];
$omega = $Parameters[7];
$gene_size = $Parameters[8];
$MatingScheme = $Parameters[9];
$FitCouple = $Parameters[10];
$FixOffSpr = $Parameters[11];
$LeftOffSpr = $Parameters[12];
$TotoOffSpr = $Parameters[13];
@fixed_mut_fit = ();
$total_fit_change = 0;
$total_fixed_mut =0;
$My_fit0 = 0;
$My_fit1 = 0;
$My_fit2 = 0;
$mutations_before = 0;
%mut_hash=();
$control = 0;
$c_repeat =0;

# Insert the Weighted Method to decide the number of OffSprings;
$Weight = 0;
#Initial Setting for $Weight = 0;


goto Normal;

print "Start a new running (A); Run from a backup (Name the backup directory except 'A').\n";
$input = <>;	chomp $input;
if ($input eq 'A'){
        print "Please choose a Selection Coefficient file.\n";
        $shu1 = <>;	chomp $shu1;
        print "Please specify a output file.\n";
        $result = <>;	chomp $result;
        mkdir "$result", 0777 unless -d "$result";
        open (OUT, ">$result/result_OUT");
        OUT -> autoflush;
	print "Please specify the parameters (using tab to separate them):\nNi\tNch\tNmut\tNrec\tomega\tgene_size\n";
	$ARGV = <>;	chomp $ARGV;
	@ARGV = split(/\t/,$ARGV);
	$Ni = $ARGV[0]; # number of individuals in the main cylcle;
	$Nch = $ARGV[1]; #number of childred per mating pair = number of gametes per individual
	$Nmut = $ARGV[2]; #number of mutations per individual
	$Nrec = $ARGV[3]; #number of meiotic recombinations for creation of one gamete
	$omega = $ARGV[4]; #determine dominance mode for gene fitness calculations
	$gene_size = $ARGV[5];
	@fixed_mut_fit = (); # total fitness of all fixed mutations in gene 0, 1, ... ($Ngene-1)
	$total_fit_change = 0; #at the stage of creating new genomic consensus the program calculates the changes in fitness of VI after removal of fix mutations under $change_fit var
	#$total_fit_change is the summ of all $change_fit
	$total_fixed_mut =0; #total number of fixed mutations
	# insert shu
	$My_fit0 = 0;
	$My_fit1 = 0;
	$My_fit2 = 0;
	$mutations_before = 0; # when $mutations_before =1 mutations will be created before recombinations, alternatively when the switch =0 mutations after recombinations
	#shu March 03/21
	%mut_hash=();
	$control = 0;
	$c_repeat =0;
	#
	goto Normal;
}else {
        $result = $input;
        mkdir "$result", 0777 unless -d "$result";
        goto Previous;
}

Previous:
$filenameA = "$result/backup_A";	$filenameB = "$result/backup_B";
print "Running from Backup\n";

$variableA = `tail -1 $result/backup_A`;	$variableB = `tail -1 $result/backup_B`;
chomp $variableA;	chomp $variableB;
@BackA = split(/\t/,$variableA);	@BackB = split(/\t/,$variableB);
#print $variableA."\t".$variableB."\n";
$succ = "SuccessfullyBackup";
if (($BackA[0] eq $succ)&&($BackB[0] ne $succ)) {
	$backup = "$result/backup_A";
	print "step1\t";
	print "run from A file\n";
}elsif  (($BackA[0] ne $succ)&&($BackB[0] eq $succ)){
	$backup = "$result/backup_B";
	print "step2\t";
	print "run from B file\n";
}elsif (($BackA[0] eq $succ)&&($BackB[0] eq $succ)){
	if ($BackA[1]> $BackB[1]){
		$backup = "$result/backup_A";
		print "step3\t";
	        print "run from A file\n";
	}else {
		$backup = "$result/backup_B";
		print "step4\t";
        	print "run from B file\n";
	}
}
#my @arrayA_file = stat("$result/backup_A");	my @arrayB_file = stat("$result/backup_B");
#if ($arrayA_file[7] > $arrayB_file[7]){
#	$backup = "$result/backup_A";
#	print "run from A file\n";
#}else {
#	$backup = "$result/backup_B";
#	print "run from B file\n";
#}

open FILEHANDLEBACKUP,$backup or die "can't open the $backup";
@Backup = <FILEHANDLEBACKUP>;	$count_gene = 0;
# Read the parameters in Backupfile into the program;
for($i = 0; $i<=$#Backup; $i++){
        $str1 = $Backup[$i];      chomp $str1;
        @tempArray = split(/\t/,$str1);
	if ($i == 0) { 
		$gbackup = $tempArray[0]+1;#when it runs from the backup, it means that the programme will start from the Next Generation;	
		$total_fixed_mut = $tempArray[1];	$Ngene = $tempArray[2];	$Ni = $tempArray[3];	$Nch = $tempArray[4];	$#Fitnesstable = $tempArray[5];
	 	$gene_size = $tempArray[6];	$Nmut = $tempArray[7];	$Nrec = $tempArray[8];	$omega = $tempArray[9];	$total_fit_change = $tempArray[10];
	 	$mutations_before = $tempArray[11];
	}
	if ($tempArray[0] =~ m/pi*/ ) {
		${$tempArray[0]}[$tempArray[1]][$tempArray[2]][$tempArray[3]] = $tempArray[4];
	}
	if ($tempArray[0] =~ m/mi*/ ) {
                ${$tempArray[0]}[$tempArray[1]][$tempArray[2]][$tempArray[3]] = $tempArray[4];
        }
	if ($#tempArray == 5){
		$meme =  $count_gene % $gene_size;
		${'matrix'.int($count_gene/$gene_size)}[$meme][0] = $tempArray[0];
		${'matrix'.int($count_gene/$gene_size)}[$meme][1] = $tempArray[1];
		${'matrix'.int($count_gene/$gene_size)}[$meme][2] = $tempArray[2];
		${'matrix'.int($count_gene/$gene_size)}[$meme][3] = $tempArray[3];
		${'matrix'.int($count_gene/$gene_size)}[$meme][4] = $tempArray[4];
		${'matrix'.int($count_gene/$gene_size)}[$meme][5] = $tempArray[5];
		$count_gene ++;
	}
	if ($#tempArray == ($Ngene + 1)){
		for ($j=0; $j< $Ngene; $j++){
			push (@{$tempArray[$Ngene].$tempArray[$Ngene+1]},$tempArray[$j]);
		}
	}

	if ($#tempArray == ($Ngene - 1)){
		for ($j=0; $j< ($Ngene-1); $j++){
                	push (@fixed_mut_fit, $tempArray[$j]);
		}
	}
}
open (OUT, ">>$result/result_OUT.new"); #write the results into a new file;
OUT -> autoflush;
$starting = $gbackup;
print "finish initiation \n";
goto BACKUP;

Normal:
print "normal starting generation.\n";
open FILEHANDLESHU,$shu1 or die "can't open the $shu1";
@Fitnesstable = <FILEHANDLESHU>; 
shift @Fitnesstable; # get rid of the first header line;
 
$Ngene = int(($#Fitnesstable)/$gene_size); 
for($i = 0; $i<=$#Fitnesstable; $i++){
        $str1 = $Fitnesstable[$i];      chomp $str1;
        @$i = split(/\t/,$str1);
        push (@{'matrix'.(int($i/$gene_size))},\@$i);
}
# v17 created a switch: $mutations_before (and after)
sub recombination {
	if ($Nrec < 1) { #addition for MAGE16 for number_of_recombinations <1
		my $chance =rand(1);
		if ($Nrec < $chance) { $rec_n =0; }
		else {$rec_n =1;}
	}
	else {$rec_n = $Nrec;}
	my @gene =(); my @rec =();
	for my $x (1..$rec_n) {push(@gene, int(rand($Ngene))); }
	my @sort_g =sort {$a <=> $b} @gene;
	for my $x (0..$#sort_g) {
		push(@rec, $sort_g[$x]);
		my $pos = int(rand($gene_size +1));
		push(@rec, $pos);
	}
	push(@rec, $Ngene);
	push(@rec, $gene_size +1);
	return @rec;
}
print OUT "ind $Ni children $Nch genes $Ngene mutations $Nmut recombinations $Nrec omega $omega genesize $gene_size  before/after $mutations_before \n"; 
## STEP 1: CREATION OF MUTATIONS IN THE GENOMES OF VI
sub mutation {
	my $pos = int(rand($gene_size));
	 SHUHAO1:
        #Shuhao         insert a label to avoid the mutation mutates into the original one;
        my $toss = int(rand(4))+1;
	my $nt = '';
        # shu
        if ($toss eq '1') {$nt = 'A';}
        elsif ($toss eq '2') {$nt = 'C';}
        elsif ($toss eq '3') {$nt = 'T';}
        else {$nt = 'G';}
        if ($nt eq ${matrix.$_[0]}[$pos][5]){
                goto SHUHAO1;
        }
# Change into No Matrix
#	if (${matrix.$_[0]}[$pos][5] eq 'A'){ $fitness = ${matrix.$_[0]}[$pos][$toss]-0;}
#	elsif (${matrix.$_[0]}[$pos][5] eq 'C'){ $fitness = ${matrix.$_[0]}[$pos][$toss]-0;}
#	elsif (${matrix.$_[0]}[$pos][5] eq 'T'){ $fitness = ${matrix.$_[0]}[$pos][$toss]-0;}
#	else { $fitness = ${matrix.$_[0]}[$pos][$toss]-0;}
	my $y = int(rand(100))+1;
        if ($y <= 81){
               $fitness = -1;
          }
        if (($y > 81)&&($y<=91)){
               $fitness = 0;
          }
        if (($y > 91)&&($y<=100)){
               $fitness = 1;
          }
	return [$pos, $nt, $fitness];
}

$starting = 1;
BACKUP:
#put the BACKUP after the $starting; then I can change $starting in front/ original setting the $starting is == 1;
for ($g=$starting; $g <= 10000; $g++) {#$g generation number
	print "gen $g \t";

@a =();
%pairs =();
%fit = ();
#$MatingScheme = 3;
if($MatingScheme == 1){
        &MatingScheme1;
} elsif ($MatingScheme == 2) {
        &MatingScheme2;
} elsif ($MatingScheme == 3){
        &MatingScheme3;
} elsif ($MatingScheme == 4){
        &MatingScheme4;
} elsif ($MatingScheme == 5) {
        &MatingScheme5;
} elsif ($MatingScheme == 6){
        &MatingScheme6;
} elsif ($MatingScheme == 7) {
        &MatingScheme7;
} else {
	&MatingScheme8;	
}


	if ($mutations_before) {
	for $k (1..$Ni) {
		$repeat = 0;
REPEAT_M:	while ( $repeat < $Nmut) {
			@mut=(); undef($block);
			my $block =int(rand($Ngene));
			my @mut = &mutation($block);
			for $c (0..$#{${'mi'.$k}[$block]}) {
                                if (${'mi'.$k}[$block][$c][0] == $mut[0][0]){
					next REPEAT_M;
				}
                        }
			push(@{${'mi'.$k}[$block]}, @mut);
			${'Fmi'.$k}[$block] += $mut[0][2];
			$repeat++;
		}
	}
	for $k (1..$Ni) {
                $repeat = 0;
REPEAT_P:       while ( $repeat < $Nmut) {
                        @mut=(); undef($block);
                        my $block =int(rand($Ngene));
                        my @mut = &mutation($block);
                        for $c (0..$#{${'pi'.$k}[$block]}) {
                                if (${'pi'.$k}[$block][$c][0] == $mut[0][0]){
					next REPEAT_P;
				}
                        }
                        push(@{${'pi'.$k}[$block]}, @mut);
                        ${'Fpi'.$k}[$block] += $mut[0][2];
                        $repeat++;
                }
        }
} #close $mutations_before switch 
# STEP 2:  CREATION OF GAMETES FROM MATERNAL AND PATERNAL CHR OF VI
for my $k (1..$Ni) {

	if($Weight == 1){
		$Nch2 = $OffSpringNumber{$k};
	}else {
		$Nch2 = $Nch;
	}
 
for $m (1..$Nch2) {
	@rec = &recombination; 
	$name = 'mi'.$k; $altern = 'pi'.$k;
	%switch = ($name => $altern, $altern => $name);
	@g =(); $start =0; $flag =1; @Fg =();
	$r01 = int(rand(2)); @keys = keys (%switch); $name = $keys[$r01]; $altern = $switch{$name};
	for ($n=0; $n < ($#rec); $n +=2) {
		@new =(); $F =0;
		if ($flag) {
			for $x ($start..($rec[$n]-1) ) {
				push (@g, ${$name}[$x]);
				push (@Fg, ${'F'.$name}[$x]);  
			}
		}
		$coord = $rec[$n+1]; 
		if ($flag) { # one recombination point in the same block (gene)
			$F = $fixed_mut_fit[$rec[$n]]; # addition to MAGEv14.pl to take into consideration fixed mutations
			for $x (0..$#{${$name}[$rec[$n]]}) {
				if (${$name}[$rec[$n]][$x][0] <= $coord) {
					push(@new, ${$name}[$rec[$n]][$x]);
					$F +=${$name}[$rec[$n]][$x][2];
				}
			}
		}
		else {  # multiple recombinations, thus, process @temp instead of @mi5 or @pi5
			$F = $fixed_mut_fit[$rec[$n]]; # addition to MAGEv14.pl to take into consideration fixed mutations
			for $x (0..$#temp) {
                		if ($temp[$x][0] <= $coord) {
                        		push(@new, $temp[$x]);
					$F +=$temp[$x][2];
                		}
        		}
		}
		for $x (0..$#{${$altern}[$rec[$n]]}) {
			if (${$altern}[$rec[$n]][$x][0] > $coord) {
				push(@new, ${$altern}[$rec[$n]][$x]);
				$F +=${$altern}[$rec[$n]][$x][2];
			}
		}
		if ($rec[$n] == $rec[$n+2]) {#check for multiple recombinations inside the same block
			$flag = 0; @temp = ();
			@temp = @{dclone(\@new)}; @new =();
			$altern = $name; $name = $switch{$name};
		}
		else{
			$flag =1; @temp =();
			$start = $rec[$n]+1;      
			$altern = $name; $name = $switch{$name};
			push (@g, [@new]); 
			push (@Fg, $F);
			@new=();
		}
	}
	pop @Fg; pop @g;
	@{'Fg'.$k. '_'.$m} = @{ dclone(\@Fg) };
	@{'g'.$k. '_'.$m} = @{ dclone(\@g) };
}
}

unless ($mutations_before) {
	for $k (1..$Ni) {

	if($Weight == 1){
                $Nch2 = $OffSpringNumber{$k};
        }else {
		$Nch2 = $Nch;
	}

	for my $m (1..$Nch2) {
		$repeat = 0;
REPEAT_M:	while ( $repeat < $Nmut) {
			@mut=(); undef($block);
			my $block =int(rand($Ngene));
			my @mut = &mutation($block);
			for $c (0..$#{${'g'.$k.'_'.$m}[$block]}) {
                                if (${'g'.$k.'_'.$m}[$block][$c][0] == $mut[0][0]){
					$c_repeat ++;
					next REPEAT_M;
				}
                        }
			push(@{${'g'.$k.'_'.$m}[$block]}, @mut);
			#shu the hash
			$mut_hash{$mut[0][2]}++;
			#
			${'Fg'.$k.'_'.$m}[$block] += $mut[0][2];
			$repeat++;
		}
	}
	}
}
	#  STEP 3 MATING PAIRS AND CREATION OF CHILDREN GENOMES
	# BLOCK FOR GENERATION OF MATING PAIRS FOR 10 individuals
	#	my @a =(); my %pairs =();
	#	for my $n (1..$Ni) {push(@a, $n);}
	#	use List::Util qw(shuffle);
	#	@ra = shuffle(@a);
	#	for my $n (1..($Ni/2)) { # number of individuals /2
	#	my $first = shift (@ra);
	#	my $second = shift (@ra);
	#	$pairs{$first} = $second;
	#	}
#SHU MatingSubroutine;i
#@a =(); 
#%pairs =();
#%fit = ();
##$MatingScheme = 3;
#if($MatingScheme == 1){
#	&MatingScheme1;
#} elsif ($MatingScheme == 2) {
#	&MatingScheme2;
#} elsif ($MatingScheme == 3){
 #       &MatingScheme3;
#} elsif ($MatingScheme == 4){
#	&MatingScheme4;
#} elsif ($MatingScheme == 5) {
#	&MatingScheme5;
#} else {
#	&MatingScheme6;
#}

# CHILDREN FITNESSES COMPARISON BLOCK
if ($MatingScheme == 1){
	foreach my $x (keys %pairs) {
		my $i1 = 'Fg'.$x;
		my $i2 = 'Fg'.$pairs{$x};
		for my $m (1..$Nch) {
			$fit{$x. '_'. $m} = &fitness(@{$i1 . '_' . $m} , @{$i2 . '_' . $m});#shu Here it stores the fitness for each block;
		}
	}
} elsif ($MatingScheme == 2){
	foreach $x (keys %pairs) {
		$i1 = 'Fg'.$x;
    	        $i2 = 'Fg'.$pairs{$x};
		$fit{$x} = &fitness(@$i1 , @$i2);
	}
} elsif ($MatingScheme == 3){
	foreach $x (keys %pairs) {
                $i1 = 'Fg'.$x;
                $i2 = 'Fg'.$pairs{$x};
                $fit{$x} = &fitness(@$i1 , @$i2);
        }
}elsif ($MatingScheme == 4){
	foreach $x (keys %pairs) {
                $i1 = 'Fg'.$x;
                $i2 = 'Fg'.$pairs{$x};
                $fit{$x} = &fitness(@$i1 , @$i2);
        }
} else {
	foreach $x (keys %pairs) {
                $i1 = 'Fg'.$x;
                $i2 = 'Fg'.$pairs{$x};
                $fit{$x} = &fitness(@$i1 , @$i2);
        }
}

while ( ($k,$v) = each %pairs ) {
#	print "$k => $v => $fit{$k}\n";
# 	print "--------------------11111111111111111111-------------------------------\n";
}
#sorting children by their fitness
@sort_fit1 = ();
@sort_fit1 = sort sort_hash (keys(%fit));

@sort_fit = ();
@sort_fit = splice (@sort_fit1,0,$Ni);

# STEP 4:  CONVERT $Ni number of fitest children to current individuals  and empty children's data
$total_fit =0;  #total fitness for the current generation (summ of fitnesses of all individuals)
for my $z (1..$Ni) {
	$total_fit +=$fit{$sort_fit[$z-1]};
	my @new_vi =split(/_/, $sort_fit[$z-1]);
	@{'old_mi'.$z} =@{dclone(\@{'mi'.$z})};
	@{'mi'.$z} = @{dclone(\@{'g'. $new_vi[0] . '_'.$new_vi[1]})};
	@{'old_pi'.$z} = @{dclone(\@{'pi'.$z})};
	@{'old_Fmi'.$z} = @{'Fmi'.$z};
	@{'Fmi'.$z} = @{'Fg'. $new_vi[0] . '_'.$new_vi[1]};
	@{'old_Fpi'.$z} = @{'Fpi'.$z};
#Shu	Different MatingSchemes Have Different Hash Pairs;
	if ($MatingScheme == 1){
		 @{'pi'.$z} = @{dclone(\@{'g'. $pairs{$new_vi[0]} . '_'.$new_vi[1]})};
		 @{'Fpi'.$z} = @{'Fg'. $pairs{$new_vi[0]} . '_'.$new_vi[1]};
	}elsif ($MatingScheme == 2) {
		my @new_Pvi = split (/_/, $pairs{$sort_fit[$z-1]});
		@{'pi'.$z} = @{dclone(\@{'g'. $new_Pvi[0] . '_'.$new_vi[1]})};
	        @{'Fpi'.$z} = @{'Fg'. $new_Pvi[0] . '_'.$new_vi[1]};
#		print "$new_Pvi[0]\t$new_vi[0]\t$new_vi[1]"."ttttttttttttttttttttt\n";
		@new_Pvi = ();
	}elsif ($MatingScheme == 3) {
		my @new_Pvi = split (/_/, $pairs{$sort_fit[$z-1]});
                @{'pi'.$z} = @{dclone(\@{'g'. $new_Pvi[0] . '_'.$new_vi[1]})};
                @{'Fpi'.$z} = @{'Fg'. $new_Pvi[0] . '_'.$new_vi[1]};
 #               print "$new_Pvi[0]\t$new_vi[0]\t$new_vi[1]"."ttttttttttttttttttttt\n";
                @new_Pvi = ();
	}elsif ($MatingScheme == 4) {
		my @new_Pvi = split (/_/, $pairs{$sort_fit[$z-1]});
                @{'pi'.$z} = @{dclone(\@{'g'. $new_Pvi[0] . '_'.$new_vi[1]})};
                @{'Fpi'.$z} = @{'Fg'. $new_Pvi[0] . '_'.$new_vi[1]};
 #               print "$new_Pvi[0]\t$new_vi[0]\t$new_vi[1]"."ttttttttttttttttttttt\n";
                @new_Pvi = ();
	}else {
		my @new_Pvi = split (/_/, $pairs{$sort_fit[$z-1]});
                @{'pi'.$z} = @{dclone(\@{'g'. $new_Pvi[0] . '_'.$new_vi[1]})};
                @{'Fpi'.$z} = @{'Fg'. $new_Pvi[0] . '_'.$new_vi[1]};
#                print "$new_Pvi[0]\t$new_vi[0]\t$new_vi[1]"."ttttttttttttttttttttt\n";
                @new_Pvi = ();
	}
}
print "total fitness is $total_fit \n";

for $k (1..$Ni) {
        for $m (1..$Nch) {
                @{'Fg'.$k . '_' .$m} = ();
                @{'g'.$k . '_' .$m} = ();
        }
}
# optional STEP 5 for creation of a new genome consensus and remove fix mutations
unless ($g%10) {
	print "generation $g \n";
	%seen =(); $new_fit=0;
	$mutation_number =0;
	for $k (1..$Ni) {
        	$nameP = 'pi'.$k;
        	$nameM = 'mi'.$k;
		$count_mut =0;
        	for $a (0..$#{$nameP}){
                	for $b (0..$#{${$nameP}[$a]}) {
				$curr =$a;
                        	for $c (0..$#{${$nameP}[$a][$b]}) {
                                	$curr .= '_'.${$nameP}[$a][$b][$c];						
                        	}
				$count_mut++;
				$seen{$curr}++;
                	}
        	}
        	for $a (0..$#{$nameM}){
                	for $b (0..$#{${$nameM}[$a]}) {
				$curr =$a;
                        	for $c (0..$#{${$nameM}[$a][$b]}) {
                                	$curr .= '_' . ${$nameM}[$a][$b][$c];
				}
				$count_mut++;
				$seen{$curr}++;
                	}
        	}
		print "IND $k total_mutations $count_mut \n";
	}
	foreach $key (sort {$seen{$b} <=> $seen{$a};} (keys(%seen))) {
		$mutation_number++;
	}
	
	$total_SNPs = 0;
	$count_fixed=0;
	$fit0 = 0; $fit1 = 0; $fit2 = 0; $control = 0;
	for $k (1..$Ni) {
                $nameP = 'pi'.$k;
                $nameM = 'mi'.$k;
		$count2_mut = 0;
		$count_fixM = $count_fixP = 0;
		@tempP = (); #major temporary file for updated individual
		for $a (0..$#{$nameP}){
			@T = (); #temporary file for updated gene
                        for $b (0..$#{${$nameP}[$a]}) {
                                $curr =$a;
                                for $c (0..$#{${$nameP}[$a][$b]}) {
                                        $curr .= '_'.${$nameP}[$a][$b][$c];
                                }
                                if ($seen{$curr} >= ($Ni*2)) {$count_fixP++;
					if($k ==1) { #conunt fitness for removed fixed mutations to control and check the algorithm	
						$fixed_mut_fit[$a] +=${$nameP}[$a][$b][2];
						$count_fixed++;
						${matrix.$a}[(${$nameP}[$a][$b][0])][5] = ${$nameP}[$a][$b][1];
					#shu insert
						if (${$nameP}[$a][$b][2] == "0"){$fit0 ++; $control ++;}
						elsif (${$nameP}[$a][$b][2] == "1"){$fit1 ++; $control ++;}
						elsif (${$nameP}[$a][$b][2] == "-1"){$fit2 ++; $control ++;} 	
						else {$control++; print "The one"."\t".${$nameP}[$a][$b][2]."\n"; }
						 #shu    Establish a hash to store the FixedMutation and Values;
                                                # $fixFit{${matrix.$a}[(${$nameP}[$a][$b][0])][0]} = ${$nameP}[$a][$b][2];


					}  #total fitness of all fixed mutations in gene $a
				}
				else {
					$count2_mut++;
					push (@T, ${$nameP}[$a][$b]);
				}
                        }
			push @tempP, [@T];
                }
		@tempM = (); #major temporary file for updated individual
                for $a (0..$#{$nameM}){
			@T = (); #temporary file for updated gene
                        for $b (0..$#{${$nameM}[$a]}) {
                                $curr =$a; 
                                for $c (0..$#{${$nameM}[$a][$b]}) {
                                        $curr .= '_' . ${$nameM}[$a][$b][$c];
                                }
				if ($seen{$curr} >= ($Ni*2)) {$count_fixM++; }
				else {
					$count2_mut++;
					push (@T, ${$nameM}[$a][$b]);
				}
			
                        }
			push @tempM, [@T];
                }
		$total_SNPs +=$count2_mut;
		print "ind $k \t countM $count_fixM \t countP $count_fixP remained SNPs $count2_mut\n";
	@{'pi'.$k} = @{dclone(\@{tempP})};
	@{'mi'.$k} = @{dclone(\@{tempM})};
	$new_fit +=&fitness(@{'Fpi' . $k} , @{'Fmi' . $k});
}
$SNP = $total_SNPs / $Ni;
$fit_change =$total_fit - $new_fit; $total_fit_change +=$fit_change;
$SNP_number = $mutation_number - $count_fixed; 
$total_fixed_mut +=$count_fixed;
# insert shu
$My_fit0 += $fit0;
$My_fit1 += $fit1;
$My_fit2 += $fit2;
$control_total += $control;
$time = localtime;

$Ave_Fit = $total_fit/$Ni;
$Ave_SNP = $SNP_number / $Ni;
print OUT $g, "\t", $Ave_Fit, "\t", $Ave_SNP, "\t", $total_fixed_mut, "\t",$My_fit0,"\t",$My_fit1,"\t",$My_fit2,"\t",$control_total,"\t",$mut_hash{0}," \t", $mut_hash{1}," \t", $mut_hash{-1},"\t",$c_repeat,"\n";
} # loop unless
#shu Create the Backup file
if ($g % 20 == 0){
	if ($g % 40 == 0){
		$tag = 'B';
		open (BACKUPB, ">$result/backup_B");
		BACKUPB -> autoflush;
	}else {
		$tag = 'A';
		open (BACKUPA, ">$result/backup_A");
		BACKUPA -> autoflush;
	}
	print {BACKUP.$tag} "$g\t$total_fixed_mut\t$Ngene\t$Ni\t$Nch\t$#Fitnesstable\t$gene_size\t$Nmut\t$Nrec\t$omega\t$total_fit_change\t$mutations_before\n";
         for $k (1..$Ni) {
                $nameP = 'pi'.$k;
                $nameM = 'mi'.$k;
                $fit_count = 0;
                 for $a (0..$#{$nameP}){
                        for $b (0..$#{${$nameP}[$a]}) {
                                for $c (0..$#{${$nameP}[$a][$b]}){
                                        print {BACKUP.$tag} $nameP."\t".$a."\t".$b."\t".$c."\t". ${$nameP}[$a][$b][$c]."\n";
                                }
                        }
                }
                 for $a (0..$#{$nameM}){
                        for $b (0..$#{${$nameM}[$a]}) {
                                for $c (0..$#{${$nameM}[$a][$b]}){
                                        print {BACKUP.$tag} $nameM."\t".$a."\t".$b."\t".$c."\t". ${$nameM}[$a][$b][$c]."\n";
                                }
                        }
                }
        }
        for($i = 0; $i<=int($#Fitnesstable/$gene_size); $i++){
                for ($x = 0; $x < $gene_size; $x++){
                        for ($y = 0; $y <6; $y++){
                                print  {BACKUP.$tag} ${'matrix'.(int($i))}[$x][$y]."\t";
                        }
                        print  {BACKUP.$tag} "\n";
                }
        }
	for $z (1..$Ni) {

		for $j (0..$#{'Fmi'.$z}){
			print {BACKUP.$tag} ${'Fmi'.$z}[$j]."\t"
		}
		print {BACKUP.$tag} "Fmi\t".$z."\n";
        	for $l (0..$#{'Fpi'.$z}){
                        print {BACKUP.$tag} ${'Fpi'.$z}[$l]."\t"
                }
                print {BACKUP.$tag} "Fpi\t".$z."\n";
	}
	
	for $u (0..$#fixed_mut_fit){
		print {BACKUP.$tag} $fixed_mut_fit[$u]."\t";
	}
	print {BACKUP.$tag} "\n";
	# Test
#	while ( ($key, $value) = each %fixFit ){
#                print {BACKUP.$tag}  $key."\t". $fixFit{$key}."\n";
#        }
	print {BACKUP.$tag} "SuccessfullyBackup\t$g\n";	

}#loop another unless
} #end loop $g

sub sort_hash {
	$fit{$b} <=> $fit{$a};
}
sub fitness {
	my $F =0;
	for my $n (0..($#_ /2)) {
		if ($_[$n] >= $_[$n+$Ngene]) {$F += $_[$n+$Ngene] + ($_[$n] - $_[$n+$Ngene])*$omega; }
		else {$F += $_[$n] + ($_[$n+$Ngene] - $_[$n])*$omega; }
	}
	return $F;
}

sub MatingScheme1{
	for $n (1..$Ni) {
		push(@a, $n);
	}
	@ra = ();
	@ra = shuffle (@a);
	for $n (1..($Ni/2)) { # number of individuals /2
        	$first = shift (@ra);
        	$second = shift (@ra);
        	$pairs{$first} = $second;
#		print "$first\t$second\n";
	}
	return %pairs;
}

sub MatingScheme2{#Mimic the process of promiscuity;	
        for $n (1..$Ni) {
                push(@a, $n);
        }
	my $Npairs = $Ni/2;
	my @Temp, @Male, @Female, @ra = ();
	my @Temp = shuffle (@a);
	my @Male = splice(@Temp,0,$Npairs);
	for $r(1..$Nch){
		@ra = shuffle (@Male);
		@Female = @Temp;
        	for $n (1..($Ni/2)) { # number of individuals /2
                	$first = shift (@ra);
                	$second = shift (@Female);
                	$pairs{$first.'_'.$r} = $second.'_'.$r;
#			print "$first.'_'.$r\t$second.'_'.$r\t$pairs{$first.'_'.$r}\n";
        	}
		@ra = ();
		@Female = ();
	}
	return %pairs;
}

sub MatingScheme3{
	for $n (1..$Ni) {
                push(@a, $n);
        }
        my $Npairs = $Ni/2;
        my @Temp, @Male, @Female, @ra = ();
        my @Temp = shuffle (@a);
        my @Male = splice(@Temp,0,$Npairs);
	for $r(1..$Nch){
		my @ra = sort @Male;
        	my @Female = sort @Temp;
		for $n (1..($Ni/2)) {
			$first = shift (@ra);
			$second = shift (@Female);
			$pairs{$first.'_'.$r} = $second.'_'.$r;
#			print "$first.'_'.$r\t$second.'_'.$r\t$pairs{$first.'_'.$r}\n";
		}
		@ra = ();
                @Female = ();
	}
	return %pairs;
}

sub MatingScheme4{
        for $n (1..$Ni) {
                push(@a, $n);
        }
        my $Npairs = $Ni/2;
        my @Temp, @Male, @Female, @ra = ();
        my @Temp = shuffle (@a);
        my @Male = splice(@Temp,0,$Npairs);
	my @ra = sort @Male;
        for $r(1..$Nch){
                my @Female = sort @Temp;
                for $n (1..($Ni/2)) {
                        $first = shift (@Female);
                        $second = $ra[0];
                        $pairs{$first.'_'.$r} = $second.'_'.$r;
                        print "$first.'_'.$r\t$second.'_'.$r\t$pairs{$first.'_'.$r}\n";
                }
                @Female = ();
        }
        return %pairs;
}

sub MatingScheme5{
        for $n (1..$Ni) {
                push(@a, $n);
        }
	my $Npairs = $Ni/2;
        my @Temp, @Male, @Female, @ra = ();
        my @Temp = shuffle (@a);
        my @Male = splice(@Temp,0,$Npairs);
        my @ra = sort @Male;
	my @Female = @Temp;
	$size = 0;
	$sizeFit = 0;	
	$sizeNotFit = 0;

	$Weight = 1;
	for $x(0..($Npairs-1)){
                if ($x <= ($FitCouple-1)){
                        $OffSpringNumber{$ra[$x]} = $FixOffSpr;
			$sizeFit += $FixOffSpr;
#			print "$ra[$x]\t$OffSpringNumber{$ra[$x]}\tFit\n";
                }else {
                        $OffSpringNumber{$ra[$x]} = $LeftOffSpr;
			$sizeNotFit += $LeftOffSpr;
#			print "$ra[$x]\t$OffSpringNumber{$ra[$x]}\tNOT Fit\n";
                }
        }
	$size = $sizeFit + $sizeNotFit;
        print "size\t$sizeFit\t$sizeNotFit\t$size\n";
        $Diffe = $size - $TotoOffSpr;
 #	print "Diffe\t$Diffe\n";
        while ($Diffe >= 1){
                $RandNum = int (rand($Npairs));
		$Diffe --;
#		print "Delete\t$RandNum\n";	
		if (($OffSpringNumber{$ra[$RandNum]}) >= 2){
#			print "$ra[$RandNum]\t$OffSpringNumber{$ra[$RandNum]}\n";
                	($OffSpringNumber{$ra[$RandNum]})--;
#			print "$ra[$RandNum]\t$OffSpringNumber{$ra[$RandNum]}\n";
		}else{
			print "WaKAKA\n";
			$Diffe ++;	
		}
        }

	for $x(0..$#ra){
		$Numb = $OffSpringNumber{$ra[$x]};
#		print "numberOfOffSpr\t$Numb\n";
		for $r(1..$Numb){# Here, the $r is the number of Children;
			$first = $ra[$x];
			$second = $Female[$x];
			$pairs{$first.'_'.$r} = $second.'_'.$r;
#			print "$first.'_'.$r\t$second.'_'.$r\t$pairs{$first.'_'.$r}\n";
		}
		$OffSpringNumber{$Female[$x]} = $OffSpringNumber{$ra[$x]};
	}
	return %OffSpringNumber;
	return %pairs;
}

sub MatingScheme6{
        for $n (1..$Ni) {
                push(@a, $n);
        }
        my $Npairs = $Ni/2;
        my @Temp, @Male, @Female, @ra = ();
        my @Temp = shuffle (@a);
        my @Male = splice(@Temp,0,$Npairs);
        my @ra = sort @Male;
        my @Female = sort @Temp;
        $size = 0;

        $Weight = 1;
        for $x(0..($Npairs-1)){
                if ($x <= ($FitCouple-1)){
                        $OffSpringNumber{$ra[$x]} = $FixOffSpr;
                        $size += $FixOffSpr;
        #               print "$ra[$x]\t$OffSpringNumber{$ra[$x]}\tFit\n";
                }else {
                        $OffSpringNumber{$ra[$x]} = $LeftOffSpr;
                        $size += $LeftOffSpr;
        #               print "$ra[$x]\t$OffSpringNumber{$ra[$x]}\tNOT Fit\n";
                }
        }

     #   print "size\t$size\n";
        $Diffe = $size - $TotoOffSpr;
 #      print "Diffe\t$Diffe\n";
        while ($Diffe >= 1){
                $RandNum = int (rand($Npairs));
                $Diffe --;
#               print "Delete\t$RandNum\n";
                if (($OffSpringNumber{$ra[$RandNum]}) >= 2){
#                       print "$ra[$RandNum]\t$OffSpringNumber{$ra[$RandNum]}\n";
                        ($OffSpringNumber{$ra[$RandNum]})--;
#                       print "$ra[$RandNum]\t$OffSpringNumber{$ra[$RandNum]}\n";
                }else{
                        print "WaKAKA\n";
                        $Diffe ++;
                }
        }

        for $x(0..$#ra){
                $Numb = $OffSpringNumber{$ra[$x]};
#               print "numberOfOffSpr\t$Numb\n";
                for $r(1..$Numb){# Here, the $r is the number of Children;
                        $first = $ra[$x];
                        $second = $Female[$x];
                        $pairs{$first.'_'.$r} = $second.'_'.$r;
#                       print "$first.'_'.$r\t$second.'_'.$r\t$pairs{$first.'_'.$r}\n";
                }
                $OffSpringNumber{$Female[$x]} = $OffSpringNumber{$ra[$x]};
        }
        return %OffSpringNumber;
        return %pairs;

}

sub MatingScheme7{
        for $n (1..$Ni) {
                push(@a, $n);
        }
        my $Npairs = $Ni/2;
        my @Temp, @Male, @Female, @ra = ();
        my @Temp = shuffle (@a);
        my @Male = splice(@Temp,0,$Npairs);
        for $r(1..$Nch){
                my @ra = sort @Male;
                my @Female = reverse (sort @Temp);
#		print "$Female[0]\t$Female[1]\t$Female[2]\t1111111111\n";
                for $n (1..($Ni/2)) {
                        $first = shift (@ra);
                        $second = shift (@Female);
                        $pairs{$first.'_'.$r} = $second.'_'.$r;
#                       print "$first.'_'.$r\t$second.'_'.$r\t$pairs{$first.'_'.$r}\n";
                }
                @ra = ();
                @Female = ();
        }
        return %pairs;
}

sub MatingScheme8{
	for $n (1..$Ni) {
                push(@a, $n);
        }
	my $Npairs = $Ni/2;
	my @Temp, @Mal, @Femal, @Male, @Female = ();
	my @Temp = shuffle (@a);
	my @Mal = splice(@Temp,0,$Npairs);
	my @Femal = @Temp;

	print "Mark1\n";

	for $n(0..$#Mal){
		print "Mark2\t$Mal[$n]\n";
		$nameP = 'pi'.{$Mal[$n]};
		$nameM = 'mi'.{$Mal[$n]};
		$count_snp = 0;
		for $a (0..$#{$nameP}){
			print "Mark3\n";	
			for $b (0..$#{${$nameP}[$a]}){
				$count_snp += $#{${$nameP}[$a][$b]};
				print "mark\t$#{${$nameP}[$a][$b]}\n";
			}
		}
		for $a (0..$#{$nameM}){
                        for $b (0..$#{${$nameM}[$a]}){
                                $count_snp += $#{${$nameM}[$a][$b]};
                        }
                }
		$SNP_hash{$Mal[$n]} = $count_snp;
	}
	
	@Male = sort {$SNP_hash{$a} cmp $SNP_hash{$b}} keys %SNP_hash;
	while (my ($k,$v)=each %SNP_hash){print "$k\t$v\tMale\n"};
	foreach (@Male) {
  		print "$_\n";
	}
	%SNP_hash = ();

	for $n(0..$#Femal){
                $namP = 'pi'.$Femal[$n];
                $namM = 'mi'.$Femal[$n];
                $count_snp = 0;
                for $a (0..$#{$namP}){
                        for $b (0..$#{${$namP}[$a]}){
                                $count_snp += $#{${$namP}[$a][$b]};
                        }
                }
                for $a (0..$#{$namM}){
                        for $b (0..$#{${$namM}[$a]}){
                                $count_snp += $#{${$namM}[$a][$b]};
                        }
                }
                $SNP_hash{$Mal[$n]} = $count_snp;

        }
	@Female = sort {$SNP_hash{$a} cmp $SNP_hash{$b}} keys %SNP_hash;
        while (my ($k,$v)=each %SNP_hash){print "$k\t$v\tFemale\n"};
        %SNP_hash = ();
	foreach (@Female) {
                print "$_\n";
        }

        for $n (1..($Ni/2)) { # number of individuals /2
                $first = shift (@Male);
                $second = shift (@Female);
                $pairs{$first} = $second;
#               print "$first\t$second\n";
        }
        return %pairs;

}