Bio/PhyloNetwork/Factory.pm

   1 #
   2 # Module for Bio::PhyloNetwork::Factory
   3 #
   4 # Please direct questions and support issues to <bioperl-l@bioperl.org>
   5 #
   6 # Cared for by Gabriel Cardona <gabriel(dot)cardona(at)uib(dot)es>
   7 #
   8 # Copyright Gabriel Cardona
   9 #
  10 # You may distribute this module under the same terms as perl itself
  11
  12 # POD documentation - main docs before the code
  13
  14 =head1 NAME
  15
  16 Bio::PhyloNetwork::Factory - Module to sequentially generate
  17 Phylogenetic Networks
  18
  19 =head1 SYNOPSIS
  20
  21  use strict;
  22  use warnings;
  23
  24  use Bio::PhyloNetwork;
  25  use Bio::PhyloNetwork::Factory;
  26
  27  # Will generate sequentially all the 4059 binary tree-child phylogenetic
  28  # networks with 4 leaves
  29
  30  my $factory=Bio::PhyloNetwork::Factory->new(-numleaves=>4);
  31
  32  my @nets;
  33
  34  while (my $net=$factory->next_network()) {
  35    push @nets,$net;
  36    print "".(scalar @nets).": ".$net->eNewick()."\n";
  37  }
  38
  39 =head1 DESCRIPTION
  40
  41 Sequentially builds a (binary tree-child) phylogenetic network each time
  42 next_network is called.
  43
  44 =head1 AUTHOR
  45
  46 Gabriel Cardona, gabriel(dot)cardona(at)uib(dot)es
  47
  48 =head1 SEE ALSO
  49
  50 L<Bio::PhyloNetwork>
  51
  52 =head1 APPENDIX
  53
  54 The rest of the documentation details each of the object methods.
  55
  56 =cut
  57
  58 package Bio::PhyloNetwork::Factory;
  59
  60 use strict;
  61 use warnings;
  62
  63 use base qw(Bio::Root::Root);
  64
  65 use Bio::PhyloNetwork;
  66 use Bio::PhyloNetwork::TreeFactory;
  67
  68 =head2 new
  69
  70  Title   : new
  71  Usage   : my $factory = new Bio::PhyloNetwork::Factory();
  72  Function: Creates a new Bio::PhyloNetwork::Factory
  73  Returns : Bio::PhyloNetwork::RandomFactory
  74  Args    : -numleaves => integer
  75             OR
  76            -leaves => reference to an array (of leaves names)
  77            -numhybrids => integer [default = numleaves -1]
  78            -recurse => boolean [optional]
  79
  80 Returns a Bio::PhyloNetwork::Factory object. Such an object will
  81 sequentially create binary tree-child phylogenetic networks
  82 each time next_network is called.
  83
  84 If the parameter -leaves=E<gt>\@leaves is given, then the set of leaves of
  85 these networks will be @leaves. If it is given the parameter
  86 -numleaves=E<gt>$numleaves, then the set of leaves will be "l1"..."l$numleaves".
  87
  88 If the parameter -numhybrids=E<gt>$numhybrids is given, then the generated
  89 networks will have at most $numhybrids hybrid nodes. Note that, necessarily,
  90 $numhybrids E<lt> $numleaves.
  91
  92 If the parameter -recurse=E<gt>1 is given, then all networks with number of hybrid
  93 nodes less or equal to $numhybrids will be given; otherwise only those with
  94 exactly $numhybrids hybrid nodes.
  95
  96 =cut
  97
  98 sub new {
  99   my ($pkg,@args)=@_;
 100
 101   my $self=$pkg->SUPER::new(@args);
 102
 103   my ($leavesR,$numleaves,$numhybrids,$recurse)=
 104     $self->_rearrange([qw(LEAVES
 105                           NUMLEAVES
 106                           NUMHYBRIDS
 107                           RECURSE)],@args);
 108   my @leaves;
 109   if ((! defined $leavesR) && (defined $numleaves)) {
 110     @leaves=map {"l$_"} (1..$numleaves);
 111     $leavesR=\@leaves;
 112   }
 113   if (! defined $leavesR) {
 114     $self->throw("No leaves set neither numleaves given");
 115   }
 116   @leaves=@$leavesR;
 117   $self->{leaves}=$leavesR;
 118   $numleaves=@leaves;
 119   $self->{numleaves}=$numleaves;
 120
 121   $recurse ||= 0;
 122   if (! defined $numhybrids) {
 123     $numhybrids=$numleaves-1;
 124     $recurse=1;
 125   }
 126   $self->{recurse}=$recurse;
 127   $self->{numhybrids}=$numhybrids;
 128   if ($numhybrids ==0) {
 129     return Bio::PhyloNetwork::TreeFactory->new(-leaves=>\@leaves);
 130   }
 131   my $parent;
 132   if ($numhybrids > 1) {
 133     $parent=new($pkg,'-leaves'=>\@leaves,
 134                 '-numhybrids'=>($numhybrids-1),
 135                 '-recurse'=>($recurse));
 136   }
 137   else {
 138     $parent=Bio::PhyloNetwork::TreeFactory->new(-leaves=>\@leaves);
 139   }
 140   $self->{parent}=$parent;
 141   my $oldnet=$parent->next_network();
 142   $self->{oldnet}=$oldnet;
 143   $self->update();
 144   $self->{found}=[];
 145   bless($self,$pkg);
 146 }
 147
 148 sub update {
 149   my ($self)=@_;
 150
 151   my @candidates=$self->{oldnet}->edges();
 152   $self->{candidates}=\@candidates;
 153   $self->{numcandidates}=(scalar @candidates);
 154   $self->{index1}=-$self->{recurse};
 155   $self->{index2}=0;
 156 }
 157
 158 =head2 next_network
 159
 160  Title   : next_network
 161  Usage   : my $net=$factory->next_network()
 162  Function: returns a network
 163  Returns : Bio::PhyloNetwork
 164  Args    : none
 165
 166 =cut
 167
 168 sub next_network {
 169   my ($self)=@_;
 170   my $numleaves=$self->{numleaves};
 171   my $numhybrids=$self->{numhybrids};
 172   START:
 173   if ($self->{index1}==-1) {
 174     $self->{index1}++;
 175     return $self->{oldnet};
 176   }
 177   if ($self->{index1} >= $self->{numcandidates}) {
 178     $self->{index2}++;
 179     $self->{index1}=0;
 180   }
 181   if ($self->{index2} >= $self->{numcandidates}) {
 182     my $oldnet=$self->{parent}->next_network();
 183     if (! $oldnet) {
 184       return 0;
 185     }
 186     $self->{oldnet}=$oldnet;
 187     $self->update();
 188     goto START;
 189   }
 190   if ((scalar $self->{oldnet}->hybrid_nodes())< $self->{numhybrids}-1) {
 191     $self->{candidates}=[];
 192     $self->{numcandidates}=0;
 193     goto START;
 194   }
 195   my $u1=$self->{candidates}->[$self->{index1}]->[0];
 196   my $v1=$self->{candidates}->[$self->{index1}]->[1];
 197   my $u2=$self->{candidates}->[$self->{index2}]->[0];
 198   my $v2=$self->{candidates}->[$self->{index2}]->[1];
 199   my $lbl=$self->{numhybrids};
 200   if ($self->{oldnet}->is_attackable($u1,$v1,$u2,$v2)) {
 201     my $net=Bio::PhyloNetwork->new(-graph=>$self->{oldnet}->graph);
 202     $net->do_attack($u1,$v1,$u2,$v2,$lbl);
 203     $self->{index1}++;
 204     my @found=@{$self->{found}};
 205     foreach my $netant (@found) {
 206       if ($net->is_mu_isomorphic($netant) ) {
 207         goto START;
 208       }
 209     }
 210     push @found,$net;
 211     $self->{found}=\@found;
 212     return $net;
 213   }
 214   else {
 215     $self->{index1}++;
 216     goto START;
 217   }
 218 }
 219
 220 1;