[Bio] / FigKernelPackages / gjonewicklib.pm Repository:
ViewVC logotype

Diff of /FigKernelPackages/gjonewicklib.pm

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.21, Sat Aug 21 17:12:51 2010 UTC revision 1.22, Fri Aug 27 23:34:33 2010 UTC
# Line 272  Line 272 
272  # \@tips    = tip_neighborhood_representative_tips( $tree, $tip, $n )  # \@tips    = tip_neighborhood_representative_tips( $tree, $tip, $n )
273  #  #
274  #===============================================================================  #===============================================================================
275    #  Random trees
276    #===============================================================================
277    #
278    #   $tree = random_equibranch_tree(  @tips, \%options )
279    #   $tree = random_equibranch_tree( \@tips, \%options )
280    #   $tree = random_equibranch_tree(  @tips )
281    #   $tree = random_equibranch_tree( \@tips )
282    #
283    #   $tree = random_ultrametric_tree(  @tips, \%options )
284    #   $tree = random_ultrametric_tree( \@tips, \%options )
285    #   $tree = random_ultrametric_tree(  @tips )
286    #   $tree = random_ultrametric_tree( \@tips )
287    #
288    #===============================================================================
289  #  Tree reading and writing:  #  Tree reading and writing:
290  #===============================================================================  #===============================================================================
291  #  Write machine-readable trees:  #  Write machine-readable trees:
# Line 391  Line 405 
405          tip_neighborhood_representative_tree          tip_neighborhood_representative_tree
406          tip_neighborhood_representative_tips          tip_neighborhood_representative_tips
407    
408            random_equibranch_tree
409            random_ultrametric_tree
410    
411          writeNewickTree          writeNewickTree
412          fwriteNewickTree          fwriteNewickTree
413          strNewickTree          strNewickTree
# Line 2983  Line 3000 
3000    
3001    
3002  #===============================================================================  #===============================================================================
3003    #  Random trees
3004    #===============================================================================
3005    #
3006    #   $tree = random_equibranch_tree(  @tips, \%options )
3007    #   $tree = random_equibranch_tree( \@tips, \%options )
3008    #   $tree = random_equibranch_tree(  @tips )
3009    #   $tree = random_equibranch_tree( \@tips )
3010    #
3011    #  Options:
3012    #
3013    #     length => $branch_length   # D = 1
3014    #
3015    #-------------------------------------------------------------------------------
3016    sub random_equibranch_tree
3017    {
3018        my $opts = $_[ 0] && ref $_[ 0] eq 'HASH' ? shift
3019                 : $_[-1] && ref $_[-1] eq 'HASH' ? pop
3020                 :                                  {};
3021        return undef if ! defined $_[0];
3022    
3023        my @tips = ref $_[0] ? @{ $_[0] } : @_;
3024        return undef if @tips < 2;
3025    
3026        my $len = $opts->{ length } ||= 1;
3027    
3028        if ( @tips == 2 )
3029        {
3030            return [ [ map { [ [], $_, $len ] } @tips ], undef, 0 ];
3031        }
3032    
3033        my $tree = [ [ ], undef, 0 ];
3034    
3035        my @links;  # \$anc_dl[i], i.e. a reference to an element in a descendent list
3036    
3037        my $anc_dl = $tree->[0];
3038        foreach my $tip ( splice( @tips, 0, 3 ) )
3039        {
3040            my $node = [ [], $tip, $len ];
3041            push @$anc_dl, $node;
3042            push @links, \$anc_dl->[-1];  #  Ref to the just added descendent list entry
3043        }
3044    
3045        foreach my $tip ( @tips )
3046        {
3047            my $link    = $links[ int( rand( scalar @links ) ) ];
3048            my $newtip  = [ [], $tip, $len ];
3049            my $new_dl  = [ $$link, $newtip ];
3050            my $newnode = [ $new_dl, undef, $len ];
3051            $$link = $newnode;
3052            push @links, \$new_dl->[0], \$new_dl->[1]
3053        }
3054    
3055        return $tree;
3056    }
3057    
3058    
3059    #-------------------------------------------------------------------------------
3060    #
3061    #   $tree = random_ultrametric_tree(  @tips, \%options )
3062    #   $tree = random_ultrametric_tree( \@tips, \%options )
3063    #   $tree = random_ultrametric_tree(  @tips )
3064    #   $tree = random_ultrametric_tree( \@tips )
3065    #
3066    #  Options:
3067    #
3068    #     depth => $root_to_tip_dist   # D = 1
3069    #
3070    #-------------------------------------------------------------------------------
3071    sub random_ultrametric_tree
3072    {
3073        my $opts = $_[ 0] && ref $_[ 0] eq 'HASH' ? shift
3074                 : $_[-1] && ref $_[-1] eq 'HASH' ? pop
3075                 :                                  {};
3076        return undef if ! defined $_[0];
3077    
3078        my @tips = ref $_[0] ? @{ $_[0] } : @_;
3079        return undef if @tips < 2;
3080    
3081        my $d2tip = $opts->{ depth } ||= 1;
3082    
3083        #  Random tip addition order (for rooted tree it matters):
3084    
3085        @tips = sort { rand() <=> 0.5 } @tips;
3086        my $tree = [ [ ], undef, 0 ];
3087    
3088        my $subtree_size = { $tree => 0 };  # total branch length of each subtree
3089    
3090        #  We start with root bifurcation:
3091    
3092        foreach my $tip ( splice( @tips, 0, 2 ) )
3093        {
3094            my $node = [ [], $tip, $d2tip ];
3095            push @{ $tree->[0] }, $node;
3096            $subtree_size->{ $node }  = $d2tip;
3097            $subtree_size->{ $tree } += $d2tip;
3098        }
3099    
3100        #  Add each remaining tip at $pos, measured along the contour length
3101        #  of the tree (with no retracing along branches).
3102    
3103        foreach my $tip ( @tips )
3104        {
3105            my $pos = rand( $subtree_size->{ $tree } );
3106            random_add_to_ultrametric_tree( $tree, $tip, $subtree_size, $pos, $d2tip );
3107        }
3108    
3109        return $tree;
3110    }
3111    
3112    
3113    sub random_add_to_ultrametric_tree
3114    {
3115        my ( $node, $tip, $subtree_size, $pos, $d2tip ) = @_;
3116        $node && $node->[0] && ref $node->[0] eq 'ARRAY' or die "Bad tree node passed to random_add_to_ultrametric_tree().\n";
3117    
3118        # Find the descendent line that it goes in:
3119    
3120        my $i;
3121        my $dl = $node->[0];
3122        my $size0 = $subtree_size->{ $dl->[0] };
3123        if ( $size0 > $pos ) { $i = 0 } else { $i = 1; $pos -= $size0 }
3124        my $desc = $dl->[$i];
3125    
3126        # Does it go within the subtree, or the branch to the subtree?
3127    
3128        my $len;
3129        my $added;
3130        if ( ( $len = $desc->[2] ) <= $pos )
3131        {
3132            $added = random_add_to_ultrametric_tree( $desc, $tip, $subtree_size, $pos - $len, $d2tip - $len );
3133        }
3134        else
3135        {
3136            # If not in subtree, then it goes in the branch to the descendent node
3137            #
3138            #     ----- node  ------------       node
3139            #       ^   /  \       ^             /  \
3140            #       |       \      | pos             \l1
3141            #       |        \     v                  \
3142            #       |      len\ ----------         newnode
3143            #       |          \                     /  \ l2
3144            # d2tip |           \                   /    \
3145            #       |           desc               /     desc
3146            #       |           /  \            l3/      /  \
3147            #       |          .    .            /      .    .
3148            #       v         .      .          /      .      .
3149            #     -----      .        .     newtip    .        .
3150    
3151            my $l1      = $pos;
3152            my $l2      = $len   - $pos;
3153            my $l3      = $d2tip - $pos;
3154            my $newtip  = [ [], $tip, $l3 ];
3155            my $newnode = [ [ $desc, $newtip ], undef, $l1 ];
3156            $dl->[$i]   = $newnode;
3157            $subtree_size->{ $newtip  } = $l3;
3158            $subtree_size->{ $newnode } = $subtree_size->{ $desc } + $l3;
3159            $desc->[2] = $l2;
3160            $subtree_size->{ $desc } -= $l1;
3161            $added = $l3;
3162        }
3163    
3164        #  New branch was inserted below this point:
3165    
3166        $subtree_size->{ $node } += $added;
3167        return $added;
3168    }
3169    
3170    
3171    
3172    #===============================================================================
3173  #  #
3174  #  Tree writing and reading  #  Tree writing and reading
3175  #  #

Legend:
Removed from v.1.21  
changed lines
  Added in v.1.22

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3