File Coverage

blib/lib/Bio/Phylo/EvolutionaryModels.pm

Criterion	Covered	Total	%
statement	396	819	48.3
branch	93	258	36.0
condition	38	143	26.5
subroutine	34	46	73.9
pod	13	21	61.9
total	574	1287	44.6

package Bio::Phylo::EvolutionaryModels;
use strict;
use base 'Exporter';
use Bio::Phylo::Forest::Tree;
use Bio::Phylo::Forest;
use Math::CDF qw'qnorm qbeta';
use List::Util qw'min max';
use POSIX qw'ceil floor';
use Config;    #Use to check whether multi-threading is available

BEGIN {

    use Bio::Phylo;
    our @EXPORT_OK =
      qw(&sample &constant_rate_birth &constant_rate_birth_death &clade_shifts);
}

=head1 NAME

Bio::Phylo::EvolutionaryModels - Evolutionary models for phylogenetic trees and methods to sample these
Klaas Hartmann, September 2007

=head1 SYNOPSIS

 #For convenience we import the sample routine (so we can write sample(...) instead of
 #Bio::Phylo::EvolutionaryModels::sample(...).
 use Bio::Phylo::EvolutionaryModels qw (sample);

 #Example#A######################################################################
 #Simulate a single tree with ten species from the constant rate birth model with parameter 0.5
 my $tree = Bio::Phylo::EvolutionaryModels::constant_rate_birth(birth_rate => .5, tree_size => 10);

 #Example#B######################################################################
 #Sample 5 trees with ten species from the constant rate birth model using the b algorithm
 my ($sample,$stats) = sample(sample_size =>5,
                              tree_size => 10,
                              algorithm => 'b',
                              algorithm_options => {rate => 1},
                              model => \&Bio::Phylo::EvolutionaryModels::constant_rate_birth,
                              model_options => {birth_rate=>.5});


 #Print a newick string for the 4th sampled tree
 print $sample->[3]->to_newick."\n";

 #Example#C######################################################################
 #Sample 5 trees with ten species from the constant rate birth and death model using
 #the bd algorithm and two threads (useful for dual core processors)
 #NB: we must specify an nstar here, an appropriate choice will depend on the birth_rate
 #    and death_rate we are giving the model

 my ($sample,$stats) = sample(sample_size =>5,
                              tree_size => 10,
                              threads => 2,
                              algorithm => 'bd',
                              algorithm_options => {rate => 1, nstar => 30},
                              model => \&Bio::Phylo::EvolutionaryModels::constant_rate_birth_death,
                              model_options => {birth_rate=>1,death_rate=>.8});

 #Example#D######################################################################
 #Sample 5 trees with ten species from the constant rate birth and death model using
 #incomplete taxon sampling
 #
 #sampling_probability is set so that the true tree has 10 species with 50% probability,
 #11 species with 30% probability and 12 species with 20% probability
 #
 #NB: we must specify an mstar here this will depend on the model parameters and the
 #    incomplete taxon sampling parameters

 my $algorithm_options = {rate => 1,
                          nstar => 30,
                          mstar => 12,
                          sampling_probability => [.5, .3, .2]};

 my ($sample,$stats) = sample(sample_size =>5,
                              tree_size => 10,
                              algorithm => 'incomplete_sampling_bd',
                              algorithm_options => $algorithm_options,
                              model => \&Bio::Phylo::EvolutionaryModels::constant_rate_birth_death,
                              model_options => {birth_rate=>1,death_rate=>.8});

 #Example#E######################################################################
 #Sample 5 trees with ten species from a Yule model using the memoryless_b algorithm

 #First we define the random function for the shortest pendant edge for a Yule model
 my $random_pendant_function = sub {
     %options = @_;
     return -log(rand)/$options{birth_rate}/$options{tree_size};
 };

 #Then we produce our sample
 my ($sample,$stats) = sample(sample_size =>5,
                              tree_size => 10,
                              algorithm => 'memoryless_b',
                              algorithm_options => {pendant_dist => $random_pendant_function},
                              model => \&Bio::Phylo::EvolutionaryModels::constant_rate_birth,
                              model_options => {birth_rate=>1});

 #Example#F#######################################################################
 #Sample 5 trees with ten species from a constant birth death rate model using the
 #constant_rate_bd algorithm
 my ($sample) = sample(sample_size => 5,
                       tree_size => 10,
                       algorithm => 'constant_rate_bd',
                       model_options => {birth_rate=>1,death_rate=>.8});

=head1 DESCRIPTION

This package contains evolutionary models for phylogenetic trees and
algorithms for sampling from these models. It is a non-OO module that
optionally exports the 'sample', 'constant_rate_birth' and
'constant_rate_birth_death' subroutines into the caller's namespace,
using the C<< use Bio::Phylo::EvolutionaryModels qw(sample constant_rate_birth constant_rate_birth_death); >>
directive. Alternatively, you can call the subroutines as class methods,
as in the synopsis.

The initial set of algorithms available in this package corresponds to those in:

Sampling trees from evolutionary models
Klaas Hartmann, Dennis Wong, Tanja Gernhard
Systematic Biology, in press

Some comments and code refers back to this paper.
Further algorithms and evolutionary are encouraged
and welcome.

To make this code as straightforward as possible to read some of the
algorithms have been implemented in a less than optimal manner. The code
also follows the structure of an earlier version of the manuscript so
there is some redundancy (eg. the birth algorithm is just a specific
instance of the birth_death algorithm)

=head1 SAMPLING

All sampling algorithms should be accessed through the generic sample
interface.

=head2 Generic sampling interface: sample()

 Type    : Interface
 Title   : sample
 Usage   : see SYNOPSIS
 Function: Samples phylogenetic trees from an evolutionary model
 Returns : A sample of phylogenetic trees and statistics from the
           sampling algorithm
 Args    : Sampling parameters in a hash

This method acts as a gateway to the various sampling algorithms. The
argument is a single hash containing the options for the sampling run.


Sampling parameters (* denotes optional parameters):

 sample_size    The number of trees to return (more trees may be returned)
 tree_size      The size that returned trees should be
 model          The evolutionary model (should be a function reference)
 model_options  A hash pointer for model options (see individual models)
 algorithm      The algorithm to use (omit the preceding sample_)
 algorithm_options A hash pointer for options for the algorithm (see individual algorithms for details)
 threads*       The number of threads to use (default is 1)
 output_format* Set to newick for newick trees (default is Bio::Phylo::Forest::Tree)
 remove_extinct Set to true to remove extinct species



Available algorithms (algorithm names in the paper are given in brackets):

 b                       For all pure birth models (simplified GSA)
 bd                      For all birth and death models (GSA)
 incomplete_sampling_bd  As above, with incomplete taxon sampling (extended GSA)
 memoryless_b            For memoryless pure birth models (PBMSA)
 constant_rate_bd        For birth and death models with constant rates (BDSA)


Model

If you create your own model it must accept an options hash as its input.
This options hash can contain any parameters you desire. Your model should
simulate a tree until it becomes extinct or the size/age limit as specified
in the options has been reached. Respectively these options are tree_size
and tree_age.


Multi-threading

Multi-thread support is very simplistic. The number of threads you specify
are created and each is assigned the task of finding sample_size/threads
samples. I had problems with using Bio::Phylo::Forest::Tree in a multi-
threaded setting. Hence the sampled trees are returned as newick strings to
the main routine where (if required) Tree objects are recreated from the
strings. For most applications this overhead seems negligible in contrast
to the sampling times.


From a code perspective this function (sample):

 Checks input arguments
 Handles multi-threading
 Calls the individual algorithms to perform sampling
 Reformats data

=cut

my @threads;

sub sample {
    my %options         = @_;
    my %methods_require = (
        b  => ['rate'],
        bd => [ 'rate', 'nstar' ],
        incomplete_sampling_bd =>
          [ 'rate', 'nstar', 'mstar', 'sampling_probability' ],
        memoryless_b     => ['pendant_dist'],
        constant_rate_bd => [],
    );

    #Default is to sample a single tree
    $options{sample_size} = 1 unless defined $options{sample_size};

    #Default is to use a single thread
    $options{threads} = 1 unless defined $options{threads};

    #Check that multiple threads are actually supported
    if ( $options{threads} > 1 && !$Config{useithreads} ) {
        Bio::Phylo::Util::Exceptions::BadArgs->throw( 'error' =>
              "your perl installation does not support multiple threads" );
    }

    #Check an algorithm was specified
    unless ( defined $options{algorithm} ) {
        Bio::Phylo::Util::Exceptions::BadArgs->throw(
            'error' => "an algorithm type must be specified" );
    }

    #Check the algorithm type is valid
    unless ( defined $methods_require{ $options{algorithm} } ) {
        Bio::Phylo::Util::Exceptions::BadFormat->throw(
            'error' => "'$options{algorithm}' is not a valid algorithm" );
    }

    #Check the algorithm options
    foreach ( @{ $methods_require{ $options{algorithm} } } ) {
        unless ( defined $options{algorithm_options}->{$_} ) {
            Bio::Phylo::Util::Exceptions::BadArgs->throw( 'error' =>
"'$_' must be specified for the '$options{algorithm}' algorithm"
            );
        }
    }

#If we are doing incomplete taxon sampling the sampling probability must be specified
    if (   defined $options{incomplete_sampling}
        && $options{incomplete_sampling}
        && !( defined $options{algorithm_options}->{sampling_probability} ) )
    {
        Bio::Phylo::Util::Exceptions::BadArgs->throw( 'error' =>
"'sampling_probability' must be specified for post hoc incomplete sampling to be applied"
        );
    }

    #Check that a model has been specified
    unless ( defined $options{model}
        || $options{algorithm} eq 'constant_rate_bd' )
    {
        Bio::Phylo::Util::Exceptions::BadArgs->throw(
            'error' => "a model must be specified" );
    }

    #Get a function pointer for the algorithm
    my $algorithm = 'sample_' . $options{algorithm};
    $algorithm = \&$algorithm;
    my @output;

    #Run the algorithm, different method for multiple threads
    if ( $options{threads} > 1 ) {
        require threads;
        require threads::shared;
        @output = ( [], [] );
        $SIG{'KILL'} = sub {
            foreach (@threads) { $_->kill('KILL')->detach(); }
            threads->exit();
        };

        #Start the threads
        for ( ( 1 .. $options{threads} ) ) {
            @_ = ();

            #Note the list context of the return argument here determines
            #the data type returned from the thread.
            ( $threads[ $_ - 1 ] ) = threads->new( \&_sample_newick, %options,
                sample_size => ceil( $options{sample_size} / $options{threads} )
            );
        }

        #Wait for them to finish and combine the data
        for ( ( 1 .. $options{threads} ) ) {
            until ( $threads[ $_ - 1 ]->is_joinable() ) { sleep(0.1); }
            my @thread_data = $threads[ $_ - 1 ]->join;
            for ( my $index = 0 ; $index < scalar @thread_data ; $index++ ) {
                $output[$index] = [] if scalar @output < $index;
                $output[$index] =
                  [ @{ $output[$index] }, @{ $thread_data[$index] } ];
            }
        }

        #Turn newick strings back into tree objects
        unless ( defined $options{output_format}
            && $options{output_format} eq 'newick' )
        {

            #Convert to newick trees
            for ( my $index = 0 ; $index < scalar @{ $output[0] } ; $index++ ) {
                $output[0]->[$index] = Bio::Phylo::IO->parse(
                    -format => 'newick',
                    -string => $output[0]->[$index]
                )->first;
            }
        }
    }
    else {

        #Get the samples
        @output = &$algorithm(%options);

        #Turn into newick trees if requested
        if ( defined $options{output_format}
            && $options{output_format} eq 'newick' )
        {

            #Convert to newick trees
            for ( my $index = 0 ; $index < scalar @{ $output[0] } ; $index++ ) {
                $output[0]->[$index] =
                  $output[0]->[$index]->to_newick( '-nodelabels' => 1 );
            }
        }
        elsif ( defined $options{output_format}
            && $options{output_format} eq 'forest' )
        {

            # save as a forest
            my $forest = Bio::Phylo::Forest->new;
            for ( my $index = 0 ; $index < scalar @{ $output[0] } ; $index++ ) {
                $forest->insert( $output[0]->[$index] );
            }
            $output[0] = $forest;
        }
    }
    return @output;
}

=begin comment

 Type    : Internal method
 Title   : _sample_newick
 Usage   : ($thread) = threads->new(\&_sample_newick, %options);
           @thread_output = $thread->join;
 Function: Wrapper for sampling routines used for multi-threading
 Returns : Output from sampling algorithms with trees replaced by newick strings
 Args    : %options to pass to sampling algorithm

=end comment

=cut

sub _sample_newick {
    my %options   = @_;
    my $algorithm = 'sample_' . $options{algorithm};

    # Thread 'cancellation' signal handler
    $SIG{'KILL'} = sub { threads->exit(); };
    $algorithm = \&$algorithm;

    #Perform the sampling
    my @output = ( &$algorithm(%options) );

    #Convert to newick trees
    for ( my $index = 0 ; $index < scalar @{ $output[0] } ; $index++ ) {
        $output[0]->[$index] =
          $output[0]->[$index]->to_newick( '-nodelabels' => 1 );
    }
    return @output;
}

=head2 Sampling algorithms

These algorithms should be accessed through the sampling interface (sample()).
Additional parameters need to be passed to these algorithms as described for
each algorithm.

=over

=item sample_b()

Sample from any birth model

 Type    : Sampling algorithm
 Title   : sample_b
 Usage   : see sample
 Function: Samples trees from a pure birth model
 Returns : see sample
 Args    : %algorithm_options requires the field:
           rate => sampling rate

=cut

sub sample_b {
    my %options = @_;

    #The sample of trees
    my @sample;

    #A list of the expected number of samples
    my @expected_summary;
    my $model = $options{model};
    my $rate  = $options{algorithm_options}->{rate};

    #While we have insufficient samples
    while ( scalar @sample < $options{sample_size} ) {

        #Generate a candidate model run
        my $candidate = &$model( %{ $options{model_options} },
            tree_size => $options{tree_size} );

        #Check that the tree has no extinctions
        unless ( $candidate->is_ultrametric(1e-6) ) {
            Bio::Phylo::Util::Exceptions::BadFormat->throw(
                'error' => "the model must be a pure birth process" );
        }

        #Get the lineage through time data
        my ( $time, $count ) = lineage_through_time($candidate);

        #The expected number of samples we want
        my $expected_samples = $rate * ( $time->[-1] - $time->[-2] );
        push( @expected_summary, $expected_samples );

        #Get the random number of samples from this candidate tree
        while ( $expected_samples > 0 ) {

            #If the number of samples remaining is greater than one
            #we take a sample, otherwise we take a sample with probability
            #according to the remaining number of samples
            if ( $expected_samples > 1 || rand(1) < $expected_samples ) {
                my $tree = copy_tree($candidate);

      #Truncate the tree at a random point during which it had tree_size species
                truncate_tree_time( $tree,
                    ( $time->[-1] - $time->[-2] ) * rand(1) + $time->[-2] );

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

 #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);
                }
            }
            $expected_samples--;
        }
    }
    return ( \@sample, \@expected_summary );
}

=item sample_bd()

Sample from any birth and death model for which nstar exists

 Type    : Sampling algorithm
 Title   : sample_bd
 Usage   : see sample
 Function: Samples trees from a birth and death model
 Returns : see sample
 Args    : %algorithm_options requires the fields:
           nstar => once a tree has nstar species there should be
           a negligible chance of returning to tree_size species
           rate => sampling rate

=cut

sub sample_bd {
    my %options = @_;

    #The sample of trees
    my @sample;

    #A list of the expected number of samples
    my @expected_summary;

    #Convenience variables
    my $model = $options{model};
    my $nstar = $options{algorithm_options}->{nstar};
    my $rate  = $options{algorithm_options}->{rate};

    #While we have insufficient samples
    while ( scalar @sample < $options{sample_size} ) {

        #Generate a candidate model run
        my $candidate =
          &$model( %{ $options{model_options} }, tree_size => $nstar );

        #Get the lineage through time data
        my ( $time, $count ) = lineage_through_time($candidate);

#Reorganise the lineage through time data
#@duration contains the length of the intervals with the right number of species
#@start contains the starting times of these intervals
#@prob contains the cumulative probability of each interval being selected
#$total_duration contains the sum of the interval lengths
        my ( @duration, @start, @prob, $total_duration );
        for ( my $index = 0 ; $index < scalar @{$time} - 1 ; $index++ ) {
            if ( $count->[$index] == $options{tree_size} ) {
                push( @duration, $time->[ $index + 1 ] - $time->[$index] );
                push( @start,    $time->[$index] );
                $total_duration += $duration[-1];
                push( @prob, $total_duration );
            }
        }
        no warnings 'uninitialized';    # FIXME
        next if $total_duration == 0;
        use warnings;
        for ( my $index = 0 ; $index < scalar @prob ; $index++ ) {
            $prob[$index] /= $total_duration;
        }

        #The expected number of samples we want
        my $expected_samples = $rate * $total_duration;
        push( @expected_summary, $expected_samples );

        #Get the random number of samples from this candidate tree
        while ( $expected_samples > 0 ) {

            #If the number of samples remaining is greater than one
            #we take a sample, otherwise we take a sample with probability
            #according to the remaining number of samples
            if ( $expected_samples > 1 || rand(1) < $expected_samples ) {
                my $tree = copy_tree($candidate);

                #Get a random interval
                my $interval_choice = rand(1);
                my $interval;
                for (
                    $interval = 0 ;
                    $interval_choice > $prob[$interval] ;
                    $interval++
                  )
                {
                }

                #Truncate the tree at a random point during this interval
                truncate_tree_time( $tree,
                    $duration[$interval] * rand(1) + $start[$interval] );

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

 #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);
                }
            }
            $expected_samples--;
        }
    }
    return ( \@sample, \@expected_summary );
}

=item sample_incomplete_sampling_bd()

Sample from any birth and death model with incomplete taxon sampling

 Type    : Sampling algorithm
 Title   : sample_incomplete_sampling_bd
 Usage   : see sample
 Function: Samples trees from a birth and death model with incomplete taxon sampling
 Returns : see sample
 Args    : %algorithm_options requires the fields:
           rate => sampling rate
           nstar => once a tree has nstar species there should be
           a negligible chance of returning to mstar species
           mstar => trees with more than mstar species form a negligible
           contribution to the final sample.
           sampling_probability => see below.

sampling_probability

 vector: must have length (mstar-tree_size+1) The ith element gives the probability
         of not sampling i species.
 scalar: the probability of sampling any individual species. Is used to calculate
         a vector as discussed in the paper.

=cut

sub sample_incomplete_sampling_bd {
    my %options = @_;

    #    %options = (%options, %{$options{algorithm_options}});
    #The sample of trees
    my @sample;

    #A list of the expected number of samples
    my @expected_summary;

    #     #Convenience variables
    my $model = $options{model};
    my $nstar = $options{algorithm_options}->{nstar};
    my $mstar = $options{algorithm_options}->{mstar};
    my $rate  = $options{algorithm_options}->{rate};
    my $sampling_probability =
      $options{algorithm_options}->{sampling_probability};

    #If sampling_probability is a list check it's length
    if ( ref $sampling_probability
        && scalar @{$sampling_probability} !=
        ( $mstar - $options{tree_size} + 1 ) )
    {
        Bio::Phylo::Util::Exceptions::BadArgs->throw( 'error' =>
"'sampling_probability' must be a scalar or a list with m_star-tree_size+1 items"
        );
    }

    #If a single sampling probability was given we calculate the
    #probability of sub-sampling given numbers of species.
    unless ( ref $sampling_probability ) {
        my $p = $sampling_probability;
        my @vec;
        my $total = 0;
        foreach ( ( $options{tree_size} .. $mstar ) ) {

#The probability of sampling tree_size species from a tree containing $_ species
            push( @vec,
                nchoosek( $_, $options{tree_size} ) *
                  ( $p**$options{tree_size} ) *
                  ( ( 1 - $p )**( $_ - $options{tree_size} ) ) );
            $total += $vec[-1];
        }
        for ( my $ii = 0 ; $ii < scalar @vec ; $ii++ ) { $vec[$ii] /= $total; }
        $sampling_probability = \@vec;
    }

    #We now normalise the sampling_probability list so that it sums to unity
    #this allows comparable sampling rates to be used here and in sample_bd
    my $total_sp = 0;
    foreach ( @{$sampling_probability} ) { $total_sp += $_; }
    no warnings;    # FIXME
    foreach ( my $index = 1 .. scalar @{$sampling_probability} ) {
        no warnings;    # FIXME
        $sampling_probability->[ $index - 1 ] /= $total_sp;
        use warnings;
    }
    use warnings;

    #While we have insufficient samples
    while ( scalar @sample < $options{sample_size} ) {

        #Generate a candidate model run
        my $candidate =
          &$model( %{ $options{model_options} }, tree_size => $nstar );

        #Get the lineage through time data
        my ( $time, $count ) = lineage_through_time($candidate);
        my @size_stats;
        for ( my $index = 0 ; $index < scalar @{$time} - 1 ; $index++ ) {
            no warnings 'uninitialized';    # FIXME
            if ( $count->[$index] >= $options{tree_size} ) {
                $size_stats[ $count->[$index] - $options{tree_size} ] +=
                  $time->[$index] *
                  $sampling_probability->[ $count->[$index] -
                  $options{tree_size} ];
            }
        }

    #Reorganise the lineage through time data
    #@duration contains the length of intervals with more than tree_size species
    #@start contains the starting times of these intervals
    #@prob contains the cumulative probability of each interval being selected
    #$total_prob contains the sum of @prob before normalisation
        my ( @duration, @start, @prob, $total_prob );
        $total_prob = 0;
        for ( my $index = 0 ; $index < scalar @{$time} - 1 ; $index++ ) {
            if ( $count->[$index] >= $options{tree_size} ) {
                push( @duration, $time->[ $index + 1 ] - $time->[$index] );
                push( @start,    $time->[$index] );
                no warnings 'uninitialized';    # FIXME
                $total_prob +=
                  $duration[-1] *
                  $sampling_probability->[ $count->[$index] -
                  $options{tree_size} ];
                push( @prob, $total_prob );
            }
        }
        next if $total_prob == 0;
        for ( my $index = 0 ; $index < scalar @prob ; $index++ ) {
            $prob[$index] /= $total_prob;
        }

        #The expected number of samples we want
        my $expected_samples = $rate * $total_prob;
        push( @expected_summary, $expected_samples );
        $expected_samples = $options{sample_size} - scalar @sample
          if $expected_samples > $options{sample_size} - scalar @sample;

        #Get the random number of samples from this candidate tree
        while ( $expected_samples > 0 ) {

            #If the number of samples remaining is greater than one
            #we take a sample, otherwise we take a sample with probability
            #according to the remaining number of samples
            if ( $expected_samples > 1 || rand(1) < $expected_samples ) {
                my $tree = copy_tree($candidate);

                #Get a random interval
                my $interval_choice = rand(1);
                my $interval;
                for (
                    $interval = 0 ;
                    $interval_choice > $prob[$interval] ;
                    $interval++
                  )
                {
                }

                #Truncate the tree at a random point during this interval
                truncate_tree_time( $tree,
                    $duration[$interval] * rand(1) + $start[$interval] );

                #Remove random species so it has the right size
                truncate_tree_size( $tree, $options{tree_size} );

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

 #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);
                }
            }
            $expected_samples--;
        }
    }
    return ( \@sample, \@expected_summary );
}

=item sample_memoryless_b()

Sample from a memoryless birth model

 Type    : Sampling algorithm
 Title   : sample_memoryless_b
 Usage   : see sample
 Function: Samples trees from a memoryless birth model
 Returns : see sample
 Args    : %algorithm_options with fields:
           pendant_dist => function reference for generating random
           shortest pendant edges

NB: The function pointed to by pendant_dist is given model_options
as it's input argument with an added field tree_size. It must return
a random value from the probability density for the shortest pendant
edges.

=cut

sub sample_memoryless_b {
    my %options = @_;

    #The sample of trees
    my @sample;

    #A list of the expected number of samples
    my @expected_summary;

    #The user specified functions
    my $model        = $options{model};
    my $pendant_dist = $options{algorithm_options}->{pendant_dist};

    #While we have insufficient samples
    while ( scalar @sample < $options{sample_size} ) {

        #Generate a tree ending just after the last speciation event
        my $tree = &$model( %{ $options{model_options} },
            tree_size => $options{tree_size} );

        #Check that the tree has no extinctions
        unless ( $tree->is_ultrametric(1e-6) ) {
            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
        foreach my $index ( 0 .. ( $n - 2 ) ) {

            #Create the parent node
            my $parent = Bio::Phylo::Forest::Node->new();
            $parent->set_name( 'ID' . ( $n + $index ) );
            push( @nodes, $parent );

            #An index for this speciation event back into the unsorted vectors
            my $spec_index = $speciation_hash{ $sorted_speciation[$index] };

            #Add the children to the parent node
            $parent->set_child( $terminals[$spec_index] );
            $terminals[$spec_index]->set_parent($parent);
            $parent->set_child( $terminals[ $spec_index + 1 ] );
            $terminals[ $spec_index + 1 ]->set_parent($parent);

            #Set the children's branch lengths
            $terminals[$spec_index]->set_branch_length(
                $sorted_speciation[$index] - $ages[$spec_index] );
            $terminals[ $spec_index + 1 ]->set_branch_length(
                $sorted_speciation[$index] - $ages[ $spec_index + 1 ] );

            #Replace the two terminals with the new one
            splice( @terminals, $spec_index, 2, $parent );
            splice( @ages,      $spec_index, 2, $sorted_speciation[$index] );

            #Update the mapping for the sorted speciation times
            foreach ( keys %speciation_hash ) {
                $speciation_hash{$_}-- if $speciation_hash{$_} > $spec_index;
            }
        }

        #Add the nodes to a tree
        my $tree = Bio::Phylo::Forest::Tree->new();
        foreach ( reverse(@nodes) ) { $tree->insert($_); }
        push( @sample, $tree );

 #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, [] );
}

=back

=head1 EVOLUTIONARY MODELS

All evolutionary models take a options hash as their input argument
and return a Bio::Phylo::Forest::Tree. This tree may contain extinct
lineages (lineages that end prior to the end of the tree).

The options hash contains any model specific parameters (see the
individual model descriptions) and one or both terminating conditions:
tree_size => the number of extant species at which to terminate the tree
tree_age => the age of the tree at which to terminate the process

Note that if the model stops due to the tree_size condition then the
tree ends immediately after the speciation event that created the last
species.

=over

=item constant_rate_birth()

A constant rate birth model (Yule/ERM)

 Type    : Evolutionary model
 Title   : constant_rate_birth
 Usage   : $tree = constant_rate_birth(%options)
 Function: Produces a tree from the model terminating at a given size/time
 Returns : Bio::Phylo::Forest::Tree
 Args    : %options with fields:
           birth_rate The birth rate parameter (default 1)
           tree_size  The size of the tree at which to terminate
           tree_age   The age of the tree at which to terminate

 NB: At least one of tree_size and tree_age must be specified

=cut

sub constant_rate_birth {
    my %options = @_;
    $options{death_rate} = 0;
    return constant_rate_birth_death(%options);
}

=item external_model()

A dummy model that takes as input a set of newick_trees and randomly samples
these.

 Type    : Evolutionary model
 Title   : external_model
 Usage   : $tree = $external_model(%options)
 Function: Returns a random tree that was given as input
 Returns : Bio::Phylo::Forest::Tree
 Args    : %options with fields:
           trees      An array of newick strings. One of these is returned at random.

 NB: The usual parameters tree_size and tree_age will be ignored. When sampling
     using this model the trees array must contain trees adhering to the requirements
     of the sampling algorithm. This is NOT checked automatically.

=cut

sub external_model {
    my %options = @_;
    my $choice  = int( rand( scalar @{ $options{trees} } ) );

    #Pick a newick string and turn it in to a Bio::Phylo::Forest::Tree object
    my $tree = Bio::Phylo::IO->parse(
        -format => 'newick',
        -string => $options{trees}->[$choice]
    )->first;
    return $tree;
}

=item constant_rate_birth_death()

A constant rate birth and death model

 Type    : Evolutionary model
 Title   : constant_rate_birth_death
 Usage   : $tree = constant_rate_birth_death(%options)
 Function: Produces a tree from the model terminating at a given size/time
 Returns : Bio::Phylo::Forest::Tree
 Args    : %options with fields:
           birth_rate The birth rate parameter (default 1)
           death_rate The death rate parameter (default no extinction)
           tree_size  The size of the tree at which to terminate
           tree_age   The age of the tree at which to terminate

 NB: At least one of tree_size and tree_age must be specified

=cut

sub constant_rate_birth_death {
    my %options = @_;

    #Check that we have a termination condition
    unless ( defined $options{tree_size} or defined $options{tree_age} ) {

        #Error here.
        return undef;
    }

    #Set the undefined condition to infinity
    $options{tree_size} = 1e6 unless defined $options{tree_size};
    $options{tree_age}  = 1e6 unless defined $options{tree_age};

    #Set default rates
    $options{birth_rate} = 1 unless defined( $options{birth_rate} );
    delete $options{death_rate}
      if defined( $options{death_rate} ) && $options{death_rate} == 0;

    #Each node gets an ID number this tracks these
    my $node_id = 0;

    #Create a new tree with a root, start the list of terminal species
    my $tree = Bio::Phylo::Forest::Tree->new();
    my $root = Bio::Phylo::Forest::Node->new();
    $root->set_branch_length(0);
    $root->set_name( 'ID' . $node_id++ );
    $tree->insert($root);
    my @terminals = ($root);
    my ( $next_extinction, $next_speciation );
    my $time      = 0;
    my $tree_size = 1;

    #Check whether we have a non-zero root edge
    if ( defined $options{root_edge} && $options{root_edge} ) {

        #Non-zero root. We set the time to the first speciation event
        $next_speciation = -log(rand) / $options{birth_rate} / $tree_size;
    }
    else {

        #Zero root, we want a speciation event straight away
        $next_speciation = 0;
    }

    #Time of the first extinction event. If no extinction we always
    #set the extinction event after the current speciation event
    if ( defined $options{death_rate} ) {
        $next_extinction = -log(rand) / $options{death_rate} / $tree_size;
    }
    else {
        $next_extinction = $next_speciation + 1;
    }

    #While the tree has not become extinct and the termination criterion
    #has not been achieved we create new speciation and extinction events
    while ($tree_size > 0
        && $tree_size < $options{tree_size}
        && $time < $options{tree_age} )
    {

        #Add the time since the last event to all terminal species
        foreach (@terminals) {
            $_->set_branch_length(
                $_->get_branch_length + min(
                    $next_extinction, $next_speciation,
                    $options{tree_age} - $time
                )
            );
        }

        #Update the time
        $time += min( $next_extinction, $next_speciation );

        #If the tree exceeds the time limit we are done
        return $tree if ( $time > $options{tree_age} );

   #Get the species effected by this event and remove it from the terminals list
        my $effected =
          splice( @terminals, int( rand( scalar @terminals ) ), 1 );

        #If we have a speciation event we add two new species
        if ( $next_speciation < $next_extinction || !defined $next_extinction )
        {
            foreach ( 1, 2 ) {

                #Create a new species
                my $child = Bio::Phylo::Forest::Node->new();
                $child->set_name( 'ID' . $node_id++ );

                #Give it a zero edge length
                $child->set_branch_length(0);

                #Add it as a child to the speciating species
                $effected->set_child($child);

                #Add it to the tree
                $tree->insert($child);

                #Add it to the terminals list
                push( @terminals, $child );
            }
        }

        #We calculate the time that the next extinction and speciation
        #events will occur (only the earliest of these will actually
        #happen). NB: this approach is only appropriate for models where
        #speciation and extinction times are exponentially distributed.
        #Windows sometimes returns 0 values for rand...
        my ( $r1, $r2 ) = ( 0, 0 );
        $r1 = rand until $r1;
        $r2 = rand until $r2;
        $tree_size = scalar @terminals;
        return $tree unless $tree_size;
        $next_speciation = -log($r1) / $options{birth_rate} / $tree_size;
        if ( defined $options{death_rate} ) {
            $next_extinction = -log($r2) / $options{death_rate} / $tree_size;
        }
        else {
            $next_extinction = $next_speciation + 1;
        }
    }
    return $tree;
}


=item diversity_dependent_speciation()

A birth and death model with speciation rate dependent on diversity as per
Etienne et. al. 2012

 Type    : Evolutionary model
 Title   : diversity_dependent_speciation
 Usage   : $tree = diversity_dependent_speciation(%options)
 Function: Produces a tree from the model terminating at a given size/time
 Returns : Bio::Phylo::Forest::Tree
 Args    : %options with fields:
           maximal_birth_rate The maximal birth rate parameter (default 1)
           death_rate The death rate parameter (default no extinction)
           K_dash     The modified carrying capacity (no default)
           tree_size  The size of the tree at which to terminate
           tree_age   The age of the tree at which to terminate

 NB: At least one of tree_size and tree_age must be specified

Reference:
Rampal S. Etienne, Bart Haegeman, Tanja Stadler, Tracy Aze, Paul N. Pearson,
Andy Purvis and Albert B. Phillimore. "Diversity-dependence brings molecular
phylogenies closer to agreement with the fossil record"
doi: 10.1098/rspb.2011.1439

=cut

sub diversity_dependent_speciation {
    my %options = @_;

    #Check that we have a termination condition
    unless ( defined $options{tree_size} or defined $options{tree_age} ) {
        #Error here.
        return undef;
    }

    #Check that we have a carrying capacity
    unless ( defined $options{K_dash} ) {
        #Error here.
        return undef;
    }

    #Set the undefined condition to infinity
    $options{tree_size} = 1e6 unless defined $options{tree_size};
    $options{tree_age}  = 1e6 unless defined $options{tree_age};

    #Set default rates
    $options{maximal_birth_rate} = 1 unless defined( $options{maximal_birth_rate} );
    delete $options{death_rate}
      if defined( $options{death_rate} ) && $options{death_rate} == 0;

    #Each node gets an ID number this tracks these
    my $node_id = 0;

    #Create a new tree with a root, start the list of terminal species
    my $tree = Bio::Phylo::Forest::Tree->new();
    my $root = Bio::Phylo::Forest::Node->new();
    $root->set_branch_length(0);
    $root->set_name( 'ID' . $node_id++ );
    $tree->insert($root);
    my @terminals = ($root);
    my ( $next_extinction, $next_speciation );
    my $time      = 0;
    my $tree_size = 1;


    $options{birth_rate} = max(0,$options{max_birth_rate}*(1-1/$options{K_dash}));


    #Check whether we have a non-zero root edge
    if ( defined $options{root_edge} && $options{root_edge} ) {

        #Non-zero root. We set the time to the first speciation event
        $next_speciation = -log(rand) / $options{birth_rate} / $tree_size;
    }
    else {

        #Zero root, we want a speciation event straight away
        $next_speciation = 0;
    }

    #Time of the first extinction event. If no extinction we always
    #set the extinction event after the current speciation event
    if ( defined $options{death_rate} ) {
        $next_extinction = -log(rand) / $options{death_rate} / $tree_size;
    }
    else {
        $next_extinction = $next_speciation + 1;
    }

    #While the tree has not become extinct and the termination criterion
    #has not been achieved we create new speciation and extinction events
    while ($tree_size > 0
        && $tree_size < $options{tree_size}
        && $time < $options{tree_age} )
    {

        #Add the time since the last event to all terminal species
        foreach (@terminals) {
            $_->set_branch_length(
                $_->get_branch_length + min(
                    $next_extinction, $next_speciation,
                    $options{tree_age} - $time
                )
            );
        }

        #Update the time
        $time += min( $next_extinction, $next_speciation );

        #If the tree exceeds the time limit we are done
        return $tree if ( $time > $options{tree_age} );

        #Get the species effected by this event and remove it from the terminals list
        my $effected =
          splice( @terminals, int( rand( scalar @terminals ) ), 1 );

        #If we have a speciation event we add two new species
        if ( $next_speciation < $next_extinction || !defined $next_extinction )
        {
            foreach ( 1, 2 ) {

                #Create a new species
                my $child = Bio::Phylo::Forest::Node->new();
                $child->set_name( 'ID' . $node_id++ );

                #Give it a zero edge length
                $child->set_branch_length(0);

                #Add it as a child to the speciating species
                $effected->set_child($child);

                #Add it to the tree
                $tree->insert($child);

                #Add it to the terminals list
                push( @terminals, $child );
            }
        }

        #We calculate the time that the next extinction and speciation
        #events will occur (only the earliest of these will actually
        #happen). NB: this approach is only appropriate for models where
        #speciation and extinction times are exponentially distributed.
        #Windows sometimes returns 0 values for rand...
        my ( $r1, $r2 ) = ( 0, 0 );
        $r1 = rand until $r1;
        $r2 = rand until $r2;
        $tree_size = scalar @terminals;
        return $tree unless $tree_size;

        $options{birth_rate} = max(0,$options{max_birth_rate}*(1-$tree_size/$options{K_dash}));
        if ($options{birth_rate}==0)
        {
        	$next_speciation = $options{tree_age};
        } else
        {
        	$next_speciation = -log($r1) / $options{birth_rate} / $tree_size;
        }
        if ( defined $options{death_rate} ) {
            $next_extinction = -log($r2) / $options{death_rate} / $tree_size;
        }
        else {
            $next_extinction = $next_speciation + 1;
        }
    }
    return $tree;
}


=item constant_rate_birth_death()

A temporal shift birth and death model

 Type    : Evolutionary model
 Title   : temporal_shift_birth_death
 Usage   : $tree = constant_rate_birth_death(%options)
 Function: Produces a tree from the model terminating at a given size/time
 Returns : Bio::Phylo::Forest::Tree
 Args    : %options with fields:
           birth_rates The birth rates
           death_rates The death rates
           rate_times  The times after which the rates apply (first element must be 0)
           tree_size  The size of the tree at which to terminate
           tree_age   The age of the tree at which to terminate

 NB: At least one of tree_size and tree_age must be specified

=cut

sub temporal_shift_birth_death {
    my %options = @_;

    #Check that we have a termination condition
    unless ( defined $options{tree_size} or defined $options{tree_age} ) {

        #Error here.
        return undef;
    }

    #Set the undefined condition to infinity
    $options{tree_size} = 1e6 unless defined $options{tree_size};
    $options{tree_age}  = 1e6 unless defined $options{tree_age};

    #Each node gets an ID number this tracks these
    my $node_id = 0;

    #Create a new tree with a root, start the list of terminal species
    my $tree = Bio::Phylo::Forest::Tree->new();
    my $root = Bio::Phylo::Forest::Node->new();
    $root->set_branch_length(0);
    $root->set_name( 'ID' . $node_id++ );
    $tree->insert($root);
    my @terminals = ($root);
    my ( $next_extinction, $next_speciation );
    my $time      = 0;
    my $tree_size = 1;

    #Load current rates
    my $birth_rate = $options{birth_rates}[0];
    my $death_rate = $options{death_rates}[0];
    my $current_rates = 0;
    my $next_rate_change = $options{rate_times}[$current_rates+1];

    #Add an additional time to the end of the rate change times to simplify checking
    push(@{$options{rate_times}},$options{tree_size}*2);

    #Check whether we have a non-zero root edge
    if ( defined $options{root_edge} && $options{root_edge} ) {

        #Non-zero root. We set the time to the first speciation event
        $next_speciation = -log(rand) / $birth_rate / $tree_size;
    }
    else {

        #Zero root, we want a speciation event straight away
        $next_speciation = 0;
    }

#        print  "RATES:".$time."|".$birth_rate."|".$death_rate."\n";
    #Time of the first extinction event. If no extinction we always
    #set the extinction event after the current speciation event
    $next_extinction = -log(rand) / $death_rate / $tree_size;

    #While the tree has not become extinct and the termination criterion
    #has not been achieved we create new speciation and extinction events
    while ($tree_size > 0
        && $tree_size < $options{tree_size}
        && $time < $options{tree_age} )
    {

#        print  "TIMES:".$next_extinction."|".$next_speciation."|".$next_rate_change."|\n";
#        print  "RATES:".$time."|".$birth_rate."|".$death_rate."\n";
        #Add the time since the last event to all terminal species
        foreach (@terminals) {
            $_->set_branch_length(
                $_->get_branch_length + min(
                    $next_extinction, $next_speciation,
                    $options{tree_age} - $time
                )
            );
        }

        #Update the time
        my $time_last = $time;
        $time += min( $next_extinction, $next_speciation, $next_rate_change-$time_last);

        #If the tree exceeds the time limit we are done
        return $tree if ( $time > $options{tree_age} );


        if ($next_rate_change-$time_last < min( $next_extinction, $next_speciation) )
        {
            $current_rates += 1;
            $birth_rate = $options{birth_rates}[$current_rates];
            $death_rate = $options{death_rates}[$current_rates];
            $next_rate_change = $options{rate_times}[$current_rates+1];

        } else
        {


            #Get the species effected by this event and remove it from the terminals list
            my $effected = splice( @terminals, int( rand( scalar @terminals ) ), 1 );

            #If we have a speciation event we add two new species
            if ( $next_speciation < $next_extinction || !defined $next_extinction )
            {
                foreach ( 1, 2 ) {

                    #Create a new species
                    my $child = Bio::Phylo::Forest::Node->new();
                    $child->set_name( 'ID' . $node_id++ );

                    #Give it a zero edge length
                    $child->set_branch_length(0);

                    #Add it as a child to the speciating species
                    $effected->set_child($child);

                    #Add it to the tree
                    $tree->insert($child);

                    #Add it to the terminals list
                    push( @terminals, $child );
                }
            }
        }

        #We calculate the time that the next extinction and speciation
        #events will occur (only the earliest of these will actually
        #happen). NB: this approach is only appropriate for models where
        #speciation and extinction times are exponentially distributed.
        #Windows sometimes returns 0 values for rand...
        my ( $r1, $r2 ) = ( 0, 0 );
        $r1 = rand until $r1;
        $r2 = rand until $r2;
        $tree_size = scalar @terminals;
        return $tree unless $tree_size;
        $next_speciation = -log($r1) / $birth_rate / $tree_size;
        $next_extinction = -log($r2) / $death_rate / $tree_size;

        if ((scalar @terminals)%100==0)
        {
                print $time."|".@terminals."|\n";
                }

    }
    return $tree;
}

=item evolving_speciation_rate()

An evolutionary model featuring evolving speciation rates. Each daughter
species is assigned its parent's speciation rate multiplied by a normally
distributed noise factor.

 Type    : Evolutionary model
 Title   : evolving_speciation_rate
 Usage   : $tree = evolving_speciation_rate(%options)
 Function: Produces a tree from the model terminating at a given size/time
 Returns : Bio::Phylo::Forest::Tree
 Args    : %options with fields:
           birth_rate The initial speciation rate (default 1)
           evolving_std The standard deviation of the normal distribution
                      from which the rate multiplier is drawn.
           tree_size  The size of the tree at which to terminate
           tree_age   The age of the tree at which to terminate

 NB: At least one of tree_size and tree_age must be specified

=cut

sub evolving_speciation_rate {
    my %options = @_;

    #Check that we have a termination condition
    unless ( defined $options{tree_size} or defined $options{tree_age} ) {

        #Error here.
        return undef;
    }

    #Set the undefined condition to infinity
    $options{tree_size} = 1e6 unless defined $options{tree_size};
    $options{tree_age}  = 1e6 unless defined $options{tree_age};

    #Set default rates
    $options{birth_rate}   = 1 unless defined( $options{birth_rate} );
    $options{evolving_std} = 1 unless defined( $options{evolving_std} );

    #Each node gets an ID number this tracks these
    my $node_id = 0;

    #Create a new tree with a root, start the list of terminal species
    my $tree = Bio::Phylo::Forest::Tree->new();
    my $root = Bio::Phylo::Forest::Node->new();
    $root->set_branch_length(0);
    $root->set_name( 'ID' . $node_id++ );
    $tree->insert($root);
    my @terminals   = ($root);
    my @birth_rates = ( $options{birth_rate} );
    my $net_rate    = $options{birth_rate};
    my $next_speciation;
    my $time      = 0;
    my $tree_size = 1;

    #Check whether we have a non-zero root edge
    if ( defined $options{root_edge} && $options{root_edge} ) {

        #Non-zero root. We set the time to the first speciation event
        $next_speciation = -log(rand) / $options{birth_rate} / $tree_size;
    }
    else {

        #Zero root, we want a speciation event straight away
        $next_speciation = 0;
    }

    #While we haven't reached termination
    while ( $tree_size < $options{tree_size} && $time < $options{tree_age} ) {

        #Add the time since the last event to all terminal species
        foreach (@terminals) {
            $_->set_branch_length( $_->get_branch_length +
                  min( $next_speciation, $options{tree_age} - $time ) );
        }

        #Update the time
        $time += $next_speciation;

        #If the tree exceeds the time limit we are done
        return $tree if ( $time > $options{tree_age} );

        #Get the species effected by this event
        my $rand_select = rand($net_rate);
        my $selected    = 0;
        for (
            ;
            $selected < scalar @terminals
            && $rand_select > $birth_rates[$selected] ;
            $selected++
          )
        {
            $rand_select -= $birth_rates[$selected];
        }

        #Remove it from the terminals list
        my $effected      = splice( @terminals,   $selected, 1 );
        my $effected_rate = splice( @birth_rates, $selected, 1 );

        #Update the net speciation rate
        $net_rate -= $effected_rate;

        #If we have a speciation event we add two new species
        foreach ( 1, 2 ) {

            #Create a new species
            my $child = Bio::Phylo::Forest::Node->new();
            $child->set_name( 'ID' . $node_id++ );

            #Give it a zero edge length
            $child->set_branch_length(0);

            #Add it as a child to the speciating species
            $effected->set_child($child);

            #Add it to the tree
            $tree->insert($child);

            #Add it to the terminals list
            push( @terminals, $child );

            #New speciation rate
            my $new_speciation_rate =
              $effected_rate * ( 1 + qnorm(rand) * $options{evolving_std} );
            if ( $new_speciation_rate < 0 ) { $new_speciation_rate = 0; }
            push( @birth_rates, $new_speciation_rate );
            $net_rate += $new_speciation_rate;
        }

        #   $net_rate = 0;
        #   foreach (@birth_rates) { $net_rate += $_; }
        #Windows sometimes returns 0 values for rand...
        my ( $r1, $r2 ) = ( 0, 0 );
        $r1 = rand until $r1;
        $tree_size = scalar @terminals;

        #If all species have stopped speciating (unlikely)
        if ( $net_rate == 0 ) {
            return $tree;
        }
        $next_speciation = -log($r1) / $net_rate / $tree_size;
        return $tree unless $tree_size;
    }
    return $tree;
}


=item clade_shifts()

A constant rate birth-death model with punctuated changes in the speciation
and extinction rates. At each change one lineage receives new pre-specified
speciation and extinction rates.

 Type    : Evolutionary model
 Title   : clade_shifts
 Usage   : $tree = clade_shifts(%options)
 Function: Produces a tree from the model terminating at a given size/time
 Returns : Bio::Phylo::Forest::Tree
 Args    : %options with fields:
           birth_rates The speciation rates
           death_rates The death rates
           rate_times  The times at which the rates are introduced to a new
             clade. The first time should be zero. The remaining must be in
             ascending order.
           tree_size  The size of the tree at which to terminate
           tree_age   The age of the tree at which to terminate

 NB: At least one of tree_size and tree_age must be specified

=cut

sub clade_shifts {
    my %options = @_;

    #Check that we have a termination condition
    unless ( defined $options{tree_size} or defined $options{tree_age} ) {

        #Error here.
        return undef;
    }

    #Set the undefined condition to infinity
    $options{tree_size} = 1e6 unless defined $options{tree_size};
    $options{tree_age}  = 1e6 unless defined $options{tree_age};

    #Each node gets an ID number this tracks these
    my $node_id = 0;

    #Create a new tree with a root, start the list of terminal species
    my $tree = Bio::Phylo::Forest::Tree->new();
    my $root = Bio::Phylo::Forest::Node->new();
    $root->set_branch_length(0);
    $root->set_name( 'ID' . $node_id++ );
    $tree->insert($root);

    #rates
    my @birth_rates_in = @{$options{birth_rates}};
    my @death_rates_in = @{$options{death_rates}};
    my @rate_times_in = @{$options{rate_times}};

    if ($rate_times_in[0] != 0)
    {
        Bio::Phylo::Util::Exceptions::BadArgs->throw( 'error' =>
              "The first rate time must be 0" );
    }
    if (scalar @birth_rates_in != scalar @death_rates_in)
    {
        Bio::Phylo::Util::Exceptions::BadArgs->throw( 'error' =>
              "birth and death rates must have the same length" );
    }
    if (scalar @birth_rates_in != scalar @rate_times_in)
    {
        Bio::Phylo::Util::Exceptions::BadArgs->throw( 'error' =>
              "birth/death rates must have the same length as rate times" );
    }
    my @birth_rates = ($birth_rates_in[0]);
    my @death_rates = ($death_rates_in[0]);

    my $net_birth_rate = $birth_rates[0];
    my $net_death_rate = $death_rates[0];

    shift(@birth_rates_in);
    shift(@death_rates_in);

    my @terminals = ($root);
    my ( $next_extinction, $next_speciation, $next_rate_change );
    my $time      = 0;
    my $tree_size = 1;
    my $inf = 9**9**9**9;

    #Check whether we have a non-zero root edge
    if ( defined $options{root_edge} && $options{root_edge} ) {

        #Non-zero root. We set the time to the first speciation event
        if ($birth_rates[0] > 0)
        {
            $next_speciation = -log(rand) / $birth_rates[0] ;
        } else
        {
            $next_speciation = $inf;
        }
    }
    else {
        #Zero root, we want a speciation event straight away
        $next_speciation = 0.0;
    }

    #Time of the first extinction event. If no extinction we always
    #set the extinction event after the current speciation event
    if ($death_rates[0] > 0)
    {
        $next_extinction = -log(rand) / $death_rates[0];
    } else
    {
        $next_extinction = $inf;
    }

    #Time of next rate change
    shift(@rate_times_in); #pop the initial 0
    $next_rate_change = shift(@rate_times_in);

    #While the tree has not become extinct and the termination criterion
    #has not been achieved we create new speciation and extinction events
    while ($tree_size > 0
        && $tree_size < $options{tree_size}
        && $time < $options{tree_age} )
    {
        #print $time."|".$tree_size."\n";
        #Update rates if a clade shift is happening
        #TODO index rates or pop off one at a time.


        #Add the time since the last event to all terminal species
        foreach (@terminals) {
            $_->set_branch_length(
                $_->get_branch_length + min(
                    $next_extinction,
                    $next_speciation,
                    $next_rate_change-$time,
                    $options{tree_age} - $time
                )
            );
        }

        #Update the time
        my $time_last = $time;
        $time += min( $next_extinction, $next_speciation, $next_rate_change-$time_last );

        #If the tree exceeds the time limit we are done
        return $tree if ( $time > $options{tree_age} );

        #We have a rate change
        if ($next_rate_change-$time_last < min($next_extinction, $next_speciation))
        {
            #Find a random species to effect
            my $effected_species = int(rand($tree_size));
            #Subtract current rates
            $net_death_rate -= $death_rates[$effected_species];
            $net_birth_rate -= $birth_rates[$effected_species];
            #Get new rates
            $death_rates[$effected_species] = shift(@death_rates_in);
            $birth_rates[$effected_species] = shift(@birth_rates_in);
            #Add new rates
            $net_death_rate += $death_rates[$effected_species];
            $net_birth_rate += $birth_rates[$effected_species];
            #Get next rate change time
            if (scalar(@rate_times_in))
            {
                $next_rate_change = shift(@rate_times_in);
            } else
            {
                $next_rate_change = $inf;
            }
        }
        #Choosing a random species to speciate
        else
        {
            my $selected    = 0;
            if ( $next_speciation < $next_extinction )
            {
                my $rand_select = rand($net_birth_rate);
                for (
                    ;
                    $selected < scalar @terminals
                    && $rand_select > $birth_rates[$selected] ;
                    $selected++
                )
                {
                    $rand_select -= $birth_rates[$selected];
                }
            } else
            {
                my $rand_select = rand($net_death_rate);
                for (
                    ;
                    $selected < scalar @terminals
                    && $rand_select > $death_rates[$selected] ;
                    $selected++
                )
                {
                    $rand_select -= $death_rates[$selected];
                }
            }
            if ($net_birth_rate == 0)
            {
                $selected = 0;
            }

            #Remove the species effected by this event and remove it from the terminals list
            my $effected      = splice( @terminals,   $selected, 1 );
            my $effected_birth_rate = splice( @birth_rates, $selected, 1 );
            my $effected_death_rate = splice( @death_rates, $selected, 1 );

            $net_birth_rate -= $effected_birth_rate;
            $net_death_rate -= $effected_death_rate;

            #If we have a speciation event we add two new species
            if ( $next_speciation < $next_extinction || !defined $next_extinction )
            {
                foreach ( 1, 2 ) {

                    #Create a new species
                    my $child = Bio::Phylo::Forest::Node->new();
                    $child->set_name( 'ID' . $node_id++ );

                    #Give it a zero edge length
                    $child->set_branch_length(0);

                    #Add it as a child to the speciating species
                    $effected->set_child($child);

                    #Add it to the tree
                    $tree->insert($child);

                    #Add it to the terminals list
                    push( @terminals, $child );

                    push( @birth_rates, $effected_birth_rate );
                    push( @death_rates, $effected_death_rate );

                    $net_death_rate += $effected_death_rate;
                    $net_birth_rate += $effected_birth_rate;

                }
            }
        }
        #We calculate the time that the next extinction and speciation
        #events will occur (only the earliest of these will actually
        #happen). NB: this approach is only appropriate for models where
        #speciation and extinction times are exponentially distributed.

        #Windows sometimes returns 0 values for rand...
        my ( $r1, $r2 ) = ( 0, 0 );
        $r1 = rand until $r1;
        $r2 = rand until $r2;

        #The current tree size
        $tree_size = scalar @terminals;

        return $tree unless $tree_size;

        if ($net_birth_rate > 0)
        {
            $next_speciation = -log($r1) / $net_birth_rate;
        } else
        {
            $next_speciation = $inf;
        }

        if ($net_death_rate > 0)
        {
            $next_extinction = -log($r2) / $net_death_rate;
        } else
        {
            $next_extinction = $inf;
        }

    }
    return $tree;
}

=item beta_binomial()

An evolutionary model featuring evolving speciation rates. From Blum2007

 Type    : Evolutionary model
 Title   : beta_binomial
 Usage   : $tree = beta_binomial(%options)
 Function: Produces a tree from the model terminating at a given size/time
 Returns : Bio::Phylo::Forest::Tree
 Args    : %options with fields:
           birth_rate The initial speciation rate (default 1)
           model_param The parameter as defined in Blum2007
           tree_size  The size of the tree at which to terminate
           tree_age   The age of the tree at which to terminate

 NB: At least one of tree_size and tree_age must be specified

=cut

sub beta_binomial {
    my %options = @_;

    #Check that we have a termination condition
    unless ( defined $options{tree_size} or defined $options{tree_age} ) {

        #Error here.
        return undef;
    }

    #Set the undefined condition to infinity
    $options{tree_size} = 1e6 unless defined $options{tree_size};
    $options{tree_age}  = 1e6 unless defined $options{tree_age};

    #Set default rates
    $options{birth_rate}  = 1 unless defined( $options{birth_rate} );
    $options{model_param} = 0 unless defined( $options{model_param} );

    #Each node gets an ID number this tracks these
    my $node_id = 0;

    #Create a new tree with a root, start the list of terminal species
    my $tree = Bio::Phylo::Forest::Tree->new();
    my $root = Bio::Phylo::Forest::Node->new();
    $root->set_branch_length(0);
    $root->set_name( 'ID' . $node_id++ );
    $tree->insert($root);
    my @terminals   = ($root);
    my @birth_rates = ( $options{birth_rate} );
    my $net_rate    = $options{birth_rate};
    my $next_speciation;
    my $time      = 0;
    my $tree_size = 1;

    #Check whether we have a non-zero root edge
    if ( defined $options{root_edge} && $options{root_edge} ) {

        #Non-zero root. We set the time to the first speciation event
        $next_speciation = -log(rand) / $options{birth_rate} / $tree_size;
    }
    else {

        #Zero root, we want a speciation event straight away
        $next_speciation = 0;
    }

    #While we haven't reached termination
    while ( $tree_size < $options{tree_size} && $time < $options{tree_age} ) {

        #Add the time since the last event to all terminal species
        foreach (@terminals) {
            $_->set_branch_length( $_->get_branch_length +
                  min( $next_speciation, $options{tree_age} - $time ) );
        }

        #Update the time
        $time += $next_speciation;

        #If the tree exceeds the time limit we are done
        return $tree if ( $time > $options{tree_age} );

        #Get the species effected by this event
        my $rand_select = rand($net_rate);
        my $selected    = 0;
        for (
            ;
            $selected < scalar @terminals
            && $rand_select > $birth_rates[$selected] ;
            $selected++
          )
        {
            $rand_select -= $birth_rates[$selected];
        }

        #Remove it from the terminals list
        my $effected      = splice( @terminals,   $selected, 1 );
        my $effected_rate = splice( @birth_rates, $selected, 1 );
        my $p =
          qbeta( rand, $options{model_param} + 1, $options{model_param} + 1 );

        #If we have a speciation event we add two new species
        foreach ( 1, 2 ) {

            #Create a new species
            my $child = Bio::Phylo::Forest::Node->new();
            $child->set_name( 'ID' . $node_id++ );

            #Give it a zero edge length
            $child->set_branch_length(0);

            #Add it as a child to the speciating species
            $effected->set_child($child);

            #Add it to the tree
            $tree->insert($child);

            #Add it to the terminals list
            push( @terminals, $child );

            #New speciation rate
            my $new_speciation_rate = $effected_rate * $p;
            $p = 1 - $p;
            if ( $new_speciation_rate < 0 ) { $new_speciation_rate = 0; }
            push( @birth_rates, $new_speciation_rate );
        }

        #Windows sometimes returns 0 values for rand...
        my ( $r1, $r2 ) = ( 0, 0 );
        $r1 = rand until $r1;
        $tree_size = scalar @terminals;

        #If all species have stopped speciating (unlikely)
        if ( $net_rate == 0 ) {
            return $tree;
        }
        $next_speciation = -log($r1) / $net_rate / $tree_size;
        return $tree unless $tree_size;
    }
    return $tree;
}

=back

=cut

=begin comment

###########################################################
#INTERNAL METHODS
#These are methods that permit additional manipulations of
#a Bio::Phylo::Tree to be easily made. As such some of
#these could easily be moved into Bio::Phylo::Forest::Tree
###########################################################

=end comment

=cut

=begin comment

 Type    : Internal method
 Title   : copy_tree
 Usage   : $tree = copy_tree($tree)
 Function: Makes a new independent copy of a tree
 Returns : the phylogenetic $tree
 Args    : the phylogenetic $tree

=end comment

=cut

sub copy_tree {
    return Bio::Phylo::IO->parse(
        -format => 'newick',
        -string => shift->to_newick( '-nodelabels' => 1 )
    )->first;
}

=begin comment

 Type    : Internal method
 Title   : truncate_tree_time
 Usage   : truncate_tree_time($tree,$age)
 Function: Truncates the tree at the specified age
 Returns : N/A
 Args    : $tree: the phylogenetic tree (which will be modified)
           $age: the age at which to cut the $tree

=end comment

=cut

sub truncate_tree_time {

    #$node and $time are used only by this function recursively
    my ( $tree, $age, $node, $time ) = @_;

    #If node and time weren't specified we are starting from the root
    $node = $tree->get_root unless defined $node;
    $time = 0 unless defined $time;

    #If we are truncating this branch
    if ( $time + $node->get_branch_length >= $age ) {

        #Collapse the node unless it is terminal
        $node->collapse unless $node->is_terminal();

        #Set the branch length appropriately
        $node->set_branch_length( $age - $time );
        return;
    }

    #If this node has no children we are done
    return if $node->is_terminal();

    #Call the function recursively on the children
    foreach ( @{ $node->get_children } ) {
        truncate_tree_time( $tree, $age, $_,
            $time + $node->get_branch_length() );
    }
}

=begin comment

 Type    : Internal method
 Title   : truncate_tree_size
 Usage   : truncate_tree_size($tree,$size)
 Function: Truncates the tree to the specified number of species
 Returns : N/A
 Args    : $tree: the phylogenetic tree (which will be modified)
           $size: random species are delete so that the tree has this many terminals

=end comment

=cut

sub truncate_tree_size {
    my ( $tree, $size ) = @_;
    my @terminals = @{ $tree->get_terminals };
    my @names;

    #Calculate the tree height and node distances from the root
    #much more efficient to do this in one hit than repeatedly
    #calling the analogous functions on the tree
    _calc_node_properties($tree);

    #Only push species that are extant and store the number of those
    #    my $tree_height = $tree->get_tallest_tip->calc_path_to_root;
    my $tree_height = $tree->get_root->get_generic('tree_height');
    foreach (@terminals) {
        if (
            abs(
                ( $_->get_generic('root_distance') - $tree_height ) /
                  $tree_height
            ) < 1e-6
          )
        {
            push( @names, $_->get_name );
        }
    }
    if ( @names < $size ) { print "Internal error\n"; }
    my %deletions;
    while ( scalar( keys %deletions ) < @names - $size ) {
        $deletions{ $names[ int( rand(@names) ) ] } = 1;
    }
    $tree = prune_tips( $tree, [ keys %deletions ] );
    return $tree;
}

sub _get_ultrametric_size {
    my ( $tree, $size ) = @_;
    my @terminals = @{ $tree->get_terminals };
    _calc_node_properties($tree);
    my @names;
    my $tree_height = $tree->get_root->get_generic('tree_height');
    foreach (@terminals) {
        if (
            abs(
                ( $_->get_generic('root_distance') - $tree_height ) /
                  $tree_height
            ) < 1e-6
          )
        {
            push( @names, $_->get_name );
        }
    }
    return scalar @names;
}

=begin comment

 Type    : Internal method
 Title   : remove_extinct_species
 Usage   : remove_extinct_species($tree)
 Function: Removes extinct species from the tree. An extinct species
           is a terminal that does not extend as far as the furthest
           terminal(s).
 Returns : N/A
 Args    : $tree: the phylogenetic tree (which will be modified)
           $age: the age at which to cut the $tree

=end comment

=cut

sub remove_extinct_species {
    my $tree = shift;

    #Calculate the tree height and node distances from the root
    #much more efficient to do this in one hit than repeatedly
    #calling the analogous functions on the tree
    _calc_node_properties($tree);
    my $height = $tree->get_root->get_generic('tree_height');
    return unless $height > 0;
    my $leaves = $tree->get_terminals;
    return unless $leaves;
    my @remove;
    foreach ( @{$leaves} ) {

        unless (
            abs( ( $_->get_generic('root_distance') - $height ) / $height ) <
            1e-6 )
        {
            push( @remove, $_->get_name );
        }
    }
    $tree = prune_tips( $tree, \@remove );
    return $tree;
}

=begin comment

 Type    : Internal method
 Title   : prune_tips
 Usage   : prune_tips($tree,$tips)
 Function: Removes named terminals from the tree
 Returns : N/A
 Args    : $tree: the phylogenetic tree (which will be modified)
           $tips: array ref of terminals to remove from the tree

NB: Available as $tree->prune_tips($tips), but had some problems with
this.

=end comment

=cut

sub prune_tips {
    my ( $self, $tips ) = @_;
    my %names_to_delete = map { $_           => 1 } @{$tips};
    my %keep            = map { $_->get_name => 1 }
      grep { not exists $names_to_delete{ $_->get_name } }
      @{ $self->get_terminals };
    $self->visit_depth_first(
        -post => sub {
            my $node = shift;
            if ( $node->is_terminal ) {
                if ( not $keep{ $node->get_name } ) {
                    $node->set_parent();
                    $self->delete($node);
                }
            }
            else {
                my $seen_tip_to_keep = 0;
                for my $tip ( @{ $node->get_terminals } ) {
                    $seen_tip_to_keep++ if $keep{ $tip->get_name };
                }
                if ( not $seen_tip_to_keep ) {
                    $node->set_parent();
                    $self->delete($node);
                }
            }
        }
    );
    $self->remove_unbranched_internals;
    return $self;
}

=begin comment

 Type    : Internal method
 Title   : lineage_through_time
 Usage   : ($time,$count) = lineage_through_time($tree)
 Function: Alternative to $tree->ltt that permits extinctions
 Returns : $time: array ref of times
           $count: array ref of species counts corresponding to the times
 Args    : the phylogenetic $tree for which to produce the ltt data

=end comment

=cut

sub lineage_through_time {
    my $tree = shift;
    my ( $speciation, $extinction ) = _recursive_ltt_helper($tree);
    my @speciation = sort { $a <=> $b } @{$speciation};
    my @extinction = sort { $a <=> $b } @{$extinction};
    my @time       = (0);
    my @count      = (1);
    my $n_species  = 1;
    my $end_time = max( @speciation, @extinction );
    return ( [], [] ) if ( $end_time == 0 );

#We remove any extinction events occurring at the very end of the tree (as they are not real extinctions)
    while ( scalar @extinction
        && ( $end_time - $extinction[-1] ) / $end_time < 1e-6 )
    {
        pop @extinction;
    }
    while ( scalar @speciation || scalar @extinction ) {
        if ( scalar @extinction == 0
            || ( scalar @speciation && $speciation[0] < $extinction[0] ) )
        {
            push( @count, ++$n_species );
            push( @time,  shift(@speciation) );
        }
        else {
            push( @count, --$n_species );
            push( @time,  shift(@extinction) );
        }
    }
    return ( \@time, \@count );
}

=begin comment

 Type    : Internal method
 Title   : _recursive_ltt_helper
 Usage   : ($speciation, $extinction) = _recursive_ltt_helper($tree)
 Function: Helper for lineage_through_time
 Returns : $speciation: array ref of speciation times
           $extinction: array ref of extinction times
 Args    : the phylogenetic $tree for which to produce the ltt data

=end comment

=cut

sub _recursive_ltt_helper {
    my ( $tree, $node, $time ) = @_;

    #If we are being invoked at the root level
    $node = $tree->get_root unless defined $node;
    $time = 0 unless defined $time;

    #The new time
    $time += $node->get_branch_length;
    return ( [], [$time] ) if ( $node->is_terminal );
    my @speciation;
    my @extinction;
    foreach ( @{ $node->get_children } ) {
        my ( $spec, $ext ) = _recursive_ltt_helper( $tree, $_, $time );
        @speciation = ( @speciation, @{$spec} );
        @extinction = ( @extinction, @{$ext} );
    }
    push( @speciation, $time );
    return ( \@speciation, \@extinction );
}

=begin comment

 Type    : Internal method.
 Title   : _calc_node_properties
 Usage   : _calc_node_properties($tree);
 Function: Calculates the distance of nodes from the root
 Returns : The maximum distance from the root
 Args    :

=end comment

=cut

sub _calc_node_properties {
    my ( $node, $root_distance );
    my $tree = shift;
    my $root = $tree->get_root;

    #Check whether we were given a node and distance
    if ( scalar @_ ) {
        $node          = shift;
        $root_distance = shift;

        #Otherwise the root is the default
    }
    else {
        $node = $root;
        $root->set_generic( tree_height => 0 );
        $root_distance = 0;
    }
    $node->set_generic( root_distance => $root_distance );
    if ( $root_distance > $root->get_generic('tree_height') ) {
        $root->set_generic( tree_height => $root_distance );
    }
    my $terminal_count = 0;
    my $children       = $node->get_children;
    if ( defined $children ) {
        foreach ( @{$children} ) {
            _calc_node_properties( $tree, $_,
                $root_distance + $_->get_branch_length() );
        }
    }
}

=begin comment

 Type    : Internal method
 Title   : nchoosek
 Usage   : $out = nchoosek($n,$k)
 Function: Returns the binomial coefficient for $n and $k
 Returns : the binomial coefficient
 Args    : $n, $k

=end comment

=cut

sub nchoosek {
    my ( $n, $k ) = @_;
    my $r = 1;
    return 0 if ( $k > $n || $k < 0 );
    for ( my $d = 1 ; $d <= $k ; $d++ ) {
        $r *= $n--;
        $r /= $d;
    }
    return $r;
}
1;