Release 1.0.0

git-svn-id: svn+ssh://imag/users/huron/danjean/svnroot/claire/altree/tags/release-1.0.0@123 cf695345-040a-0410-a956-b889e835fe2e

trunk/.cvsignore

+META.yml
+Makefile

trunk/ALTree/Base.pm

+package ALTree::Base;
+
+###########################################
+########  MAIN DATA STRUCTURES    #########
+###########################################
+
+sub _init {
+    my $self=shift;
+    if (@_) {
+        my %extra = @_;
+        @$self{keys %extra} = values %extra;
+    }
+}
+
+sub Debug {
+    my $self=shift;
+
+    #print STDERR @_;
+}
+
+1;
+

trunk/ALTree/Chi2.pm

+package ALTree::Chi2;
+
+use strict;
+use ALTree::CUtils;
+
+BEGIN {
+    use Exporter   ();
+    our ($VERSION, @ISA, @EXPORT, @EXPORT_OK, %EXPORT_TAGS);
+    
+    # set the version for version checking
+    #$VERSION     = 1.00;
+    # if using RCS/CVS, this may be preferred
+    $VERSION = do { my @r = (q$Revision$ =~ /\d+/g); 
+		    sprintf "%d."."%02d" x $#r, @r }; # must be all one line, for MakeMaker
+    
+    @ISA         = qw(Exporter);
+    @EXPORT      = qw(NON_SIGNIFICATIF SIGNIFICATIF &chi2_significatif &definition_p_chi2 &reech_chi2);
+    #%EXPORT_TAGS = ( );     # eg: TAG => [ qw!name1 name2! ],
+    
+    # your exported package globals go here,
+    # as well as any optionally exported functions
+    @EXPORT_OK   = qw();
+}
+#our @EXPORT_OK;
+
+INIT { 
+    
+}
+ 
+# exported package globals go here
+#our $Var1;
+#our %Hashit;
+
+# non-exported package globals go here
+#our @more;
+#our $stuff;
+
+# initialize package globals, first exported ones
+#$Var1   = '';
+#%Hashit = ();
+
+# then the others (which are still accessible as $Some::Module::stuff)
+#$stuff  = '';
+#@more   = ();
+
+# all file-scoped lexicals must be created before
+# the functions below that use them.
+
+# file-private lexicals go here
+#my $priv_var    = '';
+#my %secret_hash = ();
+
+# here's a file-private function as a closure,
+# callable as &$priv_func;  it cannot be prototyped.
+#my $priv_func = sub {
+#    # stuff goes here.
+#};
+
+# make all your functions, whether exported or not;
+# remember to put something interesting in the {} stubs
+#sub func1      {}    # no prototype
+#sub func2()    {}    # proto'd void
+#sub func3($$)  {}    # proto'd to 2 scalars
+
+# this one isn't exported, but could be called!
+#sub func4(\%)  {}    # proto'd to 1 hash ref
+
+END { }       # module clean-up code here (global destructor)
+
+use constant NON_SIGNIFICATIF => 0;
+use constant SIGNIFICATIF => 1;
+use constant PERM => 1000;
+##################################################################
+# Fonctions de seuil du chi2 (pr-calcul puis stockage dans un tableau)
+
+ my ($chi2_p)="chi2_p doit tre initialis !";
+#y ($chi2_p)=0.05;
+my ($test_prop_p)="test_prop_p doit tre initialis !";
+my ($chi2_seuil)=[];
+# test de significativit. Doit retourner vrai ou faux.
+sub chi2_significatif {
+    my ($ddl) = shift;
+    my ($chi2) = shift;
+
+    if ($ddl < 1) {
+	warn "chi[$ddl] asked...\n";
+    }
+
+    if (not defined($chi2_seuil->[$ddl])) {
+	#my $c=`critchi2 $chi2_p $ddl`+0; # Verif que les 2 appels sont quivalents
+	$chi2_seuil->[$ddl]=ALTree::CUtils::critchi($chi2_p, $ddl);
+	#if ($c != $$chi2_seuil[$ddl]) {
+	  # print STDERR "Critchi2 : $c != $$chi2_seuil[$ddl]\n";
+	#}
+	#print "chi2_seuil[$ddl]= $$chi2_seuil[$ddl]\n";
+	#warn "Seuil chi2 non dfini pour ddl $ddl";
+	#return 0;
+    }
+    return ($chi2 > $chi2_seuil->[$ddl]);
+}
+
+sub definition_p_chi2
+{
+    my($p)=shift;
+    my($pprop)=shift;
+    if (defined $p) {
+	$chi2_p=$p;
+    }
+    if (defined $pprop) {
+	$test_prop_p=$pprop;
+    }
+}
+
+sub chi2_fisher_significatif
+{
+    my($pvalue)=shift;
+
+    return ($pvalue < $chi2_p);
+}
+
+##################################################################
+# Rchantillonnage
+
+# Quelques variables globales pour aller plus vite (mme si c'est  viter
+# en gnral)
+my(@fils_c);
+my(@fils_m);
+my($compteur);
+my($sum_malade);
+my($sum_controle);
+my($sum_total);
+my($nb_fils);
+
+
+my(@th_c, @th_m);
+my($clades);
+sub calcule_chi2
+{
+    my($i, $chi2, $temp);
+    $chi2=0;
+    for ($i=0; $i < $nb_fils; $i++){
+	$temp=($fils_c[$i]-$th_c[$i]);
+	$chi2+=($temp)*($temp)/$th_c[$i];
+
+	$temp=($$clades[$i]-$fils_c[$i]-$th_m[$i]);
+	$chi2+=($temp)*($temp)/$th_m[$i];
+    }
+    #print "Chi2 : $chi2\n";
+    return $chi2;
+}
+
+sub reech_chi2
+{
+    $sum_malade=shift;
+    $sum_controle=shift;
+    $sum_total=$sum_malade+$sum_controle;
+    $nb_fils=shift;
+    my($chi2_reel)=shift;
+    $clades=shift;
+    my($p_val)=0;
+    my($i, $k);
+    #my($alea);
+
+# nb_fils correspond en fait a tous les groupes sur lesquelles il faut faire
+# le chi2.
+# Cet ensemble de groupe a au total: $sum_malade et $sum_controle individus
+# (respectivement malades et controles)
+# clades est une rfrence sur un tableau contenant les effectifs globaux de
+# chaque clade
+
+    
+    #print "Reechantillonage chi2 : ddl : ", ($nb_fils-1), " M: $sum_malade C: $sum_controle\n   ";
+    #print "Chi2 rel : $chi2_reel\n";
+    #print "Clades: ";
+    #for $i (@{$clades}) { print "$i "; } print "\n";
+
+    # Pr calcul des effectifs thoriques
+    for ($i=0; $i < $nb_fils; $i++){
+	$th_c[$i]=($sum_controle*$$clades[$i])/($sum_total);
+	$th_m[$i]=($sum_malade*$$clades[$i])/($sum_total);
+    }
+
+    my($clade, $alea, $malades, $controles);
+    my($nb_tests)=PERM;
+    for ($k=1;$k<=$nb_tests; $k++){
+	$malades=$sum_malade;
+	$controles=$sum_controle;
+	for($clade=0; $clade<$nb_fils; $clade++) {
+	    $fils_m[$clade]=0;
+	    $fils_c[$clade]=0;
+	    for($i=0; $i<$$clades[$clade]; $i++) {
+		$alea=rand($malades+$controles);
+		if ($alea < $malades) {
+		    $malades--;
+		    $fils_m[$clade]++;
+		} else {
+		    $controles--;
+		    $fils_c[$clade]++;
+		}
+	    }
+	}
+
+	if (calcule_chi2 >= $chi2_reel){
+	    $p_val++;
+	}
+    }
+    #DEBUG  print "CHI2=$chi2_reel\n";
+    #print"nb de chi2 non calculable (effectif nul)= $compteur\n";
+    #DEBUG print "p_val1 = ", $p_val/$nb_tests,"\n";
+    # print "chi2_p771=$chi2_p\n";
+    return ($p_val/$nb_tests);
+}
+
+sub reech_significatif 
+{
+    my($p_val)=shift;
+    my($nb_tests)=PERM;
+  #DEBUG  print "Chi2P= $chi2_p\n";
+  #DEBUG  print "p=  ", $p_val , "\n";
+  #DEBUG  print "test=", $p_val/$nb_tests<=$chi2_p, "\n";
+    return ($p_val<=$chi2_p);
+}
+
+1;

trunk/ALTree/Foret.pm

+package ALTree::Foret;
+
+################################################################
+################################################################
+####################### Foret           ########################
+################################################################
+################################################################
+
+use base qw(ALTree::Base ALTree::SiteCollection);
+
+sub New { # [classe]
+    my $class=shift;
+    my $self={};
+    bless($self, $class);
+    $self->InitSiteCollection();
+    $self->_init(@_);
+    $self->Debug("creating Foret\n");
+    return $self;
+}
+
+sub AddTree {
+    my $self=shift;
+    push @{$self->{"trees"}}, @_;
+}
+sub GetTree {
+    my $self=shift;
+    my $index=shift;
+    return $self->{"trees"}->[$index];
+}    
+sub GetTreesList {
+    my $self=shift;
+    return @{$self->{"trees"}};
+}
+
+sub ProvideSite {
+    my $self=shift;
+    my $site_nb=shift;
+
+    if (not $self->HasSiteIndex($site_nb)) {
+	$self->AddSite(ALTree::SitePerForet->New($site_nb));
+    }
+    return $self->GetSite($site_nb);
+}
+
+sub CalculVi {
+    my $self=shift;
+
+    foreach my $tree ($self->GetTreesList()) {
+	foreach my $site ($tree->GetSitesList()) {
+	    foreach my $sens ($site->GetSensList()) {
+		$self->ProvideSite($site->GetSiteNb())
+		    ->ProvideSens($sens->GetSensStruct())
+		    ->PlusVi($sens->GetVit());
+	    }
+	}
+    }
+}
+
+sub _EnsureViMax {
+    my($self)=shift;
+    if (not exists ($self->{"V_i_max"})) {
+	my @tab_trie=sort {
+	    $b->GetViMax() <=> $a->GetViMax()} 
+	   $self->GetSitesList();
+	$self->{"V_i_max"}=$tab_trie[0]->GetViMax();
+	$self->{"V_i_max_tab"}=\@tab_trie;
+    }
+}
+sub GetViMax {
+    my $self=shift;
+    $self->_EnsureViMax();
+    return $self->{"V_i_max"};
+}   
+sub GetViMaxSite {
+    my($self)=shift;
+    my($index)=shift;
+    $self->_EnsureViMax();
+    return $self->{"V_i_max_tab"}->[$index];
+}
+sub GetViMaxSiteList {
+    my($self)=shift; 
+    $self->_EnsureViMax();
+    return @{$self->{"V_i_max_tab"}}; 
+}
+sub NbViMaxSite {
+    my($self)=shift; 
+    $self->_EnsureViMax();
+    return (scalar @{$self->{"V_i_max_tab"}});
+}
+
+sub _EnsureViMaxSens {
+    my($self)=shift;
+    if (not exists ($self->{"V_i_max_sens_tab"})) {
+	my @tab_trie=sort {
+	    $b->GetVi() <=> $a->GetVi()
+	    } map { $_->GetSensList(); } $self->GetSitesList();
+	#if ($tab_trie[0]->GetVi() != $self->GetViMax()) {
+	#    die "Arghh\n";
+	#}
+	$self->{"V_i_max_sens_tab"}=\@tab_trie;
+    }
+}
+sub GetViMaxSens {
+    my($self)=shift;
+    my($index)=shift;
+    $self->_EnsureViMaxSens();
+    return $self->{"V_i_max_sens_tab"}->[$index];
+}
+sub GetViMaxSensList {
+    my($self)=shift; 
+    $self->_EnsureViMaxSens();
+    return @{$self->{"V_i_max_sens_tab"}}; 
+}
+sub NbViMaxSens {
+    my($self)=shift; 
+    $self->_EnsureViMaxSens();
+    return (scalar @{$self->{"V_i_max_sens_tab"}});
+}
+
+1;

trunk/ALTree/Import.pm

+package ALTree::Import;
+
+use strict;
+use ALTree::Utils;
+use ALTree::Input;
+use ALTree::Sens;
+use ALTree::Tree;
+use ALTree::Foret;
+use ALTree::Node;
+use ALTree::SitePerTree;
+use ALTree::SitePerForet;
+use ALTree::SiteSensPerForet;
+
+BEGIN {
+    use Exporter   ();
+    our ($VERSION, @ISA, @EXPORT, @EXPORT_OK, %EXPORT_TAGS);
+    
+    # set the version for version checking
+    #$VERSION     = 1.00;
+    # if using RCS/CVS, this may be preferred
+    $VERSION = do { my @r = (q$Revision$ =~ /\d+/g); 
+		    sprintf "%d."."%02d" x $#r, @r }; # must be all one line, for MakeMaker
+    
+    @ISA         = qw(Exporter);
+    @EXPORT      = qw();
+    #%EXPORT_TAGS = ( );     # eg: TAG => [ qw!name1 name2! ],
+    
+    # your exported package globals go here,
+    # as well as any optionally exported functions
+    @EXPORT_OK   = qw();
+}
+#our @EXPORT_OK;
+
+INIT { 
+    
+}
+ 
+# exported package globals go here
+#our $Var1;
+#our %Hashit;
+
+# non-exported package globals go here
+#our @more;
+#our $stuff;
+
+# initialize package globals, first exported ones
+#$Var1   = '';
+#%Hashit = ();
+
+# then the others (which are still accessible as $Some::Module::stuff)
+#$stuff  = '';
+#@more   = ();
+
+# all file-scoped lexicals must be created before
+# the functions below that use them.
+
+# file-private lexicals go here
+#my $priv_var    = '';
+#my %secret_hash = ();
+
+# here's a file-private function as a closure,
+# callable as &$priv_func;  it cannot be prototyped.
+#my $priv_func = sub {
+#    # stuff goes here.
+#};
+
+# make all your functions, whether exported or not;
+# remember to put something interesting in the {} stubs
+#sub func1      {}    # no prototype
+#sub func2()    {}    # proto'd void
+#sub func3($$)  {}    # proto'd to 2 scalars
+
+# this one isn't exported, but could be called!
+#sub func4(\%)  {}    # proto'd to 1 hash ref
+
+END { }       # module clean-up code here (global destructor)
+
+##################################################################
+# Fonctions de seuil du chi2 (pr-calcul puis stockage dans un tableau)
+
+1;

trunk/ALTree/Node.pm

+package ALTree::Node;
+
+################################################################
+################################################################
+####################### Node          ##########################
+################################################################
+################################################################
+
+use base 'ALTree::Base';
+
+# Structure Node
+#   "id" -> String
+#   "children" -> Array of (Node)
+#   "father" -> Node
+#   "apo" -> Hash of ('num_apo' => SiteSens)
+#   "br_len" -> Integer
+#   "case" -> Integer
+#   "control" -> Integer
+#   "level" -> Integer 
+#   "height" -> Integer
+#   "sequence" ->  string # used only with phylip
+#   "oldfather" -> Node  # aprs fusion des branches nulles (anciennes relations)
+#   "oldchildren" -> Array of (Node)  # aprs fusion des branches nulles (anciennes relations de parent)
+#   "label" -> string # nom noeuds apres fusion 3+(4) par exemple 
+#   "res_test" -> Array of (TestResults)
+
+#Creation d'une structure
+sub New { # [classe] id  
+    my $class=shift;
+    my $self={};
+    my $id=shift;
+    bless($self, $class);
+    $self->_init("id" => $id, "apo" => {}, 
+		 "children" => [], @_);
+    $self->Debug("creating Node $id\n");
+    return $self;
+}
+
+sub GetId {
+    my $self=shift;
+    return $self->{"id"};
+}
+
+sub SetFather {
+    my $self=shift;
+    my $father=shift;
+    $self->{"father"}=$father;
+}
+sub GetFather {
+    my $self=shift;
+    return $self->{"father"};
+}
+sub HasFather {
+    my $self=shift;
+    return defined($self->{"father"});
+}
+sub Father {
+    my $self=shift;
+    my $newfather=shift;
+    if (defined($newfather)) {
+	$self->SetFather($newfather);
+    }
+    return $self->GetFather();
+}
+sub RecordOrigFather {
+    my $self=shift;
+    my $father=shift;
+    if ($self->HasFather()) {
+    	$self->{"orig_father"}=$self->GetFather();
+    }
+}
+sub GetOrigFather {
+    my $self=shift;
+    if (exists($self->{"orig_father"})) {
+    	return $self->{"orig_father"};
+    } else {
+    	return $self->GetFather();
+    }
+}
+
+sub SetCase {
+    my $self=shift;
+    my $value=shift;
+    $self->{"case"}=$value;
+}
+sub GetCase {
+    my $self=shift;
+    return $self->{"case"};
+}
+sub EraseCase {
+    my $self=shift;
+    delete($self->{"case"});
+}
+
+sub SetControl {
+    my $self=shift;
+    my $value=shift;
+    $self->{"control"}=$value;
+}
+sub GetControl {
+    my $self=shift;
+    return $self->{"control"};
+}
+sub EraseControl {
+    my $self=shift;
+    delete($self->{"control"});
+}
+
+sub SetBrLen {
+    my $self=shift;
+    my $br_len=shift;
+    $self->{"br_len"}=$br_len;
+}
+sub GetBrLen {
+    my $self=shift;
+    return $self->{"br_len"};
+}
+sub HasBrLen {
+    my $self=shift;
+    return exists($self->{"br_len"});
+}
+
+sub AddOldChild {
+    my $self=shift;
+    push @{$self->{"oldchildren"}}, @_;
+}
+sub GetOldChildrenList {
+    my $self=shift;
+    return @{$self->{"oldchildren"}};
+}
+
+sub AddChild {
+    my $self=shift;
+    push @{$self->{"children"}}, @_;
+}
+sub GetChildrenList {
+    my $self=shift;
+    return @{$self->{"children"}};
+}
+sub DeleteChild {
+    my $self=shift;
+    my $todelete=shift;
+    my @newchidren=grep {$_ != $todelete} $self->GetChildrenList();
+    $self->{"children"}=\@newchidren;
+    return;
+#    my $children=$self->{"children"};
+#    my($i);
+#    for ($i=0; $i<=$#$children; $i++) { 
+#	if ($children->[$i]==$todelete) {
+#	    splice(@{$children}, $i, 1);
+#	    return;
+#	}
+#    }
+}
+sub NbChildren {
+    my $self=shift;
+    return scalar(@{$self->{"children"}});
+}
+sub GetChild {
+    my $self=shift;
+    my $index=shift;
+
+    return $self->{"children"}->[$index];
+}
+sub ForgetChildren {
+    my $self=shift;
+    $self->{"children"}=[];
+}
+
+sub SetSequence {
+    my($self)=shift;
+    my($sequence)=shift;
+    $self->{"sequence"}=$sequence;
+}
+sub GetSequence {
+    my($self)=shift;
+    return $self->{"sequence"};
+}
+
+sub AddApo {
+    my $self=shift;
+    my $site_sens=shift;
+    my $apo=$site_sens->GetSiteNb();
+    if (exists($self->{"apo"}->{$apo})) {
+	die "Adding aleardy present apo $apo in ", $self->Name(),"\n";
+    }
+    $self->Debug("Adding apo $apo in ", $self->Name(),"\n");
+    $self->{"apo"}->{$apo}=$site_sens;
+    $site_sens->_AddNode($self);
+}
+sub NbApo {
+    my $self=shift;
+    return scalar (keys(%{$self->{"apo"}}));
+}
+sub GetApoIndexList {
+    my $self=shift;
+    return keys(%{$self->{"apo"}});
+}
+sub GetApoList {
+    my $self=shift;
+    return values(%{$self->{"apo"}});
+}
+sub GetApo {
+    my $self=shift;
+    my $apo=shift;
+    return $self->{"apo"}->{$apo};
+}
+sub DeleteAllApo {
+    my $self=shift;
+    foreach my $apo ($self->GetApoList()) {
+	$apo->_DeleteNode($self);
+    }
+    $self->{"apo"}={};
+}
+
+##Modifie Claire
+sub NbApoStep {
+    my $self=shift;
+    my($compteur)=0;
+    my($apo);
+    foreach $apo ($self->GetApoList()) {
+	$compteur+=$apo->GetStep();
+    }
+    return $compteur;
+}
+
+sub Name {
+    my $self=shift;
+    if (defined ($self->{"label"})) {
+	return $self->{"label"};
+    } else {
+	return $self->GetId();
+    }
+}
+
+#Faire meme chose pour case, control et br_len (HasBrLen), level, high, sequence, label, oldfather
+#Faire oldchildren (GetOldChildren, SetOldChildren, ) 
+
+
+1;

trunk/ALTree/Sens.pm

+package ALTree::Sens;
+
+################################################################
+################################################################
+####################### Sens          ##########################
+################################################################
+################################################################
+
+use base 'ALTree::Base';
+
+#Creation d'une structure
+sub New { # [classe] sens  
+    my $class=shift;
+    my $self={};
+    my $sens=shift;
+    bless($self, $class);
+    
+    if ($sens !~ /^\s*([0-9a-zA-Z?]+)\s*[-=_]*[>]\s*([0-9a-zA-Z?]+)\s*$/) {
+	die "Invalid sens : $sens\n";
+    }
+
+    $self->_init("start" => $1, "end" => $2, "switch" => 0, @_);
+    return $self;
+}
+
+#Creation d'une structure
+sub NewRev { # [classe] sens  
+    my $class=shift;
+    my $self={};
+    my $sens=shift;
+    bless($self, $class);
+    
+    $self->_init("start" => $sens->{"end"}, 
+		 "end" => $sens->{"start"}, "switch" => 0, @_);
+    return $self;
+}
+
+sub _makeLabel {
+    my $self=shift;
+    my $start=shift;
+    my $end=shift;
+
+    if ($self->{'switch'}) {
+	return $end."-->".$start;
+    } else {
+	return $start."-->".$end;
+    }
+}
+
+sub GetLabel {
+    my $self=shift;
+    return $self->_makeLabel($self->{"start"}, $self->{"end"});
+}
+
+sub GetLabelRev {
+    my $self=shift;
+    return $self->_makeLabel($self->{"end"}, $self->{"start"});
+}
+
+sub Switch {
+    my $self=shift;
+
+    $self->{'switch'} = 1-$self->{'switch'};
+}
+
+1;

trunk/ALTree/Site.pm

+package ALTree::Site;
+
+################################################################
+################################################################
+####################### Site          ##########################
+################################################################
+################################################################
+
+use base 'ALTree::Base';
+use ALTree::Sens;
+
+# Structure Site
+#   "site_nb" -> Integer
+#   "sens_struct" -> Hash of ('sens_label' -> SiteSens)
+
+sub InitSite { # [obj] $site
+    my $self=shift;
+    my $site_nb=shift;
+    $self->_init("site_nb" => $site_nb, "sens_struct" => {});
+    $self->Debug("creating Site $site_nb\n");
+}
+
+sub GetSiteNb {
+    my $self=shift;
+    return $self->{"site_nb"};
+}
+
+sub HasSensIndex {
+    my $self=shift;
+    my $sens=shift;
+
+    return exists($self->{"sens_struct"}->{$sens->GetLabel()});
+}
+sub AddSens {
+    my $self=shift;
+    my $sens=shift;
+
+    my($ref_site_sens)=$self->NewSens($sens);
+    $self->{"sens_struct"}->{$sens->GetLabel()}=$ref_site_sens;
+    my $sensRev=ALTree::Sens->NewRev($sens);
+    my($ref_site_sens_rev)=$self->NewSens($sensRev);
+    $self->{"sens_struct"}->{$sensRev->GetLabel()}=$ref_site_sens_rev;
+    ALTree::SiteSens::LinkRev($ref_site_sens, $ref_site_sens_rev);
+}
+sub GetSens {
+    my $self=shift;
+    my $sens=shift;
+
+    return $self->{"sens_struct"}->{$sens->GetLabel()};
+}
+sub ProvideSens {
+    my $self=shift;
+    my $sens=shift;
+    if (not $self->HasSensIndex($sens)) {
+	$self->AddSens($sens);
+    } # creation du hash ref_site_sens et d'une ref dessus
+    return $self->GetSens($sens);
+}
+sub GetSensIndexList {
+    my $self=shift;
+    return keys(%{$self->{"sens_struct"}});
+}
+sub GetSensList {
+    my $self=shift;
+    return values(%{$self->{"sens_struct"}});
+}
+
+sub NewSens {
+    my $self=shift;
+    my $sens=shift;
+    
+    die "This method needs to be overwriten\n";
+}
+
+1;

trunk/ALTree/SiteCollection.pm

+package ALTree::SiteCollection;
+
+################################################################
+################################################################
+####################### SiteCollection  ########################
+################################################################
+################################################################
+
+use base 'ALTree::Base';
+
+sub InitSiteCollection {
+    my $self=shift;
+    $self->_init("sites" => {});
+}
+
+sub AddSite {
+    my $self=shift;
+    my $site=shift; 
+    my $site_nb;
+    $site_nb=$site->GetSiteNb();
+    $self->{"sites"}->{$site_nb}=$site;
+}
+sub GetSite {
+    my $self=shift;
+    my $site_nb=shift;
+    #my $site=$self->{"sites"}->{$site_nb};
+    #if (not defined($site)) {
+    #die "The site number $site_nb does not exist";
+    #}
+    return $self->{"sites"}->{$site_nb};
+}
+sub HasSiteIndex {
+    my $self=shift;
+    my $site_nb=shift;
+    return exists($self->{"sites"}->{$site_nb});
+}
+sub GetSitesList {
+    my $self=shift;
+    return values(%{$self->{"sites"}});
+}
+1;

trunk/ALTree/SitePerForet.pm

+package ALTree::SitePerForet;
+
+################################################################
+################################################################
+####################### SitePerForet  ##########################
+################################################################
+################################################################
+
+use base qw(ALTree::Base ALTree::Site);
+use sort '_quicksort';
+
+sub New { # [classe] site_nb
+    my $class=shift;
+    my $self={};
+    my $site_nb=shift;
+    bless($self, $class);
+    $self->InitSite($site_nb);
+    $self->_init("nb_mut" => 0, @_);
+    return $self;
+}
+
+sub NewSens {
+    my $self=shift;
+    my $sens=shift;
+    return ALTree::SiteSensPerForet->New($sens, $self);
+}
+
+sub mysort {
+    my $func=shift;
+    my @tab=@_;
+
+    my @tab2=();
+
+    if (scalar(@tab) == 0) {
+	return @tab2;
+    }
+ #   print "Site ", $tab[0]->GetSiteNb(), "\n";
+#    foreach my $sens (@tab) {
+#	print $sens->GetVi()," ", $sens->GetSensLabel(), "\n";
+#    }
+    my $size;
+    while (($size=scalar(@tab)) > 1) {
+	my $min=0;
+	for(my $i=1; $i<$size; $i++) {
+	    if ($func->($tab[$i], $tab[$min]) < 0) {
+		$min=$i;
+	    }
+	}
+	push @tab2, $tab[$min];
+#	print "min=$min adding ",$tab[$min]->GetSensLabel(),"\n";
+#	print "Old tab\n";
+#	foreach my $sens (@tab) {
+#	    print $sens->GetVi()," ", $sens->GetSensLabel(), "\n";
+#	}
+	#my @tab3=splice(@tab, $min, 1);
+	my @tab3;
+	for (my $i=0; $i<=$#tab; $i++) {
+	    if ($i != $min) {
+		push @tab3, $tab[$i];
+	    }
+	}
+#	print "NewTab\n";
+#	foreach my $sens (@tab3) {
+#	    print $sens->GetVi()," ", $sens->GetSensLabel(), "\n";
+#	}
+	@tab=@tab3;
+    }
+    push @tab2, $tab[0];
+#    print "Site ", $tab[0]->GetSiteNb(), " result\n";
+#    foreach my $sens (@tab2) {
+#	print $sens->GetVi()," ", $sens->GetSensLabel(), "\n";
+#    }
+    return @tab2;
+
+}
+
+sub functrie {
+    my $a=shift;
+    my $b=shift;
+    #print "a=$a, b=$b\n";
+    #if (not defined($a)) {
+	#return 0;
+	#print "a=$a b=$b\n";
+	#die "beauk\n";
+    #}
+   # my $value=($b->GetVi() <=> $a->GetVi());
+   # print "Comparing ", $b->GetVi(), " and ", $a->GetVi(), " = $value\n";
+   # return $value;#
+    return ($b->GetVi() <=> $a->GetVi());
+}
+
+sub _EnsureViMax {
+    my($self)=shift;
+    if (not exists ($self->{"V_i_max"})) {
+	my @toto=$self->GetSensList();
+	#print "site: ", $self->GetSiteNb(), "\n";
+	#for my $ref (@toto) {
+	    #print "Ref: $ref\n";
+	#}
+	my @tab_trie=mysort \&functrie, @toto;
+	#for my $ref (@toto) {
+	    #print "Ref: $ref\n";
+	#}
+
+	#my @tab3;
+	#for (my $i=0; $i<=$#toto; $i++) {
+	#    print "value $i = $toto[$i]\n";
+	#    push @tab3, $toto[$i];
+	#}
+
+	#my @tab_trie=sort {
+	#    $b->GetVi() <=> $a->GetVi();
+	#} @tab3;
+	
+	my($Vimax)=$self->{"V_i_max"}=$tab_trie[0]->GetVi();
+	for (my $i=1; $i<=$#tab_trie; $i++) {
+	    if ($tab_trie[$i]->GetVi() !=  $Vimax) {
+		splice(@tab_trie, $i); # Elimine tout  partir de indice $i
+		last;
+	    }
+	}
+	#print "Vi max : ", $Vimax, "\n";
+	$self->{"V_i_max_tab"}=\@tab_trie;
+    }
+}
+sub GetViMax {
+    my $self=shift;
+    $self->_EnsureViMax();
+    return $self->{"V_i_max"};
+}   
+sub GetViMaxSens {
+    my($self)=shift;
+    my($index)=shift;
+    $self->_EnsureViMax();
+    return $self->{"V_i_max_tab"}->[$index];
+}
+#NbViMaxSens() => nombre sitesensperforet
+sub NbViMaxSens {
+    my($self)=shift; 
+    $self->_EnsureViMax();
+    return (scalar @{$self->{"V_i_max_tab"}});
+}
+#GetViMaxSensList() => array de tous les sitessensperforet (avec Vi==ViMax)
+sub GetViMaxSensList {
+    my($self)=shift; 
+    $self->_EnsureViMax();
+    return @{$self->{"V_i_max_tab"}}; # Comment a marche ce truc?
+}
+ 
+1;

trunk/ALTree/SitePerTree.pm

+package ALTree::SitePerTree;
+
+################################################################
+################################################################
+####################### SitePerTree   ##########################
+################################################################
+################################################################
+
+use base qw(ALTree::Base ALTree::Site);
+use ALTree::SiteSensPerTree;
+
+# Structure SitePerTree
+#   "site_nb" -> Integer
+#   "sens_struct" -> Hash of ('sens_label' -> SiteSens)
+#   "nb_mut" -> Integer
+
+sub New { # [classe] site_nb
+    my $class=shift;
+    my $self={};
+    my $site_nb=shift;
+    bless($self, $class);
+    $self->InitSite($site_nb);
+    $self->_init("nb_mut" => 0, @_);
+    return $self;
+}
+
+sub NewSens {
+    my $self=shift;
+    my $sens=shift;
+    return ALTree::SiteSensPerTree->New($sens, $self);
+}
+
+1;

trunk/ALTree/SiteSens.pm

+package ALTree::SiteSens;
+
+################################################################
+################################################################
+####################### SiteSens        ########################
+################################################################
+################################################################
+
+use base 'ALTree::Base';
+
+# Structure SiteSens
+#   "site_struct" -> Site
+#   "sens" -> Sens
+#   "rev" -> SiteSens
+
+sub InitSiteSens { # [Obj] Sens_label SiteSens
+    my $self=shift;
+    my $sens=shift;
+    my $site=shift;
+
+    $self->_init("sens" => $sens, "site_struct" => $site, @_);
+    $self->Debug("creating SiteSens $sens\n");
+}
+
+sub GetSensLabel {
+    my $self=shift;
+    return $self->{"sens"}->GetLabel();
+}
+
+sub GetSensStruct {
+    my $self=shift;
+    return $self->{"sens"};
+}
+
+sub GetSite {
+    my $self=shift;
+    return $self->{"site_struct"};
+}
+
+sub GetSiteNb {
+    my $self=shift;
+    return $self->GetSite()->GetSiteNb();
+}
+
+sub LinkRev {
+    my $siteSens=shift;
+    my $siteSensRev=shift;
+
+    $siteSens->{'rev'}=$siteSensRev;
+    $siteSensRev->{'rev'}=$siteSens;
+}
+
+sub GetSensRev {
+    my $self=shift;
+
+    if (not defined($self->{'rev'})) {
+	die "LinkRev not called\n";
+    }
+    return $self->{'rev'};
+}    
+
+1;

trunk/ALTree/SiteSensPerForet.pm

+package ALTree::SiteSensPerForet;
+
+################################################################
+################################################################
+####################### SiteSensPerForet #######################
+################################################################
+################################################################
+
+use base qw(ALTree::Base ALTree::SiteSens);
+
+# Structure SiteSensPerForet
+#   "site_struct" -> Site
+#   "sens_label" -> String
+#   "V_i" -> Float
+
+sub New { # [classe] sens_label site_struct
+    my $class=shift;
+    my $self={};
+    my $sens=shift;
+    my $site=shift;
+    bless($self, $class);
+    $self->InitSiteSens($sens, $site);
+    $self->_init("V_i" => 0, @_);
+    return $self;
+}
+
+sub PlusVi {
+    my $self=shift;
+    my $value=shift;
+    $self->{"V_i"} += $value;
+}
+sub GetVi {
+    my $self=shift;
+    return $self->{"V_i"};
+}    
+
+1;

trunk/ALTree/SiteSensPerTree.pm

+package ALTree::SiteSensPerTree;
+
+################################################################
+################################################################
+####################### SiteSensPerTree ########################
+################################################################
+################################################################
+
+use base qw(ALTree::Base ALTree::SiteSens);
+
+# Structure SiteSens
+#   "site_struct" -> Site
+#   "sens_label" -> String
+#   "node_list" -> Array of (Node)
+#    "m_it" -> Interger # nb mutation of this chage in the tree
+#    "R_it" -> Interger # nb co-mutations of this change with character S
+#    "V_it" -> Integer  # (R_it-E_it)/sqrt(E_it)
+
+sub New { # [classe] sens_label site_struct
+    my $class=shift;
+    my $self={};
+    my $sens=shift;
+    my $site=shift;
+    bless($self, $class);
+    $self->InitSiteSens($sens, $site);
+    $self->_init("node_list" => [], "R_it" => 0, @_);
+    return $self;
+}
+
+# Appel par Node->AddApo
+sub _AddNode {
+    my $self=shift;
+    my $node=shift;
+    if (exists($self->{"m_it"})) {
+	die "_AddNode called after GetMit";	
+    }
+    push @{$self->{"node_list"}}, $node;
+}
+# Appel par Node->DeleteAllApo
+sub _DeleteNode {
+    my $self=shift;
+    my $node=shift;
+    my @new_node_list=grep { $_ != $node } @{$self->{"node_list"}};
+    if (scalar(@new_node_list)+1 != $self->NbNodes()) {
+	die "Error while removing a node from a SiteSensPerTree";
+    }
+    $self->{"node_list"}=\@new_node_list;
+}
+sub NbNodes {
+    my $self=shift;
+
+    return scalar (@{$self->{"node_list"}});
+}
+sub GetNode {
+    my $self=shift;
+    my $index=shift;
+    return $self->{"node_list"}->[$index];
+}
+sub GetNodesList {
+    my $self=shift;
+    return @{$self->{"node_list"}};
+}
+
+sub SetStep {
+    my $self=shift;
+    my $step=shift;
+    $self->{"Step"}=$step;
+}
+
+sub GetStep {
+    my $self=shift;
+    my $step=shift;
+    return $self->{"Step"};
+}
+
+sub GetMit {
+    my($self)=shift;
+    if (not exists($self->{"m_it"})) {
+	$self->{"m_it"}=$self->NbNodes();
+    }
+    return ($self->{"m_it"});
+}
+
+sub IncRit {
+    my($self)=shift;
+    if (exists($self->{"V_it"})) {
+	die "IncRit called after GetVit";	
+    }    
+    $self->{"R_it"}++;
+}
+sub GetRit {
+    my($self)=shift;
+    return ($self->{"R_it"});
+}
+
+sub GetVit {
+    my $self=shift;
+
+    my $Rit=$self->GetRit();
+    my $Eit=$self->GetEit();
+
+    if (not exists($self->{"V_it"})) {
+	if ($Eit == 0) {
+	    $self->{"V_it"}=0;
+	} else {
+	    $self->{"V_it"}=($Rit-$Eit)/sqrt($Eit);
+	}
+    }
+    return ($self->{"V_it"});
+}
+
+sub SetEit {
+    my($self)=shift;
+    my($value)=shift;
+    if (exists($self->{"V_it"})) {
+	die "SetEit called after GetVit";	
+    }    
+    $self->{"E_it"}=$value;
+}
+
+sub GetEit {
+    my($self)=shift;
+    return ($self->{"E_it"});
+}
+1;

trunk/ALTree/Tree.pm

+package ALTree::Tree;
+
+################################################################
+################################################################
+####################### Tree            ########################
+################################################################
+################################################################
+
+use base qw(ALTree::Base ALTree::SiteCollection);
+
+# "nodes" -> Hash of ('id' -> Node)
+# "sites" -> Hash of ('site_nb' -> SitePerTree)
+
+sub New { # [classe]
+    my $class=shift;
+    my $self={};
+    bless($self, $class);
+    $self->InitSiteCollection();
+    $self->_init("nodes" => {}, @_);
+    $self->Debug("creating Tree\n");
+    return $self;
+}
+
+sub AddNode {
+    my $self=shift;
+    my $node=shift;
+
+    my $id=$node->GetId();
+    $self->{"nodes"}->{$id}=$node;
+}
+sub GetNode {
+    my $self=shift;
+    my $id=shift;
+    if (defined ($self->{"nodes"}->{$id})) {
+	return $self->{"nodes"}->{$id};
+    } else {
+	return undef;
+    }
+}    
+sub HasNodeIndex {
+    my $self=shift;
+    my $id=shift;
+    return exists($self->{"nodes"}->{$id});
+}
+sub GetNodesIndexList {
+    my $self=shift;
+    return keys(%{$self->{"nodes"}});
+}
+sub GetNodesList {
+    my $self=shift;
+    return values(%{$self->{"nodes"}});
+}
+
+#sub _SetOutgroup {
+#    my $self=shift;
+#    my($outgroup)=shift;
+#    $self->{"outgroup"}=$outgroup;
+#}
+sub GetOutgroup {
+     my $self=shift;
+     if (not exists($self->{"outgroup"})) {
+       FIND: {
+	   foreach my $clef ($self->GetNodesIndexList()) {
+	       if ($clef eq "H000") {
+		   $self->{"outgroup"}=$self->GetNode($clef);
+		   last FIND;
+	       }
+	   }
+	   die "No outgroup corresponding to id=H000 identified in the tree\n";
+       }
+     }
+     return ($self->{"outgroup"});
+}
+
+sub SetNbBrNonNulle {
+    my $self=shift;
+    my($value)=shift;
+    $self->{"nb_br_non_nulle"}=$value;
+} 
+sub GetNbBrNonNulle {
+    my $self=shift;
+    return ($self->{"nb_br_non_nulle"});
+}
+
+sub _SetRoot {
+    my $self=shift;
+
+    my(@roots);
+    #print STDERR "nodes: \n";
+    foreach my $node ($self->GetNodesList()) {
+	if (not $node->HasFather()) {
+	    push @roots, $node;
+	    #print STDERR "root: ", $node->Name(), "\n";
+	}	    
+	#else { print STDERR "node: ", $node->Name(), " father ", $node->GetFather()->Name(), "\n"; }
+    }
+    #Verification that there is only one root
+    if (scalar(@roots) !=1) {
+	die "error in the tree: ", scalar(@roots), " roots for the tree!\n";
+    }
+    $self->{"root"}=$roots[0];
+} 
+sub GetRoot {
+    my $self=shift;
+    if (not exists($self->{"root"})) {
+	$self->_SetRoot();
+    }
+    return ($self->{"root"});
+}
+
+sub ChangeRoot {
+    my $self=shift;
+    my $newroot=shift;
+
+    if ($newroot == $self->GetRoot()) {
+	return;
+    }
+    my $oldfather=$newroot->GetFather();
+    if (not defined($oldfather)) {
+	die ("Non root node of the tree has no father\n"
+	     ."Do you use a node of the correct tree ?");
+    }
+    $self->ChangeRoot($oldfather);
+    
+    # Some integrity tests...
+    if ($oldfather->NbApo() != 0) {
+	die "Root has apomorphies !";
+    }
+    if ($oldfather->GetBrLen() != 0) {
+	die "Root has non null BrLen !";
+    }
+    
+    foreach my $apo ($newroot->GetApoList()) {
+	$oldfather->AddApo($apo->GetSensRev());
+    }
+    $newroot->DeleteAllApo();
+    #print "Setting BRLen to ",$newroot->GetBrLen()," for ", $oldfather->Name()," from ", $newroot->Name(),"\n";
+    $oldfather->SetBrLen($newroot->GetBrLen());
+    $newroot->SetBrLen(0);
+
+    $oldfather->SetFather($newroot);
+    $newroot->SetFather(undef);
+    $newroot->AddChild($oldfather);
+    $oldfather->DeleteChild($newroot);
+
+    $self->{"root"}=$newroot;
+}
+
+1;

trunk/ALTree/Utils.pm

+package ALTree::Utils;
+use Pod::Usage;
+
+use strict;
+
+BEGIN {
+    use Exporter   ();
+    use vars       qw($VERSION @ISA @EXPORT @EXPORT_OK %EXPORT_TAGS);
+    
+    # set the version for version checking
+    #$VERSION     = 1.00;
+    # if using RCS/CVS, this may be preferred
+    #$VERSION = do { my @r = (q$Revision$ =~ /\d+/g); sprintf "%d."."%02d" x $#r, @r }; # must be all one line, for MakeMaker
+    
+    @ISA         = qw(Exporter);
+    @EXPORT      = qw(); #(&func1 &func2 &func4);
+    %EXPORT_TAGS = ( );     # eg: TAG => [ qw!name1 name2! ],
+    
+    # your exported package globals go here,
+    # as well as any optionally exported functions
+    @EXPORT_OK   = qw(&erreur &internal_error);
+}
+use vars      @EXPORT_OK;
+
+
+sub erreur
+{
+    my $msg=shift;
+    my $use_pod=shift;
+
+    if (not defined($use_pod) or $use_pod) {
+	pod2usage("Error: ".$msg);
+    } else {
+	print STDERR "Error: ".$msg;
+	exit 1;
+    }
+}
+
+sub internal_error
+{
+    my $msg=shift;
+
+    die("Internal error: $msg\n".
+	"Please, report this bug (with all is needed to reproduce it) to:\n".
+	"bardel\@vjf.inserm.fr\n");
+}
+
+1;
+
+package DataType;
+use constant SNP => 0;
+use constant DNA => 1;
+
+package PhylProg;
+use constant PHYLIP => 0;
+use constant PAUP => 1;
+use constant PAML => 2;
+
+package ANC;
+use constant Rooted => 0;
+use constant Unrooted => 1;
+use constant OutGroup => 2;
+
+no strict;
+
+@Name=("rooted using an ancestral sequence",
+       "unrooted", "rooted using an outgroup");

trunk/CUtils/.cvsignore

+Makefile

trunk/CUtils/CUtils.xs

+#include "EXTERN.h"
+#include "perl.h"
+#include "XSUB.h"
+
+#include "ppport.h"
+
+#include <c_sources/fisher.h>
+#include <c_sources/chisq.h>
+#include <c_sources/double_permutation.h>
+
+#include "const-c.inc"
+
+MODULE = ALTree::CUtils		PACKAGE = ALTree::CUtils		
+
+INCLUDE: const-xs.inc
+
+double
+bilateral(a, b, c, d)
+	double	a
+	double	b
+	double	c
+	double	d
+
+double
+critchi(p, df)
+	double	p
+	int	df
+
+double
+left(a, b, c, d)
+	double	a
+	double	b
+	double	c
+	double	d
+
+double
+pochisq(x, df)
+	double	x
+	int	df
+
+double
+right(a, b, c, d)
+	double	a
+	double	b
+	double	c
+	double	d
+
+SV *
+double_permutation(nb_sample, nb_chi2, data)
+	int nb_sample
+	int nb_chi2
+	SV * data
+    INIT:
+        matrice_t mat;
+        replicat_t rep;
+        int i,j,n;
+        datatype_t min;
+	HV * rh;
+        AV * ra;
+    CODE:
+	//fprintf(stderr, "\nSTART(%i, %i)\n", nb_sample, nb_chi2);
+        if ((nb_sample <= 0)
+	    || (nb_chi2 <= 0)
+	    || (!SvROK(data))
+            || (SvTYPE(SvRV(data)) != SVt_PVAV)
+            || (av_len((AV *)SvRV(data)) != nb_sample*nb_chi2-1))
+        {
+		XSRETURN_UNDEF;
+	}
+        rh = (HV *)sv_2mortal((SV *)newHV());
+	
+        mat=alloc_matrice(nb_sample, nb_chi2);
+	rep=alloc_replicat(nb_chi2);
+	
+        n=0;
+	for (i=0; i<nb_sample; i++) {
+		for (j=0; j<nb_chi2; j++) {
+			//fprintf(stderr, "[%i][%i](%i)...\n", i, j, n);
+                        /* Attention: on place un rplicat par colonne
+                         * (et pas par ligne) */
+			mat[j][i]=SvNV(*av_fetch((AV *)SvRV(data), n, 0));
+			//fprintf(stderr, "[%i][%i](%i)=%lf\n", i, j, n,mat[j][i]);
+			n++;
+		}
+	}
+	min=calcul(nb_sample, nb_chi2, mat, rep);
+
+	hv_store(rh, "pmin", 4, newSVnv(min), 0);
+
+	ra = (AV *)sv_2mortal((SV *)newAV());
+	for (j=0; j<nb_chi2; j++) {
+		av_push(ra, newSVnv(rep[j]));
+		//fprintf(stderr, "chi2[%i]=%lf\n", j ,rep[j]);
+	}
+	
+	hv_store(rh, "chi2", 4, newRV((SV *)ra), 0);
+
+        free_matrice(mat, nb_sample, nb_chi2);
+        free_replicat(rep);
+
+	RETVAL = newRV((SV *)rh);
+    OUTPUT:
+        RETVAL

trunk/CUtils/Makefile.PL

+use 5.008;
+use ExtUtils::MakeMaker;
+# See lib/ExtUtils/MakeMaker.pm for details of how to influence
+# the contents of the Makefile that is written.
+WriteMakefile(
+    NAME              => 'ALTree::CUtils',
+    VERSION_FROM      => 'lib/ALTree/CUtils.pm', # finds $VERSION
+    PREREQ_PM         => {}, # e.g., Module::Name => 1.1
+    ($] >= 5.005 ?     ## Add these new keywords supported since 5.005
+      (ABSTRACT_FROM  => 'lib/ALTree/CUtils.pm', # retrieve abstract from module
+       AUTHOR         => 'Claire Bardel <bardel@vjf.inserm.fr>') : ()),
+    LIBS              => ['-lm '], # e.g., '-lm'
+    DEFINE            => '', # e.g., '-DHAVE_SOMETHING'
+    INC               => '-I.', # e.g., '-I. -I/usr/include/other'
+	# Un-comment this if you add C files to link with later:
+    # OBJECT            => '$(O_FILES)', # link all the C files too
+    MYEXTLIB          => 'c_sources/libaltree-cutils$(LIB_EXT)',
+);
+if  (eval {require ExtUtils::Constant; 1}) {
+  # If you edit these definitions to change the constants used by this module,
+  # you will need to use the generated const-c.inc and const-xs.inc
+  # files to replace their "fallback" counterparts before distributing your
+  # changes.
+  my @names = (qw());
+  ExtUtils::Constant::WriteConstants(
+                                     NAME         => 'ALTree::CUtils',
+                                     NAMES        => \@names,
+                                     DEFAULT_TYPE => 'IV',
+                                     C_FILE       => 'const-c.inc',
+                                     XS_FILE      => 'const-xs.inc',
+                                  );
+
+}
+else {
+  use File::Copy;
+  use File::Spec;
+  foreach my $file ('const-c.inc', 'const-xs.inc') {
+    my $fallback = File::Spec->catfile('fallback', $file);
+    copy ($fallback, $file) or die "Can't copy $fallback to $file: $!";
+  }
+}
+sub MY::postamble {
+'
+$(MYEXTLIB): c_sources/Makefile
+	$(MAKE) -C c_sources $(PASSTHRU)
+';
+}
+
+