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 ()

virtual bool	setStatRecorders (std::string &token)

void	setFreqRecorders (age_t AGE)

void	setFreqRecordersPerPatch (age_t AGE)

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)

void	setDxyRecorders (age_t AGE, bool patchwise)

Allele and genotype frequencies:
void	setAdultAlleleFreq ()

void	setOffspringAlleleFreq ()

void	setHeterozygosity (age_t AGE)

void	setAdultHeterozygosity ()

void	setOffspringHeterozygosity ()

double	getGlobalAlleleFreq (unsigned int loc, unsigned int all)

double	getHeterozygosity (unsigned int loc)

F-stats:
void	setAlleleTables (age_t AGE)

void	setHeteroTable (age_t AGE)

void	allocateTables (unsigned int loci, unsigned int all)

DataTable< double > *	getAlleleFreqTable ()
	Accessor to the table of allele frequencies, per patch. More...

DataTable< unsigned int > *	getAlleleCountTable ()

DataTable< double > *	getHeteroTable ()

TMatrix *	getGlobalFreqs ()
	Accessor to the table of allele frequencies in the whole population. More...

void	setFstMatrix (age_t AGE, unsigned char dim)
	Computes the weighted within and between patch Fst's as well as the overall Fst (Theta). More...

void	setAdultsFstMatrix ()

void	setAdultsFstWithin ()

void	setAdultsFstBetween ()

void	setOffsprgFstMatrix ()

void	setOffsprgFstWithin ()

void	setOffsprgFstBetween ()

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

double	getFst_ij (unsigned int i)
	Accessor to the Fst matrix as set by setFstMatrix(). More...

void	setFst_li (unsigned int N, unsigned int L, double **array)
	Computes the per-locus per-patch Fst values using Weir&Hill 2002 approach. More...

void	setFstat (age_t AGE)
	Computes the F-statistics following Nei & Chesser (1983). More...

void	setFstat_bitstring (age_t AGE)
	Streaming F-stat computation for diallelic bitstring traits. More...

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 ()

deque< double >	setHo2 (age_idx age_pos)
	New version of Nei & Chesser. More...

deque< double >	setHs2 (age_idx age_pos)

deque< double >	setHt2 (age_idx age_pos)

void	setFstatWeirCockerham (age_t AGE)
	Computes the Weir & Cockerham (1984) Fstat values (Theta, F, and f). More...

void	setFstatWeirCockerham_bitstring (age_t AGE)
	Streaming W&C Fstat for diallelic bitstring traits. More...

void	setFstatWeirCockerham_MS (age_t AGE)

void	setOffspringFstatWeirCockerham ()

void	setAdultsFstatWeirCockerham ()

double	getFstWC ()

double	getFisWC ()

double	getFitWC ()

void	setLociDivCounter (age_t AGE)
	Sets the allelic diversity counters. More...

double	getNbAllLocal ()

double	getNbAllGlobal ()

double	getFixLocLocal ()

double	getFixLocGlobal ()

Coancestries
double	Coancestry (const TTrait ind1, const TTrait ind2, unsigned int nb_locus)
	Gives the coancestry (probability of identity by state) of two gene sequences. More...

void	setCoaMatrix (age_idx age_pos, unsigned char dim)
	Computes the within and between patches coancestry coefficients. More...

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. More...

double	getMeanTheta ()

double	getMeanAlpha ()

double	getTheta_FF ()
	Gives the mean within females coancestry coefficient. More...

double	getTheta_MM ()
	Gives the mean within males coancestry coefficient. More...

double	getTheta_FM ()
	Gives the mean between males and females coancestry coefficient. More...

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 ()

Sequence divergence: Dxy (Nei & Li 1979; Nei 1987)
double	getDxyOffspringPerPatch (unsigned int patch1, unsigned patch2)

double	getDxyAdultPerPatch (unsigned int patch1, unsigned patch2)

double	getDxyPerPatch (age_idx age, unsigned int patch1, unsigned patch2)

double	getDxy (unsigned int age_class)

Public Member Functions inherited from TraitStatHandler< TProtoNeutralGenes, TTNeutralGenesSH >
	TraitStatHandler (TProtoNeutralGenes *trait_proto)

virtual	~TraitStatHandler ()

Public Member Functions inherited from StatHandler< SH >
	StatHandler ()

virtual	~StatHandler ()

virtual void	clear ()
	Empties the _recorders list, they are destroyed in StatHandlerBase::reset(). More...

virtual StatRecorder< SH > *	add (std::string Title, std::string Name, age_t AGE, unsigned int ARG1, unsigned int ARG2, double(SH::getStatNoArg)(void), double(SH::getStatOneArg)(unsigned int), double(SH::getStatTwoArg)(unsigned int, unsigned int), void(SH::setStat)(void))
	Adds a StatRecorder to the list, it is also added to the StatHandlerBase::_stats list. More...

Public Member Functions inherited from StatHandlerBase
	StatHandlerBase ()

virtual	~StatHandlerBase ()

virtual void	reset ()
	Empties the _stats list and calls clear() (defined in the derived class). More...

Metapop *	get_pop_ptr ()

void	set_service (StatServices *srv)

StatServices *	get_service ()

unsigned int	getOccurrence ()

unsigned int	getNumOccurrences ()

unsigned int	getCurrentOccurrence ()

unsigned int	getNbRecorders ()

std::list< StatRecBase * > &	getStats ()

virtual void	add (StatRecBase *rec)

virtual void	update ()
	This function is left empty as the StatServices calls StatRecorder::setVal directly. More...

Public Member Functions inherited from Handler
virtual	~Handler ()

Private Attributes
DataTable< unsigned int >	_alleleCountTable

DataTable< double >	_alleleFreqTable

DataTable< double >	_heteroTable

TMatrix	_globalAlleleFreq

unsigned int	_table_set_gen

unsigned int	_table_set_age

unsigned int	_table_set_repl

bool	_is_diallelic_bitstring

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. More...

double	_sib_coa [4]

double	_ho
	F-statistics. More...

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
	Weir & Hill (2002) F-stat estimates. More...

double	_fst_WC
	Weir & Cockerham (1984) F-stat estimates. More...

double	_fis_WC

double	_fit_WC

double *	_fst_WC_loc
	Per-locus F-stats (Weir&Cockerham). More...

double *	_fis_WC_loc

double *	_fit_WC_loc

double	_fst_W1

double	_fst_W2

TMatrix *	_fst_matrix
	Pairwise Fst matrix. More...

TMatrix *	_D

double	_meanD

Additional Inherited Members
Protected Types inherited from StatHandler< SH >
typedef std::list< StatRecorder< SH > * >::iterator	REC_IT

Protected Attributes inherited from TraitStatHandler< TProtoNeutralGenes, TTNeutralGenesSH >
TProtoNeutralGenes *	_SHLinkedTrait
	Pointer to a TraitProtoype object. More...

int	_SHLinkedTraitIndex
	Index of the trait in the Individual::Traits table. More...

Protected Attributes inherited from StatHandler< SH >
std::list< StatRecorder< SH > * >	_recorders
	The list of stat recorders. More...

Protected Attributes inherited from StatHandlerBase
Metapop *	_pop
	Link to the current population, set through the link to the StatService. More...

Detailed Description

The stat handler for neutral markers.

Constructor & Destructor Documentation

◆ TTNeutralGenesSH()

TTNeutralGenesSH::TTNeutralGenesSH ( TProtoNeutralGenes * TP )

inline

     : TraitStatHandler<TProtoNeutralGenes, TTNeutralGenesSH>(TP), _table_set_gen(0), _table_set_age(0),
     _table_set_repl(0), _is_diallelic_bitstring(false), _coa_matrix(0), _fst_WC_loc(0), _fis_WC_loc(0),
     _fit_WC_loc(0), _fst_matrix(0), _D(0)
         { }

◆ ~TTNeutralGenesSH()

virtual TTNeutralGenesSH::~TTNeutralGenesSH ( )

inlinevirtual

   {
     if(_coa_matrix != NULL) delete _coa_matrix;
     if(_fst_matrix != NULL) delete _fst_matrix;
     if(_fst_WC_loc) delete[]_fst_WC_loc;
     if(_fis_WC_loc) delete[]_fis_WC_loc;
     if(_fit_WC_loc) delete[]_fit_WC_loc;
     if(_D != NULL) delete _D;
   }

References _coa_matrix, _D, _fis_WC_loc, _fit_WC_loc, _fst_matrix, and _fst_WC_loc.

Member Function Documentation

◆ allocateTables()

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

 {
   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;
   
   //reset the time info, to force recalculation after allocate
   _table_set_age = NONE;
   _table_set_gen = 0;
   _table_set_repl = 0;
 }

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

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

Here is the caller graph for this function:

◆ Coancestry()

double TTNeutralGenesSH::Coancestry	(	const TTrait *	ind1,
		const TTrait *	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

 {
   unsigned int p = 0;
  
   for (unsigned int k = 0; k < nb_locus; ++k)
     p += (ind1->get_allele(k,0)==ind2->get_allele(k,0))
        + (ind1->get_allele(k,0)==ind2->get_allele(k,1))
        + (ind1->get_allele(k,1)==ind2->get_allele(k,0))
        + (ind1->get_allele(k,1)==ind2->get_allele(k,1));
  
   return (double)p/(4.0*nb_locus);
 }

References TTrait::get_allele().

◆ getAlleleCountTable()

DataTable< unsigned int >* TTNeutralGenesSH::getAlleleCountTable ( )

inline

416 {return &_alleleCountTable;}

References _alleleCountTable.

Referenced by TTNeutralGenesFH::write_varcompWC().

Here is the caller graph for this function:

◆ getAlleleFreqTable()

DataTable<double>* TTNeutralGenesSH::getAlleleFreqTable ( )

inline

Accessor to the table of allele frequencies, per patch.

414 {return &_alleleFreqTable;}

References _alleleFreqTable.

Referenced by TTNeutralGenesFH::write_varcompWC().

Here is the caller graph for this function:

◆ getCoa()

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.

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

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

◆ getDxy()

double TTNeutralGenesSH::getDxy ( unsigned int age_class )

 {
   double D = 0;
   unsigned int p = 0;
   
   age_t AGE = (static_cast<age_idx>(age_class) == OFFSx ? OFFSPRG : ADULTS);
   
   for (unsigned int p1 = 0; p1 < _pop->getPatchNbr(); ++p1) {
     
     if( !_pop->size( AGE , p1) ) continue;
       
     for (unsigned int p2 = p1 + 1; p2 < _pop->getPatchNbr(); ++p2) {
       
       if( !_pop->size( AGE , p2) ) continue;
  
       D += getDxyPerPatch(static_cast<age_idx>(age_class), p1, p2);
     
       p++;
     }
   }
   
   return (p != 0 ? D/p : nanf("NULL"));
 }

References StatHandlerBase::_pop, ADULTS, getDxyPerPatch(), Metapop::getPatchNbr(), OFFSPRG, OFFSx, and Metapop::size().

Referenced by setDxyRecorders().

Here is the caller graph for this function:

◆ getDxyAdultPerPatch()

double TTNeutralGenesSH::getDxyAdultPerPatch	(	unsigned int	patch1,
		unsigned	patch2
	)

inline

562 {return getDxyPerPatch(ADLTx, patch1, patch2);}

ADLTx

@ ADLTx

Definition: types.h:41

References ADLTx, and getDxyPerPatch().

Referenced by setDxyRecorders().

Here is the caller graph for this function:

◆ getDxyOffspringPerPatch()

double TTNeutralGenesSH::getDxyOffspringPerPatch	(	unsigned int	patch1,
		unsigned	patch2
	)

inline

561 {return getDxyPerPatch(OFFSx, patch1, patch2);}

References getDxyPerPatch(), and OFFSx.

Referenced by setDxyRecorders().

Here is the caller graph for this function:

◆ getDxyPerPatch()

double TTNeutralGenesSH::getDxyPerPatch	(	age_idx	age,
		unsigned int	patch1,
		unsigned	patch2
	)

 {
   double D = 0;
   Patch* patch_1 = _pop->getPatchPtr(patch1);
   Patch* patch_2 = _pop->getPatchPtr(patch2);
   unsigned int size_1 = patch_1->size(FEM, age);
   unsigned int size_2 = patch_2->size(FEM, age);
   unsigned int N = 0;//4*size_1 * size_2;// (2N)^2, tot num of haplotype comparisons
   
   
   unsigned int num_loc = _SHLinkedTrait->get_locus_num();
  
   TTrait *trait_1, *trait_2;
  
   //females:
   for (unsigned int i = 0; i < size_1; ++i) {
  
     trait_1 = patch_1->get(FEM, age, i)->getTrait(_SHLinkedTraitIndex);
  
     for (unsigned int j = 0; j < size_2; ++j) {
  
       trait_2 = patch_2->get(FEM, age, j)->getTrait(_SHLinkedTraitIndex);
  
       for (unsigned int l = 0; l < num_loc; ++l) {
         D += (trait_1->get_allele(l, 0) != trait_2->get_allele(l, 0));
         D += (trait_1->get_allele(l, 1) != trait_2->get_allele(l, 1));
         D += (trait_1->get_allele(l, 0) != trait_2->get_allele(l, 1));
         D += (trait_1->get_allele(l, 1) != trait_2->get_allele(l, 0));
       }
       N += 4;
     }
  
     for (unsigned int j = 0; j < patch_2->size(MAL, age); ++j) {
  
       trait_2 = patch_2->get(MAL, age, j)->getTrait(_SHLinkedTraitIndex);
  
       for (unsigned int l = 0; l < num_loc; ++l) {
         D += (trait_1->get_allele(l, 0) != trait_2->get_allele(l, 0));
         D += (trait_1->get_allele(l, 1) != trait_2->get_allele(l, 1));
         D += (trait_1->get_allele(l, 0) != trait_2->get_allele(l, 1));
         D += (trait_1->get_allele(l, 1) != trait_2->get_allele(l, 0));
       }
       N += 4;
  
     }
   }
  
   //same for the males:
   size_1 = patch_1->size(MAL, age);
   size_2 = patch_2->size(MAL, age);
  
 //  N += 4 * size_1 * size_2;
  
 //  if( N == 0) return 0; //needless to go further
  
   for (unsigned int i = 0; i < size_1; ++i) {
  
     trait_1 = patch_1->get(MAL, age, i)->getTrait(_SHLinkedTraitIndex);
  
     for (unsigned int j = 0; j < size_2; ++j) {
  
       trait_2 = patch_2->get(MAL, age, j)->getTrait(_SHLinkedTraitIndex);
  
       for (unsigned int l = 0; l < num_loc; ++l) {
         D += (trait_1->get_allele(l, 0) != trait_2->get_allele(l, 0));
         D += (trait_1->get_allele(l, 1) != trait_2->get_allele(l, 1));
         D += (trait_1->get_allele(l, 0) != trait_2->get_allele(l, 1));
         D += (trait_1->get_allele(l, 1) != trait_2->get_allele(l, 0));
       }
       N += 4;
     }
  
     for (unsigned int j = 0; j < patch_2->size(FEM, age); ++j) {
  
       trait_2 = patch_2->get(FEM, age, j)->getTrait(_SHLinkedTraitIndex);
  
       for (unsigned int l = 0; l < num_loc; ++l) {
         D += (trait_1->get_allele(l, 0) != trait_2->get_allele(l, 0));
         D += (trait_1->get_allele(l, 1) != trait_2->get_allele(l, 1));
         D += (trait_1->get_allele(l, 0) != trait_2->get_allele(l, 1));
         D += (trait_1->get_allele(l, 1) != trait_2->get_allele(l, 0));
       }
       N += 4;
  
     }
   }
   
   if( N == 0) return 0;
  
   return D/N;
 }

References StatHandlerBase::_pop, TraitStatHandler< TProtoNeutralGenes, TTNeutralGenesSH >::_SHLinkedTrait, TraitStatHandler< TProtoNeutralGenes, TTNeutralGenesSH >::_SHLinkedTraitIndex, FEM, Patch::get(), TTrait::get_allele(), TProtoNeutralGenes::get_locus_num(), Metapop::getPatchPtr(), Individual::getTrait(), MAL, and Patch::size().

Referenced by getDxy(), getDxyAdultPerPatch(), and getDxyOffspringPerPatch().

Here is the caller graph for this function:

◆ getFis()

double TTNeutralGenesSH::getFis ( )

inline

470 {return _fis;}

TTNeutralGenesSH::_fis

double _fis

Definition: ttneutralgenes.h:339

References _fis.

Referenced by setFstatRecorders().

Here is the caller graph for this function:

◆ getFisWC()

double TTNeutralGenesSH::getFisWC ( )

inline

488 {return _fis_WC;}

TTNeutralGenesSH::_fis_WC

double _fis_WC

Definition: ttneutralgenes.h:343

References _fis_WC.

Referenced by setFstatWCRecorders().

Here is the caller graph for this function:

◆ getFit()

double TTNeutralGenesSH::getFit ( )

inline

471 {return _fit;}

TTNeutralGenesSH::_fit

double _fit

Definition: ttneutralgenes.h:339

References _fit.

Referenced by setFstatRecorders().

Here is the caller graph for this function:

◆ getFitWC()

double TTNeutralGenesSH::getFitWC ( )

inline

489 {return _fit_WC;}

TTNeutralGenesSH::_fit_WC

double _fit_WC

Definition: ttneutralgenes.h:343

References _fit_WC.

Referenced by setFstatWCRecorders().

Here is the caller graph for this function:

◆ getFixLocGlobal()

double TTNeutralGenesSH::getFixLocGlobal ( )

inline

495 {return _fix_loc_global;}

TTNeutralGenesSH::_fix_loc_global

double _fix_loc_global

Definition: ttneutralgenes.h:339

References _fix_loc_global.

Referenced by setFstatRecorders().

Here is the caller graph for this function:

◆ getFixLocLocal()

double TTNeutralGenesSH::getFixLocLocal ( )

inline

494 {return _fix_loc_local;}

TTNeutralGenesSH::_fix_loc_local

double _fix_loc_local

Definition: ttneutralgenes.h:339

References _fix_loc_local.

Referenced by setFstatRecorders().

Here is the caller graph for this function:

◆ getFst()

double TTNeutralGenesSH::getFst ( )

inline

469 {return _fst;}

TTNeutralGenesSH::_fst

double _fst

Definition: ttneutralgenes.h:339

References _fst.

Referenced by setFstatRecorders().

Here is the caller graph for this function:

◆ getFst_ij()

double TTNeutralGenesSH::getFst_ij ( unsigned int i )

inline

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

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

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

Referenced by setFstMatrixRecorders().

Here is the caller graph for this function:

◆ getFstWC()

double TTNeutralGenesSH::getFstWC ( )

inline

487 {return _fst_WC;}

TTNeutralGenesSH::_fst_WC

double _fst_WC

Weir & Cockerham (1984) F-stat estimates.

Definition: ttneutralgenes.h:343

References _fst_WC.

Referenced by setFstatWCRecorders().

Here is the caller graph for this function:

◆ getGlobalAlleleFreq()

double TTNeutralGenesSH::getGlobalAlleleFreq	(	unsigned int	loc,
		unsigned int	all
	)

inline

                                                                   {
     return _globalAlleleFreq.get(loc, all);
   }

References _globalAlleleFreq, and TMatrix::get().

Referenced by setFreqRecorders().

Here is the caller graph for this function:

◆ getGlobalFreqs()

TMatrix* TTNeutralGenesSH::getGlobalFreqs ( )

inline

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

421 {return &_globalAlleleFreq;}

References _globalAlleleFreq.

Referenced by TTNeutralGenesFH::write_varcompWC().

Here is the caller graph for this function:

◆ getHeteroTable()

DataTable<double>* TTNeutralGenesSH::getHeteroTable ( )

inline

418 {return &_heteroTable;}

TTNeutralGenesSH::_heteroTable

DataTable< double > _heteroTable

Definition: ttneutralgenes.h:324

References _heteroTable.

Referenced by TTNeutralGenesFH::write_varcompWC().

Here is the caller graph for this function:

◆ getHeterozygosity()

double TTNeutralGenesSH::getHeterozygosity ( unsigned int loc )

inline

                                               {
     double het = 0;
     for(unsigned int i = 0; i < _heteroTable.getNumGroups(); ++i )
       het += _heteroTable.get(i, loc, 0);
     return het/_heteroTable.getNumGroups(); //mean per patch heterozygosity
   }

References _heteroTable, DataTable< T >::get(), and DataTable< T >::getNumGroups().

Referenced by setFreqRecorders().

Here is the caller graph for this function:

◆ getHo()

double TTNeutralGenesSH::getHo ( )

inline

466 {return _ho;}

TTNeutralGenesSH::_ho

double _ho

F-statistics.

Definition: ttneutralgenes.h:338

References _ho.

Referenced by setFstatRecorders().

Here is the caller graph for this function:

◆ getHs()

double TTNeutralGenesSH::getHs ( )

inline

467 {return _hs;}

TTNeutralGenesSH::_hs

double _hs

Definition: ttneutralgenes.h:338

References _hs.

◆ getHsnei()

double TTNeutralGenesSH::getHsnei ( )

inline

464 {return _hsnei;}

TTNeutralGenesSH::_hsnei

double _hsnei

Definition: ttneutralgenes.h:338

References _hsnei.

Referenced by setFstatRecorders().

Here is the caller graph for this function:

◆ getHt()

double TTNeutralGenesSH::getHt ( )

inline

468 {return _ht;}

TTNeutralGenesSH::_ht

double _ht

Definition: ttneutralgenes.h:338

References _ht.

◆ getHtnei()

double TTNeutralGenesSH::getHtnei ( )

inline

465 {return _htnei;}

TTNeutralGenesSH::_htnei

double _htnei

Definition: ttneutralgenes.h:338

References _htnei.

Referenced by setFstatRecorders().

Here is the caller graph for this function:

◆ getMeanAlpha()

double TTNeutralGenesSH::getMeanAlpha ( )

inline

533 {return _mean_alpha;}

TTNeutralGenesSH::_mean_alpha

double _mean_alpha

Definition: ttneutralgenes.h:330

References _mean_alpha.

◆ getMeanNeiGeneticDistance()

double TTNeutralGenesSH::getMeanNeiGeneticDistance ( )

inline

556 {return _meanD;}

TTNeutralGenesSH::_meanD

double _meanD

Definition: ttneutralgenes.h:353

References _meanD.

Referenced by setNeiGeneticDistanceRecorders().

Here is the caller graph for this function:

◆ getMeanTheta()

double TTNeutralGenesSH::getMeanTheta ( )

inline

532 {return _mean_theta;}

TTNeutralGenesSH::_mean_theta

double _mean_theta

Definition: ttneutralgenes.h:330

References _mean_theta.

◆ getNbAllGlobal()

double TTNeutralGenesSH::getNbAllGlobal ( )

inline

493 {return _nb_all_global;}

TTNeutralGenesSH::_nb_all_global

double _nb_all_global

Definition: ttneutralgenes.h:338

References _nb_all_global.

Referenced by setFstatRecorders().

Here is the caller graph for this function:

◆ getNbAllLocal()

double TTNeutralGenesSH::getNbAllLocal ( )

inline

492 {return _nb_all_local;}

TTNeutralGenesSH::_nb_all_local

double _nb_all_local

Definition: ttneutralgenes.h:338

References _nb_all_local.

Referenced by setFstatRecorders().

Here is the caller graph for this function:

◆ getNeiGeneticDistance()

double TTNeutralGenesSH::getNeiGeneticDistance ( unsigned int i )

inline

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

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

Referenced by setNeiGeneticDistanceRecorders().

Here is the caller graph for this function:

◆ getSibCoaMeans()

double TTNeutralGenesSH::getSibCoaMeans ( unsigned int i )

inline

543 {return _sib_coa[i];}

TTNeutralGenesSH::_sib_coa

double _sib_coa[4]

Definition: ttneutralgenes.h:335

References _sib_coa.

◆ getSibProportions()

double TTNeutralGenesSH::getSibProportions ( unsigned int i )

inline

542 {return _sib_prop[i];}

TTNeutralGenesSH::_sib_prop

double _sib_prop[4]

Kinship classes proportions.

Definition: ttneutralgenes.h:334

References _sib_prop.

◆ getTheta_FF()

double TTNeutralGenesSH::getTheta_FF ( )

inline

Gives the mean within females coancestry coefficient.

535 {return Theta_FF;}

TTNeutralGenesSH::Theta_FF

double Theta_FF

Definition: ttneutralgenes.h:329

References Theta_FF.

◆ getTheta_FM()

double TTNeutralGenesSH::getTheta_FM ( )

inline

Gives the mean between males and females coancestry coefficient.

539 {return Theta_FM;}

TTNeutralGenesSH::Theta_FM

double Theta_FM

Definition: ttneutralgenes.h:329

References Theta_FM.

◆ getTheta_MM()

double TTNeutralGenesSH::getTheta_MM ( )

inline

Gives the mean within males coancestry coefficient.

537 {return Theta_MM;}

TTNeutralGenesSH::Theta_MM

double Theta_MM

Definition: ttneutralgenes.h:329

References Theta_MM.

◆ getWeightedFst()

double TTNeutralGenesSH::getWeightedFst ( )

inline

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

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

441 {return _fst_WH;}

TTNeutralGenesSH::_fst_WH

double _fst_WH

Weir & Hill (2002) F-stat estimates.

Definition: ttneutralgenes.h:341

References _fst_WH.

Referenced by setFstMatrixRecorders().

Here is the caller graph for this function:

◆ init()

void TTNeutralGenesSH::init ( )

virtual

Reimplemented from StatHandlerBase.

 {
   StatHandler<TTNeutralGenesSH>::init();
  
   _is_diallelic_bitstring = (_SHLinkedTrait->get_allele_num() == 2);
  
   // For diallelic bitstring, defer DataTable allocation.
   // setFstat_bitstring() computes F-stats in streaming mode without DataTables.
   // Other stat methods (fstat2, fstWC, etc.) will trigger lazy allocation if needed.
   if (!_is_diallelic_bitstring)
     allocateTables(_SHLinkedTrait->get_locus_num(),_SHLinkedTrait->get_allele_num());
 }

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

◆ setAdltNeiGeneticDistance()

void TTNeutralGenesSH::setAdltNeiGeneticDistance ( )

inline

548 {setNeiGeneticDistance(ADULTS);}

TTNeutralGenesSH::setNeiGeneticDistance

void setNeiGeneticDistance(age_t AGE)

Definition: stats_fstat.cc:1458

References ADULTS, and setNeiGeneticDistance().

Referenced by setNeiGeneticDistanceRecorders().

Here is the caller graph for this function:

◆ setAdultAlleleFreq()

void TTNeutralGenesSH::setAdultAlleleFreq ( )

inline

389 {setAlleleTables(ADULTS);}

TTNeutralGenesSH::setAlleleTables

void setAlleleTables(age_t AGE)

Definition: stats_fstat.cc:75

References ADULTS, and setAlleleTables().

Referenced by setFreqRecorders().

Here is the caller graph for this function:

◆ setAdultHeterozygosity()

void TTNeutralGenesSH::setAdultHeterozygosity ( )

inline

392 {setHeterozygosity(ADULTS);}

TTNeutralGenesSH::setHeterozygosity

void setHeterozygosity(age_t AGE)

Definition: stats_fstat.cc:252

References ADULTS, and setHeterozygosity().

Referenced by setFreqRecorders().

Here is the caller graph for this function:

◆ setAdults_Theta()

void TTNeutralGenesSH::setAdults_Theta ( )

 {
   unsigned int Fsize,Msize,i,j,k, 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);
     
     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),
                                P1->get(FEM, ADLTx, k)->getTrait(_SHLinkedTraitIndex),
                                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),
                                P1->get(MAL, ADLTx, k)->getTrait(_SHLinkedTraitIndex),
                                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),
                                P1->get(MAL, ADLTx, k)->getTrait(_SHLinkedTraitIndex),
                                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;
 }

References ADLTx, FEM, Patch::get(), Individual::getTrait(), MAL, and Patch::size().

◆ setAdultsCoaBetween()

void TTNeutralGenesSH::setAdultsCoaBetween ( )

inline

519 {setCoaMatrix(ADLTx, 2);}

TTNeutralGenesSH::setCoaMatrix

void setCoaMatrix(age_idx age_pos, unsigned char dim)

Computes the within and between patches coancestry coefficients.

Definition: stats_coa.cc:57

References ADLTx, and setCoaMatrix().

◆ setAdultsCoaMatrix()

void TTNeutralGenesSH::setAdultsCoaMatrix ( )

inline

515 {setCoaMatrix(ADLTx, 3);}

References ADLTx, and setCoaMatrix().

◆ setAdultsCoaWithin()

void TTNeutralGenesSH::setAdultsCoaWithin ( )

inline

517 {setCoaMatrix(ADLTx, 1);}

References ADLTx, and setCoaMatrix().

◆ setAdultsFstat()

void TTNeutralGenesSH::setAdultsFstat ( )

inline

460 {setFstat(ADULTS);}

TTNeutralGenesSH::setFstat

void setFstat(age_t AGE)

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

Definition: stats_fstat.cc:281

References ADULTS, and setFstat().

Referenced by setFstatRecorders().

Here is the caller graph for this function:

◆ setAdultsFstatWeirCockerham()

void TTNeutralGenesSH::setAdultsFstatWeirCockerham ( )

inline

486 {setFstatWeirCockerham(ADULTS);}

TTNeutralGenesSH::setFstatWeirCockerham

void setFstatWeirCockerham(age_t AGE)

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

Definition: stats_fstat.cc:975

References ADULTS, and setFstatWeirCockerham().

Referenced by setFstatWCRecorders().

Here is the caller graph for this function:

◆ setAdultsFstBetween()

void TTNeutralGenesSH::setAdultsFstBetween ( )

inline

435 {setFstMatrix(ADULTS, 2);}

TTNeutralGenesSH::setFstMatrix

void setFstMatrix(age_t AGE, unsigned char dim)

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

Definition: stats_fstat.cc:869

References ADULTS, and setFstMatrix().

Referenced by setFstMatrixRecorders().

Here is the caller graph for this function:

◆ setAdultsFstMatrix()

void TTNeutralGenesSH::setAdultsFstMatrix ( )

inline

433 {setFstMatrix(ADULTS, 3);}

References ADULTS, and setFstMatrix().

Referenced by setFstMatrixRecorders().

Here is the caller graph for this function:

◆ setAdultsFstWithin()

void TTNeutralGenesSH::setAdultsFstWithin ( )

inline

434 {setFstMatrix(ADULTS, 1);}

References ADULTS, and setFstMatrix().

Referenced by setFstMatrixRecorders().

Here is the caller graph for this function:

◆ setAlleleTables()

void TTNeutralGenesSH::setAlleleTables ( age_t AGE )

 {  
   
   if(_table_set_age == AGE 
      && _table_set_gen == _pop->getCurrentGeneration()
      && _table_set_repl == _pop->getCurrentReplicate())
     return;
   
   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;
   TTrait* trait;
  
   //check if the population size has changed:
   if(_alleleCountTable.getNumGroups() != patchNbr)
     allocateTables(_SHLinkedTrait->get_locus_num(), _SHLinkedTrait->get_allele_num());
  
   _alleleCountTable.init(0);
  
   if (_is_diallelic_bitstring) {
  
     unsigned int* count1 = new unsigned int[nb_locus]();
  
     for(unsigned int k = 0; k < patchNbr; ++k) {
  
       crnt_patch = _pop->getPatch(k);
       memset(count1, 0, nb_locus * sizeof(unsigned int));
  
       for (int sx = 0; sx < 2; ++sx) {
  
         for(unsigned int i = 0, size = crnt_patch->size(sex_t(sx), age_pos); i < size; ++i) {
  
           const TTNeutralGenes_bitstring* bs = static_cast<const TTNeutralGenes_bitstring*>(
               crnt_patch->get(sex_t(sx), age_pos, i)->getTrait(_SHLinkedTraitIndex));
  
           for (int chr = 0; chr < 2; ++chr) {
             const bitstring& seq = bs->get_bit_sequence(chr);
             size_t nwords = seq.nb_words();
             for (size_t w = 0; w < nwords; ++w) {
               bitstring::_ul word = *seq.getword_atIdx(w);
               while (word) {
                 count1[w * BITS_PER_WORD + __builtin_ctzl(word)]++;
                 word &= word - 1;
               }
             }
           }
         }
       }
  
       for (unsigned int j = 0; j < nb_locus; ++j) {
         _alleleCountTable.set(k, j, 1, count1[j]);
         _alleleCountTable.set(k, j, 0,
           crnt_patch->size(age_pos) * 2 - count1[j]);
       }
     }
  
     delete[] count1;
  
   } else {
  
     //counting the copies of each allele present in each patch:
     for(unsigned int k = 0; k < patchNbr; ++k) {
  
       crnt_patch = _pop->getPatch(k);
  
       for (int sx = 0; sx < 2; ++sx) {
  
         for(unsigned int i = 0, size = crnt_patch->size(sex_t(sx), age_pos); i < size; ++i) {
  
           trait = crnt_patch->get(sex_t(sx), age_pos, i)->getTrait(_SHLinkedTraitIndex);
  
           for (unsigned int j = 0; j < nb_locus; ++j)
             _alleleCountTable.increment(k, j, trait->get_allele(j, 0));
  
           for (unsigned int j = 0; j < nb_locus; ++j)
             _alleleCountTable.increment(k, j, trait->get_allele(j, 1));
  
         }
       }
     }
   }
   
   //allelic frequencies:
   for (unsigned int i = 0; i < patchNbr; ++i) {
     
     patch_size = _pop->getPatch(i)->size(age_pos) * 2;
     
     if (patch_size) {
       
     for (unsigned int l = 0; l < nb_locus; ++l)
       for (unsigned int u = 0; u < nb_allele; ++u)
         _alleleFreqTable.set(i, l, u, (double)_alleleCountTable.get(i, l, u) / (double)patch_size);
     
     } else {
       for (unsigned int l = 0; l < nb_locus; ++l)
         for (unsigned int u = 0; u < nb_allele; ++u)
           _alleleFreqTable.set(i, l, u, 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();
 }

References _alleleCountTable, _alleleFreqTable, _globalAlleleFreq, _is_diallelic_bitstring, StatHandlerBase::_pop, TraitStatHandler< TProtoNeutralGenes, TTNeutralGenesSH >::_SHLinkedTrait, TraitStatHandler< TProtoNeutralGenes, TTNeutralGenesSH >::_SHLinkedTraitIndex, _table_set_age, _table_set_gen, _table_set_repl, ADLTx, ADULTS, allocateTables(), TMatrix::assign(), BITS_PER_WORD, TMatrix::divide(), Patch::get(), DataTable< T >::get(), TTrait::get_allele(), TProtoNeutralGenes::get_allele_num(), TTNeutralGenes_bitstring::get_bit_sequence(), TProtoNeutralGenes::get_locus_num(), Metapop::getCurrentGeneration(), Metapop::getCurrentReplicate(), DataTable< T >::getNumGroups(), Metapop::getPatch(), Metapop::getPatchNbr(), Individual::getTrait(), bitstring::getword_atIdx(), DataTable< T >::increment(), DataTable< T >::init(), bitstring::nb_words(), OFFSx, TMatrix::plus(), DataTable< T >::set(), Metapop::size(), and Patch::size().

Referenced by setAdultAlleleFreq(), setFst_li(), setFstat(), setFstatWeirCockerham(), setFstatWeirCockerham_MS(), setFstMatrix(), setNeiGeneticDistance(), setOffspringAlleleFreq(), and TTNeutralGenesFH::write_varcompWC().

Here is the caller graph for this function:

◆ setCoaMatrix()

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

 {
   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 (theta)
     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; ++j)
             for(k = j+1; k < Fsize; ++k)
               _coa_matrix->plus(i, i, Coancestry(P1->get(FEM, age_pos, j)->getTrait(_SHLinkedTraitIndex),
                                                  P1->get(FEM, age_pos, k)->getTrait(_SHLinkedTraitIndex),
                                                  nb_locus));
         //mal-mal coa
         for(j = 0; j < Msize; ++j)
           for(k = j+1; k < Msize; ++k)
             _coa_matrix->plus(i, i, Coancestry(P1->get(MAL, age_pos, j)->getTrait(_SHLinkedTraitIndex),
                                                P1->get(MAL, age_pos, k)->getTrait(_SHLinkedTraitIndex),
                                                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),
                                                P1->get(MAL, age_pos, k)->getTrait(_SHLinkedTraitIndex),
                                                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),
                                                  P2->get(FEM, age_pos, k)->getTrait(_SHLinkedTraitIndex),
                                                  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),
                                                  P2->get(MAL, age_pos, k)->getTrait(_SHLinkedTraitIndex),
                                                  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),
                                                  P2->get(MAL, age_pos, k)->getTrait(_SHLinkedTraitIndex),
                                                  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),
                                                  P2->get(FEM, age_pos, k)->getTrait(_SHLinkedTraitIndex),
                                                  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
   
 }

References FEM, Patch::get(), Individual::getTrait(), MAL, and Patch::size().

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

Here is the caller graph for this function:

◆ setCoaMatrixRecorders()

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

◆ setDxyRecorders()

void TTNeutralGenesSH::setDxyRecorders	(	age_t	AGE,
		bool	patchwise
	)

 {
   string prefix = (AGE == OFFSPRG ? "off." : "adlt.");
   age_idx age = (AGE == OFFSPRG ? OFFSx : ADLTx );
   string name = "Sequence divergence - Dxy";
   string sub_name = "Dxy";
  
  
   if (!patchwise) {
  
     add(name, prefix + sub_name, AGE, (unsigned int)age, 0, 0, &TTNeutralGenesSH::getDxy, 0, 0);
  
   } else {
  
     for (unsigned int p1 = 0; p1 < _pop->getPatchNbr(); ++p1) {
       for (unsigned int p2 = p1 + 1; p2 < _pop->getPatchNbr(); ++p2) {
         add(name, prefix + sub_name + ".p" + tstring::int2str(p1+1) + "p" + tstring::int2str(p2+1),
             AGE, p1, p2,
             0, 0, (AGE == OFFSPRG ? &TTNeutralGenesSH::getDxyOffspringPerPatch :
                 &TTNeutralGenesSH::getDxyAdultPerPatch),
                 0);
       }
     }
  
   }
 }

References StatHandlerBase::_pop, StatHandler< SH >::add(), ADLTx, getDxy(), getDxyAdultPerPatch(), getDxyOffspringPerPatch(), Metapop::getPatchNbr(), tstring::int2str(), OFFSPRG, and OFFSx.

Referenced by setStatRecorders().

Here is the caller graph for this function:

◆ setFreqRecorders()

void TTNeutralGenesSH::setFreqRecorders ( age_t AGE )

 {
   string prefix = (AGE == OFFSPRG ? "off." : "adlt.");
  
   void (TTNeutralGenesSH::* setter) () = (AGE == ADULTS ? &TTNeutralGenesSH::setAdultAlleleFreq :
       &TTNeutralGenesSH::setOffspringAlleleFreq );
  
   unsigned int nb_allele = _SHLinkedTrait->get_allele_num();
   unsigned int nb_locus = _SHLinkedTrait->get_locus_num();
  
   for (unsigned int l = 0; l < nb_locus; ++l) {
     for (unsigned int u = 0; u < nb_allele-1; ++u) {
       add("", prefix + "ntrl.l" + tstring::int2str(l+1) + ".a" + tstring::int2str(u+1), AGE, l, u,
           0, 0, &TTNeutralGenesSH::getGlobalAlleleFreq, setter);
     }
   }
  
   setter = (AGE == ADULTS ? &TTNeutralGenesSH::setAdultHeterozygosity :
       &TTNeutralGenesSH::setOffspringHeterozygosity );
  
   for (unsigned int l = 0; l < nb_locus; ++l) {
     add("", prefix + "ntrl.l" + tstring::int2str(l+1) + ".Het", AGE, l, 0, 0,
         &TTNeutralGenesSH::getHeterozygosity, 0, setter);
     //    add("", prefix + "ntrl.l" + tstring::int2str(l) + ".Hom", AGE, l, 0, 0,
     //        &TTNeutralGenesSH::getHomozygosity, 0, 0);
   }
  
 }

References TraitStatHandler< TProtoNeutralGenes, TTNeutralGenesSH >::_SHLinkedTrait, StatHandler< SH >::add(), ADULTS, TProtoNeutralGenes::get_allele_num(), TProtoNeutralGenes::get_locus_num(), getGlobalAlleleFreq(), getHeterozygosity(), tstring::int2str(), OFFSPRG, setAdultAlleleFreq(), setAdultHeterozygosity(), setOffspringAlleleFreq(), and setOffspringHeterozygosity().

Referenced by setStatRecorders().

Here is the caller graph for this function:

◆ setFreqRecordersPerPatch()

void TTNeutralGenesSH::setFreqRecordersPerPatch ( age_t AGE )

◆ setFst_li()

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.

 {
   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;
 }

References _alleleFreqTable, _globalAlleleFreq, StatHandlerBase::_pop, TraitStatHandler< TProtoNeutralGenes, TTNeutralGenesSH >::_SHLinkedTrait, ADULTS, DataTable< T >::get(), TMatrix::get(), TProtoNeutralGenes::get_allele_num(), TProtoNeutralGenes::get_locus_num(), Metapop::getPatchNbr(), setAlleleTables(), and Metapop::size().

Referenced by TTNeutralGenesFH::write_Fst_i().

Here is the caller graph for this function:

◆ setFstat()

void TTNeutralGenesSH::setFstat ( age_t AGE )

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

 {
   if (_is_diallelic_bitstring) {
     setFstat_bitstring(AGE);
     return;
   }
  
   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") );
 }

References _alleleCountTable, _fis, _fit, _fst, _ho, _hs, _hsnei, _ht, _htnei, _is_diallelic_bitstring, StatHandlerBase::_pop, TraitStatHandler< TProtoNeutralGenes, TTNeutralGenesSH >::_SHLinkedTrait, ADLTx, ADULTS, allocateTables(), TProtoNeutralGenes::get_allele_num(), TProtoNeutralGenes::get_locus_num(), DataTable< T >::getNumGroups(), Metapop::getPatchNbr(), OFFSx, setAlleleTables(), setFstat_bitstring(), setHo(), setHs(), setHt(), setLociDivCounter(), and Metapop::size().

Referenced by setAdultsFstat(), and setOffsprgFstat().

Here is the caller graph for this function:

◆ setFstat_bitstring()

void TTNeutralGenesSH::setFstat_bitstring ( age_t AGE )

Streaming F-stat computation for diallelic bitstring traits.

Processes loci in word-sized chunks (64 loci) to avoid allocating per-patch x per-locus DataTables.

 {
   unsigned int patchNbr = _pop->getPatchNbr();
   unsigned int nloc = _SHLinkedTrait->get_locus_num();
   age_idx age_pos = (AGE == ADULTS ? ADLTx : OFFSx);
  
   // harmonic mean of non-empty patch sizes
   double harmonic = 0, nbpatch = 0;
   unsigned int* patch_2N = new unsigned int[patchNbr];
  
   for (unsigned int i = 0; i < patchNbr; ++i) {
     unsigned int psize = _pop->size(AGE, i);
     patch_2N[i] = psize * 2;
     if (psize != 0) {
       nbpatch++;
       harmonic += 1.0 / psize;
     }
   }
   harmonic = (nbpatch > 0) ? nbpatch / harmonic : 0;
  
   unsigned int total_2N = _pop->size(AGE) * 2;
  
   // Ho: use existing bitstring-optimized method (no DataTables needed)
   _ho = setHo(age_pos);
  
   // Streaming computation of Hs, Ht, and loci diversity counters.
   // Process BITS_PER_WORD loci at a time to avoid per-patch x per-locus DataTables.
   size_t nwords = BITSET_WORDS(nloc);
  
   double hs_sum = 0, ht_sum = 0;
   unsigned int nb_all_local_sum = 0;
   unsigned int fix_loc_local_count = 0, fix_loc_global_count = 0;
   bool pop_has_allele1, pop_has_allele0;
   double nb_all_global_sum = 0;
  
   // per-patch allele-1 count for each bit in current word
   unsigned int* count1 = new unsigned int[patchNbr * BITS_PER_WORD];
  
   for (size_t w = 0; w < nwords; ++w) {
  
     unsigned int loci_in_word = (unsigned int)BITS_PER_WORD;
     if (w == nwords - 1 && nloc % BITS_PER_WORD != 0)
       loci_in_word = nloc % BITS_PER_WORD;
  
     memset(count1, 0, patchNbr * BITS_PER_WORD * sizeof(unsigned int));
  
     // accumulate allele "1" counts per patch for this word
     for (unsigned int k = 0; k < patchNbr; ++k) {
       Patch* patch = _pop->getPatch(k);
       unsigned int base = k * BITS_PER_WORD;
  
       for (int sx = 0; sx < 2; ++sx) {
         for (unsigned int ind = 0, sz = patch->size(sex_t(sx), age_pos);
              ind < sz; ++ind)
         {
           const TTNeutralGenes_bitstring* bs =
             static_cast<const TTNeutralGenes_bitstring*>(
               patch->get(sex_t(sx), age_pos, ind)->getTrait(_SHLinkedTraitIndex));
  
           for (int chr = 0; chr < 2; ++chr) {
             bitstring::_ul word = *bs->get_bit_sequence(chr).getword_atIdx(w);
             while (word) {
               count1[base + __builtin_ctzl(word)]++;
               word &= word - 1;
             }
           }
         }
       }
     }
  
     // compute stats for each locus in this word
     for (unsigned int bit = 0; bit < loci_in_word; ++bit) {
  
       unsigned int global_count1 = 0;
       pop_has_allele1 = false;
       pop_has_allele0 = false;
  
       for (unsigned int k = 0; k < patchNbr; ++k) {
         if (patch_2N[k] == 0) continue;
  
         unsigned int c1 = count1[k * BITS_PER_WORD + bit];
         unsigned int c0 = patch_2N[k] - c1;
         global_count1 += c1;
  
         // Hs: within-patch expected heterozygosity 1 - p^2 - q^2 = 2pq
         double p = (double)c1 / patch_2N[k];
         hs_sum += 2.0 * p * (1.0 - p);
  
         // allelic diversity per patch: count alleles present
         nb_all_local_sum += (c0 > 0) + (c1 > 0);
  
         // fixed loci per patch
         fix_loc_local_count += (c1 == 0 || c1 == patch_2N[k]);
  
         pop_has_allele1 |= (c1 > 0);
         pop_has_allele0 |= (c0 > 0);
       }
  
       // Ht: total expected heterozygosity
       if (total_2N > 0) {
         double p_global = (double)global_count1 / total_2N;
         ht_sum += 2.0 * p_global * (1.0 - p_global);
       }
  
       // global allelic diversity
       nb_all_global_sum += (pop_has_allele0 ? 1.0 : 0.0)
                          + (pop_has_allele1 ? 1.0 : 0.0);
  
       // fixed loci globally
       fix_loc_global_count += (global_count1 == 0 || global_count1 == total_2N);
     }
   }
  
   delete[] count1;
   delete[] patch_2N;
  
   _hs = (nbpatch > 0) ? hs_sum / (nloc * nbpatch) : 0;
   _ht = (nloc > 0) ? ht_sum / nloc : 0;
  
   _nb_all_local  = (nbpatch > 0) ? nb_all_local_sum / (double)(nloc * nbpatch)
                                   : nanf("NULL");
   _nb_all_global = nb_all_global_sum / nloc;
   _fix_loc_local = (nbpatch > 0) ? (double)fix_loc_local_count / nbpatch : 0;
   _fix_loc_global = fix_loc_global_count;
  
   // 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") );
 }

References _fis, _fit, _fix_loc_global, _fix_loc_local, _fst, _ho, _hs, _hsnei, _ht, _htnei, _nb_all_global, _nb_all_local, StatHandlerBase::_pop, TraitStatHandler< TProtoNeutralGenes, TTNeutralGenesSH >::_SHLinkedTrait, TraitStatHandler< TProtoNeutralGenes, TTNeutralGenesSH >::_SHLinkedTraitIndex, ADLTx, ADULTS, BITS_PER_WORD, BITSET_WORDS, Patch::get(), TTNeutralGenes_bitstring::get_bit_sequence(), TProtoNeutralGenes::get_locus_num(), Metapop::getPatch(), Metapop::getPatchNbr(), Individual::getTrait(), bitstring::getword_atIdx(), OFFSx, setHo(), Metapop::size(), and Patch::size().

Referenced by setFstat().

Here is the caller graph for this function:

◆ setFstatRecorders()

void TTNeutralGenesSH::setFstatRecorders ( age_t AGE )

 {
   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);
   }
 }

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

Referenced by setStatRecorders().

Here is the caller graph for this function:

◆ setFstatWCRecorders()

void TTNeutralGenesSH::setFstatWCRecorders ( age_t AGE )

 {
   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);
  
   }
 }

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

Referenced by setStatRecorders().

Here is the caller graph for this function:

◆ setFstatWeirCockerham()

void TTNeutralGenesSH::setFstatWeirCockerham ( age_t AGE )

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

 {
   if (_is_diallelic_bitstring) {
     setFstatWeirCockerham_bitstring(AGE);
     return;
   }
   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);
   setHeteroTable(AGE);
   
   //init
   double *pop_sizes = new double [patchNbr];
   double tot_size, inv_ntot;
   double sum_weights = 0;
   double nbar, nc, inv_nbar;
   unsigned int extantPs = 0;
   
   tot_size = _pop->size(AGE);
   
   for(unsigned int i = 0; i < patchNbr; i++) {
     pop_sizes[i] = _pop->size(AGE, i); 
     if(pop_sizes[i]) {
       extantPs++;
       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);
   
   delete [] pop_sizes;   
 }

References _alleleFreqTable, _fis_WC, _fit_WC, _fst_WC, _globalAlleleFreq, _heteroTable, _is_diallelic_bitstring, StatHandlerBase::_pop, TraitStatHandler< TProtoNeutralGenes, TTNeutralGenesSH >::_SHLinkedTrait, DataTable< T >::get(), TMatrix::get(), TProtoNeutralGenes::get_allele_num(), TProtoNeutralGenes::get_locus_num(), Metapop::getPatchNbr(), setAlleleTables(), setFstatWeirCockerham_bitstring(), setHeteroTable(), and Metapop::size().

Referenced by setAdultsFstatWeirCockerham(), and setOffspringFstatWeirCockerham().

Here is the caller graph for this function:

◆ setFstatWeirCockerham_bitstring()

void TTNeutralGenesSH::setFstatWeirCockerham_bitstring ( age_t AGE )

Streaming W&C Fstat for diallelic bitstring traits.

 {
   unsigned int patchNbr = _pop->getPatchNbr();
   unsigned int nloc = _SHLinkedTrait->get_locus_num();
   age_idx age_pos = (AGE == ADULTS ? ADLTx : OFFSx);
  
   // patch sizes and W&C constants
   double *pop_sizes = new double[patchNbr];
   double tot_size = _pop->size(AGE);
   double inv_ntot = 1.0 / tot_size;
   double sum_weights = 0;
   unsigned int extantPs = 0;
  
   for (unsigned int i = 0; i < patchNbr; i++) {
     pop_sizes[i] = _pop->size(AGE, i);
     if (pop_sizes[i]) {
       extantPs++;
       sum_weights += (pop_sizes[i] * pop_sizes[i] / tot_size);
     }
   }
  
   double nbar = tot_size / extantPs;
   double nc = (tot_size - sum_weights) / (extantPs - 1);
   double inv_nbar = 1.0 / (nbar - 1);
   double s2_denom = 1.0 / ((extantPs - 1) * nbar);
   double r = (double)(extantPs - 1) / extantPs;
   double hbar_factor = (2 * nbar - 1) / (4 * nbar);
  
   double a = 0, b = 0, c = 0;
  
   size_t nwords = BITSET_WORDS(nloc);
  
   // per-patch allele-1 count and heterozygote count for each bit in current word
   unsigned int* count1 = new unsigned int[patchNbr * BITS_PER_WORD];
   unsigned int* het    = new unsigned int[patchNbr * BITS_PER_WORD];
  
   for (size_t w = 0; w < nwords; ++w) {
  
     unsigned int loci_in_word = (unsigned int)BITS_PER_WORD;
     if (w == nwords - 1 && nloc % BITS_PER_WORD != 0)
       loci_in_word = nloc % BITS_PER_WORD;
  
     memset(count1, 0, patchNbr * BITS_PER_WORD * sizeof(unsigned int));
     memset(het,    0, patchNbr * BITS_PER_WORD * sizeof(unsigned int));
  
     for (unsigned int k = 0; k < patchNbr; ++k) {
       Patch* patch = _pop->getPatch(k);
       unsigned int base = k * BITS_PER_WORD;
  
       for (int sx = 0; sx < 2; ++sx) {
         for (unsigned int ind = 0, sz = patch->size(sex_t(sx), age_pos);
              ind < sz; ++ind)
         {
           const TTNeutralGenes_bitstring* bs =
             static_cast<const TTNeutralGenes_bitstring*>(
               patch->get(sex_t(sx), age_pos, ind)->getTrait(_SHLinkedTraitIndex));
  
           bitstring::_ul w0 = *bs->get_bit_sequence(0).getword_atIdx(w);
           bitstring::_ul w1 = *bs->get_bit_sequence(1).getword_atIdx(w);
  
           // allele 1 count from both chromosomes
           bitstring::_ul word = w0;
           while (word) { count1[base + __builtin_ctzl(word)]++; word &= word - 1; }
           word = w1;
           while (word) { count1[base + __builtin_ctzl(word)]++; word &= word - 1; }
  
           // heterozygote count from XOR
           word = w0 ^ w1;
           while (word) { het[base + __builtin_ctzl(word)]++; word &= word - 1; }
         }
       }
     }
  
     // W&C variance components for each locus in this word.
     // For 2 alleles the contributions from allele 0 and allele 1 are identical,
     // so computing for one allele gives the same Fst ratios (factor of 2 cancels).
     for (unsigned int bit = 0; bit < loci_in_word; ++bit) {
  
       unsigned int global_count1 = 0;
       double hbar = 0;
  
       for (unsigned int k = 0; k < patchNbr; ++k) {
         if (pop_sizes[k] == 0) continue;
         global_count1 += count1[k * BITS_PER_WORD + bit];
         hbar += het[k * BITS_PER_WORD + bit];
       }
  
       double pbar = (double)global_count1 / (2.0 * tot_size);
       hbar *= inv_ntot;
  
       double s2 = 0;
       for (unsigned int k = 0; k < patchNbr; ++k) {
         if (pop_sizes[k] == 0) continue;
         double p_i = (double)count1[k * BITS_PER_WORD + bit] / (2.0 * pop_sizes[k]);
         double var = p_i - pbar;
         s2 += var * var * pop_sizes[k];
       }
       s2 *= s2_denom;
  
       double x = pbar * (1 - pbar) - r * s2;
       a += s2 - inv_nbar * (x - 0.25 * hbar);
       b += x - hbar_factor * hbar;
       c += hbar;
     }
   }
  
   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);
  
   delete[] count1;
   delete[] het;
   delete[] pop_sizes;
 }

References _fis_WC, _fit_WC, _fst_WC, StatHandlerBase::_pop, TraitStatHandler< TProtoNeutralGenes, TTNeutralGenesSH >::_SHLinkedTrait, TraitStatHandler< TProtoNeutralGenes, TTNeutralGenesSH >::_SHLinkedTraitIndex, ADLTx, ADULTS, BITS_PER_WORD, BITSET_WORDS, Patch::get(), TTNeutralGenes_bitstring::get_bit_sequence(), TProtoNeutralGenes::get_locus_num(), Metapop::getPatch(), Metapop::getPatchNbr(), Individual::getTrait(), bitstring::getword_atIdx(), OFFSx, Metapop::size(), and Patch::size().

Referenced by setFstatWeirCockerham().

Here is the caller graph for this function:

◆ setFstatWeirCockerham_MS()

void TTNeutralGenesSH::setFstatWeirCockerham_MS ( age_t AGE )

Computes W&C F-stats using the Mean Squares approach, similar to the implementation in Hierfstat. This code gives the exact same results as the other version.

 { 
         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);
         setHeteroTable(AGE);
  
         //init
         double *pop_sizes = new double [patchNbr];
         double tot_size;
         double sum_weights = 0;
         double nc;
         unsigned int extantPs = 0;
  
         tot_size = _pop->size(AGE);
  
         for(unsigned int i = 0; i < patchNbr; i++) {
                 pop_sizes[i] = _pop->size(AGE, i);
                 if(pop_sizes[i]) {
                         extantPs++;
                         sum_weights += (pop_sizes[i] * pop_sizes[i] / tot_size);
                 }
         }
  
         nc = (tot_size - sum_weights)/(extantPs-1);
  
 //      unsigned int np = extantPs;
         unsigned int npl = extantPs; //all loci typed in all patches
  
         //p = _alleleFreqTable
         //pb = _globalAlleleFreq
  
         unsigned int *alploc = new unsigned int [nb_locus];
  
         unsigned int **alploc_table = new unsigned int* [nb_locus];
  
         for(unsigned int i = 0; i < nb_locus; ++i)
                 alploc_table[i] = new unsigned int[nb_allele];
  
         unsigned int tot_num_allele = 0;
  
         for(unsigned int l = 0; l < nb_locus; ++l){
  
                 alploc[l] = 0;
  
                 for(unsigned int cnt,  a = 0; a < nb_allele; ++a) {
  
                         cnt=0;
  
                         for(unsigned int i = 0; i < patchNbr; i++) {
  
                                  cnt += _alleleCountTable.get(i,l,a);
  
                         }
                         alploc_table[l][a] = (cnt != 0);
                         alploc[l] += (cnt != 0);
                 }
  
                 tot_num_allele += alploc[l];
         }
  
         //n, and nal are given by pop_sizes, same num ind typed at all loci in each patch
         //nc is the same for each locus
         //nt is given by tot_size, same tot num of ind typed for all loci
  
         //SSG: het/2 for each allele
         double *SSG = new double[tot_num_allele];
         double *SSP = new double[tot_num_allele];
         double *SSi = new double[tot_num_allele];
  
         unsigned int all_cntr = 0;
  
         double het, freq, var;
  
         for(unsigned int l = 0; l < nb_locus; ++l) {
  
                 for(unsigned int a = 0; a < nb_allele & all_cntr < tot_num_allele; ++a) {
  
                         if(alploc_table[l][a] == 0) continue; //do not consider alleles not present in the pop
  
                         SSG[all_cntr] = 0;
                         SSi[all_cntr] = 0;
                         SSP[all_cntr] = 0;
  
                         for(unsigned int p = 0; p < patchNbr; ++p){
  
                                 if(!_pop->size(AGE, p)) continue; //skip empty patches
  
                                 het = _heteroTable.get(p, l, a);
  
                                 freq = _alleleFreqTable.get(p, l, a);
  
                                 var = freq - _globalAlleleFreq.get(l, a); //(p_liu - pbar_u)^2
  
                                 var *= var;
  
                                 SSG[all_cntr] += het;
  
                                 SSi[all_cntr] += 2*pop_sizes[p]*freq*(1-freq) - het/2;
  
                                 SSP[all_cntr] += 2*pop_sizes[p]*var;
                         }
  
                         all_cntr++;
                 }
  
         }
  
  
         assert(all_cntr == tot_num_allele);
  
         double *MSG = new double[tot_num_allele];
         double *MSP = new double[tot_num_allele];
         double *MSI = new double[tot_num_allele];
         double *sigw = new double[tot_num_allele];
         double *siga = new double[tot_num_allele];
         double *sigb = new double[tot_num_allele];
  
 //      double *FST_pal = new double[tot_num_allele];
 //      double *FIS_pal = new double[tot_num_allele];
  
         double SIGA = 0, SIGB = 0, SIGW = 0;
  
         for(unsigned int i = 0; i < tot_num_allele; ++i){
  
                 MSG[i] = SSG[i] / (2 * tot_size);
                 sigw[i] = MSG[i]; //wasted!
  
                 MSP[i] = SSP[i] / (npl-1);
  
                 MSI[i] = SSi[i]/ (tot_size - npl);
  
                 sigb[i] = 0.5*(MSI[i] - MSG[i]);
  
                 siga[i] = (MSP[i] - MSI[i])/(2*nc);
  
 //              FST_pal[i] = siga[i]/(siga[i]+sigb[i]+sigw[i]);
 //              FIS_pal[i] = sigb[i]/(sigb[i]+sigw[i]);
  
                 SIGA += siga[i];
                 SIGB += sigb[i];
                 SIGW += sigw[i];
         }
  
         //per locus stats:
         if(_fst_WC_loc) delete [] _fst_WC_loc; _fst_WC_loc = new double [nb_locus];
         if(_fis_WC_loc) delete [] _fis_WC_loc; _fis_WC_loc = new double [nb_locus];
         if(_fit_WC_loc) delete [] _fit_WC_loc; _fit_WC_loc = new double [nb_locus];
  
         double lsiga, lsigb, lsigw;
  
 //      cout<<"  computing sigma per locus\n";
  
         for(unsigned int allcntr = 0, i = 0; i < nb_locus; ++i) {
  
                 lsiga = lsigb = lsigw = 0;
  
                 for(unsigned int l = 0; l < alploc[i]; ++l) {
  
                         lsiga += siga[allcntr];
                         lsigb += sigb[allcntr];
                         lsigw += sigw[allcntr];
  
                         allcntr++;
  
                 }
  
                 _fst_WC_loc[i] = lsiga /(lsiga + lsigb + lsigw);
                 _fis_WC_loc[i] = lsigb /(lsigb + lsigw);
                 _fit_WC_loc[i] = (lsiga +lsigb) /(lsiga + lsigb + lsigw);
  
         }
  
  
         // Total F-stats
         _fst_WC = SIGA / (SIGA + SIGB + SIGW);
         _fit_WC = (SIGA + SIGB) / (SIGA + SIGB + SIGW);
         _fis_WC = SIGB / (SIGB + SIGW);
  
         delete[]pop_sizes;
         delete[]alploc;
         for(unsigned int i = 0; i < nb_locus; ++i)
                 delete[]alploc_table[i];
         delete[]alploc_table;
         delete[]SSG;
         delete[]SSi;
         delete[]SSP;
         delete[]MSG;
         delete[]MSI;
         delete[]MSP;
         delete[]sigw;
         delete[]siga;
         delete[]sigb;
  
 }

References _alleleCountTable, _alleleFreqTable, _fis_WC, _fis_WC_loc, _fit_WC, _fit_WC_loc, _fst_WC, _fst_WC_loc, _globalAlleleFreq, _heteroTable, StatHandlerBase::_pop, TraitStatHandler< TProtoNeutralGenes, TTNeutralGenesSH >::_SHLinkedTrait, DataTable< T >::get(), TMatrix::get(), TProtoNeutralGenes::get_allele_num(), TProtoNeutralGenes::get_locus_num(), Metapop::getPatchNbr(), setAlleleTables(), setHeteroTable(), and Metapop::size().

◆ setFstMatrix()

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

 {
   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;
 }

References _alleleFreqTable, _fst_matrix, _fst_WH, _globalAlleleFreq, StatHandlerBase::_pop, TraitStatHandler< TProtoNeutralGenes, TTNeutralGenesSH >::_SHLinkedTrait, TMatrix::assign(), DataTable< T >::get(), TMatrix::get(), TProtoNeutralGenes::get_allele_num(), TProtoNeutralGenes::get_locus_num(), Metapop::getPatchNbr(), TMatrix::length(), TMatrix::reset(), TMatrix::set(), setAlleleTables(), and Metapop::size().

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

Here is the caller graph for this function:

◆ setFstMatrixRecorders()

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

 {
   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("");
       }
     }
   }
  
 }

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

Referenced by setStatRecorders().

Here is the caller graph for this function:

◆ setHeteroTable()

void TTNeutralGenesSH::setHeteroTable ( age_t AGE )

 {
   unsigned int patchNbr = this->_pop->getPatchNbr();
   unsigned int nb_locus = this->_SHLinkedTrait->get_locus_num();
   unsigned int nb_allele = this->_SHLinkedTrait->get_allele_num();
   Patch* patch;
   unsigned int isHetero, a0, a1;
   TTrait* trait;
   age_idx age_pos = (AGE == ADULTS ? ADLTx : OFFSx);
  
   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.update(patchNbr, nb_locus, sizes);
   _heteroTable.init(0);
  
  
   for (unsigned int i = 0; i < patchNbr; ++i) {
  
     patch = _pop->getPatch(i);
  
     if(!patch->size(age_pos)) continue;
  
     for (int sx = 0; sx < 2; ++sx) {
  
       for(unsigned int j = 0, size = patch->size(sex_t(sx), age_pos); j < size; ++j) {
  
         trait = patch->get(sex_t(sx), age_pos, j)->getTrait(_SHLinkedTraitIndex);
  
         for (unsigned int l = 0; l < nb_locus; ++l) {
           a0 = trait->get_allele(l, 0);
           a1 = trait->get_allele(l, 1);
           isHetero = (a0 != a1);
           _heteroTable.plus(i, l, a0, isHetero);
           _heteroTable.plus(i, l, a1, isHetero);
         }
       }
     }
   }  
   for(unsigned int i = 0; i < patchNbr; ++i) delete [] sizes[i];
   delete [] sizes;
  
 }

References _heteroTable, StatHandlerBase::_pop, TraitStatHandler< TProtoNeutralGenes, TTNeutralGenesSH >::_SHLinkedTrait, TraitStatHandler< TProtoNeutralGenes, TTNeutralGenesSH >::_SHLinkedTraitIndex, ADLTx, ADULTS, Patch::get(), TTrait::get_allele(), TProtoNeutralGenes::get_allele_num(), TProtoNeutralGenes::get_locus_num(), Metapop::getPatch(), Metapop::getPatchNbr(), DataTable< T >::init(), OFFSx, DataTable< T >::plus(), Patch::size(), and DataTable< T >::update().

Referenced by setFstatWeirCockerham(), setFstatWeirCockerham_MS(), setHeterozygosity(), and TTNeutralGenesFH::write_varcompWC().

Here is the caller graph for this function:

◆ setHeterozygosity()

void TTNeutralGenesSH::setHeterozygosity ( age_t AGE )

 {
   unsigned int patchNbr  = this->_pop->getPatchNbr();
   unsigned int nb_locus  = this->_SHLinkedTrait->get_locus_num();
   unsigned int nb_allele = this->_SHLinkedTrait->get_allele_num();
   Patch* patch;
   age_idx age_pos = (AGE == ADULTS ? ADLTx : OFFSx);
  
   setHeteroTable(AGE);
   
   for (unsigned int i = 0, size; i < patchNbr; ++i) {
     
     patch = _pop->getPatch(i);
     size = patch->size(age_pos);
     
     if(!size) continue;
  
     for (unsigned int l = 0; l < nb_locus; ++l) {
       for (unsigned int u = 0; u < nb_allele; ++u) {
         _heteroTable.divide(i, l , u, size);
       }
     }
   }
   
   
 }

References _heteroTable, StatHandlerBase::_pop, TraitStatHandler< TProtoNeutralGenes, TTNeutralGenesSH >::_SHLinkedTrait, ADLTx, ADULTS, DataTable< T >::divide(), TProtoNeutralGenes::get_allele_num(), TProtoNeutralGenes::get_locus_num(), Metapop::getPatch(), Metapop::getPatchNbr(), OFFSx, setHeteroTable(), and Patch::size().

Referenced by setAdultHeterozygosity(), and setOffspringHeterozygosity().

Here is the caller graph for this function:

◆ setHo()

double TTNeutralGenesSH::setHo ( age_idx age_pos )

 {
   unsigned int nloc = _SHLinkedTrait->get_locus_num(), nbpatch = _pop->getPatchNbr(), psize;
   double indloc = 0, hetero = 0;
   Patch* current_patch;
  
   if (_is_diallelic_bitstring) {
  
     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) {
         const TTNeutralGenes_bitstring* bs = static_cast<const TTNeutralGenes_bitstring*>(
             current_patch->get(FEM, age_pos, j)->getTrait(_SHLinkedTraitIndex));
         hetero += bs->get_bit_sequence(0).count_xor(bs->get_bit_sequence(1));
       }
  
       psize = current_patch->size(MAL, age_pos);
       indloc += psize;
       for(unsigned int j = 0; j < psize; ++j) {
         const TTNeutralGenes_bitstring* bs = static_cast<const TTNeutralGenes_bitstring*>(
             current_patch->get(MAL, age_pos, j)->getTrait(_SHLinkedTraitIndex));
         hetero += bs->get_bit_sequence(0).count_xor(bs->get_bit_sequence(1));
       }
     }
  
   } else {
  
     TTrait* trait;
  
     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) {
         trait = current_patch->get(FEM, age_pos, j)->getTrait(_SHLinkedTraitIndex);
         for (unsigned int k = 0; k < nloc; ++k)  hetero += (trait->get_allele(k, 0) != trait->get_allele(k, 1));
       }
     }
     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) {
         trait = current_patch->get(MAL, age_pos, j)->getTrait(_SHLinkedTraitIndex);
         for (unsigned int k = 0; k < nloc; ++k)  hetero += (trait->get_allele(k, 0) != trait->get_allele(k, 1));
       }
     }
   }
  
   indloc *= nloc;
  
   return (indloc != 0 ? hetero/indloc : 0.0);
 }

References _is_diallelic_bitstring, StatHandlerBase::_pop, TraitStatHandler< TProtoNeutralGenes, TTNeutralGenesSH >::_SHLinkedTrait, TraitStatHandler< TProtoNeutralGenes, TTNeutralGenesSH >::_SHLinkedTraitIndex, bitstring::count_xor(), FEM, Patch::get(), TTrait::get_allele(), TTNeutralGenes_bitstring::get_bit_sequence(), TProtoNeutralGenes::get_locus_num(), Metapop::getPatch(), Metapop::getPatchNbr(), Individual::getTrait(), MAL, and Patch::size().

Referenced by setFstat(), and setFstat_bitstring().

Here is the caller graph for this function:

◆ setHo2()

deque< double > TTNeutralGenesSH::setHo2 ( age_idx age_pos )

New version of Nei & Chesser.

 {
   unsigned int nloc = _SHLinkedTrait->get_locus_num(), nbpatch = _pop->getPatchNbr();
   unsigned int psize;
   deque<double> hetero(nloc,0);
   Patch* current_patch;
  
   if (_is_diallelic_bitstring) {
  
     for (unsigned int i = 0; i < nbpatch; ++i) {
       current_patch = _pop->getPatch(i);
  
       for (int sx = 0; sx < 2; ++sx) {
         sex_t sex = (sx == 0) ? FEM : MAL;
         for(unsigned int j = 0, size = current_patch->size(sex, age_pos); j < size; ++j) {
           const TTNeutralGenes_bitstring* bs = static_cast<const TTNeutralGenes_bitstring*>(
               current_patch->get(sex, age_pos, j)->getTrait(_SHLinkedTraitIndex));
           const bitstring& s0 = bs->get_bit_sequence(0);
           const bitstring& s1 = bs->get_bit_sequence(1);
           size_t nwords = s0.nb_words();
           for (size_t w = 0; w < nwords; ++w) {
             bitstring::_ul xword = *s0.getword_atIdx(w) ^ *s1.getword_atIdx(w);
             while (xword) {
               hetero[w * BITS_PER_WORD + __builtin_ctzl(xword)]++;
               xword &= xword - 1;
             }
           }
         }
       }
     }
  
   } else {
  
     TTrait* trait;
  
     for (unsigned int i = 0; i < nbpatch; ++i) {
       current_patch = _pop->getPatch(i);
  
       for(unsigned int j = 0, fsize = current_patch->size(FEM, age_pos); j < fsize; ++j) {
         trait = current_patch->get(FEM, age_pos, j)->getTrait(_SHLinkedTraitIndex);
         for (unsigned int k = 0; k < nloc; ++k)
           hetero[k] += (trait->get_allele(k, 0) != trait->get_allele(k, 1));
       }
  
       for(unsigned int j = 0, msize = current_patch->size(MAL, age_pos); j < msize; ++j) {
         trait = current_patch->get(MAL, age_pos, j)->getTrait(_SHLinkedTraitIndex);
         for (unsigned int k = 0; k < nloc; ++k)
           hetero[k] += (trait->get_allele(k, 0) != trait->get_allele(k, 1));
       }
     }
   }
  
   psize = _pop->size(age_pos);
  
   if(psize != 0)
     for (unsigned int k = 0; k < nloc; ++k)
       hetero[k] /= psize;
  
   return hetero;
 }

References _is_diallelic_bitstring, StatHandlerBase::_pop, TraitStatHandler< TProtoNeutralGenes, TTNeutralGenesSH >::_SHLinkedTrait, TraitStatHandler< TProtoNeutralGenes, TTNeutralGenesSH >::_SHLinkedTraitIndex, BITS_PER_WORD, FEM, Patch::get(), TTrait::get_allele(), TTNeutralGenes_bitstring::get_bit_sequence(), TProtoNeutralGenes::get_locus_num(), Metapop::getPatch(), Metapop::getPatchNbr(), Individual::getTrait(), bitstring::getword_atIdx(), MAL, bitstring::nb_words(), Metapop::size(), and Patch::size().

Referenced by TTNeutralGenesFH::write_varcompWC().

Here is the caller graph for this function:

◆ setHs()

double TTNeutralGenesSH::setHs ( age_idx age_pos )

 {
   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);
 }

References _alleleFreqTable, StatHandlerBase::_pop, TraitStatHandler< TProtoNeutralGenes, TTNeutralGenesSH >::_SHLinkedTrait, Patch::assign(), DataTable< T >::get(), TProtoNeutralGenes::get_allele_num(), TProtoNeutralGenes::get_locus_num(), Metapop::getPatch(), Metapop::getPatchNbr(), and Patch::size().

Referenced by setFstat().

Here is the caller graph for this function:

◆ setHs2()

deque< double > TTNeutralGenesSH::setHs2 ( age_idx age_pos )

 {
   unsigned int i, k, x, nb_patch = 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>hsloc(nbLoc,0);
   double freq;
   Patch* current_patch;
   
   for (i = 0; i < patchNbr; ++i) {
     if(_pop->size(age_pos, i) != 0)
       nb_patch++;
   }
   
   //compute the expected heterozygosity for each patch and return average:
   for (k = 0; k < nbLoc; ++k) {    
     
     for (i = 0; i < patchNbr; ++i) {
       
       current_patch = _pop->getPatch(i);
       
       genehet[k] = 1;
       
       if(current_patch->size(age_pos) != 0) {
         
         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
         }
         
       }//end_if size!=0
       
       hsloc[k] += genehet[k];
       
     }//end patch loop
     
     if(nb_patch !=0) hsloc[k] /= nb_patch;
     
   }//end loci loop
   
   return  hsloc;
 }

References _alleleFreqTable, StatHandlerBase::_pop, TraitStatHandler< TProtoNeutralGenes, TTNeutralGenesSH >::_SHLinkedTrait, DataTable< T >::get(), TProtoNeutralGenes::get_allele_num(), TProtoNeutralGenes::get_locus_num(), Metapop::getPatch(), Metapop::getPatchNbr(), Metapop::size(), and Patch::size().

◆ setHt()

double TTNeutralGenesSH::setHt ( age_idx age_pos )

 {
   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);
 }

References _globalAlleleFreq, TraitStatHandler< TProtoNeutralGenes, TTNeutralGenesSH >::_SHLinkedTrait, TMatrix::get(), TProtoNeutralGenes::get_allele_num(), and TProtoNeutralGenes::get_locus_num().

Referenced by setFstat().

Here is the caller graph for this function:

◆ setHt2()

deque< double > TTNeutralGenesSH::setHt2 ( age_idx age_pos )

 {
   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
     }
   }
   
   return genehet;
 }

References _globalAlleleFreq, TraitStatHandler< TProtoNeutralGenes, TTNeutralGenesSH >::_SHLinkedTrait, TMatrix::get(), TProtoNeutralGenes::get_allele_num(), and TProtoNeutralGenes::get_locus_num().

◆ setLociDivCounter()

void TTNeutralGenesSH::setLociDivCounter ( age_t AGE )

Sets the allelic diversity counters.

 {
   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 );
  
 }

References _alleleCountTable, _alleleFreqTable, _fix_loc_global, _fix_loc_local, _globalAlleleFreq, _nb_all_global, _nb_all_local, StatHandlerBase::_pop, TraitStatHandler< TProtoNeutralGenes, TTNeutralGenesSH >::_SHLinkedTrait, DataTable< T >::get(), TMatrix::get(), TProtoNeutralGenes::get_allele_num(), TProtoNeutralGenes::get_locus_num(), Metapop::getPatchNbr(), and Metapop::size().

Referenced by setFstat().

Here is the caller graph for this function:

◆ setNeiGeneticDistance()

void TTNeutralGenesSH::setNeiGeneticDistance ( age_t AGE )

 {
   //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, equ. 6 in Nei 1978 (Genetics)
     }
   }
   
   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;
 }

References _alleleFreqTable, _D, _meanD, StatHandlerBase::_pop, TraitStatHandler< TProtoNeutralGenes, TTNeutralGenesSH >::_SHLinkedTrait, TMatrix::assign(), DataTable< T >::get(), TProtoNeutralGenes::get_allele_num(), TProtoNeutralGenes::get_locus_num(), Metapop::getPatchNbr(), TMatrix::length(), TMatrix::reset(), TMatrix::set(), setAlleleTables(), and Metapop::size().

Referenced by setAdltNeiGeneticDistance(), and setOffsprgNeiGeneticDistance().

Here is the caller graph for this function:

◆ setNeiGeneticDistanceRecorders()

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

 {
   string name, sub_name;
   unsigned int nbpatch =  _pop->getPatchNbr();
   string 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,  AGE, 0, 0,
       &TTNeutralGenesSH::getMeanNeiGeneticDistance, 0, 0, setter);
  
   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" + tstring::int2str(i+1) + " and pop" + tstring::int2str(j+1);
         sub_name = prefix + "D.p" + tstring::int2str(i+1) + "p" + tstring::int2str(j+1);
         add(name, sub_name,  AGE, i*scale + j, 0, 0, &TTNeutralGenesSH::getNeiGeneticDistance, 0, 0);
       }
     }
   }
 }

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

Referenced by setStatRecorders().

Here is the caller graph for this function:

◆ setOffsprgCoaBetween()

void TTNeutralGenesSH::setOffsprgCoaBetween ( )

inline

520 {setCoaMatrix(OFFSx, 2);}

References OFFSx, and setCoaMatrix().

◆ setOffsprgCoaMatrix()

void TTNeutralGenesSH::setOffsprgCoaMatrix ( )

inline

516 {setCoaMatrix(OFFSx, 3);}

References OFFSx, and setCoaMatrix().

◆ setOffsprgCoaWithin()

void TTNeutralGenesSH::setOffsprgCoaWithin ( )

inline

518 {setCoaMatrix(OFFSx, 1);}

References OFFSx, and setCoaMatrix().

◆ setOffsprgFstat()

void TTNeutralGenesSH::setOffsprgFstat ( )

inline

459 {setFstat(OFFSPRG);}

References OFFSPRG, and setFstat().

Referenced by setFstatRecorders().

Here is the caller graph for this function:

◆ setOffsprgFstBetween()

void TTNeutralGenesSH::setOffsprgFstBetween ( )

inline

438 {setFstMatrix(OFFSPRG, 2);}

References OFFSPRG, and setFstMatrix().

Referenced by setFstMatrixRecorders().

Here is the caller graph for this function:

◆ setOffsprgFstMatrix()

void TTNeutralGenesSH::setOffsprgFstMatrix ( )

inline

436 {setFstMatrix(OFFSPRG, 3);}

References OFFSPRG, and setFstMatrix().

Referenced by setFstMatrixRecorders().

Here is the caller graph for this function:

◆ setOffsprgFstWithin()

void TTNeutralGenesSH::setOffsprgFstWithin ( )

inline

437 {setFstMatrix(OFFSPRG, 1);}

References OFFSPRG, and setFstMatrix().

Referenced by setFstMatrixRecorders().

Here is the caller graph for this function:

◆ setOffsprgNeiGeneticDistance()

void TTNeutralGenesSH::setOffsprgNeiGeneticDistance ( )

inline

549 {setNeiGeneticDistance(OFFSPRG);}

References OFFSPRG, and setNeiGeneticDistance().

Referenced by setNeiGeneticDistanceRecorders().

Here is the caller graph for this function:

◆ setOffspringAlleleFreq()

void TTNeutralGenesSH::setOffspringAlleleFreq ( )

inline

390 {setAlleleTables(OFFSPRG);}\

References OFFSPRG, and setAlleleTables().

Referenced by setFreqRecorders().

Here is the caller graph for this function:

◆ setOffspringFstatWeirCockerham()

void TTNeutralGenesSH::setOffspringFstatWeirCockerham ( )

inline

485 {setFstatWeirCockerham(OFFSPRG);}

References OFFSPRG, and setFstatWeirCockerham().

Referenced by setFstatWCRecorders().

Here is the caller graph for this function:

◆ setOffspringHeterozygosity()

void TTNeutralGenesSH::setOffspringHeterozygosity ( )

inline

393 {setHeterozygosity(OFFSPRG);}

References OFFSPRG, and setHeterozygosity().

Referenced by setFreqRecorders().

Here is the caller graph for this function:

◆ setSibCoa()

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

 {
   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),
                               I2->getTrait(_SHLinkedTraitIndex),
                               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),
                               I2->getTrait(_SHLinkedTraitIndex),
                               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),
                               I2->getTrait(_SHLinkedTraitIndex),
                               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),
                               I2->getTrait(_SHLinkedTraitIndex),
                               nb_locus);
   }
 }  

References Individual::getFatherID(), Individual::getMotherID(), and Individual::getTrait().

◆ setSibStats()

void TTNeutralGenesSH::setSibStats ( )

 {
   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;
   }
 }

References FEM, Patch::get(), MAL, OFFSx, and Patch::size().

◆ setStatRecorders()

bool TTNeutralGenesSH::setStatRecorders ( std::string & token )

virtual

Implements StatHandlerBase.

 {
 #ifdef _DEBUG_
   message("-TTNeutralGenesSH::setStatRecorders ");
 #endif
 //  if(token == "coa") {
 //
 //    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 == "adlt.coa") {
 //
 //    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 == "off.coa") {
 //
 //    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 == "adlt.coa.persex") {
 //
 //    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 == "adlt.coa.within") {
 //
 //    add("Wtn Patch Coancestry   (adult)","adlt.theta",ADULTS,0,0,
 //        &TTNeutralGenesSH::getMeanTheta,0,0,&TTNeutralGenesSH::setAdultsCoaWithin);
 //
 //  } else if(token == "off.coa.within") {
 //
 //    add("Wtn Patch Coancestry (offsprg)","off.theta",OFFSPRG,0,0,
 //        &TTNeutralGenesSH::getMeanTheta,0,0,&TTNeutralGenesSH::setOffsprgCoaWithin);
 //
 //  } else if(token == "adlt.coa.between") {
 //
 //    add("Btn Patch Coancestry   (adult)","adlt.alpha",ADULTS,0,0,
 //        &TTNeutralGenesSH::getMeanAlpha,0,0,&TTNeutralGenesSH::setAdultsCoaBetween);
 //
 //  } else if(token == "off.coa.between") {
 //
 //    add("Btn Patch Coancestry (offsprg)","off.alpha",OFFSPRG,0,0,
 //        &TTNeutralGenesSH::getMeanAlpha,0,0,&TTNeutralGenesSH::setOffsprgCoaBetween);
 //
 //  } else if(token == "coa.matrix") {
 //
 //    setCoaMatrixRecorders(OFFSPRG, 3);
 //    setCoaMatrixRecorders(ADULTS, 3);
 //
 //  } else if(token == "off.coa.matrix") {
 //
 //    setCoaMatrixRecorders(OFFSPRG, 3);
 //
 //  } else if(token == "adlt.coa.matrix") {
 //
 //    setCoaMatrixRecorders(ADULTS, 3);
 //
 //  } else if(token == "coa.matrix.within") {
 //
 //    setCoaMatrixRecorders(OFFSPRG, 1);
 //    setCoaMatrixRecorders(ADULTS, 1);
 //
 //  } else if(token == "off.coa.matrix.within") {
 //
 //    setCoaMatrixRecorders(OFFSPRG, 1);
 //
 //  } else if(token == "adlt.coa.matrix.within") {
 //
 //    setCoaMatrixRecorders(ADULTS, 1);
 //
 //  } else if(token == "sibcoa") {
 //
 //    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 == "ntrl.freq") {
  
     setFreqRecorders(OFFSPRG);
     setFreqRecorders(ADULTS);
  
   } else if (token == "off.ntrl.freq") {
  
     setFreqRecorders(OFFSPRG);
  
   } else if (token == "adlt.ntrl.freq") {
  
     setFreqRecorders(ADULTS);
  
     //  } else if (token == "ntrl.freq.patch") {
     //
     //    setFreqRecordersPerPatch(ALL);
     //
     //  } else if (token == "off.ntrl.freq.patch") {
     //
     //    setFreqRecordersPerPatch(OFFSPRG);
     //
     //  } else if (token == "adlt.ntrl.freq.patch") {
     //
     //    setFreqRecordersPerPatch(ADULTS);
  
   } else if(token == "off.fstat") {
  
     setFstatRecorders(OFFSPRG);
  
   } else if(token == "adlt.fstat") {
  
     setFstatRecorders(ADULTS);
  
   } else if(token == "fstat") {
  
     setFstatRecorders(ALL);
  
 //  } else if(token == "off.fstat2") {
 //
 //    setFstat2Recorders(OFFSPRG);
 //
 //  } else if(token == "adlt.fstat2") {
 //
 //    setFstat2Recorders(ADULTS);
 //
 //  } else if(token == "fstat2") {
 //
 //    setFstat2Recorders(ALL);
  
   } else if(token == "fstWC" || token == "fstatWC") {
  
     setFstatWCRecorders(ALL);
  
   } else if(token == "off.fstWC" || token == "off.fstatWC") {
  
     setFstatWCRecorders(OFFSPRG);
  
   } else if(token == "adlt.fstWC" || token == "adlt.fstatWC") {
  
     setFstatWCRecorders(ADULTS);
  
   } else if(token == "weighted.fst") {
  
     setFstMatrixRecorders(OFFSPRG, 0);
     setFstMatrixRecorders(ADULTS, 0);
  
   } else if(token == "off.weighted.fst") {
  
     setFstMatrixRecorders(OFFSPRG, 0);
  
   } else if(token == "adlt.weighted.fst") {
  
     setFstMatrixRecorders(ADULTS, 0);
  
   } else if(token == "weighted.fst.matrix") {
  
     setFstMatrixRecorders(OFFSPRG, 3);
     setFstMatrixRecorders(ADULTS, 3);
  
   } else if(token == "off.weighted.fst.matrix") {
  
     setFstMatrixRecorders(OFFSPRG, 3);
  
   } else if(token == "adlt.weighted.fst.matrix") {
  
     setFstMatrixRecorders(ADULTS, 3);
  
   } else if(token == "weighted.fst.within") {
  
     setFstMatrixRecorders(OFFSPRG, 1);
     setFstMatrixRecorders(ADULTS, 1);
  
   } else if(token == "off.weighted.fst.within") {
  
     setFstMatrixRecorders(OFFSPRG, 1);
  
   } else if(token == "adlt.weighted.fst.within") {
  
     setFstMatrixRecorders(ADULTS, 1);
  
   } else if(token == "adlt.NeiDistance") {
  
     setNeiGeneticDistanceRecorders(ADULTS, true);
  
   } else if(token == "off.NeiDistance") {
  
     setNeiGeneticDistanceRecorders(OFFSPRG, true);
  
   } else if(token == "NeiDistance") {
  
     setNeiGeneticDistanceRecorders(ADULTS, true);
     setNeiGeneticDistanceRecorders(OFFSPRG, true);
  
   } else if(token == "mean.NeiDistance") {
  
     setNeiGeneticDistanceRecorders(ADULTS, false);
     setNeiGeneticDistanceRecorders(OFFSPRG, false);
  
   } else if(token == "adlt.mean.NeiDistance") {
  
     setNeiGeneticDistanceRecorders(ADULTS, false);
  
   } else if(token == "off.mean.NeiDistance") {
  
     setNeiGeneticDistanceRecorders(OFFSPRG, false);
  
   } else if(token == "Dxy") {
  
     setDxyRecorders(ADULTS, false);
     setDxyRecorders(OFFSPRG, false);
  
   } else if(token == "off.Dxy") {
  
     setDxyRecorders(OFFSPRG, false);
  
   } else if(token == "adlt.Dxy") {
  
     setDxyRecorders(ADULTS, false);
  
   } else if(token == "Dxy.patch") {
  
     setDxyRecorders(ADULTS, true);
     setDxyRecorders(OFFSPRG, true);
  
   } else if(token == "off.Dxy.patch") {
  
     setDxyRecorders(OFFSPRG, true);
  
   } else if(token == "adlt.Dxy.patch") {
  
     setDxyRecorders(ADULTS, true);
  
   } else
     return false;
  
   return true;
 }