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

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

  my $tree = $sim->next_tree;

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

  my $stats = new Bio::PopGen::Statistics();
  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 = new Bio::PopGen::IO(-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 = new Bio::AlignIO(-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);

     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;
    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 tajima_D");
	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 $D = ( $pi - ( $seg_sites / $a1 ) ) / 
	sqrt ( ($e1 * $seg_sites) + (( $e2 * $seg_sites) * ( $seg_sites - 1)));

    return $D;

    my ($self,$individuals,$numsites) = @_;
    my (%data,@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
	my %marker_total;
	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};
	}
    } else { 
	$self->throw("expected an array reference of a list of Bio::PopGen::IndividualI to pi");
    }
    # doing all pairwise combinations

    # For now we assume that all individuals have the same markers
    my ($diffcount,$totalcompare) = (0,0);
    my $pi = 0;
    foreach my $markerdat ( values %data ) {
	my $totalalleles; # this will only be different among markers
	                  # when there is missing data
	my @alleles = keys %$markerdat;
	foreach my $al ( @alleles ) { $totalalleles += $markerdat->{$al} }
	for( my $i =0; $i < scalar @alleles -1; $i++ ) {
	    my ($a1,$a2) = ( $alleles[$i], $alleles[$i+1]);
	    $pi += $self->heterozygosity($n, 
					 $markerdat->{$a1} / $totalalleles,
					 $markerdat->{$a2} / $totalalleles);
	}
    }
    $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;
    }
    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;

   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}++;
	       }
	   }
       }
   } elsif( $otype->isa('Bio::PopGen::PopulationI') ) { 
       $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 = new Bio::PopGen::Population(-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];
	my (%genotypes, %total_genotype_count,
    	%total_pairwisegeno_count,%pairwise_genotypes);
	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{"$site1,$site2"}->{"$genostr1,$genostr2"}++;
		# some individuals 
		$total_pairwisegeno_count{"$site1,$site2"}++;
	    }
	    for my $site ( %allele_freqs ) {
		for my $al ( keys %{ $allele_freqs{$site} } ) {
		    $allele_freqs{$site}->{$al} /= $allele_count{$site};
		}
	    }
	    my $n = $total_pairwisegeno_count{"$site1,$site2"};	# number of inds
	    # '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{"$site1,$site2"}->{$N1genostr} || 0;
	    # count of AABb in 
	    my $n2 = $pairwise_genotypes{"$site1,$site2"}->{$N2genostr} || 0;
	    # count of AaBB in 
	    my $n4 = $pairwise_genotypes{"$site1,$site2"}->{$N4genostr} || 0;
	    # count of AaBb in 
	    my $n5 = $pairwise_genotypes{"$site1,$site2"}->{$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;

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

  http://bugzilla.open-bio.org/

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