#######################################################################
# This software has been created by Genome Research Limited (GRL).    # 
# GRL hereby grants permission to use, copy, modify and distribute    # 
# this software and its documentation for non-commercial purposes     # 
# without fee at the user's own risk on the basis set out below.      #
# GRL neither undertakes nor accepts any duty whether contractual or  # 
# otherwise in connection with the software, its use or the use of    # 
# any derivative, and makes no representations or warranties, express #
# or implied, concerning the software, its suitability, fitness for   #
# a particular purpose or non-infringement.                           #
# In no event shall the authors of the software or GRL be responsible # 
# or liable for any loss or damage whatsoever arising in any way      # 
# directly or indirectly out of the use of this software or its       # 
# derivatives, even if advised of the possibility of such damage.     #
# Our software can be freely distributed under the conditions set out # 
# above, and must contain this copyright notice.                      #
#######################################################################
# $Id: AutoOrder2.pm,v 1.20 1997/08/04 12:33:41 rmd Exp $

# $Log: AutoOrder2.pm,v $
# Revision 1.20  1997/08/04  12:33:41  rmd
# Added new order by tree search method
#

package AutoOrder2;

###########################################################################################################
# Simulated annealing contig ordering ...
#
#
#
use IPC::Open2;

my @Data;
my @Matrix;
my %Uniques;
my %Pass;
my %ContOff;
##############################################################################################################

sub FirstPass{

    my $i =0;
    my $line =0;
    my $fs;
    my $sf;
    my %fsfs;

    foreach $pair (@Data)
    {
	my @p;
	if (($$pair[5] == 1)) # && ($$pair[8] == 0)) # Has fwd rev and not bad
	{
	    $Uniques{$$pair[13]}++;
	    $Uniques{$$pair[14]}++;
	    push(@{$Pass{$$pair[13]}},$line); 
	    push(@{$Pass{$$pair[14]}},$line); 
	}      
	$line++;	
    }
   
}

sub newAutoOrder {
    my (@Assem) = @_;
    my %notOrderedYet = map { $_, 1 } @Assem;
    @Data = &ByTempl::Get_Assem;
    &FirstPass;
    %ContOff = &ByTempl::Get_CLR;
    my @newOrder;
    my @firstset = sort keys(%Uniques);
    if ($#firstset <= 1) { return @Assem; }
    my %cmap;                                                               # Contig map
    my $first;
    foreach $first (@firstset) {                                            # Find which contigs are joined by read pairs
	my $second;
	foreach $second (@firstset) {
	    if ($first ne $second) {
		my $inCommon = &newInCommon($first, $second);
		if ($inCommon) { $cmap{$first}->{$second} = $inCommon; }
	    }
	}
    }
    my $contig;
    foreach $contig (@Assem) {
	next unless ($notOrderedYet{$contig});                             # Don't try to do a contig twice
	if ($cmap{$contig}) {                                              # See if contig is in a group
	    my $subset = &getSubset(\%cmap, $contig);                      # Find the interesting subset of the contig map
	    my @subsetKeys;
	    while (@subsetKeys = keys %{$subset}) {
		my $leaves = &findLeafNodes($subset);                          # Find a good starting point : either a leaf node,
		unless (scalar(@{$leaves})) { push(@{$leaves}, $contig); }     # or a random place if there are no leaf nodes.
		my $leaf = pop(@{$leaves});
		my $bestScore = 0;
		my @bestList;
		my %visited;
		&searchNet($subset, $leaf, 0, [], \$bestScore, \@bestList, \%visited);    # Find the longest list from the starting point.
		if (scalar(@{$leaves}) > 1) {
		    $leaf = pop(@bestList);
		    $bestScore = 0;
		    @bestList = ();
		    &searchNet($subset, $leaf, 0, [], \$bestScore, \@bestList, \%visited);  # If there were more than 2 leaf nodes, try the most
		}                                                                           # distant one, as it could give a better result
		push(@newOrder, @bestList);                                   # Push list on to the new contig order
		map { $notOrderedYet{$_} = 0; } @bestList;                    # Make sure we don't include any of the contigs on the list again.
		map { &delCmapEntry(\%cmap, $_); } @bestList;
		map { &delCmapEntry($subset, $_); } @bestList;
	    }
	} else {                                                          # If contig wasn't in a group, put it in on its own.
	    push(@newOrder, $contig);
	    $notOrderedYet{$contig} = 0;
	}
    }
    return @newOrder;
}

sub delCmapEntry {
    my ($cmap, $entry) = @_;
    my $children = $cmap->{$entry};
    map { delete $cmap->{$_}->{$entry}; } keys %{$children};
    map { unless (keys %{$cmap->{$_}}) { delete $cmap->{$_}; } } keys %{$children};
    delete $cmap->{$entry};
}

sub getSubset {
    my ($cmap, $contig) = @_;
    my %visited;
    my %subSet;
    $visited{$contig} = 1;
    my @toInclude = keys %{$cmap->{$contig}};
    $subSet{$contig} = $cmap->{$contig};
    my $nextContig;
    while ($nextContig = pop(@toInclude)) {
	unless ($visited{$nextContig}) {
	    $visited{$nextContig} = 1;
	    $subSet{$nextContig} = $cmap->{$nextContig};
	    push(@toInclude, (keys %{$cmap->{$nextContig}}));
	}
    }
    return \%subSet;
}

sub searchNet {
    my ($cmap, $contig, $score, $currentList, $bestScore, $bestList, $visited) = @_;
    push(@{$currentList}, $contig);
    $visited->{$contig} = 1;
    my @children = grep { (!$visited->{$_}) && defined($cmap->{$_}); } (keys %{$cmap->{$contig}});
    if (scalar(@children) == 0) {
	if ($$bestScore < $score) {
	    $$bestScore = $score;
	    @{$bestList} = @{$currentList};
	}
	$visited->{$contig} = 0;
	pop(@{$currentList});
	return;
    }
    my $child;
    foreach $child (@children) {
	unless ($visited->{$child}) {
	    &searchNet($cmap, $child, $score + $cmap->{$contig}->{$child}, $currentList, $bestScore, $bestList, $visited);
	}
    }
    $visited->{$contig} = 0;
    pop(@{$currentList});
    return;
}

sub findLeafNodes {
    my ($cmap) = @_;
    my @leaves = grep { scalar(keys %{$cmap->{$_}}) == 1; } (keys %{$cmap});
    return \@leaves;
}

sub newInCommon{

    my ($f,$s) = @_;
    my @both;
    my @fd = @{$Pass{$f}};
    my @sd = @{$Pass{$s}};
    my $gcount =0;

    my $fn;
    foreach $fn (@fd) {
	my $sn;
	foreach $sn (@sd) {
	    if ($fn == $sn) { push(@both,$fn); }
	}
    }
   
    if ($#both == -1) { return (0); }

    foreach (@both)
    {
	my $pair = $Data[$_];
	my $pt = &PointTo($$pair[9],$$pair[10],$$pair[11],$$pair[12]);
	if ($pt == 1) { $gcount++; }   # good stuff
    }
    return $gcount;
}

##############################################################################################################
sub AutoOrder{

    my (@Assem) = @_;

    my $f = $#Assem;
#    if($f > 999){print "Crow : Failed : Too Many Contigs !!!!\n"; return 0;}
    @Data = &ByTempl::Get_Assem;
    &FirstPass;
    %ContOff = &ByTempl::Get_CLR;
    my @NewOrder;
    
 #   print "Crow : Warning .... this feature has not been tested and may be incorrect ";
 #   &main::statusbar("Warning .... this feature has not been tested and may be incorrect");

    my @firstset = keys(%Uniques);
    if($#firstset <= 1 ) {return (@Assem);}
     
    foreach $f (@firstset)	# 
    {
	my @vals;		
	push(@vals, $f);
       
	foreach $s (@firstset)
	{
	    if($s ne $f){push(@vals,&InCommon($f,$s));}else{push(@vals,0);} 
	}
	push(@Matrix, \@vals);
    } 


    open2(NEWORDER,MATRIX,"Anneal");

 #   print "\n";        
    foreach $line (@Matrix)
    {
	foreach (@{$line}){print MATRIX "$_\t";}
	print MATRIX "\n";
#	print "\n";  
    }

    close(MATRIX);
    $i =0;  

    while ($line = <NEWORDER>)
    {
#	print "Got Back $line ";
	chop($line);

	$NewOrder[$i] = $line;
 	$i++;
    }	
    close(NEWORDER);
    my $pos =0;
#    print "\nNew Order \n";
    foreach $no (@NewOrder)
    {
#	print "$no \n";
	$have{$no} = $pos;
	$pos++;
    }

    foreach(@Assem){if(!defined($have{$_})){push(@NewOrder,$_);}}
        
    undef(@Matrix);
    undef(%Pass);
    undef(@Data);
    undef(%Uniques);

    @NewOrder;   
}

 		       
    


###############################################################################################################
# num from 0 to 3 ...
#

sub InCommon{

    my ($f,$s) = @_;
    my @both;
    my %temp;
    my @fd = @{$Pass{$f}};
    my @sd = @{$Pass{$s}};
    my $pair;
    my $gcount =0;
    my $bcount =0;
    my $total =0;
    my $pt =0;
    my $d;
    my $dist;


    foreach $fn (@fd)
    {
	foreach $sn(@sd)
	{
          push(@both,$fn) if($fn == $sn);
	}
    }
   
    return (0) if($#both == -1);
    foreach  (@both)
    {
	$pair = $Data[$_];
	
	$pt = &PointTo($$pair[9],$$pair[10],$$pair[11],$$pair[12]);
#	$d = &InDist($$pair[9],$$pair[10],$$pair[11],$$pair[12],$$pair[13],$$pair[14]);

	$gcount++ if($pt == 1);   # good stuff
#	$bcount++ if($pt == 0);  # bad stuff 	    
	
    }
    $total = $gcount;

    
    $total /= ($#both +1);
    $total = int($total *10);
    $total = 10 if($total > 10);
    $total;
}
###############################################################################################################
sub PointTo{


    my ($fs1,$fs2,$ss1,$ss2) =@_;

    # has to be near an end
  
    if((($fs1 < $fs2) && ($ss1 > $ss2)) ||(($fs1 > $fs2) && ($ss1 < $ss2))){return (1)};
    return(0);

}
###############################################################################################################
sub InDist{

 my ($fs1,$fs2,$ss1,$ss2,$lcont,$rcont) =@_;
 if ((abs($fs1 - $ContOff{$lcont}->[0]) > 2000) || (abs($fs1 - $ContOff{$lcont}->[1]) > 2000)) {return (0)};
 if ((abs($ss1 - $ContOff{$rcont}->[0]) > 2000) || (abs($fs1 - $ContOff{$rcont}->[1]) > 2000)) {return (0)};
 return (1);
}

1;



