view release on metacpan or  search on metacpan

lib/Bio/Phylo/EvolutionaryModels.pm  view on Meta::CPAN

            Bio::Phylo::Util::Exceptions::BadFormat->throw(
                'error' => "the model must be a pure birth process" );
        }

        #Get the random length to add after the last speciation event
        my $pendant_add = &$pendant_dist( %{ $options{model_options} },
            tree_size => $options{tree_size} );

        #Add the final length
        foreach ( @{ $tree->get_terminals } ) {
            $_->set_branch_length( $_->get_branch_length + $pendant_add );
        }

        #Add to the sample
        push( @sample,           $tree );
        push( @expected_summary, 1 );

 #Update the sample counter (for GUI or other applications that want to know how
 #many samples have been obtained)
        if ( defined $options{counter} ) {
            my $counter = $options{counter};
            &$counter(1);
        }
    }
    return ( \@sample, \@expected_summary );
}

=item sample_constant_rate_bd()

Sample from a constant rate birth and death model

 Type    : Sampling algorithm
 Title   : sample_constant_rate_bd
 Usage   : see sample
 Function: Samples trees from a memoryless birth model
 Returns : see sample
 Args    : no specific algorithm options but see below

NB: This algorithm only applies to constant rate birth and death 
processes. Consequently a model does not need to be specified (and
will be ignored if it is). But birth_rate and death_rate model 
options must be given. 

=cut

sub sample_constant_rate_bd {
    my %options = @_;

    #Store parameters in shorter variables (for clarity)
    my ( $br, $dr, $n ) = (
        $options{model_options}->{birth_rate},
        $options{model_options}->{death_rate},
        $options{tree_size}
    );
    my @sample;

    #Loop for sampling each tree
    while ( scalar @sample < $options{sample_size} ) {
        my @nodes;

       #Compute the random tree age from the inverse CDF (different formulas for
       #birth rate == death rate and otherwise)
        my $tree_age;

        #The uniform random variable
        my $r = rand;
        if ( $br == $dr ) {
            $tree_age = 1 / ( $br * ( $r**( -1 / $n ) - 1 ) );
        }
        else {
            $tree_age =
              1 /
              ( $br - $dr ) *
              log(
                ( 1 - $dr / $br * $r**( 1 / $n ) ) / ( 1 - $r**( 1 / $n ) ) );
        }

        #Find the random speciation times
        my @speciation;
        foreach ( 0 .. ( $n - 2 ) ) {
            if ( $br == $dr ) {
                my $r = rand;
                $speciation[$_] =
                  $r * $tree_age / ( 1 + $br * $tree_age * ( 1 - $r ) );
            }
            else {

                #Two repeated parts of the inverse CDF for clarity
                my $a = $br - $dr * exp( ( $dr - $br ) * $tree_age );
                my $b = ( 1 - exp( ( $dr - $br ) * $tree_age ) ) * rand;

                #The random speciation time from the inverse CDF
                $speciation[$_] =
                  1 /
                  ( $br - $dr ) *
                  log( ( $a - $dr * $b ) / ( $a - $br * $b ) );
            }
        }

        #Create the initial terminals and a vector for their ages
        my @terminals;
        my @ages;
        foreach ( 0 .. ( $n - 1 ) ) {

            #Add a new terminal
            $terminals[$_] = Bio::Phylo::Forest::Node->new();
            $terminals[$_]->set_name( 'ID' . $_ );
            $ages[$_] = 0;
        }
        @nodes = @terminals;

        #Sort the speciation times
        my @sorted_speciation = sort { $a <=> $b } @speciation;

        #Make a hash for easily finding the index of a given speciation event
        my %speciation_hash;
        foreach ( 0 .. ( $n - 2 ) ) {
            $speciation_hash{ $speciation[$_] } = $_;
        }

        #Construct the tree