use Bio::PopGen::Statistics;
  use Bio::AlignIO;
  use Bio::PopGen::IO;
  use Bio::PopGen::Simulation::Coalescent;

  my $sim = Bio::PopGen::Simulation::Coalescent->new( -sample_size => 12);

  my $tree = $sim->next_tree;

  $sim->add_Mutations($tree,20);

  my $stats = Bio::PopGen::Statistics->new();
  my $individuals = [ $tree->get_leaf_nodes];
  my $pi = $stats->pi($individuals);
  my $D  = $stats->tajima_D($individuals);

  # Alternatively to do this on input data from
  # See the tests in t/PopGen.t for more examples
  my $parser = Bio::PopGen::IO->new(-format => 'prettybase',
                                   -file   => 't/data/popstats.prettybase');
  my $pop = $parser->next_population;
  # Note that you can also call the stats as a class method if you like
  # the only reason to instantiate it (as above) is if you want
  # to set the verbosity for debugging
  $pi     = Bio::PopGen::Statistics->pi($pop);
  $theta  = Bio::PopGen::Statistics->theta($pop);

  # Pi and Theta also take additional arguments,
  # see the documentation for more information

  use Bio::PopGen::Utilities;
  use Bio::AlignIO;

  my $in = Bio::AlignIO->new(-file   => 't/data/t7.aln',
                            -format => 'clustalw');
  my $aln = $in->next_aln;
  # get a population, each sequence is an individual and 
  # for the default case, every site which is not monomorphic
  # is a 'marker'.  Each individual will have a 'genotype' for the
  # site which will be the specific base in the alignment at that
  # site

  my $pop = Bio::PopGen::Utilities->aln_to_population(-alignment => $aln);

    eval { require Text::NSP::Measures::2D::Fisher2::twotailed };
    if( $@ ) { $has_twotailed = 0; }
    else { $has_twotailed = 1;

     my ($self,$ingroup,$outgroup) = @_;

    my ($seg_sites,$n,$ancestral,$derived) = (0,0,0,0);
    if( ref($ingroup) =~ /ARRAY/i ) {
	$n = scalar @$ingroup;
	# pi - all pairwise differences 
	$seg_sites   = $self->segregating_sites_count($ingroup);
    } elsif( ref($ingroup) && 
	     $ingroup->isa('Bio::PopGen::PopulationI')) {
	$n = $ingroup->get_number_individuals;
	$seg_sites   = $self->segregating_sites_count($ingroup);
    } else { 
	$self->throw("expected an array reference of a list of Bio::PopGen::IndividualI OR a Bio::PopGen::PopulationI object to fu_and_li_D");
	return 0;
    }
    
    if( $seg_sites <= 0 ) { 
	$self->warn("mutation total was not > 0, cannot calculate a Fu and Li D");
	return 0;
    }

    if( ! defined $outgroup ) {
	$self->warn("Need to provide either an array ref to the outgroup individuals or the number of external mutations");
	return 0;
    } elsif( ref($outgroup) ) {
	($ancestral,$derived) = $self->derived_mutations($ingroup,$outgroup);
	$ancestral = 0 unless defined $ancestral;
    } else { 
	$ancestral = $outgroup;
    }
   
    return $self->fu_and_li_D_counts($n,$seg_sites,
				     $ancestral,$derived);

    my ($self,$n,$seg_sites, $external_mut) = @_;
    my $a_n = 0;
    if( $n <= 3 ) {
	$self->warn("n is $n, too small, must be > 3\n");
	return;
    }
    for(my $k= 1; $k < $n; $k++ ) {
	$a_n += ( 1 / $k );
    }
    my $b = 0;
    for(my $k= 1; $k < $n; $k++ ) {
        $b += ( 1 / $k**2 );
    }

    my $c = 2 * ( ( ( $n * $a_n ) - (2 * ( $n -1 ))) /
                  ( ( $n - 1) * ( $n - 2 ) ) );

    my $v = 1 + ( ( $a_n**2 / ( $b + $a_n**2 ) ) * 
		  ( $c - ( ( $n + 1) /
			   ( $n - 1) ) ));
    
    my $u = $a_n - 1 - $v;

    ($seg_sites - $a_n * $external_mut) / 
	sqrt( ($u * $seg_sites) + ($v * $seg_sites*$seg_sites));

    my ($self,$individuals) = @_;

    my ($seg_sites,$n,$singletons);
    if( ref($individuals) =~ /ARRAY/i ) {
	$n = scalar @$individuals;
	$seg_sites   = $self->segregating_sites_count($individuals);
	$singletons  = $self->singleton_count($individuals);
    } elsif( ref($individuals) && 
	     $individuals->isa('Bio::PopGen::PopulationI')) {
	my $pop = $individuals;
	$n = $pop->get_number_individuals;
	$seg_sites   = $self->segregating_sites_count($pop);
	$singletons  = $self->singleton_count($pop);
    } else { 
	$self->throw("expected an array reference of a list of Bio::PopGen::IndividualI OR a Bio::PopGen::PopulationI object to fu_and_li_D_star");
	return 0;
    }

    return $self->fu_and_li_D_star_counts($n,$seg_sites, $singletons);

    my ($self,$n,$seg_sites, $singletons) = @_;
    my $a_n;
    for(my $k = 1; $k < $n; $k++ ) {
	$a_n += ( 1 / $k );
    }

    my $a1 = $a_n + 1 / $n;

    my $b = 0;
    for(my $k= 1; $k < $n; $k++ ) {
        $b += ( 1 / $k**2 );
    }

    my $c = 2 * ( ( ( $n * $a_n ) - (2 * ( $n -1 ))) /
                  ( ( $n - 1) * ( $n - 2 ) ) );

    my $d = $c + ($n -2) / ($n - 1)**2 +
	2 / ($n -1) * 
	( 1.5 - ( (2*$a1 - 3) / ($n -2) ) - 
	  1 / $n ); 
    
    my $v_star = ( ( ($n/($n-1) )**2)*$b + (($a_n**2)*$d) -
		 (2*( ($n*$a_n*($a_n+1)) )/(($n-1)**2)) )  /
		   (($a_n**2) + $b);

    my $u_star = ( ($n/($n-1))*
		   ($a_n - ($n/
			  ($n-1)))) - $v_star;


    return (($n / ($n - 1)) * $seg_sites - 
	    $a_n * $singletons) / 
	    sqrt( ($u_star * $seg_sites) + ($v_star * $seg_sites*$seg_sites));

    my ($self,$ingroup,$outgroup) = @_;
    my ($seg_sites,$pi,$n,$external,$internal);
    if( ref($ingroup) =~ /ARRAY/i ) {
	$n = scalar @$ingroup;
	# pi - all pairwise differences 
	$pi          = $self->pi($ingroup);  
	$seg_sites   = $self->segregating_sites_count($ingroup);
    } elsif( ref($ingroup) && 
	     $ingroup->isa('Bio::PopGen::PopulationI')) {
	$n = $ingroup->get_number_individuals;
	$pi          = $self->pi($ingroup);
	$seg_sites   = $self->segregating_sites_count($ingroup);
    } else { 
	$self->throw("expected an array reference of a list of Bio::PopGen::IndividualI OR a Bio::PopGen::PopulationI object to Fu and Li's F");
	return 0;
    }
    
    if( ! defined $outgroup ) {
	$self->warn("Need to provide either an array ref to the outgroup individuals or the number of external mutations");
	return 0;
    } elsif( ref($outgroup) ) {
	($external,$internal) = $self->derived_mutations($ingroup,$outgroup);
    } else { 
	$external = $outgroup;
    }
    $self->fu_and_li_F_counts($n,$pi,$seg_sites,$external);

    my ($self,$n,$pi,$seg_sites, $external) = @_;
    my $a_n = 0;
    for(my $k= 1; $k < $n; $k++ ) {
	$a_n += ( 1 / $k );
    }

    my $a1 = $a_n + (1 / $n );

    my $b = 0;
    for(my $k= 1; $k < $n; $k++ ) {
	$b += ( 1 / $k**2 );
    }

    my $c = 2 * ( ( ( $n * $a_n ) - (2 * ( $n -1 ))) / 
		  ( ( $n - 1) * ( $n - 2 ) ) );

    my $v_F = ( $c + ( (2*(($n**2)+$n+3)) / 
		       ( (9*$n)*($n-1) ) ) -
		(2/($n-1)) ) / ( ($a_n**2)+$b );

    my $u_F = ( 1 + ( ($n+1)/(3*($n-1)) )-
		( 4*( ($n+1)/(($n-1)**2) ))*
		($a1 - ((2*$n)/($n+1))) ) /
		$a_n - $v_F;

    # warn("$v_F vf $u_F uf n = $n\n");
    my $F = ($pi - $external) / ( sqrt( ($u_F*$seg_sites) +
					($v_F*($seg_sites**2)) ) );

    return $F;

    my ($self,$individuals) = @_;

    my ($seg_sites,$pi,$n,$singletons);
    if( ref($individuals) =~ /ARRAY/i ) {
	$n = scalar @$individuals;
	# pi - all pairwise differences 
	$pi          = $self->pi($individuals);  
	$seg_sites   = $self->segregating_sites_count($individuals);
	$singletons  = $self->singleton_count($individuals);
    } elsif( ref($individuals) && 
	     $individuals->isa('Bio::PopGen::PopulationI')) {
	my $pop = $individuals;
	$n = $pop->get_number_individuals;
	$pi          = $self->pi($pop);
	$seg_sites   = $self->segregating_sites_count($pop);
	$singletons  = $self->singleton_count($pop);
    } else { 
	$self->throw("expected an array reference of a list of Bio::PopGen::IndividualI OR a Bio::PopGen::PopulationI object to fu_and_li_F_star");
	return 0;
    }
    return $self->fu_and_li_F_star_counts($n,
					  $pi,
					  $seg_sites,
					  $singletons);

    my ($self,$n,$pi,$seg_sites, $singletons) = @_;
    if( $n <= 1 ) {
	$self->warn("N must be > 1\n");
	return;
    }
    if( $n == 2) { 
	return 0;
    } 

    my $a_n = 0;
    

    my $b = 0;
    for(my $k= 1; $k < $n; $k++ ) {
	$b += (1 / ($k**2));
	$a_n += ( 1 / $k );     # Eq (2)
    }
    my $a1 = $a_n + (1 / $n );

    # warn("a_n is $a_n a1 is $a1 n is $n b is $b\n");

    # From Simonsen et al (1995) instead of Fu and Li 1993
    my $v_F_star = ( (( 2 * $n ** 3 + 110 * $n**2 - (255 * $n) + 153)/
		      (9 * ($n ** 2) * ( $n - 1))) +
		     ((2 * ($n - 1) * $a_n ) / $n ** 2) -
		     (8 * $b / $n) ) / 
		     ( ($a_n ** 2) + $b );
    
    my $u_F_star = ((( (4* ($n**2)) + (19 * $n) + 3 - (12 * ($n + 1)* $a1)) /
		    (3 * $n * ( $n - 1))) / $a_n) - $v_F_star;

    # warn("vf* = $v_F_star uf* = $u_F_star n = $n\n");
    my $F_star = ( $pi - ($singletons*( ( $n-1) / $n)) ) /
	sqrt ( $u_F_star*$seg_sites + $v_F_star*$seg_sites**2);
    return $F_star;

    my ($self,$individuals) = @_;
    my ($seg_sites,$pi,$n);

    if( ref($individuals) =~ /ARRAY/i ) {
	$n = scalar @$individuals;
	# pi - all pairwise differences 
	$pi          = $self->pi($individuals);  
	$seg_sites = $self->segregating_sites_count($individuals);

    } elsif( ref($individuals) && 
	     $individuals->isa('Bio::PopGen::PopulationI')) {
	my $pop = $individuals;
	$n = $pop->get_number_individuals;
	$pi          = $self->pi($pop);
	$seg_sites = $self->segregating_sites_count($pop);
    } else { 
	$self->throw("expected an array reference of a list of Bio::PopGen::IndividualI OR a Bio::PopGen::PopulationI object to tajima_D");
	return 0;
    }
    $self->tajima_D_counts($n,$seg_sites,$pi);

    my ($self,$n,$seg_sites,$pi) = @_;
    my $a1 = 0; 
    for(my $k= 1; $k < $n; $k++ ) {
	$a1 += ( 1 / $k );
    }

     my $a2 = 0;
     for(my $k= 1; $k < $n; $k++ ) {
	 $a2 += ( 1 / $k**2 );
     }
    
    my $b1 = ( $n + 1 ) / ( 3* ( $n - 1) );
    my $b2 = ( 2 * ( $n ** 2 + $n + 3) ) / 
	     ( ( 9 * $n) * ( $n - 1) );
    my $c1 = $b1 - ( 1 / $a1 );
    my $c2 = $b2 - ( ( $n + 2 ) /
		     ( $a1 * $n))+( $a2 / $a1 ** 2);
    my $e1 = $c1 / $a1;
    my $e2 = $c2 / ( $a1**2 + $a2 );
    
    my $denom = sqrt ( ($e1 * $seg_sites) + (( $e2 * $seg_sites) * ( $seg_sites - 1)));
    return if $denom == 0;
    my $D = ( $pi - ( $seg_sites / $a1 ) ) / $denom;
    return $D;

    my ($self,$individuals,$numsites) = @_;
    my (%data,%marker_total,@marker_names,$n);

    if( ref($individuals) =~ /ARRAY/i ) {
	# one possible argument is an arrayref of Bio::PopGen::IndividualI objs
	@marker_names = $individuals->[0]->get_marker_names;
	$n = scalar @$individuals;

	# Here we are calculating the allele frequencies
	foreach my $ind ( @$individuals ) {
	    if( ! $ind->isa('Bio::PopGen::IndividualI') ) {
		$self->warn("Expected an arrayref of Bio::PopGen::IndividualI objects, this is a ".ref($ind)."\n");
		return 0;
	    }
	    foreach my $m ( @marker_names ) {
		foreach my $allele (map { $_->get_Alleles} 
				    $ind->get_Genotypes($m) ) {
		    $data{$m}->{$allele}++;
		    $marker_total{$m}++;
		}
	    }
	}
#	while( my ($marker,$count) =  each %marker_total ) {
#	    foreach my $c ( values %{$data{$marker}} ) {
#		$c /= $count;
#	    }
#	}
	# %data will contain allele frequencies for each marker, allele
    } elsif( ref($individuals) &&
	     $individuals->isa('Bio::PopGen::PopulationI') ) {
	my $pop = $individuals;
	$n = $pop->get_number_individuals;
	foreach my $marker( $pop->get_Markers ) {
	    push @marker_names, $marker->name;
	    #$data{$marker->name} = {$marker->get_Allele_Frequencies};
	    my @genotypes = $pop->get_Genotypes(-marker => $marker->name);
	    for my $al ( map { $_->get_Alleles} @genotypes ) {
	      $data{$marker->name}->{$al}++;
	      $marker_total{$marker->name}++;
	   }
	}
    } else {
	$self->throw("expected an array reference of a list of Bio::PopGen::IndividualI to pi");
    }
    # based on Kevin Thornton's code:
    # http://molpopgen.org/software/libsequence/doc/html/PolySNP_8cc-source.html#l00152
    # For now we assume that all individuals have the same markers
    my ($diffcount,$totalcompare) = (0,0);
    my $pi = 0;
    while ( my ($marker,$markerdat) = each %data ) {
      my $sampsize = $marker_total{$marker};
      my $ssh = 0;
      my @alleles = keys %$markerdat;
      if ( $sampsize > 1 ) {
	my $denom = $sampsize * ($sampsize - 1.0);
	foreach my $al ( @alleles ) {
	  $ssh += ($markerdat->{$al} * ($markerdat->{$al} - 1)) / $denom;
	}
	$pi += 1.0 - $ssh;
      }
    }
    $self->debug( "pi=$pi\n");
    if( $numsites ) {
	return $pi / $numsites;
    } else {
	return $pi;
    }

    my $self = shift;
    my ( $n, $seg_sites,$totalsites) = @_;
    if( ref($n) =~ /ARRAY/i ) {
	my $samps = $n;
	$totalsites = $seg_sites; # only 2 arguments if one is an array
	my %data;
	my @marker_names = $samps->[0]->get_marker_names;
	# we need to calculate number of polymorphic sites
	$seg_sites = $self->segregating_sites_count($samps);
	$n = scalar @$samps;

    } elsif(ref($n) &&
	    $n->isa('Bio::PopGen::PopulationI') ) {
	# This will handle the case when we pass in a PopulationI object
	my $pop = $n;
	$totalsites = $seg_sites; # shift the arguments over by one
	$n = $pop->haploid_population->get_number_individuals;
	$seg_sites = $self->segregating_sites_count($pop);
    }
    my $a1 = 0; 
    for(my $k= 1; $k < $n; $k++ ) {
	$a1 += ( 1 / $k );
    }    
    if( $totalsites ) { # 0 and undef are the same can't divide by them
	$seg_sites /= $totalsites;
    }
    if( $a1 == 0 ) { 
	return 0;
    } 
    return $seg_sites / $a1;

    my ($self,$individuals) = @_;

    my @inds;
    if( ref($individuals) =~ /ARRAY/ ) {
	@inds = @$individuals;
    } elsif( ref($individuals) && 
	     $individuals->isa('Bio::PopGen::PopulationI') ) {
	my $pop = $individuals;
	@inds = $pop->get_Individuals();
	unless( @inds ) { 
	    $self->warn("Need to provide a population which has individuals loaded, not just a population with allele frequencies");
	    return 0;
	}
    } else {
	$self->warn("Expected either a PopulationI object or an arrayref of IndividualI objects");
	return 0;
    }
    # find number of sites where a particular allele is only seen once

    my ($singleton_allele_ct,%sites) = (0);
    # first collect all the alleles into a hash structure
    
    foreach my $n ( @inds ) {
	if( ! $n->isa('Bio::PopGen::IndividualI') ) {
	    $self->warn("Expected an arrayref of Bio::PopGen::IndividualI objects, this is a ".ref($n)."\n");
	    return 0;
	}
	foreach my $g ( $n->get_Genotypes ) {
	    my ($nm,@alleles) = ($g->marker_name, $g->get_Alleles);
	    foreach my $allele (@alleles ) {
		$sites{$nm}->{$allele}++;
	    }
	}
    }
    foreach my $site ( values %sites ) { # don't really care what the name is
	foreach my $allelect ( values %$site ) { # 
            # find the sites which have an allele with only 1 copy
 	    $singleton_allele_ct++ if( $allelect == 1 );
	}
    }
    return $singleton_allele_ct;
}

# Yes I know that singleton_count and segregating_sites_count are
# basically processing the same data so calling them both is
# redundant, something I want to fix later but want to make things
# correct and simple first

   my ($self,$individuals) = @_;
   my $type = ref($individuals);
   my $seg_sites = 0;
   if( $type =~ /ARRAY/i ) {
       my %sites;
       foreach my $n ( @$individuals ) {
	   if( ! $n->isa('Bio::PopGen::IndividualI') ) {
	       $self->warn("Expected an arrayref of Bio::PopGen::IndividualI objects, this is a ".ref($n)."\n");
	       return 0;
	   }
	   foreach my $g ( $n->get_Genotypes ) {
	       my ($nm,@alleles) = ($g->marker_name, $g->get_Alleles);
	       foreach my $allele (@alleles ) {
		   $sites{$nm}->{$allele}++;
	       }
	   }
       }
       foreach my $site ( values %sites ) { # use values b/c we don't 
	                                    # really care what the name is
	   # find the sites which >1 allele
	   $seg_sites++ if( keys %$site > 1 );
       }
   } elsif( $type && $individuals->isa('Bio::PopGen::PopulationI') ) {
       foreach my $marker ( $individuals->haploid_population->get_Markers ) {  
	   my @alleles = $marker->get_Alleles;	    
	   $seg_sites++ if ( scalar @alleles > 1 );
       }
   } else { 
       $self->warn("segregating_sites_count expects either a PopulationI object or a list of IndividualI objects");
       return 0;
   } 
   return $seg_sites;

    my ($self,$samp_size, $freq1,$freq2) = @_;
    if( ! $freq2 ) { $freq2 = 1 - $freq1 }
    if( $freq1 > 1 || $freq2 > 1 ) { 
	$self->warn("heterozygosity expects frequencies to be less than 1");
    }
    my $sum = ($freq1**2) + (($freq2)**2);
    my $h = ( $samp_size*(1- $sum) ) / ($samp_size - 1) ;
    return $h;

   my ($self,$ingroup,$outgroup) = @_;
   my (%indata,%outdata,@marker_names);

   # basically we have to do some type checking
   # if that perl were typed...
   my ($itype,$otype) = (ref($ingroup),ref($outgroup));

   return $outgroup unless( $otype ); # we expect arrayrefs or objects, nums
                                      # are already the value we 
                                      # are searching for
   # pick apart the ingroup
   # get the data
   if( ref($ingroup) =~ /ARRAY/i ) {
       if( ! ref($ingroup->[0]) ||
	   ! $ingroup->[0]->isa('Bio::PopGen::IndividualI') ) {
	   $self->warn("Expected an arrayref of Bio::PopGen::IndividualI objects or a Population for ingroup in external_mutations");
	   return 0;
       }
       # we assume that all individuals have the same markers 
       # i.e. that they are aligned
       @marker_names = $ingroup->[0]->get_marker_names;
       for my $ind ( @$ingroup ) {
	   for my $m ( @marker_names ) {
	       for my $allele ( map { $_->get_Alleles }
				    $ind->get_Genotypes($m) ) {
		   $indata{$m}->{$allele}++;
	       }
	   }
       }	   
   } elsif( ref($ingroup) && $ingroup->isa('Bio::PopGen::PopulationI') ) {
       @marker_names = $ingroup->get_marker_names;
       for my $ind ( $ingroup->haploid_population->get_Individuals() ) {
	   for my $m ( @marker_names ) {
	       for my $allele ( map { $_->get_Alleles} 
				    $ind->get_Genotypes($m) ) {
		   $indata{$m}->{$allele}++;
	       }
	   }
       }
   } else { 
       $self->warn("Need an arrayref of Bio::PopGen::IndividualI objs or a Bio::PopGen::Population for ingroup in external_mutations");
       return 0;
   }
    
   if( $otype =~ /ARRAY/i ) {
       if( ! ref($outgroup->[0]) ||
	   ! $outgroup->[0]->isa('Bio::PopGen::IndividualI') ) {
	   $self->warn("Expected an arrayref of Bio::PopGen::IndividualI objects or a Population for outgroup in external_mutations");
	   return 0;
       }
       for my $ind ( @$outgroup ) {
	   for my $m ( @marker_names ) {
	       for my $allele ( map { $_->get_Alleles }
				$ind->get_Genotypes($m) ) {
		   $outdata{$m}->{$allele}++;
	       }
	   }
       }
   
   } elsif( $otype->isa('Bio::PopGen::PopulationI') ) {
       for my $ind ( $outgroup->haploid_population->get_Individuals() ) {
	   for my $m ( @marker_names ) {
	       for my $allele ( map { $_->get_Alleles} 
				    $ind->get_Genotypes($m) ) {
		   $outdata{$m}->{$allele}++;
	       }
	   }
       }
   } else {
       $self->warn("Need an arrayref of Bio::PopGen::IndividualI objs or a Bio::PopGen::Population for outgroup in external_mutations");
       return 0;
   }
   
   # derived mutations are defined as 
   # 
   # ingroup  (G A T)
   # outgroup (A)
   # derived mutations are G and T, A is the external mutation
   
   # ingroup  (A T)
   # outgroup (C)
   # derived mutations A,T no external/ancestral mutations
   
   # ingroup  (G A T)
   # outgroup (A T)
   # cannot determine
  
   my ($internal,$external);
   foreach my $marker ( @marker_names ) {
       my @outalleles = keys %{$outdata{$marker}};
       my @in_alleles = keys %{$indata{$marker}};
       next if( @outalleles > 1 || @in_alleles == 1);
       for my $allele ( @in_alleles ) {
	   if( ! exists $outdata{$marker}->{$allele} ) { 
	       if( $indata{$marker}->{$allele} == 1 ) { 
		   $external++;
	       } else { 
		   $internal++;
	       }
	   }
       }
   }
   return ($external, $internal);

    my ($self,$pop) = @_;
    if( ref($pop) =~ /ARRAY/i ) {
	if( ref($pop->[0]) && $pop->[0]->isa('Bio::PopGen::IndividualI') ) {
	    $pop = Bio::PopGen::Population->new(-individuals => @$pop);
	} else { 
	    $self->warn("composite_LD expects a Bio::PopGen::PopulationI or an arrayref of Bio::PopGen::IndividualI objects");
	    return ();
	}
    } elsif( ! ref($pop) || ! $pop->isa('Bio::PopGen::PopulationI') ) {
	$self->warn("composite_LD expects a Bio::PopGen::PopulationI or an arrayref of Bio::PopGen::IndividualI objects");
	return ();
    }

    my @marker_names = $pop->get_marker_names;
    my @inds = $pop->get_Individuals;
    my $num_inds = scalar @inds;
    my (%lookup);
    # calculate allele frequencies for each marker from the population
    # use the built-in get_Marker to get the allele freqs
    # we still need to calculate the genotype frequencies
    foreach my $marker_name ( @marker_names ) {	
	my(%allelef);

	foreach my $ind ( @inds ) {
	    my ($genotype) = $ind->get_Genotypes(-marker => $marker_name);
	    if( ! defined $genotype ) { 
		$self->warn("no genotype for marker $marker_name for individual ". $ind->unique_id. "\n");
		next;
	    }
	    my @alleles  = sort $genotype->get_Alleles;
	    next if( scalar @alleles != 2);
	    my $genostr  = join(',', @alleles);
            $allelef{$alleles[0]}++;
            $allelef{$alleles[1]}++;
	}

	# we should check for cases where there > 2 alleles or
	# only 1 allele and throw out those markers.
	my @alleles      = sort keys %allelef;
	my $allele_count = scalar @alleles;
	# test if site is polymorphic
	if( $allele_count != 2) { 
	    # only really warn if we're seeing multi-allele
	    $self->warn("Skipping $marker_name because it has $allele_count alleles (".join(',',@alleles)."),\n composite_LD will currently only work for biallelic markers") if $allele_count > 2;
	    next;		# skip this marker
	}

	# Need to do something here to detect alleles which aren't 
	# a single character
	if( length($alleles[0]) != 1 ||
	    length($alleles[1]) != 1 ) {
	    $self->warn("An individual has an allele which is not a single base, this is currently not supported in composite_LD - consider recoding the allele as a single character");
	    next;
	}

	# fix the call for allele 1 (A or B) and 
	# allele 2 (a or b) in terms of how we'll do the 
	# N square from Weir p.126
	$self->debug( "$alleles[0] is 1, $alleles[1] is 2 for $marker_name\n");
	$lookup{$marker_name}->{'1'} = $alleles[0];
	$lookup{$marker_name}->{'2'} = $alleles[1];
    }

    @marker_names = sort keys %lookup;
    my $site_count   = scalar @marker_names;
    # where the final data will be stored
    my %stats_for_sites;

    # standard way of generating pairwise combos
    # LD is done by comparing all the pairwise site (marker)
    # combinations and keeping track of the genotype and 
    # pairwise genotype (ie genotypes of the 2 sites) frequencies
    for( my $i = 0; $i < $site_count - 1; $i++ ) {
	my $site1 = $marker_names[$i];

	for( my $j = $i+1; $j < $site_count ; $j++) { 	 
	    my (%genotypes, %total_genotype_count,$total_pairwisegeno_count,
		%pairwise_genotypes);
	 
	    my $site2 = $marker_names[$j];
	    my (%allele_count,%allele_freqs) = (0,0);
	    foreach my $ind ( @inds ) {
		# build string of genotype at site 1
		my ($genotype1) = $ind->get_Genotypes(-marker => $site1);
		my @alleles1  = sort $genotype1->get_Alleles;

                # if an individual has only one available allele
		# (has a blank or N for one of the chromosomes)
		# we don't want to use it in our calculation

		next unless( scalar @alleles1 == 2);
		my $genostr1  = join(',', @alleles1);

		# build string of genotype at site 2
		my ($genotype2) = $ind->get_Genotypes(-marker => $site2);
		my @alleles2  = sort $genotype2->get_Alleles;
		my $genostr2  = join(',', @alleles2);
		
		next unless( scalar @alleles2 == 2);
		for (@alleles1) {
		    $allele_count{$site1}++;
		    $allele_freqs{$site1}->{$_}++;
		}
		$genotypes{$site1}->{$genostr1}++;
		$total_genotype_count{$site1}++;

		for (@alleles2) {
		    $allele_count{$site2}++;
		    $allele_freqs{$site2}->{$_}++;
		}
		$genotypes{$site2}->{$genostr2}++;
		$total_genotype_count{$site2}++;

		# We are using the $site1,$site2 to signify
		# a unique key
		$pairwise_genotypes{"$genostr1,$genostr2"}++;
		# some individuals 
		$total_pairwisegeno_count++;
	    }
	    for my $site ( %allele_freqs ) {
		for my $al ( keys %{ $allele_freqs{$site} } ) {
		    $allele_freqs{$site}->{$al} /= $allele_count{$site};
		}
	    }
	    my $n = $total_pairwisegeno_count;	# number of pairs of comparisons
	    # 'A' and 'B' are two loci or in our case site1 and site2  
	    my $allele1_site1 = $lookup{$site1}->{'1'};	# this is the BigA allele
	    my $allele1_site2 = $lookup{$site2}->{'1'};	# this is the BigB allele
	    my $allele2_site1 = $lookup{$site1}->{'2'};	# this is the LittleA allele
	    my $allele2_site2 = $lookup{$site2}->{'2'};	# this is the LittleB allele
	    # AABB
	    my $N1genostr = join(",",( $allele1_site1, $allele1_site1,
				       $allele1_site2, $allele1_site2));
	    $self->debug(" [$site1,$site2](AABB) N1genostr=$N1genostr\n");
	    # AABb
	    my $N2genostr = join(",",( $allele1_site1, $allele1_site1,
				       $allele1_site2, $allele2_site2));
	    $self->debug(" [$site1,$site2](AABb) N2genostr=$N2genostr\n");
	    # AaBB
	    my $N4genostr = join(",",( $allele1_site1, $allele2_site1,
				       $allele1_site2, $allele1_site2));
	    $self->debug(" [$site1,$site2](AaBB) N4genostr=$N4genostr\n");
	    # AaBb
	    my $N5genostr = join(",",( $allele1_site1, $allele2_site1,
				       $allele1_site2, $allele2_site2));
	    $self->debug(" [$site1,$site2](AaBb) N5genostr=$N5genostr\n");
	    # count of AABB in 
	    my $n1 = $pairwise_genotypes{$N1genostr} || 0;
	    # count of AABb in 
	    my $n2 = $pairwise_genotypes{$N2genostr} || 0;
	    # count of AaBB in 
	    my $n4 = $pairwise_genotypes{$N4genostr} || 0;
	    # count of AaBb in 
	    my $n5 = $pairwise_genotypes{$N5genostr} || 0;

	    my $homozA_site1 = join(",", ($allele1_site1,$allele1_site1));
	    my $homozB_site2 = join(",", ($allele1_site2,$allele1_site2));
	    my $p_AA = ($genotypes{$site1}->{$homozA_site1} || 0) / $n;
	    my $p_BB = ($genotypes{$site2}->{$homozB_site2} || 0) / $n;
	    my $p_A  = $allele_freqs{$site1}->{$allele1_site1} || 0;	# an individual allele freq
	    my $p_a  =  1 - $p_A;

	    my $p_B  = $allele_freqs{$site2}->{$allele1_site2} || 0;	# an individual allele freq
	    my $p_b  =  1 - $p_B;

	    # variance of allele frequencies
	    my $pi_A = $p_A * $p_a;
	    my $pi_B = $p_B * $p_b;

	    # hardy weinberg
	    my $D_A  = $p_AA - $p_A**2;
	    my $D_B  = $p_BB - $p_B**2;
	    my $n_AB = 2*$n1 + $n2 + $n4 + 0.5 * $n5;
	    $self->debug("n_AB=$n_AB -- n1=$n1, n2=$n2 n4=$n4 n5=$n5\n");

	    my $delta_AB = (1 / $n ) * ( $n_AB ) - ( 2 * $p_A * $p_B );
	    $self->debug("delta_AB=$delta_AB -- n=$n, n_AB=$n_AB p_A=$p_A, p_B=$p_B\n");
	    $self->debug(sprintf(" (%d * %.4f) / ( %.2f + %.2f) * ( %.2f + %.2f)\n ",
				 $n,$delta_AB**2, $pi_A, $D_A, $pi_B, $D_B));
	    
	    my $chisquared;
	    eval { $chisquared = ( $n * ($delta_AB**2) ) / 
		       ( ( $pi_A + $D_A) * ( $pi_B + $D_B) );
	       };
	    if( $@ ) {
		$self->debug("Skipping the site because the denom is 0.\nsite1=$site1, site2=$site2 : pi_A=$pi_A, pi_B=$pi_B D_A=$D_A, D_B=$D_B\n");
		next;
	    }
	    # this will be an upper triangular matrix
	    $stats_for_sites{$site1}->{$site2} = [$delta_AB,$chisquared];
	}
    }
    return %stats_for_sites;

    my ($self,@args) = @_;
    my ($ingroup, $outgroup,$polarized) = 
	$self->_rearrange([qw(INGROUP OUTGROUP POLARIZED)],@args);
    my $verbose = $self->verbose;
    my $outgroup_count;
    my $gapchar = '\-';
    if( ref($outgroup) =~ /ARRAY/i ) {
	$outgroup_count = scalar @$outgroup;
    } elsif( UNIVERSAL::isa($outgroup,'Bio::PopGen::PopulationI') ) {
	$outgroup_count = $outgroup->get_number_individuals;
    } else {
	$self->throw("Expected an ArrayRef of Individuals OR a Bio::PopGen::PopulationI");
    }
	
    if( $polarized ) {
	if( $outgroup_count < 2 ) {
	    $self->throw("Need 2 outgroups with polarized option\n");
	}
    } elsif( $outgroup_count > 1 ) {
	$self->warn(sprintf("%s outgroup sequences provided, but only first will be used",$outgroup_count ));
    } elsif( $outgroup_count == 0 ) {
	$self->throw("No outgroup sequence provided");
    }
    
    my $codon_path = Bio::MolEvol::CodonModel->codon_path;
    
    my (%marker_names,%unique,@inds);
    for my $p ( $ingroup, $outgroup)  {
	if( ref($p) =~ /ARRAY/i ) {
	    push @inds, @$p;
	} else {
	    push @inds, $p->get_Individuals;
	}
    }
    for my $i ( @inds ) {
	if( $unique{$i->unique_id}++ ) {
	    $self->warn("Individual ". $i->unique_id. " is seen more than once in the ingroup or outgroup set\n");
	}
	for my $n ( $i->get_marker_names ) {
	    $marker_names{$n}++;
	}
    }

    my @marker_names = keys %marker_names;
    if( $marker_names[0] =~ /^(Site|Codon)/ ) {
	# sort by site or codon number and do it in 
	# a schwartzian transformation baby!
	@marker_names = map { $_->[1] } 
	sort { $a->[0] <=> $b->[0] }
	map { [$_ =~ /^(?:Codon|Site)-(\d+)/, $_] } @marker_names;
    }


    my $num_inds = scalar @inds;
    my %vals = ( 'ingroup'  => $ingroup,
		 'outgroup' => $outgroup,		 
		 );

    # Make the Codon Table type a parameter!
    my $table = Bio::Tools::CodonTable->new(-id => $codon_table);
    my @vt = qw(outgroup ingroup);
    my %changes;
    my %status;
    my %two_by_two = ( 'fixed_N' => 0,
		       'fixed_S' => 0,
		       'poly_N'  => 0,
		       'poly_S'  => 0);

    for my $codon ( @marker_names ) {
	my (%codonvals);
	my %all_alleles;
	for my $t ( @vt ) {
	    my $outcount = 1;
	    for my $ind ( @{$vals{$t}} ) {
		my @alleles = $ind->get_Genotypes($codon)->get_Alleles;
		if( @alleles > 2 ) {
		    warn("Codon $codon saw ", scalar @alleles, " alleles for ind ", $ind->unique_id, "\n");
		    die;
		} else {
		    my ($allele) = shift @alleles;
		    $all_alleles{$ind->unique_id} = $allele;
		    my $AA = $table->translate($allele);
		    next if( $AA eq 'X' || $AA eq '*' || $allele =~ /N/i);

		    my $label = $t;
		    if( $t eq 'outgroup' ) {
			$label = $t.$outcount++;
		    }
		    $codonvals{$label}->{$allele}++;
		    $codonvals{all}->{$allele}++;
		}
	    }
	}
	my $total = sum ( values %{$codonvals{'ingroup'}} );
	next if( $total && $total < 2 ); # skip sites with < alleles
	# process all the seen alleles (codons) 
	# this is a vertical slide through the alignment
	if( keys %{$codonvals{all}} <= 1 ) {
	    # no changes or no VALID codons - monomorphic
	} else { 
	    # grab only the first outgroup codon (what to do with rest?)
	    my ($outcodon) = keys %{$codonvals{'outgroup1'}};
            if( ! $outcodon ) { 
		$status{"no outgroup codon $codon"}++;
		next;
	    }
	    my $out_AA = $table->translate($outcodon);
	    my ($outcodon2) = keys %{$codonvals{'outgroup2'}};
	    if( ($polarized && ($outcodon ne $outcodon2)) ||
		$out_AA eq 'X' || $out_AA eq '*' ) {
		# skip if outgroup codons are different 
		# (when polarized option is on)
		# or skip if the outcodon is STOP or 'NNN'
		if( $verbose > 0 ) {
		    $self->debug("skipping $out_AA and $outcodon $outcodon2\n");
		}
		$status{'outgroup codons different'}++;
		next;
	    }

	    # check if ingroup is actually different from outgroup -
	    # if there are the same number of alleles when considering
	    # ALL or just the ingroup, then there is nothing new seen
	    # in the outgroup so it must be a shared allele (codon)

	    # so we just count how many total alleles were seen
	    # if this is the same as the number of alleles seen for just 
	    # the ingroup then the outgroup presents no new information

	    my @ingroup_codons = keys %{$codonvals{'ingroup'}};
	    my $diff_from_out = ! exists $codonvals{'ingroup'}->{$outcodon};

	    if( $verbose > 0 ) {
		$self->debug("alleles are in: ", join(",", @ingroup_codons),
			     " out: ", join(",", keys %{$codonvals{outgroup1}}),
			     " diff_from_out=$diff_from_out\n");

		for my $ind ( sort keys %all_alleles ) {
		    $self->debug( "$ind\t$all_alleles{$ind}\n");
		}
	    }
	    # are all the ingroup alleles the same and diferent from outgroup?
	    # fixed differences between species
	    if( $diff_from_out ) {
		if( scalar @ingroup_codons == 1 ) { 
		    # fixed differences
		    if( $outcodon =~ /^$gapchar/ ) {
			$status{'outgroup codons with gaps'}++;
			next;
		    } elsif( $ingroup_codons[0] =~ /$gapchar/) {
			$status{'ingroup codons with gaps'}++;
			next;
		    }
		    my $path = $codon_path->{uc $ingroup_codons[0].$outcodon};
		    $two_by_two{fixed_N} += $path->[0];
		    $two_by_two{fixed_S} += $path->[1];
		    if( $verbose > 0 ) {
			$self->debug("ingroup is @ingroup_codons outcodon is $outcodon\n");
			$self->debug("path is ",join(",",@$path),"\n");
			$self->debug
			    (sprintf("%-15s fixeddiff - %s;%s(%s) %d,%d\tNfix=%d Sfix=%d Npoly=%d Spoly=%s\n",$codon,$ingroup_codons[0], $outcodon,$out_AA,
				     @$path, map { $two_by_two{$_} } 
				     qw(fixed_N fixed_S poly_N poly_S)));
		    }
		} else { 
		    # polymorphic and all are different from outgroup
		    # Here we find the minimum number of NS subst
		    my ($Ndiff,$Sdiff) = (3,0);	# most different path
		    for my $c ( @ingroup_codons ) {
			next if( $c =~ /$gapchar/ || $outcodon =~ /$gapchar/);
			my $path = $codon_path->{uc $c.$outcodon};
			my ($tNdiff,$tSdiff) = @$path;
			if( $path->[0] < $Ndiff ||
			    ($tNdiff == $Ndiff &&
			     $tSdiff  <= $Sdiff)) {
			    ($Ndiff,$Sdiff) = ($tNdiff,$tSdiff);
			}
		    }
		    $two_by_two{fixed_N} += $Ndiff;
		    $two_by_two{fixed_S} += $Sdiff;
	            if( @ingroup_codons > 2 ) { 
			$status{"more than 2 ingroup codons $codon"}++;
			warn("more than 2 ingroup codons (@ingroup_codons)\n");	
		    } else {
		    	my $path = $codon_path->{uc join('',@ingroup_codons)};

		    	$two_by_two{poly_N} += $path->[0];
		    	$two_by_two{poly_S} += $path->[1];
		    	if( $verbose > 0 ) {
			    $self->debug(sprintf("%-15s polysite_all - %s;%s(%s) %d,%d\tNfix=%d Sfix=%d Npoly=%d Spoly=%s\n",$codon,join(',',@ingroup_codons), $outcodon,$out_AA,@$path, map { $two_by_two{$_} } qw(fixed_N fixed_S poly_N poly_S)));
			}
		    } 
		} 
	    } else {
		my %unq = map { $_ => 1 } @ingroup_codons;
		delete $unq{$outcodon};
		my @unique_codons = keys %unq;

		# calc path for diff add to poly
		# Here we find the minimum number of subst bw
		# codons
		my ($Ndiff,$Sdiff) = (3,0); # most different path
		for my $c ( @unique_codons ) {
		    my $path = $codon_path->{uc $c.$outcodon };
		    if( ! defined $path ) {
			die " cannot get path for ", $c.$outcodon, "\n";
		    }
		    my ($tNdiff,$tSdiff) = @$path;
		    if( $path->[0] < $Ndiff ||
			($tNdiff == $Ndiff &&
			 $tSdiff  <= $Sdiff)) {
			($Ndiff,$Sdiff) = ($tNdiff,$tSdiff);
		    }
		}

		if( @unique_codons == 2 ) {
		    my $path = $codon_path->{uc join('',@unique_codons)};
		    if( ! defined $path ) {
			$self->throw("no path for @unique_codons\n");
		    }
		    $Ndiff += $path->[0];
		    $Sdiff += $path->[1];
		}
		$two_by_two{poly_N} += $Ndiff;
		$two_by_two{poly_S} += $Sdiff;
		if( $verbose > 0 ) {
		    $self->debug(sprintf("%-15s polysite - %s;%s(%s) %d,%d\tNfix=%d Sfix=%d Npoly=%d Spoly=%s\n",$codon,join(',',@ingroup_codons), $outcodon,$out_AA,
					 $Ndiff, $Sdiff, map { $two_by_two{$_} } 
					 qw(fixed_N fixed_S poly_N poly_S)));
		}
	    }
	}	    
    }
    return ( $two_by_two{'poly_N'},
	     $two_by_two{'fixed_N'},
	     $two_by_two{'poly_S'},
	     $two_by_two{'fixed_S'},
	     {%status});
    
}

*MK =\& mcdonald_kreitman;

    my ($self,$Npoly,$Nfix,$Spoly,$Sfix) = @_;
    if( $has_twotailed ) {
	return &Text::NSP::Measures::2D::Fisher2::twotailed::calculateStatistic 
	    (n11=>$Npoly,
	     n1p=>$Npoly+$Spoly,
	     np1=>$Npoly+$Nfix,
	     npp=>$Npoly+$Nfix+$Spoly+$Sfix);
    } else {
	$self->warn("cannot call mcdonald_kreitman_counts because no Fisher's exact is available - install Text::NSP::Measures::2D::Fisher2::twotailed");
	return 0;
    }
}


1;

  bioperl-l@bioperl.org                  - General discussion
  http://bioperl.org/wiki/Mailing_lists  - About the mailing lists

  https://redmine.open-bio.org/projects/bioperl/

 Title   : new
 Usage   : my $obj = Bio::PopGen::Statistics->new();
 Function: Builds a new Bio::PopGen::Statistics object 
 Returns : an instance of Bio::PopGen::Statistics
 Args    : none