The stat handler for neutral markers. More...

#include <ttneutralgenes.h>

Inheritance diagram for TTNeutralGenesSH:

Collaboration diagram for TTNeutralGenesSH:

Public Member Functions
	TTNeutralGenesSH (TProtoNeutralGenes *TP)
virtual	~TTNeutralGenesSH ()
virtual void	init ()
	Inits state.
virtual bool	setStatRecorders (std::string &token)
void	setFstatRecorders (age_t AGE)
void	setFstatWCRecorders (age_t AGE)
void	setCoaMatrixRecorders (age_t AGE, unsigned char dim)
void	setFstMatrixRecorders (age_t AGE, unsigned char dim)
void	setNeiGeneticDistanceRecorders (age_t AGE, bool pairwise)
F-stats:
void	setAlleleTables (age_t AGE)
void	allocateTables (unsigned int loci, unsigned int all)
void	setFstMatrix (age_t AGE, unsigned char dim)
	Computes the weighted within and between patch Fst's as well as the overall Fst (Theta).
void	setAdultsFstMatrix ()
void	setAdultsFstWithin ()
void	setAdultsFstBetween ()
void	setOffsprgFstMatrix ()
void	setOffsprgFstWithin ()
void	setOffsprgFstBetween ()
double	getWeightedFst ()
	Returns the weighted Fst using Weir & Hill (2002) method.
double	getFst_ij (unsigned int i)
	Accessor to the Fst matrix as set by setFstMatrix().
void	setFst_li (unsigned int N, unsigned int L, double **array)
	Computes the per-locus per-patch Fst values using Weir&Hill 2002 approach.
DataTable< double > *	getAlleleFreqTable ()
	Computes raw Fst following the original definition (=var(p)/p_bar(1 - p_bar)).
TMatrix *	getGlobalFreqs ()
	Accessor to the table of allele frequencies in the whole population.
void	setFstat (age_t AGE)
	Computes the F-statistics following Nei & Chesser (1983).
void	setOffsprgFstat ()
void	setAdultsFstat ()
double	setHo (age_idx age_pos)
double	setHs (age_idx age_pos)
double	setHt (age_idx age_pos)
double	getHsnei ()
double	getHtnei ()
double	getHo ()
double	getHs ()
double	getHt ()
double	getFst ()
double	getFis ()
double	getFit ()
void	setFstatWeirCockerham (age_t AGE)
	Computes the Weir & Cockerham (1984) Fstat values (Theta, F, and f).
void	setOffspringFstatWeirCockerham ()
void	setAdultsFstatWeirCockerham ()
double	getFstWC ()
double	getFisWC ()
double	getFitWC ()
void	setLociDivCounter (age_t AGE)
	Sets the allelic diversity counters.
double	getNbAllLocal ()
double	getNbAllGlobal ()
double	getFixLocLocal ()
double	getFixLocGlobal ()
Coancestries
double	Coancestry (void ind1, void ind2, unsigned int nb_locus)
	Gives the coancestry (probability of identity by state) of two gene sequences.
void	setCoaMatrix (age_idx age_pos, unsigned char dim)
	Computes the within and between patches coancestry coefficients.
void	setAdultsCoaMatrix ()
void	setOffsprgCoaMatrix ()
void	setAdultsCoaWithin ()
void	setOffsprgCoaWithin ()
void	setAdultsCoaBetween ()
void	setOffsprgCoaBetween ()
void	setAdults_Theta ()
double	getCoa (unsigned int i)
	Gets the given coancestry coefficient from the coancestry matrix.
double	getMeanTheta ()
double	getMeanAlpha ()
double	getTheta_FF ()
	Gives the mean within females coancestry coefficient.
double	getTheta_MM ()
	Gives the mean within males coancestry coefficient.
double	getTheta_FM ()
	Gives the mean between males and females coancestry coefficient.
void	setSibStats ()
void	setSibCoa (Individual I1, Individual I2)
double	getSibProportions (unsigned int i)
double	getSibCoaMeans (unsigned int i)
Nei's genetic distance:
void	setAdltNeiGeneticDistance ()
void	setOffsprgNeiGeneticDistance ()
void	setNeiGeneticDistance (age_t AGE)
double	getNeiGeneticDistance (unsigned int i)
double	getMeanNeiGeneticDistance ()
Private Attributes
DataTable< unsigned int >	_alleleCountTable
DataTable< double >	_alleleFreqTable
TMatrix	_globalAlleleFreq
unsigned int	_table_set_gen
unsigned int	_table_set_age
unsigned int	_table_set_repl
double	Theta_FF
double	Theta_MM
double	Theta_FM
double	_mean_theta
double	_mean_alpha
TMatrix *	_coa_matrix
double	_sib_prop [4]
	Kinship classes proportions.
double	_sib_coa [4]
double	_ho
double	_hs
double	_ht
double	_hsnei
double	_htnei
double	_nb_all_local
double	_nb_all_global
double	_fst
double	_fis
double	_fit
double	_fix_loc_local
double	_fix_loc_global
double	_fst_WH
double	_fst_WC
double	_fis_WC
double	_fit_WC
double	_fst_W1
double	_fst_W2
TMatrix *	_fst_matrix
TMatrix *	_D
double	_meanD

Detailed Description

The stat handler for neutral markers.

Definition at line 229 of file ttneutralgenes.h.

Constructor & Destructor Documentation

TTNeutralGenesSH::TTNeutralGenesSH ( TProtoNeutralGenes * TP ) [inline]

Definition at line 257 of file ttneutralgenes.h.

    : TraitStatHandler<TProtoNeutralGenes, TTNeutralGenesSH>(TP), _table_set_gen(0), _table_set_age(0),
    _table_set_repl(0), _coa_matrix(0), _fst_matrix(0), _D(0) { }

virtual TTNeutralGenesSH::~TTNeutralGenesSH ( ) [inline, virtual]

Definition at line 261 of file ttneutralgenes.h.

References _coa_matrix, _D, and _fst_matrix.

  {
    if(_coa_matrix != NULL) delete _coa_matrix;
    if(_fst_matrix != NULL) delete _fst_matrix;
    if(_D != NULL) delete _D;
  }

Member Function Documentation

void TTNeutralGenesSH::allocateTables	(	unsigned int	loci,
		unsigned int	all
	)

Definition at line 1077 of file ttneutralgenes.cc.

References _alleleCountTable, _alleleFreqTable, _globalAlleleFreq, StatHandlerBase::_pop, DataTable< T >::allocate(), Metapop::getPatchNbr(), and TMatrix::reset().

Referenced by init(), TTNeutralGenesFH::write_freq(), TTNeutralGenesFH::write_Fst_i(), and TTNeutralGenesFH::write_varcompWC().

{
  unsigned int nb_patch = _pop->getPatchNbr();
  unsigned int **sizes;

  sizes = new unsigned int * [nb_patch];
  
  for(unsigned int i = 0; i < nb_patch; ++i) {
    sizes[i] = new unsigned int [loci];
    for(unsigned int j = 0; j < loci; ++j)
      sizes[i][j] = all;
  }

  _alleleCountTable.allocate(nb_patch, loci, sizes);
  
  _alleleFreqTable.allocate(nb_patch, loci, sizes);
  
  _globalAlleleFreq.reset(loci, all);

  for(unsigned int i = 0; i < nb_patch; ++i)
    delete [] sizes[i];
  delete [] sizes;
}

double TTNeutralGenesSH::Coancestry	(	void **	ind1,
		void **	ind2,
		unsigned int	nb_locus
	)

Gives the coancestry (probability of identity by state) of two gene sequences.

The probability returned is the average probability of having two identical alleles at a locus between the two sequences.

Parameters:

ind1	first sequence, treated as of type (unsigned char**)
ind2	second sequence, treated as of type (unsigned char**)
nb_locus	number of loci present in each sequence

Definition at line 43 of file stats_coa.cc.

{
  unsigned int p = 0;
  unsigned char **seq1, **seq2;
  
  seq1 = (unsigned char**)ind1;
  seq2 = (unsigned char**)ind2;
  
  for (unsigned int k = 0; k < nb_locus; ++k)
        p += !(seq1[0][k]^seq2[0][k]) + !(seq1[0][k]^seq2[1][k]) + !(seq1[1][k]^seq2[0][k]) + !(seq1[1][k]^seq2[1][k]);

  return (double)p/(4.0*nb_locus);
}

DataTable<double>* TTNeutralGenesSH::getAlleleFreqTable ( ) [inline]

Computes raw Fst following the original definition (=var(p)/p_bar(1 - p_bar)).

Accessor to the table of allele frequencies, per patch.

Definition at line 314 of file ttneutralgenes.h.

References _alleleFreqTable.

Referenced by TTNeutralGenesFH::write_freq(), and TTNeutralGenesFH::write_varcompWC().

{return &_alleleFreqTable;}

double TTNeutralGenesSH::getCoa ( unsigned int i ) [inline]

Gets the given coancestry coefficient from the coancestry matrix.

Parameters:

i	combination of the row and column indexes (see setCoaMatrixRecorders()).

Note:: the upper half and the diagonal of the matrix are filled, other positions are set to 0.

Definition at line 377 of file ttneutralgenes.h.

References _coa_matrix, TMatrix::get(), and TMatrix::getNbCols().

Referenced by setCoaMatrixRecorders().

  {
    unsigned int scale = (unsigned int)pow( 10.0, (int)log10((float)_coa_matrix->getNbCols()) + 1 );
    return _coa_matrix->get(i/scale, i%scale);
  }

double TTNeutralGenesSH::getFis ( ) [inline]

Definition at line 331 of file ttneutralgenes.h.

References _fis.

Referenced by setFstatRecorders().

{return _fis;}

double TTNeutralGenesSH::getFisWC ( ) [inline]

Definition at line 338 of file ttneutralgenes.h.

References _fis_WC.

Referenced by setFstatWCRecorders().

{return _fis_WC;}

double TTNeutralGenesSH::getFit ( ) [inline]

Definition at line 332 of file ttneutralgenes.h.

References _fit.

Referenced by setFstatRecorders().

{return _fit;}

double TTNeutralGenesSH::getFitWC ( ) [inline]

Definition at line 339 of file ttneutralgenes.h.

References _fit_WC.

Referenced by setFstatWCRecorders().

{return _fit_WC;}

double TTNeutralGenesSH::getFixLocGlobal ( ) [inline]

Definition at line 345 of file ttneutralgenes.h.

References _fix_loc_global.

Referenced by setFstatRecorders().

{return _fix_loc_global;}

double TTNeutralGenesSH::getFixLocLocal ( ) [inline]

Definition at line 344 of file ttneutralgenes.h.

References _fix_loc_local.

Referenced by setFstatRecorders().

{return _fix_loc_local;}

double TTNeutralGenesSH::getFst ( ) [inline]

Definition at line 330 of file ttneutralgenes.h.

References _fst.

Referenced by setFstatRecorders().

{return _fst;}

double TTNeutralGenesSH::getFst_ij ( unsigned int i ) [inline]

Accessor to the Fst matrix as set by setFstMatrix().

Definition at line 304 of file ttneutralgenes.h.

References _fst_matrix, TMatrix::get(), and TMatrix::getNbCols().

Referenced by setFstMatrixRecorders().

  {
    unsigned int scale = (unsigned int)pow( 10.0, (int)log10((float)_fst_matrix->getNbCols()) + 1 );
    return _fst_matrix->get(i/scale, i%scale);
  }

double TTNeutralGenesSH::getFstWC ( ) [inline]

Definition at line 337 of file ttneutralgenes.h.

References _fst_WC.

Referenced by setFstatWCRecorders().

{return _fst_WC;}

TMatrix* TTNeutralGenesSH::getGlobalFreqs ( ) [inline]

Accessor to the table of allele frequencies in the whole population.

Definition at line 316 of file ttneutralgenes.h.

References _globalAlleleFreq.

Referenced by TTNeutralGenesFH::write_varcompWC().

{return &_globalAlleleFreq;}

double TTNeutralGenesSH::getHo ( ) [inline]

Definition at line 327 of file ttneutralgenes.h.

References _ho.

Referenced by setFstatRecorders().

{return _ho;}

double TTNeutralGenesSH::getHs ( ) [inline]

Definition at line 328 of file ttneutralgenes.h.

References _hs.

{return _hs;}

double TTNeutralGenesSH::getHsnei ( ) [inline]

Definition at line 325 of file ttneutralgenes.h.

References _hsnei.

Referenced by setFstatRecorders().

{return _hsnei;}

double TTNeutralGenesSH::getHt ( ) [inline]

Definition at line 329 of file ttneutralgenes.h.

References _ht.

{return _ht;}

double TTNeutralGenesSH::getHtnei ( ) [inline]

Definition at line 326 of file ttneutralgenes.h.

References _htnei.

Referenced by setFstatRecorders().

{return _htnei;}

double TTNeutralGenesSH::getMeanAlpha ( ) [inline]

Definition at line 383 of file ttneutralgenes.h.

References _mean_alpha.

Referenced by setCoaMatrixRecorders(), and setStatRecorders().

{return _mean_alpha;}

double TTNeutralGenesSH::getMeanNeiGeneticDistance ( ) [inline]

Definition at line 406 of file ttneutralgenes.h.

References _meanD.

Referenced by setNeiGeneticDistanceRecorders().

{return _meanD;}

double TTNeutralGenesSH::getMeanTheta ( ) [inline]

Definition at line 382 of file ttneutralgenes.h.

References _mean_theta.

Referenced by setCoaMatrixRecorders(), and setStatRecorders().

{return _mean_theta;}

double TTNeutralGenesSH::getNbAllGlobal ( ) [inline]

Definition at line 343 of file ttneutralgenes.h.

References _nb_all_global.

Referenced by setFstatRecorders().

{return _nb_all_global;}

double TTNeutralGenesSH::getNbAllLocal ( ) [inline]

Definition at line 342 of file ttneutralgenes.h.

References _nb_all_local.

Referenced by setFstatRecorders().

{return _nb_all_local;}

double TTNeutralGenesSH::getNeiGeneticDistance ( unsigned int i ) [inline]

Definition at line 401 of file ttneutralgenes.h.

References _D, TMatrix::get(), and TMatrix::getNbCols().

Referenced by setNeiGeneticDistanceRecorders().

  { 
    unsigned int scale = (unsigned int)pow( 10.0, (int)log10((float)_D->getNbCols()) + 1 );
    return _D->get(i/scale,i%scale);
  }

double TTNeutralGenesSH::getSibCoaMeans ( unsigned int i ) [inline]

Definition at line 393 of file ttneutralgenes.h.

References _sib_coa.

Referenced by setStatRecorders().

{return _sib_coa[i];}

double TTNeutralGenesSH::getSibProportions ( unsigned int i ) [inline]

Definition at line 392 of file ttneutralgenes.h.

References _sib_prop.

Referenced by setStatRecorders().

{return _sib_prop[i];}

double TTNeutralGenesSH::getTheta_FF ( ) [inline]

Gives the mean within females coancestry coefficient.

Definition at line 385 of file ttneutralgenes.h.

References Theta_FF.

Referenced by setStatRecorders().

{return Theta_FF;}

double TTNeutralGenesSH::getTheta_FM ( ) [inline]

Gives the mean between males and females coancestry coefficient.

Definition at line 389 of file ttneutralgenes.h.

References Theta_FM.

Referenced by setStatRecorders().

{return Theta_FM;}

double TTNeutralGenesSH::getTheta_MM ( ) [inline]

Gives the mean within males coancestry coefficient.

Definition at line 387 of file ttneutralgenes.h.

References Theta_MM.

Referenced by setStatRecorders().

{return Theta_MM;}

double TTNeutralGenesSH::getWeightedFst ( ) [inline]

Returns the weighted Fst using Weir & Hill (2002) method.

This Fst is set by a previous call to setFstMatrix().

Definition at line 302 of file ttneutralgenes.h.

References _fst_WH.

Referenced by setFstMatrixRecorders().

{return _fst_WH;}

void TTNeutralGenesSH::init ( ) [virtual]

Inits state.

Reimplemented from StatHandlerBase.

Definition at line 1068 of file ttneutralgenes.cc.

References TraitStatHandler< TProtoNeutralGenes, TTNeutralGenesSH >::_SHLinkedTrait, allocateTables(), TProtoNeutralGenes::get_allele_num(), and TProtoNeutralGenes::get_locus_num().

{
  StatHandler<TTNeutralGenesSH>::init();

  allocateTables(_SHLinkedTrait->get_locus_num(),_SHLinkedTrait->get_allele_num());
}

void TTNeutralGenesSH::setAdltNeiGeneticDistance ( ) [inline]

Definition at line 398 of file ttneutralgenes.h.

References ADULTS, and setNeiGeneticDistance().

Referenced by setNeiGeneticDistanceRecorders().

{setNeiGeneticDistance(ADULTS);}

void TTNeutralGenesSH::setAdults_Theta ( )

Definition at line 203 of file stats_coa.cc.

Referenced by setStatRecorders().

{
  unsigned int Fsize,Msize,i,j,k, tot_size, FFsize, MMsize, FMsize;
  unsigned int patchNbr = this->_pop->getPatchNbr();
  unsigned int nb_locus = this->_SHLinkedTrait->get_locus_num();
  Patch *P1;
  double mean = 0, grand_mean = 0;
  
  Theta_FF = 0;
  Theta_MM = 0;
  Theta_FM = 0;
  
  _mean_theta = 0;
  
  for(i = 0; i < patchNbr; ++i) {
    
    P1 = _pop->getPatch(i);
        Fsize = P1->size(FEM, ADLTx);
        Msize = P1->size(MAL, ADLTx);
    tot_size = Fsize + Msize;
    
    FFsize = Fsize*(Fsize-1)/2;
    MMsize = Msize*(Msize-1)/2;
    FMsize = Fsize*Msize;
    
    grand_mean = 0;
    
        if(Fsize != 0 || Msize != 0) {
      
      if(Fsize != 0) {
        mean = 0;
        for(j = 0; j < Fsize-1;++j)
          for(k = j+1; k < Fsize;++k)
            mean += Coancestry(P1->get(FEM, ADLTx, j)->getTrait(_SHLinkedTraitIndex)->get_sequence(),
                                   P1->get(FEM, ADLTx, k)->getTrait(_SHLinkedTraitIndex)->get_sequence(),
                                   nb_locus);
        Theta_FF += mean/FFsize;
      }
      grand_mean += mean;
      if(Msize != 0) {
        mean = 0;
        for(j = 0; j < Msize-1;++j)
          for(k = j+1; k < Msize;++k)
            mean += Coancestry(P1->get(MAL, ADLTx, j)->getTrait(_SHLinkedTraitIndex)->get_sequence(),
                                   P1->get(MAL, ADLTx, k)->getTrait(_SHLinkedTraitIndex)->get_sequence(),
                                   nb_locus);
        Theta_MM += mean/MMsize;
      }
      grand_mean += mean;
      if(Fsize != 0 && Msize != 0) {
        mean = 0;
        for(j = 0; j < Fsize;++j)
          for(k = 0; k < Msize;++k)
            mean += Coancestry(P1->get(FEM, ADLTx, j)->getTrait(_SHLinkedTraitIndex)->get_sequence(),
                                   P1->get(MAL, ADLTx, k)->getTrait(_SHLinkedTraitIndex)->get_sequence(),
                                   nb_locus);
        Theta_FM += mean/FMsize;
      }
      grand_mean += mean;
      _mean_theta += grand_mean / (FFsize + MMsize + FMsize);
    }
  }
  
  _mean_theta /= patchNbr;
  Theta_FF /= patchNbr;
  Theta_MM /= patchNbr;
  Theta_FM /= patchNbr;
}

void TTNeutralGenesSH::setAdultsCoaBetween ( ) [inline]

Definition at line 369 of file ttneutralgenes.h.

References ADLTx, and setCoaMatrix().

Referenced by setCoaMatrixRecorders(), and setStatRecorders().

{setCoaMatrix(ADLTx, 2);}

void TTNeutralGenesSH::setAdultsCoaMatrix ( ) [inline]

Definition at line 365 of file ttneutralgenes.h.

References ADLTx, and setCoaMatrix().

Referenced by setCoaMatrixRecorders(), and setStatRecorders().

{setCoaMatrix(ADLTx, 3);}

void TTNeutralGenesSH::setAdultsCoaWithin ( ) [inline]

Definition at line 367 of file ttneutralgenes.h.

References ADLTx, and setCoaMatrix().

Referenced by setCoaMatrixRecorders(), and setStatRecorders().

{setCoaMatrix(ADLTx, 1);}

void TTNeutralGenesSH::setAdultsFstat ( ) [inline]

Definition at line 321 of file ttneutralgenes.h.

References ADULTS, and setFstat().

Referenced by setFstatRecorders().

{setFstat(ADULTS);}

void TTNeutralGenesSH::setAdultsFstatWeirCockerham ( ) [inline]

Definition at line 336 of file ttneutralgenes.h.

References ADULTS, and setFstatWeirCockerham().

Referenced by setFstatWCRecorders().

{setFstatWeirCockerham(ADULTS);}

void TTNeutralGenesSH::setAdultsFstBetween ( ) [inline]

Definition at line 296 of file ttneutralgenes.h.

References ADULTS, and setFstMatrix().

Referenced by setFstMatrixRecorders().

{setFstMatrix(ADULTS, 2);}

void TTNeutralGenesSH::setAdultsFstMatrix ( ) [inline]

Definition at line 294 of file ttneutralgenes.h.

References ADULTS, and setFstMatrix().

Referenced by setFstMatrixRecorders().

{setFstMatrix(ADULTS, 3);}

void TTNeutralGenesSH::setAdultsFstWithin ( ) [inline]

Definition at line 295 of file ttneutralgenes.h.

References ADULTS, and setFstMatrix().

Referenced by setFstMatrixRecorders().

{setFstMatrix(ADULTS, 1);}

void TTNeutralGenesSH::setAlleleTables ( age_t AGE )

Definition at line 41 of file stats_fstat.cc.

Referenced by TTNeutralGenesFH::write_freq(), and TTNeutralGenesFH::write_varcompWC().

{  
  
  if(_table_set_age == AGE 
     && _table_set_gen == _pop->getCurrentGeneration()
     && _table_set_repl == _pop->getCurrentReplicate())
    return;
  
  unsigned char** genes;
  unsigned int nb_locus = _SHLinkedTrait->get_locus_num(), nb_allele = _SHLinkedTrait->get_allele_num();
  unsigned int patch_size, patchNbr = _pop->getPatchNbr();
  age_idx age_pos = (AGE == ADULTS ? ADLTx : OFFSx);
  Patch* crnt_patch;
  
  _alleleCountTable.init(0);

  //counting the copies of each allele present in the pop:
  for(unsigned int k = 0; k < patchNbr; ++k) {

    crnt_patch = _pop->getPatch(k);

    for(unsigned int i = 0, size = crnt_patch->size(FEM, age_pos); i < size; ++i) {

      genes = (unsigned char**)crnt_patch->get(FEM, age_pos, i)->getTrait(_SHLinkedTraitIndex)->get_sequence();
      
      for (unsigned int j = 0; j < nb_locus; ++j)
        _alleleCountTable.increment(k, j, genes[0][j]);
      
      for (unsigned int j = 0; j < nb_locus; ++j)
        _alleleCountTable.increment(k, j, genes[1][j]);
    } 
    
    for(unsigned int i = 0, size = crnt_patch->size(MAL, age_pos); i < size; ++i) {
      
      genes = (unsigned char**)crnt_patch->get(MAL, age_pos, i)->getTrait(_SHLinkedTraitIndex)->get_sequence();
      
      for (unsigned int j = 0; j < nb_locus; ++j)
        _alleleCountTable.increment(k, j, genes[0][j]);
      
      for (unsigned int j = 0; j < nb_locus; ++j)
        _alleleCountTable.increment(k, j, genes[1][j]);
    }
  }

  //allelic frequencies:
  for (unsigned int i = 0; i < patchNbr; ++i) {
    
    patch_size = _pop->getPatch(i)->size(age_pos) * 2;
    
    for (unsigned int l = 0; l < nb_locus; ++l)
      
      for (unsigned int u = 0; u < nb_allele; ++u)
        
        _alleleFreqTable.set(i, l, u, (patch_size != 0 ? (double)_alleleCountTable.get(i, l, u) / patch_size : 0) );
  }
  
  //population-wide allelic frequencies:  
  unsigned int tot_size = _pop->size(AGE) * 2;

  _globalAlleleFreq.assign(0);
  
  for(unsigned int i = 0; i < patchNbr; i++)
    for (unsigned int l = 0; l < nb_locus; ++l)
      for (unsigned int u = 0; u < nb_allele; ++u)
        _globalAlleleFreq.plus(l, u, _alleleCountTable.get(i, l, u));
  
  for (unsigned int l = 0; l < nb_locus; ++l)
    for (unsigned int u = 0; u < nb_allele; ++u)
      _globalAlleleFreq.divide(l, u, tot_size);
  
  _table_set_age = AGE;
  
  _table_set_gen = _pop->getCurrentGeneration();
  
  _table_set_repl = _pop->getCurrentReplicate();
}

void TTNeutralGenesSH::setCoaMatrix	(	age_idx	age_pos,
		unsigned char	dim
	)

Computes the within and between patches coancestry coefficients.

Parameters:

age_pos	the age class index
dim	the dimension of the matrix to fill: 1 = the diagonal (i.e. the wihtin patch coancestries or theta's) 2 = the upper half (i.e. the between patch coancestries or alpha's) 3 = both

Definition at line 59 of file stats_coa.cc.

Referenced by setAdultsCoaBetween(), setAdultsCoaMatrix(), setAdultsCoaWithin(), setOffsprgCoaBetween(), setOffsprgCoaMatrix(), and setOffsprgCoaWithin().

{
  unsigned int Fsize,Msize,tot_size,size_i,size_l, wt;
  unsigned int nb_coeff = 0;
  unsigned int i, j, k, l;
  unsigned int patchNbr = this->_pop->getPatchNbr();
  unsigned int nb_locus = this->_SHLinkedTrait->get_locus_num();
  Patch *P1, *P2;
  
  if(_coa_matrix == NULL)
    
    _coa_matrix = new TMatrix(patchNbr, patchNbr);
  
  else if( _coa_matrix->length() != patchNbr * patchNbr)
  
    _coa_matrix->reset(patchNbr, patchNbr);
  
  _coa_matrix->assign(0);
  
  if(dim & 1) {
    
    _mean_theta = 0;
    
    wt = 0;
    
    //first fill the diagonale: within deme coancestry (betha)
    for(i = 0; i < patchNbr; ++i) {
      
      P1 = _pop->getPatch(i);
      Fsize = P1->size(FEM, age_pos);
      Msize = P1->size(MAL, age_pos);
      
      tot_size = Fsize + Msize;
      
      if(tot_size != 0) {
        
        //fem-fem coa
        for(j = 0; j < Fsize-1; ++j)
          for(k = j+1; k < Fsize; ++k)
            _coa_matrix->plus(i, i, Coancestry(P1->get(FEM, age_pos, j)->getTrait(_SHLinkedTraitIndex)->get_sequence(),
                                               P1->get(FEM, age_pos, k)->getTrait(_SHLinkedTraitIndex)->get_sequence(),
                                               nb_locus));
        //mal-mal coa
        for(j = 0; j < Msize-1; ++j)
          for(k = j+1; k < Msize; ++k)
            _coa_matrix->plus(i, i, Coancestry(P1->get(MAL, age_pos, j)->getTrait(_SHLinkedTraitIndex)->get_sequence(),
                                               P1->get(MAL, age_pos, k)->getTrait(_SHLinkedTraitIndex)->get_sequence(),
                                               nb_locus));
        //fem-mal coa
        for(j = 0; j < Fsize; ++j)
          for(k = 0; k < Msize; ++k)
            _coa_matrix->plus(i, i, Coancestry(P1->get(FEM, age_pos, j)->getTrait(_SHLinkedTraitIndex)->get_sequence(),
                                               P1->get(MAL, age_pos, k)->getTrait(_SHLinkedTraitIndex)->get_sequence(),
                                               nb_locus));
        
        _coa_matrix->divide(i, i, tot_size*(tot_size -1)/2.0); 
      }//end if
      
      _mean_theta += tot_size * _coa_matrix->get(i, i);
    
      wt += tot_size;
    
    }//end for patchNbr
    
    //weighted average:
    _mean_theta /= wt;
    
  } //end if diag
  
  if(dim & 2) {
    //fill the first upper half of the matrix: between deme coancestry (alpha)
   
    _mean_alpha = 0;
    
    wt = 0;

    for(i = 0; i < patchNbr-1; ++i) {
      
      P1 = _pop->getPatch(i);
      
      for(l = i+1; l < patchNbr; ++l) {
        
        P2 = _pop->getPatch(l);
                
        tot_size = P1->size(age_pos) * P2->size(age_pos);
        
        if(tot_size != 0) {
          //females i vs. females l
          size_i = P1->size(FEM, age_pos);
          size_l = P2->size(FEM, age_pos);
          nb_coeff = size_i * size_l;
          
          for(j = 0; j < size_i; ++j)
            for(k = 0; k < size_l; ++k)
              _coa_matrix->plus(i, l, Coancestry(P1->get(FEM, age_pos, j)->getTrait(_SHLinkedTraitIndex)->get_sequence(),
                                                 P2->get(FEM, age_pos, k)->getTrait(_SHLinkedTraitIndex)->get_sequence(),
                                                 nb_locus));
          //females i vs. males l
          size_l = P2->size(MAL, age_pos);
          nb_coeff += size_i * size_l;
          
          for(j = 0; j < size_i; ++j)
            for(k = 0; k < size_l; ++k)
              _coa_matrix->plus(i, l, Coancestry(P1->get(FEM, age_pos, j)->getTrait(_SHLinkedTraitIndex)->get_sequence(),
                                                 P2->get(MAL, age_pos, k)->getTrait(_SHLinkedTraitIndex)->get_sequence(),
                                                 nb_locus));
          //males i vs. males l
          size_i = P1->size(MAL, age_pos);
          size_l = P2->size(MAL, age_pos);
          nb_coeff += size_i * size_l;
          
          for(j = 0; j < size_i; ++j)
            for(k = 0; k < size_l; ++k)
              _coa_matrix->plus(i, l, Coancestry(P1->get(MAL, age_pos, j)->getTrait(_SHLinkedTraitIndex)->get_sequence(),
                                                 P2->get(MAL, age_pos, k)->getTrait(_SHLinkedTraitIndex)->get_sequence(),
                                                 nb_locus));
          //males i vs. females l
          size_l = P2->size(FEM, age_pos);
          nb_coeff += size_i * size_l;
          
          for(j = 0; j < size_i; ++j)
            for(k = 0; k < size_l; ++k)
              _coa_matrix->plus(i, l, Coancestry(P1->get(MAL, age_pos, j)->getTrait(_SHLinkedTraitIndex)->get_sequence(),
                                                 P2->get(FEM, age_pos, k)->getTrait(_SHLinkedTraitIndex)->get_sequence(),
                                                 nb_locus));
          _coa_matrix->divide(i, l, nb_coeff);
        }//endif
        
        _mean_alpha += tot_size * _coa_matrix->get(i, l);
        
        wt += tot_size;

      }//end for P2
    }//end for P1 
    
    //weighted average:
    _mean_alpha /= wt;
    
  }//end if upper half

}

void TTNeutralGenesSH::setCoaMatrixRecorders	(	age_t	AGE,
		unsigned char	dim
	)

Definition at line 1293 of file ttneutralgenes.cc.

References StatHandlerBase::_pop, StatHandler< TTNeutralGenesSH >::add(), ADULTS, getCoa(), getMeanAlpha(), getMeanTheta(), Metapop::getPatchNbr(), setAdultsCoaBetween(), setAdultsCoaMatrix(), setAdultsCoaWithin(), setOffsprgCoaBetween(), setOffsprgCoaMatrix(), and setOffsprgCoaWithin().

Referenced by setStatRecorders().

{
  std::ostringstream name, sub_name;
  
  void (TTNeutralGenesSH::* setter) () = (AGE == ADULTS ?
                                          ( dim & 3 ? &TTNeutralGenesSH::setAdultsCoaMatrix : 
                                            dim & 2 ? &TTNeutralGenesSH::setAdultsCoaBetween : 
                                            &TTNeutralGenesSH::setAdultsCoaWithin ) :
                                          ( dim & 3 ? &TTNeutralGenesSH::setOffsprgCoaMatrix :
                                            dim & 2 ? &TTNeutralGenesSH::setOffsprgCoaBetween :
                                            &TTNeutralGenesSH::setOffsprgCoaWithin) );
  
  const char *prefix = (AGE == ADULTS ? "adlt." : "off.");
  
  unsigned int nbpatch =  _pop->getPatchNbr();
  unsigned int scale = (unsigned int)pow(10.0, (int)log10((float)nbpatch) + 1);
  
  if(dim & 1) {
    name<<"Wtn Patch Coancestry ("<<prefix<<")";
    sub_name<< prefix << "theta";
    add(name.str(), sub_name.str(),  AGE,0, 0, &TTNeutralGenesSH::getMeanTheta, 0, 0, setter);
    name.str("");
    sub_name.str("");
  }
  if(dim & 2){
    name<<"Btn Patch Coancestry ("<<prefix<<")";
    sub_name<< prefix << "alpha";
    add(name.str(), sub_name.str(),  AGE, 0, 0, &TTNeutralGenesSH::getMeanAlpha, 0, 0, (!(dim&1) ? setter : 0) );
    name.str("");
    sub_name.str("");
  }  
  if(dim & 1) {
    for(unsigned int i = 0; i < nbpatch; ++i) {
      name<<"coancestry "<<i+1<<"."<<i+1;
      sub_name<< prefix << "coa" << i+1 << "." << i+1;
      add(name.str(), sub_name.str(),  AGE, i*scale + i, 0, 0, &TTNeutralGenesSH::getCoa, 0, 0);
      name.str("");
      sub_name.str("");
    }
  }
  if(dim & 2){
    for(unsigned int i = 0; i < nbpatch; ++i) {
      for(unsigned int j = i+1; j < nbpatch; ++j) {
        name<<"coancestry "<<i+1<< "." <<j+1;
        sub_name<< prefix << "coa" << i+1 << "." << j+1;
        add(name.str(), sub_name.str(),  AGE, i*scale + j, 0, 0, &TTNeutralGenesSH::getCoa, 0, 0);
        name.str("");
        sub_name.str("");
      }
    }
  }
  
}

void TTNeutralGenesSH::setFst_li	(	unsigned int	N,
		unsigned int	L,
		double **	array
	)

Computes the per-locus per-patch Fst values using Weir&Hill 2002 approach.

Definition at line 584 of file stats_fstat.cc.

Referenced by TTNeutralGenesFH::write_Fst_i().

{
  unsigned int patchNbr = this->_pop->getPatchNbr();
  unsigned int nb_locus = this->_SHLinkedTrait->get_locus_num();
  unsigned int nb_allele = this->_SHLinkedTrait->get_allele_num();
  
  assert(N == patchNbr);
  assert(L == nb_locus);
  
  setAlleleTables(ADULTS);
  
  double *pop_weights = new double[patchNbr];
  double *pop_sizes = new double[patchNbr];
  double *denominator = new double[nb_locus];
  double sum_weights = 0;
  
  double tot_size = _pop->size(ADULTS) * 2; //diploids!
  
  for(unsigned int i = 0; i < patchNbr; ++i) {
    pop_sizes[i] = _pop->size(ADULTS, i) * 2; 
    pop_weights[i] = pop_sizes[i] - (pop_sizes[i] * pop_sizes[i] / tot_size); //n_ic in Weir & Hill 2002
    sum_weights += pop_weights[i];
    for(unsigned int j = 0; j < nb_locus; ++j)
      array[i][j] = nanf("NULL");
  }
  
  for(unsigned int l = 0; l < nb_locus; ++l)
    denominator[l] = 0;
  
  double p, pq, var;
  
  for (unsigned int i = 0; i < patchNbr; ++i) {
    
    if( !pop_sizes[i] ) continue;
    
    for (unsigned int l = 0; l < nb_locus; ++l) {
      
      array[i][l] = 0;
      
      for (unsigned int u = 0; u < nb_allele; ++u) {
        
        p = _alleleFreqTable.get(i, l, u); //p_liu
        
        pq = p * (1 - p);
        
        var = p - _globalAlleleFreq.get(l, u); //(p_liu - pbar_u)
        
        var *= var;
        
        array[i][l] += pq * pop_sizes[i] / (pop_sizes[i] -1);
        
        denominator[l] += pop_sizes[i] * var + pop_weights[i] * pq;
        
      } // end for allele
    }// end for locus
  }//end for pop
  
  
  for (unsigned int i = 0; i < patchNbr; ++i) {
    if( !pop_sizes[i] ) continue;
    for (unsigned int l = 0; l < nb_locus; ++l) 
      array[i][l] = 1 - ( array[i][l] * sum_weights / denominator[l]);
  }
 
  delete [] pop_weights;
  delete [] pop_sizes;
  delete [] denominator;
}

void TTNeutralGenesSH::setFstat ( age_t AGE )

Computes the F-statistics following Nei & Chesser (1983).

Definition at line 120 of file stats_fstat.cc.

Referenced by setAdultsFstat(), and setOffsprgFstat().

{
  double harmonic = 0, nbpatch = 0, nbind;
  unsigned int patchNbr = _pop->getPatchNbr();
  age_idx age_pos = (AGE == ADULTS ? ADLTx : OFFSx);
  
  //check if the population size has changed:
  if(_alleleCountTable.getNumGroups() != patchNbr)
    allocateTables(_SHLinkedTrait->get_locus_num(), _SHLinkedTrait->get_allele_num());
  
  setAlleleTables(AGE);
  
  setLociDivCounter(AGE);

  //harmonic mean of patch sizes:
  for (unsigned int i = 0; i < patchNbr; ++i){
        
    nbind = _pop->size(AGE, i);
        if(nbind != 0){
      nbpatch++;
          harmonic += 1.0/nbind;
    }
  }

  harmonic = nbpatch / harmonic;
  
  _ho = setHo(age_pos);
  _hs = setHs(age_pos);
  _ht = setHt(age_pos);
  //Nei's corrections:
  _hsnei = (nbpatch != 0 ? harmonic/(harmonic-1.0)*(_hs-(_ho/(2.0*harmonic))) : nanf("NULL") );
  _htnei = (nbpatch != 0 ? _ht + (_hsnei/(harmonic*nbpatch))
            -(_ho/(2.0*harmonic*nbpatch)) : nanf("NULL") );
  
  _fis = ( _hsnei ? 1.0-(_ho/_hsnei) : nanf("NULL") );
  _fit = ( _htnei ? 1.0-(_ho/_htnei) : nanf("NULL") );
  _fst = ( _htnei ? 1.0-(_hsnei/_htnei) : nanf("NULL") );
}

void TTNeutralGenesSH::setFstatRecorders ( age_t AGE )

Definition at line 1349 of file ttneutralgenes.cc.

References StatHandler< TTNeutralGenesSH >::add(), ADULTS, getFis(), getFit(), getFixLocGlobal(), getFixLocLocal(), getFst(), getHo(), getHsnei(), getHtnei(), getNbAllGlobal(), getNbAllLocal(), OFFSPRG, setAdultsFstat(), and setOffsprgFstat().

Referenced by setStatRecorders().

{
  if(AGE & OFFSPRG) {
    add("Nbr of Alleles per Locus - local  (offsprg)","off.allnbp",OFFSPRG,0,0,
        &TTNeutralGenesSH::getNbAllLocal,0,0,&TTNeutralGenesSH::setOffsprgFstat);
    add("Nbr of Alleles per Locus - global (offsprg)","off.allnb",OFFSPRG,0,0,&TTNeutralGenesSH::getNbAllGlobal,0,0,0);
    add("Nbr of Fixed Loci per Patch  (offsprg)","off.fixlocp",OFFSPRG,0,0,&TTNeutralGenesSH::getFixLocLocal,0,0,0);
    add("Nbr of Fixed Loci in the Pop (offsprg)","off.fixloc",OFFSPRG,0,0,&TTNeutralGenesSH::getFixLocGlobal,0,0,0);
    add("Ho     (offsprg)","off.ho",OFFSPRG,0,0,&TTNeutralGenesSH::getHo,0,0,0);
    add("Hs-Nei (offsprg)","off.hsnei",OFFSPRG,0,0,&TTNeutralGenesSH::getHsnei,0,0,0);
    add("Ht-Nei (offsprg)","off.htnei",OFFSPRG,0,0,&TTNeutralGenesSH::getHtnei,0,0,0);
    add("Fis    (offsprg)","off.fis",OFFSPRG,0,0,&TTNeutralGenesSH::getFis,0,0,0);
    add("Fst    (offsprg)","off.fst",OFFSPRG,0,0,&TTNeutralGenesSH::getFst,0,0,0);
    add("Fit    (offsprg)","off.fit",OFFSPRG,0,0,&TTNeutralGenesSH::getFit,0,0,0);
  }
  
  if(AGE & ADULTS) {
    add("Nbr of Alleles per Locus - local  (adult)","adlt.allnbp",ADULTS,0,0,
        &TTNeutralGenesSH::getNbAllLocal,0,0,&TTNeutralGenesSH::setAdultsFstat);
    add("Nbr of Alleles per Locus - global (adult)","adlt.allnb",ADULTS,0,0,&TTNeutralGenesSH::getNbAllGlobal,0,0,0);
    add("Nbr of Fixed Loci per Patch  (adult)","adlt.fixlocp",ADULTS,0,0,&TTNeutralGenesSH::getFixLocLocal,0,0,0);
    add("Nbr of Fixed Loci in the Pop (adult)","adlt.fixloc",ADULTS,0,0,&TTNeutralGenesSH::getFixLocGlobal,0,0,0);
    add("Ho       (adult)","adlt.ho",ADULTS,0,0,&TTNeutralGenesSH::getHo,0,0,0);
    add("Hs-Nei   (adult)","adlt.hsnei",ADULTS,0,0,&TTNeutralGenesSH::getHsnei,0,0,0);
    add("Ht-Nei   (adult)","adlt.htnei",ADULTS,0,0,&TTNeutralGenesSH::getHtnei,0,0,0);
    add("Fis      (adult)","adlt.fis",ADULTS,0,0,&TTNeutralGenesSH::getFis,0,0,0);
    add("Fst      (adult)","adlt.fst",ADULTS,0,0,&TTNeutralGenesSH::getFst,0,0,0);
    add("Fit      (adult)","adlt.fit",ADULTS,0,0,&TTNeutralGenesSH::getFit,0,0,0);
  }
}

void TTNeutralGenesSH::setFstatWCRecorders ( age_t AGE )

Definition at line 1382 of file ttneutralgenes.cc.

References StatHandler< TTNeutralGenesSH >::add(), ADULTS, getFisWC(), getFitWC(), getFstWC(), OFFSPRG, setAdultsFstatWeirCockerham(), and setOffspringFstatWeirCockerham().

Referenced by setStatRecorders().

{
  if(AGE & OFFSPRG) {
    
    add("Fis Weir & Cockerham","off.fis.WC",OFFSPRG,0,0,&TTNeutralGenesSH::getFisWC,0,0,&TTNeutralGenesSH::setOffspringFstatWeirCockerham);
    add("Fst Weir & Cockerham","off.fst.WC",OFFSPRG,0,0,&TTNeutralGenesSH::getFstWC,0,0,0);
    add("Fit Weir & Cockerham","off.fit.WC",OFFSPRG,0,0,&TTNeutralGenesSH::getFitWC,0,0,0);
  }
  
  if(AGE & ADULTS) {
    
    add("Fis Weir & Cockerham","adlt.fis.WC",ADULTS,0,0,&TTNeutralGenesSH::getFisWC,0,0,&TTNeutralGenesSH::setAdultsFstatWeirCockerham);
    add("Fst Weir & Cockerham","adlt.fst.WC",ADULTS,0,0,&TTNeutralGenesSH::getFstWC,0,0,0);
    add("Fit Weir & Cockerham","adlt.fit.WC",ADULTS,0,0,&TTNeutralGenesSH::getFitWC,0,0,0);
    
  }
}

void TTNeutralGenesSH::setFstatWeirCockerham ( age_t AGE )

Computes the Weir & Cockerham (1984) Fstat values (Theta, F, and f).

Definition at line 456 of file stats_fstat.cc.

Referenced by setAdultsFstatWeirCockerham(), and setOffspringFstatWeirCockerham().

{
  unsigned int patchNbr = this->_pop->getPatchNbr();
  unsigned int nb_locus = this->_SHLinkedTrait->get_locus_num();
  unsigned int nb_allele = this->_SHLinkedTrait->get_allele_num();
  age_idx age_pos = (AGE == ADULTS ? ADLTx : OFFSx);
  
  setAlleleTables(AGE);
  
  DataTable< unsigned int > heteroTable;
  unsigned int **sizes;
  
  sizes = new unsigned int * [patchNbr];
  
  for(unsigned int i = 0; i < patchNbr; ++i) {
    sizes[i] = new unsigned int [nb_locus];
    for(unsigned int j = 0; j < nb_locus; ++j)
      sizes[i][j] = nb_allele;
  }
  
  heteroTable.allocate(patchNbr,nb_locus, sizes);
  heteroTable.init(0);
  
  Patch* patch;
  unsigned int isHetero;
  unsigned int extantPs = 0;
  unsigned char** genes;
  
  for (unsigned int i = 0; i < patchNbr; ++i) {
        
    patch = _pop->getPatch(i);
    
    if(patch->size(age_pos)) extantPs++;
    else continue;
    
        for(unsigned int j = 0, size = patch->size(FEM, age_pos); j < size; ++j) {
          
          genes = (unsigned char**)patch->get(FEM, age_pos, j)->getTrait(_SHLinkedTraitIndex)->get_sequence();
          
          for (unsigned int l = 0; l < nb_locus; ++l) {
        isHetero = (genes[0][l] != genes[1][l]);
        heteroTable.plus(i, l, genes[0][l], isHetero);
        heteroTable.plus(i, l, genes[1][l], isHetero);
      }
        }
    
    for(unsigned int j = 0, size = patch->size(MAL, age_pos); j < size; ++j) {
      
      genes = (unsigned char**)patch->get(MAL, age_pos, j)->getTrait(_SHLinkedTraitIndex)->get_sequence();
      
      for (unsigned int l = 0; l < nb_locus; ++l) {
        isHetero = (genes[0][l] != genes[1][l]);
        heteroTable.plus(i, l, genes[0][l], isHetero);
        heteroTable.plus(i, l, genes[1][l], isHetero);
      }
    }
  }
  
  //init
  double *pop_sizes = new double [patchNbr];
  double tot_size, inv_ntot;
  double sum_weights = 0;
  double nbar, nc, inv_nbar;
  
  tot_size = _pop->size(AGE);
  
  for(unsigned int i = 0; i < patchNbr; i++) {
    pop_sizes[i] = _pop->size(AGE, i); 
    sum_weights += (pop_sizes[i] * pop_sizes[i] / tot_size);
  }
  nbar = tot_size/extantPs;
  nc = (tot_size - sum_weights)/(extantPs-1);
  inv_nbar = 1/(nbar - 1);
  inv_ntot = 1/tot_size;
  
  double var;
  double s2, pbar, hbar;
  double s2_denom = 1.0/((extantPs-1)*nbar), r = (double)(extantPs-1)/extantPs, hbar_factor=(2*nbar-1)/(4*nbar);
  double a = 0, b = 0, c = 0, x;
  
  for (unsigned int l = 0; l < nb_locus; ++l) {
    
    for (unsigned int u = 0; u < nb_allele; ++u) {
      
      s2 = pbar = hbar = 0;
      
      for (unsigned int i = 0; i < patchNbr; ++i) {
        
        var = _alleleFreqTable.get(i, l, u) - _globalAlleleFreq.get(l, u); //(p_liu - pbar_u)^2 
        
        var *= var;
        
        s2 += var * pop_sizes[i];
        
        hbar += heteroTable.get(i, l, u);
        
      }//end for pop
      
      s2 *= s2_denom;
      pbar = _globalAlleleFreq.get(l, u);
      hbar *= inv_ntot;
      
      x = pbar * (1 - pbar) - r * s2;
      a += s2 - inv_nbar*( x - 0.25 * hbar);
      b += x - hbar_factor * hbar;
      c += hbar;
      
      
    } // end for allele
    
  }// end for locus
    
  a *= nbar/nc;
  b *= nbar/(nbar - 1);
  c *= 0.5;
  
  _fst_WC = a / (a + b + c);
  _fit_WC = (a + b) / (a + b + c);
  _fis_WC = b / (b + c);
  
  
  for(unsigned int i = 0; i < patchNbr; ++i) delete [] sizes[i];
  delete [] sizes;
  delete [] pop_sizes;   
}

void TTNeutralGenesSH::setFstMatrix	(	age_t	AGE,
		unsigned char	dim
	)

Computes the weighted within and between patch Fst's as well as the overall Fst (Theta).

The method used here is that of Weir & Hill 2002, Ann. Rev. Genet. 36:721-750. The weighting is done for samples (patches) of unequal sizes.

Parameters:

AGE	the age class
dim	the dimension of the matrix to fill: 1 = the diagonal (i.e. the wihtin patch Fst or theta_ii) 2 = the upper half (i.e. the between patch Fst or theta_ii') 3 = both

Definition at line 349 of file stats_fstat.cc.

Referenced by setAdultsFstBetween(), setAdultsFstMatrix(), setAdultsFstWithin(), setOffsprgFstBetween(), setOffsprgFstMatrix(), and setOffsprgFstWithin().

{
  unsigned int patchNbr = this->_pop->getPatchNbr();
  unsigned int nb_locus = this->_SHLinkedTrait->get_locus_num();
  unsigned int nb_allele = this->_SHLinkedTrait->get_allele_num();

  setAlleleTables(AGE);
  
  if(_fst_matrix == NULL)
    
    _fst_matrix = new TMatrix(patchNbr, patchNbr);
  
  else if( _fst_matrix->length() != patchNbr * patchNbr)
    
    _fst_matrix->reset(patchNbr, patchNbr);
  
  _fst_matrix->assign(nanf("NULL"));
  
  //init
  double *pop_weights = new double[patchNbr];
  double *pop_sizes = new double[patchNbr];
  double **numerator = new double*[patchNbr];
  for(unsigned int i = 0; i < patchNbr; i++) numerator[i] = new double [patchNbr];
  double tot_size;
  double numerator_W = 0;
  double denominator = 0;
  double sum_weights = 0;
  
  tot_size = _pop->size(AGE) * 2; 
  
  for(unsigned int i = 0; i < patchNbr; ++i) {
    pop_sizes[i] = _pop->size(AGE, i) * 2; 
    pop_weights[i] = pop_sizes[i] - (pop_sizes[i] * pop_sizes[i] / tot_size); //n_ic in Weir & Hill 2002
    sum_weights += pop_weights[i];
    for(unsigned int j = 0; j < patchNbr; j++)
      numerator[i][j] = 0;
  }
  
   double p, pq, var, num;

   for (unsigned int i = 0; i < patchNbr; ++i) {
     
     if( !pop_sizes[i] ) continue;
     
     for (unsigned int l = 0; l < nb_locus; ++l) {
       
       for (unsigned int u = 0; u < nb_allele; ++u) {
         
         p = _alleleFreqTable.get(i, l, u); //p_liu
         
         pq = p * (1 - p);
         
         var = p - _globalAlleleFreq.get(l, u); //(p_liu - pbar_u)^2 
         
         var *= var;
         
         num = pq * pop_sizes[i] / (pop_sizes[i] -1);

         numerator[i][i] += num;
         
         numerator_W += num * pop_sizes[i]; //see equ. 9, Weir & Hill 2002
         
         denominator += pop_sizes[i] * var + pop_weights[i] * pq; //common denominator
         
       } // end for allele
     }// end for locus
   }//end for pop

  for (unsigned int i = 0; i < patchNbr; ++i) {
    if( !pop_sizes[i] ) continue;
    _fst_matrix->set(i, i, 1 - (numerator[i][i] * sum_weights / denominator) );
  }   
   _fst_WH = 1 - ((numerator_W * sum_weights) / (denominator * tot_size)); //equ. 9 Weir & Hill 2002
   
   //pairwise Fst:
   if(dim & 2) {
     double pi, pj;
     for (unsigned int l = 0; l < nb_locus; ++l)
       for (unsigned int u = 0; u < nb_allele; ++u)
         for (unsigned int i = 0; i < patchNbr - 1; ++i) {
           if( !pop_sizes[i] ) continue;
           for (unsigned int j = i + 1; j < patchNbr; ++j) {
             if( !pop_sizes[j] ) continue;
             pi = _alleleFreqTable.get(i, l, u);
             pj = _alleleFreqTable.get(j, l, u);
             numerator[i][j] += pi * (1 - pj) + pj * (1 - pi); //equ. 7 of Weir & Hill 2002
           }
         }
     for (unsigned int i = 0; i < patchNbr - 1; ++i){
       if( !pop_sizes[i] ) continue;
       for (unsigned int j = i + 1; j < patchNbr; ++j){
         if( !pop_sizes[j] ) continue;
         _fst_matrix->set(i, j, 1 - ( (numerator[i][j] * sum_weights) / (2 * denominator)) );
       }
     }
   }
 
  
  delete [] pop_weights;
  delete [] pop_sizes;
  for(unsigned int i = 0; i < patchNbr; i++) delete [] numerator[i];
  delete [] numerator;
}

void TTNeutralGenesSH::setFstMatrixRecorders	(	age_t	AGE,
		unsigned char	dim
	)

Definition at line 1402 of file ttneutralgenes.cc.

References StatHandlerBase::_pop, StatHandler< TTNeutralGenesSH >::add(), ADULTS, getFst_ij(), Metapop::getPatchNbr(), getWeightedFst(), setAdultsFstBetween(), setAdultsFstMatrix(), setAdultsFstWithin(), setOffsprgFstBetween(), setOffsprgFstMatrix(), and setOffsprgFstWithin().

Referenced by setStatRecorders().

{
  std::ostringstream name, sub_name;
  
  void (TTNeutralGenesSH::* setter) () = (AGE == ADULTS ?
                                          ( dim & 3 ? &TTNeutralGenesSH::setAdultsFstMatrix : 
                                            dim & 2 ? &TTNeutralGenesSH::setAdultsFstBetween : 
                                            &TTNeutralGenesSH::setAdultsFstWithin ) :
                                          ( dim & 3 ? &TTNeutralGenesSH::setOffsprgFstMatrix :
                                            dim & 2 ? &TTNeutralGenesSH::setOffsprgFstBetween :
                                            &TTNeutralGenesSH::setOffsprgFstWithin) );
  
  const char *prefix = (AGE == ADULTS ? "adlt." : "off.");
  
  unsigned int nbpatch =  _pop->getPatchNbr();
  unsigned int scale = (unsigned int)pow(10.0, (int)log10((float)nbpatch) + 1);
  
  name<<"Weir&Hill weighted Fst ("<<prefix<<")";
  sub_name<< prefix << "fst.WH";
  add(name.str(), sub_name.str(),  AGE, 0, 0, &TTNeutralGenesSH::getWeightedFst, 0, 0, setter);
  name.str("");
  sub_name.str("");
  
  if(dim & 1) {
    for(unsigned int i = 0; i < nbpatch; ++i) {
      name<<"Weighted Fst "<<i+1<<"."<<i+1;
      sub_name<< prefix << "fst" << i+1 << "." << i+1;
      add(name.str(), sub_name.str(),  AGE, i*scale + i, 0, 0, &TTNeutralGenesSH::getFst_ij, 0,  0);
      name.str("");
      sub_name.str("");
    }
  }
  if(dim & 2){
    for(unsigned int i = 0; i < nbpatch; ++i) {
      for(unsigned int j = i+1; j < nbpatch; ++j) {
        name<<"Weighted Fst "<<i+1<< "." <<j+1;
        sub_name<< prefix << "fst" << i+1 << "." << j+1;
        add(name.str(), sub_name.str(),  AGE, i*scale + j, 0, 0, &TTNeutralGenesSH::getFst_ij, 0, 0);
        name.str("");
        sub_name.str("");
      }
    }
  }
  
}

double TTNeutralGenesSH::setHo ( age_idx age_pos )

Definition at line 229 of file stats_fstat.cc.

{
  unsigned int nloc = _SHLinkedTrait->get_locus_num(), nbpatch = _pop->getPatchNbr(), psize;
  double indloc = 0, hetero = 0;
  unsigned char** genes;
  Patch* current_patch;
  
  indloc = _pop->size(age_pos) * _SHLinkedTrait->get_locus_num();
  
  for (unsigned int i = 0; i < nbpatch; ++i) {
        
    current_patch = _pop->getPatch(i);
    
        psize = current_patch->size(FEM, age_pos);
    
    indloc += psize;
    
        for(unsigned int j = 0; j < psize; ++j) {
          
          genes = (unsigned char**)current_patch->get(FEM, age_pos, j)->getTrait(_SHLinkedTraitIndex)->get_sequence();
          
          for (unsigned int k = 0; k < nloc; ++k)  hetero += (genes[0][k] != genes[1][k]);
        }
  }
  //doing some unrolling, loops for the males here:
  for (unsigned int i = 0; i < nbpatch; ++i) {
        
    current_patch = _pop->getPatch(i);
        
    psize = current_patch->size(MAL, age_pos);
    
    indloc += psize;
    
        for(unsigned int j = 0; j < psize; ++j) {
          
          genes = (unsigned char**)current_patch->get(MAL, age_pos, j)->getTrait(_SHLinkedTraitIndex)->get_sequence();
          
          for (unsigned int k = 0; k < nloc; ++k)  hetero += (genes[0][k] != genes[1][k]);
        }
  }
  
  indloc *= nloc;

  return (indloc != 0 ? hetero/indloc : 0.0);
}

double TTNeutralGenesSH::setHs ( age_idx age_pos )

Definition at line 278 of file stats_fstat.cc.

{
  unsigned int i, k, x, nb_patch=0;
  double hs = 0;
  unsigned int nbLoc = _SHLinkedTrait->get_locus_num(), nbAll = _SHLinkedTrait->get_allele_num();
  unsigned int patchNbr = _pop->getPatchNbr();
  deque<double>genehet(nbLoc,1);
  deque<double>hspop(patchNbr,0);
  double freq;
  Patch* current_patch;
  
  //compute the expected heterozygosity for each patch and return average:
  for (i = 0; i < patchNbr; ++i) {
    
    current_patch = _pop->getPatch(i);
                
        genehet.assign(nbLoc, 1);
    
        if(current_patch->size(age_pos) != 0) {
      
      nb_patch++;
      
      for (k = 0; k < nbLoc; ++k) {
        
                for (x = 0; x < nbAll; ++x) {
          
          freq = _alleleFreqTable.get(i, k, x);
          
                  freq *= freq; //squared frequencies (expected _homozygosity)
          
                  genehet[k] -= freq; //1 - sum of p2 = expected heterozygosity
                }
                //expected heterozygosity:
        hspop[i] += genehet[k];
      }
    }//end_if size!=0
    
        hs += hspop[i];
  }//end patch loop
  
  return (nb_patch !=0 ? hs/(nbLoc*nb_patch) : 0.0);
}

double TTNeutralGenesSH::setHt ( age_idx age_pos )

Definition at line 323 of file stats_fstat.cc.

{
  double ht = 0;
  unsigned int nbLoc = _SHLinkedTrait->get_locus_num(), nbAll = _SHLinkedTrait->get_allele_num();
  deque<double> genehet(nbLoc, 1);
  double freq;
  
  //get global allele frequencies per locus
  for (unsigned int l = 0; l < nbLoc; ++l) {
    for (unsigned int u = 0; u < nbAll; ++u) {
      
      freq = _globalAlleleFreq.get(l, u);
      
      freq *= freq; //squared frequencies
      
          genehet[l] -= freq;  //1 - sum of p2 = expected heterozygosity
    }
    
    ht += genehet[l];
  }
  
  return (ht/nbLoc);
}

void TTNeutralGenesSH::setLociDivCounter ( age_t AGE )

Sets the allelic diversity counters.

Definition at line 161 of file stats_fstat.cc.

{
  unsigned int i, j, k;
  unsigned int nbLoc = _SHLinkedTrait->get_locus_num(), nbAll = _SHLinkedTrait->get_allele_num();
  unsigned int nbpatch = 0, patchNbr = _pop->getPatchNbr();
  double patch_mean, pop_mean = 0;
  bool **pop_div;
  
  //number of alleles per locus, Patch and pop counters:
  pop_div = new bool * [nbLoc];
  
  for(i = 0; i < nbLoc; ++i) {
    pop_div[i] = new bool [nbAll];
    for(j = 0; j < nbAll;++j)
      pop_div[i][j] = 0;
  }
  
  for(k = 0; k < patchNbr; ++k) 
    nbpatch += (_pop->size(AGE, k) != 0);
  
  for(k = 0; k < patchNbr; ++k) {
        
    patch_mean = 0;
        
    for(i = 0; i < nbLoc; ++i)    
          for(j = 0; j < nbAll;++j) {
                patch_mean += (_alleleCountTable.get(k,i,j) != 0);
        pop_div[i][j] |= (_alleleCountTable.get(k,i,j) != 0);
      }
        //add mean nb of alleles per locus for Patch k to the pop mean
        pop_mean += patch_mean/nbLoc;
  }
  
  _nb_all_local = (nbpatch ? pop_mean/nbpatch : nanf("NULL"));

  _nb_all_global = 0;
  
  for(i = 0; i < nbLoc; ++i)      
    for(j = 0; j < nbAll;++j)
      _nb_all_global += pop_div[i][j];
  
  _nb_all_global /= nbLoc;
  
  for(i = 0; i < nbLoc; ++i)
    delete [] pop_div[i];
  delete [] pop_div;
      
  //number of fixed loci, local and global counters:
  _fix_loc_local = 0;
  
  for (k = 0; k < patchNbr; ++k)
    for (i = 0; i < nbLoc; ++i)
      for (j = 0; j < nbAll; ++j)        
        _fix_loc_local += ( _alleleFreqTable.get(k, i, j) == 1 );
  
  _fix_loc_local /= nbpatch;
  
  _fix_loc_global = 0;
      
  //globally:  
  for (i = 0; i < nbLoc; ++i)
    for (j = 0; j < nbAll; ++j)  
      _fix_loc_global += ( _globalAlleleFreq.get(i, j) == 1 );

}

void TTNeutralGenesSH::setNeiGeneticDistance ( age_t AGE )

Definition at line 656 of file stats_fstat.cc.

Referenced by setAdltNeiGeneticDistance(), and setOffsprgNeiGeneticDistance().

{
  //see: Nei, M. 1978. Genetics 89:583-590
  unsigned int patchNbr = this->_pop->getPatchNbr();
  unsigned int nb_locus = this->_SHLinkedTrait->get_locus_num();
  unsigned int nb_allele = this->_SHLinkedTrait->get_allele_num();
  
  setAlleleTables(AGE);
  
  if(_D == NULL) 
    
    _D = new TMatrix(patchNbr, patchNbr);
  
  else if( _D->length() != patchNbr*patchNbr )
    
    _D->reset(patchNbr, patchNbr);
  
  _D->assign(nanf("NULL"));
  
  double num, denom, p1, p2, sum;
  double *sum_square = new double [patchNbr];
  double *sample_size = new double [patchNbr];
  
  for(unsigned int i = 0; i < patchNbr; ++i) {
    
    sample_size[i] = 2 * _pop->size(AGE, i);
    
    sum_square[i] = 0;
    
    if( !sample_size[i] ) continue;
    
    for(unsigned int l = 0; l < nb_locus; ++l) {
      
      sum = 0;
      
      for(unsigned int u = 0; u < nb_allele ; ++u) {
        p1 = _alleleFreqTable.get(i, l, u);
        sum += p1*p1;
      }
      
      sum_square[i] += (sample_size[i] * sum -1) / (sample_size[i] -1); //unbiased estimate
    }
  }
  
  unsigned int pairs = 0;
  double Dpair;
  _meanD = 0;
  
  for(unsigned int i = 0; i < patchNbr-1; ++i) {
    
    if( !sample_size[i] ) continue;
    
    for(unsigned int j = i + 1; j < patchNbr; ++j) {

      if( !sample_size[j] ) continue;
      
      num = 0;
      pairs++;
      
      for(unsigned int l = 0; l < nb_locus; ++l) {
        for(unsigned int u = 0; u < nb_allele ; ++u) {
          p1 = _alleleFreqTable.get(i, l, u);
          p2 = _alleleFreqTable.get(j, l, u);
          num += p1*p2;
        }
      }      
      
      denom = sqrt(sum_square[i] * sum_square[j]);
      Dpair = -log(num/denom);
      _meanD += Dpair;
      _D->set(i, j, Dpair);
    }
  }
  _meanD /= pairs;
  
  delete [] sum_square;
  delete [] sample_size;
}

void TTNeutralGenesSH::setNeiGeneticDistanceRecorders	(	age_t	AGE,
		bool	pairwise
	)

Definition at line 1450 of file ttneutralgenes.cc.

References StatHandlerBase::_pop, StatHandler< TTNeutralGenesSH >::add(), ADULTS, getMeanNeiGeneticDistance(), getNeiGeneticDistance(), Metapop::getPatchNbr(), setAdltNeiGeneticDistance(), and setOffsprgNeiGeneticDistance().

Referenced by setStatRecorders().

{
  std::ostringstream name, sub_name;
  unsigned int nbpatch =  _pop->getPatchNbr();
  const char *prefix = (AGE == ADULTS ? "adlt." : "off.");
  unsigned int scale = (unsigned int)pow(10.0, (int)log10((float)nbpatch) + 1);

  void (TTNeutralGenesSH::* setter) () = (AGE == ADULTS ?
                                          &TTNeutralGenesSH::setAdltNeiGeneticDistance 
                                          : &TTNeutralGenesSH::setOffsprgNeiGeneticDistance );

  sub_name<< prefix << "D";
  add("Average between pop Nei's D", sub_name.str(),  AGE, 0, 0, 
      &TTNeutralGenesSH::getMeanNeiGeneticDistance, 0, 0, setter);
  sub_name.str("");
  
  if(pairwise) {
    
    for(unsigned int i = 0; i < nbpatch -1; i++){
      for(unsigned int j = i+1; j < nbpatch; j++) {
        name<<"Nei's D between pop"<<i+1<<" and pop"<<j+1;
        sub_name<< prefix << "D" << i+1 << "." << j+1;
        add(name.str(), sub_name.str(),  AGE, i*scale + j, 0, 0, &TTNeutralGenesSH::getNeiGeneticDistance, 0, 0);
        name.str("");
        sub_name.str("");
      }
    }
  }
}

void TTNeutralGenesSH::setOffsprgCoaBetween ( ) [inline]

Definition at line 370 of file ttneutralgenes.h.

References OFFSx, and setCoaMatrix().

Referenced by setCoaMatrixRecorders(), and setStatRecorders().

{setCoaMatrix(OFFSx, 2);}

void TTNeutralGenesSH::setOffsprgCoaMatrix ( ) [inline]

Definition at line 366 of file ttneutralgenes.h.

References OFFSx, and setCoaMatrix().

Referenced by setCoaMatrixRecorders(), and setStatRecorders().

{setCoaMatrix(OFFSx, 3);}

void TTNeutralGenesSH::setOffsprgCoaWithin ( ) [inline]

Definition at line 368 of file ttneutralgenes.h.

References OFFSx, and setCoaMatrix().

Referenced by setCoaMatrixRecorders(), and setStatRecorders().

{setCoaMatrix(OFFSx, 1);}

void TTNeutralGenesSH::setOffsprgFstat ( ) [inline]

Definition at line 320 of file ttneutralgenes.h.

References OFFSPRG, and setFstat().

Referenced by setFstatRecorders().

{setFstat(OFFSPRG);}

void TTNeutralGenesSH::setOffsprgFstBetween ( ) [inline]

Definition at line 299 of file ttneutralgenes.h.

References OFFSPRG, and setFstMatrix().

Referenced by setFstMatrixRecorders().

{setFstMatrix(OFFSPRG, 2);}

void TTNeutralGenesSH::setOffsprgFstMatrix ( ) [inline]

Definition at line 297 of file ttneutralgenes.h.

References OFFSPRG, and setFstMatrix().

Referenced by setFstMatrixRecorders().

{setFstMatrix(OFFSPRG, 3);}

void TTNeutralGenesSH::setOffsprgFstWithin ( ) [inline]

Definition at line 298 of file ttneutralgenes.h.

References OFFSPRG, and setFstMatrix().

Referenced by setFstMatrixRecorders().

{setFstMatrix(OFFSPRG, 1);}

void TTNeutralGenesSH::setOffsprgNeiGeneticDistance ( ) [inline]

Definition at line 399 of file ttneutralgenes.h.

References OFFSPRG, and setNeiGeneticDistance().

Referenced by setNeiGeneticDistanceRecorders().

{setNeiGeneticDistance(OFFSPRG);}

void TTNeutralGenesSH::setOffspringFstatWeirCockerham ( ) [inline]

Definition at line 335 of file ttneutralgenes.h.

References OFFSPRG, and setFstatWeirCockerham().

Referenced by setFstatWCRecorders().

{setFstatWeirCockerham(OFFSPRG);}

void TTNeutralGenesSH::setSibCoa	(	Individual *	I1,
		Individual *	I2
	)

Definition at line 348 of file stats_coa.cc.

{
  unsigned int nb_locus = this->_SHLinkedTrait->get_locus_num();
  //non sibs [3]
  if((I1->getMotherID() != I2->getMotherID()) && (I1->getFatherID() != I2->getFatherID())) {
    _sib_prop[3]++;
    _sib_coa[3] += Coancestry(I1->getTrait(_SHLinkedTraitIndex)->get_sequence(),
                              I2->getTrait(_SHLinkedTraitIndex)->get_sequence(),
                              nb_locus);
  }
  //maternal half sibs [2]
  else if((I1->getMotherID() == I2->getMotherID()) && (I1->getFatherID() != I2->getFatherID())) {
    _sib_prop[2]++;
    _sib_coa[2] += Coancestry(I1->getTrait(_SHLinkedTraitIndex)->get_sequence(),
                              I2->getTrait(_SHLinkedTraitIndex)->get_sequence(),
                              nb_locus);
  }
  //paternal half sibs [1]
  else if((I1->getMotherID() != I2->getMotherID()) && (I1->getFatherID() == I2->getFatherID())) {
    _sib_prop[1]++;
    _sib_coa[1] += Coancestry(I1->getTrait(_SHLinkedTraitIndex)->get_sequence(),
                              I2->getTrait(_SHLinkedTraitIndex)->get_sequence(),
                              nb_locus);
  }
  //full sibs [0]
  else if((I1->getMotherID() == I2->getMotherID()) && (I1->getFatherID() == I2->getFatherID())) {
    _sib_prop[0]++;
    _sib_coa[0] += Coancestry(I1->getTrait(_SHLinkedTraitIndex)->get_sequence(),
                              I2->getTrait(_SHLinkedTraitIndex)->get_sequence(),
                              nb_locus);
  }
}

void TTNeutralGenesSH::setSibStats ( )

Definition at line 277 of file stats_coa.cc.

Referenced by setStatRecorders().

{
  unsigned int i,j,k;
  unsigned int patchNbr = this->_pop->getPatchNbr(), Fsize, Msize;
  Patch* current_patch;
  Individual *I1, *I2;

  for(i = 0;i < 4; ++i) {
        _sib_prop[i] = 0.0;
    _sib_coa[i] = 0.0;
  }
  
  for( i = 0; i < patchNbr; ++i) {
    
    current_patch = _pop->getPatch(i);
    
        if ( (Fsize = current_patch->size(FEM, OFFSx)) != 0) {
      
          for(j = 0; j < Fsize -1; ++j) {
        
        I1 = current_patch->get(FEM, OFFSx, j);
        
                for(k = j+1; k < Fsize; ++k) {
          
          I2 = current_patch->get(FEM, OFFSx, k);
          
          setSibCoa(I1, I2);
                }
          }
    }
    
    if ( (Msize = current_patch->size(MAL, OFFSx)) != 0) {
      
      for(j = 0; j < Msize -1; ++j) {
        
        I1 = current_patch->get(MAL, OFFSx, j);
        
        for(k = j+1; k < Msize; ++k) {
          
          I2 = current_patch->get(MAL, OFFSx, k);
          
          setSibCoa(I1, I2);
        }
      }
    }
    
    //male-female
    for(j = 0; j < Msize; ++j) {
      
      I1 = current_patch->get(MAL, OFFSx, j);
      
      for(k = 0; k < Fsize; ++k) {
        
        I2 = current_patch->get(FEM, OFFSx, k);
        
        setSibCoa(I1, I2);
      }
    }
    
  }
  
  double tot = _sib_prop[0] + _sib_prop[1] + _sib_prop[2] + _sib_prop[3];
  
  for(i = 0 ; i < 4; ++i) {
        _sib_coa[i] = ( (_sib_prop[i] != 0) ? _sib_coa[i]/_sib_prop[i] : nanf("NULL"));
        _sib_prop[i] /= tot;
  }
}

bool TTNeutralGenesSH::setStatRecorders ( std::string & token ) [virtual]

Implements StatHandlerBase.

Definition at line 1103 of file ttneutralgenes.cc.

References StatHandler< TTNeutralGenesSH >::add(), ADULTS, ALL, getMeanAlpha(), getMeanTheta(), getSibCoaMeans(), getSibProportions(), getTheta_FF(), getTheta_FM(), getTheta_MM(), message(), OFFSPRG, setAdults_Theta(), setAdultsCoaBetween(), setAdultsCoaMatrix(), setAdultsCoaWithin(), setCoaMatrixRecorders(), setFstatRecorders(), setFstatWCRecorders(), setFstMatrixRecorders(), setNeiGeneticDistanceRecorders(), setOffsprgCoaBetween(), setOffsprgCoaMatrix(), setOffsprgCoaWithin(), and setSibStats().

{
#ifdef _DEBUG_
  message("-TTNeutralGenesSH::setStatRecorders ");
#endif
  if(token.compare("coa") == 0) {
        
        add("Wtn Patch Coancestry (offsprg)","off.theta",OFFSPRG,0,0,
        &TTNeutralGenesSH::getMeanTheta,0,0,&TTNeutralGenesSH::setOffsprgCoaMatrix);
        add("Btn Patch Coancestry (offsprg)","off.alpha",OFFSPRG,0,0,&TTNeutralGenesSH::getMeanAlpha,0,0,0);
    
        add("Wtn Patch Coancestry   (adult)","adlt.theta",ADULTS,0,0,
        &TTNeutralGenesSH::getMeanTheta,0,0,&TTNeutralGenesSH::setAdultsCoaMatrix);
        add("Btn Patch Coancestry   (adult)","adlt.alpha",ADULTS,0,0,&TTNeutralGenesSH::getMeanAlpha,0,0,0);
    
  } else if(token.compare("adlt.coa") == 0) {
    
        add("Wtn Patch Coancestry   (adult)","adlt.theta",ADULTS,0,0,
        &TTNeutralGenesSH::getMeanTheta,0,0,&TTNeutralGenesSH::setAdultsCoaMatrix);
        add("Btn Patch Coancestry   (adult)","adlt.alpha",ADULTS,0,0,&TTNeutralGenesSH::getMeanAlpha,0,0,0);
    
  } else if(token.compare("off.coa") == 0) {
    
    add("Wtn Patch Coancestry (offsprg)","off.theta",OFFSPRG,0,0,
        &TTNeutralGenesSH::getMeanTheta,0,0,&TTNeutralGenesSH::setOffsprgCoaMatrix);
        add("Btn Patch Coancestry (offsprg)","off.alpha",OFFSPRG,0,0,&TTNeutralGenesSH::getMeanAlpha,0,0,0);    
        
  } else if(token.compare("adlt.coa.persex") == 0) {
        
        add("Female Theta      (adult)","adlt.theta",ADULTS,0,0,
        &TTNeutralGenesSH::getMeanTheta, 0, 0, &TTNeutralGenesSH::setAdults_Theta);
        add("Female Theta      (adult)","adlt.thetaFF",ADULTS,0,0,&TTNeutralGenesSH::getTheta_FF,0,0,0);
        add("Male Theta        (adult)","adlt.thetaMM",ADULTS,0,0,&TTNeutralGenesSH::getTheta_MM,0,0,0);
        add("Female-Male Theta (adult)","adlt.thetaFM",ADULTS,0,0,&TTNeutralGenesSH::getTheta_FM,0,0,0);
    
  } else if(token.compare("adlt.coa.within") == 0) {
    
        add("Wtn Patch Coancestry   (adult)","adlt.theta",ADULTS,0,0,
        &TTNeutralGenesSH::getMeanTheta,0,0,&TTNeutralGenesSH::setAdultsCoaWithin);
    
  } else if(token.compare("off.coa.within") == 0) {
    
    add("Wtn Patch Coancestry (offsprg)","off.theta",OFFSPRG,0,0,
        &TTNeutralGenesSH::getMeanTheta,0,0,&TTNeutralGenesSH::setOffsprgCoaWithin);   
    
  } else if(token.compare("adlt.coa.between") == 0) {
    
        add("Btn Patch Coancestry   (adult)","adlt.alpha",ADULTS,0,0,
        &TTNeutralGenesSH::getMeanAlpha,0,0,&TTNeutralGenesSH::setAdultsCoaBetween);
    
  } else if(token.compare("off.coa.between") == 0) {
    
    add("Btn Patch Coancestry (offsprg)","off.alpha",OFFSPRG,0,0,
        &TTNeutralGenesSH::getMeanAlpha,0,0,&TTNeutralGenesSH::setOffsprgCoaBetween);   
    
  } else if(token.compare("coa.matrix") == 0) {
    
    setCoaMatrixRecorders(OFFSPRG, 3);
    setCoaMatrixRecorders(ADULTS, 3);
    
  } else if(token.compare("off.coa.matrix") == 0) {
        
    setCoaMatrixRecorders(OFFSPRG, 3);
    
  } else if(token.compare("adlt.coa.matrix") == 0) {
    
    setCoaMatrixRecorders(ADULTS, 3);
    
  } else if(token.compare("coa.matrix.within") == 0) {
    
    setCoaMatrixRecorders(OFFSPRG, 1);
    setCoaMatrixRecorders(ADULTS, 1);
    
  } else if(token.compare("off.coa.matrix.within") == 0) {
    
    setCoaMatrixRecorders(OFFSPRG, 1);
    
  } else if(token.compare("adlt.coa.matrix.within") == 0) {
    
    setCoaMatrixRecorders(ADULTS, 1);
    
  } else if(token.compare("sibcoa") == 0) {

    add("Proportion of full-sib offspring","prop.fsib",OFFSPRG,0,0,
        0,&TTNeutralGenesSH::getSibProportions,0,&TTNeutralGenesSH::setSibStats);
    add("Proportion of paternal half-sib ","prop.phsib",OFFSPRG,1,0,0,&TTNeutralGenesSH::getSibProportions,0,0);
    add("Proportion of maternal half-sib ","prop.mhsib",OFFSPRG,2,0,0, &TTNeutralGenesSH::getSibProportions,0,0);
    add("Proportion of non-sib offspring ","prop.nsib",OFFSPRG,3,0,0,&TTNeutralGenesSH::getSibProportions,0,0);
    add("Coancestry of full-sib offspring","coa.fsib",OFFSPRG,0,0,0,&TTNeutralGenesSH::getSibCoaMeans,0,0);
    add("Coancestry of paternal half-sib ","coa.phsib",OFFSPRG,1,0,0,&TTNeutralGenesSH::getSibCoaMeans,0,0);
    add("Coancestry of maternal half-sib ","coa.mhsib",OFFSPRG,2,0,0,&TTNeutralGenesSH::getSibCoaMeans,0,0);
    add("Coancestry of non-sib offspring ","coa.nsib",OFFSPRG,3,0,0,&TTNeutralGenesSH::getSibCoaMeans,0,0);
    
  } else if(token.compare("offsprgfstat") == 0 || token.compare("off.fstat") == 0) {
        
    setFstatRecorders(OFFSPRG);
    
  } else if(token.compare("adultfstat") == 0 || token.compare("adlt.fstat") == 0) {

    setFstatRecorders(ADULTS);
    
  } else if(token.compare("fstat") == 0) {
    
    setFstatRecorders(ALL);
    
  } else if(token.compare("fstWC") == 0 || token.compare("fstatWC") == 0) {
    
    setFstatWCRecorders(ALL);
     
  } else if(token.compare("off.fstWC") == 0 || token.compare("off.fstatWC") == 0) {
    
    setFstatWCRecorders(OFFSPRG);
    
  } else if(token.compare("adlt.fstWC") == 0 || token.compare("adlt.fstatWC") == 0) {
    
    setFstatWCRecorders(ADULTS);
    
  } else if(token.compare("weighted.fst") == 0) {
    
    setFstMatrixRecorders(OFFSPRG, 0);
    setFstMatrixRecorders(ADULTS, 0);
    
  } else if(token.compare("off.weighted.fst") == 0) {
    
    setFstMatrixRecorders(OFFSPRG, 0);
    
  } else if(token.compare("adlt.weighted.fst") == 0) {
    
    setFstMatrixRecorders(ADULTS, 0);
    
  } else if(token.compare("weighted.fst.matrix") == 0) {
    
    setFstMatrixRecorders(OFFSPRG, 3);
    setFstMatrixRecorders(ADULTS, 3);
    
  } else if(token.compare("off.weighted.fst.matrix") == 0) {
    
    setFstMatrixRecorders(OFFSPRG, 3);
    
  } else if(token.compare("adlt.weighted.fst.matrix") == 0) {
    
    setFstMatrixRecorders(ADULTS, 3);
    
  } else if(token.compare("weighted.fst.within") == 0) {
    
    setFstMatrixRecorders(OFFSPRG, 1);
    setFstMatrixRecorders(ADULTS, 1);
    
  } else if(token.compare("off.weighted.fst.within") == 0) {
    
    setFstMatrixRecorders(OFFSPRG, 1);
    
  } else if(token.compare("adlt.weighted.fst.within") == 0) {
    
    setFstMatrixRecorders(ADULTS, 1);
    
  } else if(token.compare("adlt.NeiDistance") == 0) {

    setNeiGeneticDistanceRecorders(ADULTS, true);
    
  } else if(token.compare("off.NeiDistance") == 0) {
    
    setNeiGeneticDistanceRecorders(OFFSPRG, true);
    
  } else if(token.compare("NeiDistance") == 0) {
    
    setNeiGeneticDistanceRecorders(ADULTS, true);
    setNeiGeneticDistanceRecorders(OFFSPRG, true);
    
  } else if(token.compare("mean.NeiDistance") == 0) {
   
    setNeiGeneticDistanceRecorders(ADULTS, false);
    setNeiGeneticDistanceRecorders(OFFSPRG, false);
    
  } else if(token.compare("adlt.mean.NeiDistance") == 0) {
    
    setNeiGeneticDistanceRecorders(ADULTS, false);
    
  } else if(token.compare("off.mean.NeiDistance") == 0) {
    
    setNeiGeneticDistanceRecorders(OFFSPRG, false);
        
  } else
        return false;
  
  return true;
}