Main Page Namespace List Class Hierarchy Alphabetical List Compound List File List Namespace Members Compound Members File Members

matvec::Population Class Reference

#include <population.h>

List of all members.

Detailed Description

a genetic population consisting of individuals

See also:: Individual Pedigree

Definition at line 68 of file population.h.

Public Methods

void mark_descendant_recursive (const Individual *ind, const unsigned counter=1)

void mark_ancestor_recursive (const Individual *ind, const unsigned counter=1)

void mark_relative_recursive (Individual *ind, const unsigned counter=1)

void partial_iterative_peeling (Individual *I, doubleMatrix &wspace)

void setPenetrance (penetranceType p)

Population (void)

Population (GeneticDist *D)

Population (const int maxsize, GeneticDist *D)

Population (const Population &A)

~Population (void)

const Population & operator= (const Population &A)

void copyfrom (const Population &A)

void initialize (void)

void inbcoef_quaas (void)

void inbcoef_meuwissen (void)

void confirm_sex (const int col_sex)

doubleMatrix rela (void)

doubleMatrix relv (const double er)

doubleMatrix relv_nomissing (const double er)

double cprob_op (const unsigned goff[], const unsigned gparent[], const int ori)

double llhood_genotype (const unsigned c_id, const unsigned l_id)

double llhood_phenotype (const unsigned nrep)

doubleMatrix relv_inv (const double er)

void allele_freq (const char type[], const unsigned nal,...)

void penetrance_ml (const char fname[])

void release (void)

void resize (const unsigned maxn, GeneticDist *D)

void build_offs_info (void)

void build_spouse_info (void)

void build_trans_mat (void)

void renum (void)

void display (const char key[])

unsigned input_ped (const char fname[], const char recfmt[], GeneticDist *D)

unsigned input_ped0 (matvec::Pedigree &P, GeneticDist *D)

unsigned input_ped (matvec::Pedigree &P, GeneticDist *D)

unsigned input_data (Data *D)

unsigned input_data (const char fname[], GeneticDist *G)

unsigned size (void) const

unsigned nbase (void) const

unsigned get_id (const char name[])

unsigned input_markerData (Data *D)

unsigned input_descentGraph (char *dgfile)

void output_descentGraph (char *dgfile)

const char * ind_name (const unsigned id) const

void compute_pdm (const Individual *ind, Matrix< double > &pdm)

Vector< double > inbcoef (void)

void setup_m_ww (SparseMatrix &mme, Vector< double > &rhs, Vector< int > &start_addr, const double ev_e)

void add_Ga_inv (SparseMatrix &mme, const int start_addr, const double ev_a)

void add_Gv_inv (SparseMatrix &mme, const int start_addr, const double ev_v, const double er)

SparseMatrix * setup_m_mme1 (Vector< double > &rhs, const double ev_v, const double ev_u, const double ev_e, const double er)

SparseMatrix * setup_m_mme (Vector< double > &rhs, const double ev_v, const double ev_u, const double ev_e, const double er)

Individual * member (unsigned k)

Vector< double > * mblup1 (const double ev_v, const double ev_u, const double ev_e, const double er)

Vector< double > * mblup (const double ev_v, const double ev_u, const double ev_e, const double er)

void remove_mark (void)

unsigned mark_descendant_of (const unsigned id, const unsigned counter=1)

unsigned mark_ancestor_of (const unsigned id, const unsigned counter=1)

unsigned mark_relative_of (const unsigned id, const unsigned counter=1)

unsigned mark_families (const unsigned start_family_id=1)

Population * sub (const unsigned subsize)

Population * ancestor_of (const unsigned id)

Population * descendant_of (const unsigned id)

Population * relative_of (const unsigned id)

Population * family_of (const unsigned id)

Population * ancestor_of (const char name[])

Population * descendant_of (const char name[])

Population * relative_of (const char name[])

Population * family_of (const char name[])

unsigned fetch_families (const char filename[])

unsigned fetch_families (std::ostream &stream)

unsigned n_nufamily (void) const

unsigned ind_gamete (const Individual *I, const unsigned jc, unsigned g_id[], double fq[])

double cprob_children (const Individual *I, const unsigned jc, unsigned sg[], unsigned dg[], double sgf[], double dgf[])

unsigned full_cdist (Individual *I, const unsigned jc, unsigned **cdist_value, double *cdist_prob, unsigned **gid_mat, double **freq_mat)

unsigned partial_cdist (Individual *I, const unsigned jc, unsigned **cdist_value, double *cdist_prob, unsigned **gid_mat, double **freq_mat)

int peeling_sequence (void)

int detect_loop (void)

void maxant_maxpost (void)

void break_loop (void)

void cond_genotype_config ()

void genotype_config (const char type[])

void build_pop_gamete ()

void gibbs_iterate (const int count_genotype=0)

void count_genotype (Individual *I)

void release_genotype_counter (void)

void resize_genotype_counter (void)

void build_nufamily (void)

void anterior (Individual *I, doubleMatrix &wspace)

void posterior (Individual *I, Individual *J, doubleMatrix &wspace, const unsigned pj)

void anterior_posterior (Individual *I, doubleMatrix &wspace)

void genotype_dist_peeling (const int prtlevel=1, const int estifreq=1)

double log_likelihood_peeling (const unsigned maxit=3)

double get_posterior (Individual *I, Individual *exclJ, Vector< double > &vec)

double fullsibs_prob (Individual *dam, Individual *sire, Individual *excludeI, doubleMatrix &wspace)

Vector< double > get_genotype_freq (void)

void maxant_maxpost_old (void)

NuFamily ** nufamily (void)

const DataNode * record_ind (const unsigned i) const

void set_switches (void)

double multipoint_likelihood (int maxit)

double multipoint_init_parm (Fpmm &Farg)

void multi_geno_dist_peeling (const unsigned prtlevel)

double multi_m_log_likelihood_peeling (int maxit)

void multi_m_geno_dist_peeling (const unsigned prtlevel)

void build_connected_groups (unsigned)

int update_allele_vectors (unsigned, Individual *ind)

void build_founder_allele_neighbors (unsigned lcs)

int build_allele_vector (unsigned lcs, unsigned connected_group)

void process_alleles_neighbors (unsigned lcs, unsigned founder_allele)

void show_descent_graph_stuff (void)

double calc_prior_descent_graph (unsigned lcs)

void descent_graph_init_parm (void)

double descent_graph_log_lhood ()

void descent_graph_setup (Data *)

double M_H_sample_ext (double l_hood, unsigned &option)

double MH_ibd_sample (unsigned l, double log_lhood, unsigned &option)

double calc_log_q_ratio (void)

void show_change (const int, const int)

void sum_descentState (unsigned)

void sum_genotype_freq1 (unsigned)

void sum_genotype_freq2 (unsigned)

void sum_genotype_freq3 (unsigned)

void output_pdq (int n_samples, char *fname)

void Gibbs_sample (void)

void update_ibdMatrix (unsigned locus)

void output_ibdMatrix (unsigned n_samples, char *fname)

void countFounders (void)

void initJointAlleleNodeList (unsigned howManyLoci)

void initAlleleNodeList (unsigned whichLocus)

creates the allele node list, and the founder, penetrance, and transmission sets for each allele node

void initGenotypeNodeList (unsigned whichLocus)

creates the genotype node list, and the founder, penetrance, and transmission sets for each genotype node

void getInitialGNodeListSample (unsigned maxsize, unsigned startBlock, unsigned stopBlock, unsigned sizeBlock, string samplerType)

void getInitialGNodeListSample (unsigned maxsize, unsigned numLoci, string samplerType)

void getOldGNodeListProbability (unsigned maxsize, unsigned startBlock, unsigned stopBlock, unsigned sizeBlock, string samplerType)

recompute the needed quantities for the existing (old) sample (for the MH step)

void getGNodeListSample (unsigned maxsize, unsigned startBlock, unsigned stopBlock, unsigned sizeBlock, string samplerType)

obtain a new candidate sample using peeling and cutting (if necessary)

void getGNodeListSample (unsigned maxsize, unsigned numLoci, string samplerType)

void setSegregationIndex (unsigned atLocus, string samplerType)

set the segregation indicators when these can be infered from the ordered genotypes that have been sampled

void setupRSampler (Pedigree &P, Data *D, GeneticDist *G)

method to setup the structures needed by the R-sampler

void SimpleGenotypeElimination (unsigned locus)

eliminates genotypes using the individuals phenotypes

void LGGenotypeElimination (unsigned locus)

method to do Lange and Goradia genotype elimination

void setAlleleStateVectors (unsigned locus)

void buildRecombinationMatrix (void)

void finishUpAlleleStateInit (unsigned locus)

void switchStateBack (string samplerType)

void copyGNodeStatesToCandidateStates (string samplerType)

saves the GNode states (alleles or genotypes) for candidate sample

void copyCandidateToAccepted (string samplerType)

save candidate state as accepted states; allele origins for accepted states are already in (p)m_gameteOld

void storeSampledGametes (void)

store sampled gametes to be able to retreive them if we retreive the gametes sampled by Gibbs if MH rejects the sample

void retreiveSampledGametes (void)

void initGenotypeFreq (void)

void countHaplotypes (string samplerType)

void countGenotypes (string samplerType)

void displayHaplotypeFrequencies (unsigned numSamples, std::ostream &os=cout)

void displayGenotypeFrequencies (unsigned numSamples, std::ostream &outfile=cout)

void displaySegregationIndicators (std::ostream &outfile=cout)

unsigned calcTotalNumberOfGenotypes (unsigned locus)

method to do count total number of possible genotypes across the whole pedigree

void sampleSegregationIndicators (void)

sample the unknown segregation indicators of non-founders based on the ordered parental genotypes

void PrintSegregationIndicators (void)

print the segregation indicators of the non-founders

void resizeSegregationIndicators (void)

double MH_ibd_sample_map (double log_lhood, unsigned option)

void sum_descentState_map (void)

void setHaploSegregationIndex (unsigned atLocus, unsigned atMarker)

void ListAlleleFounders (void)

void SetPossibleHaplotypes (void)

void SetFreqHaploFounders (void)

void UpdateFreqHaploFounders (void)

void CalcFreqHaploFounders (unsigned numOfSamples)

void DisplayFreqHaploFounders (void)

Static Public Methods

void gtindex (const unsigned num, const unsigned ni, const Vector< unsigned > &ngt, Vector< unsigned > &gtvec)

Public Attributes

Matrix< float > ibdMatrix

int stdid

int mark_need_remove

int sex_confirmed

int offs_info_built

int spouse_info_built

int fdone

bool nuFamiliesDone

double scaling_val

NuFamily ** nufamily_vec

unsigned num_nufamily

unsigned num_t_nufamily

unsigned nonloop_t_nufamily

unsigned tn_genotype

unsigned tn_gamete

unsigned maxnamelen

unsigned tn_qtl

unsigned n_markerLoci

doubleMatrix * trans_mat

unsigned popsize

unsigned maxpopsize

unsigned numchrom

unsigned numtrait

Individual * pm_storage

Individual ** popmember

Vector< double > blupsol

std::string evaluate_method

Genome * pop_gamete

Genome * gametebase

GeneticDist * prior

int nmarker

std::string * markersymbol

penetranceType penetrance_f

SparseMatrix hmmec

HashTable hashtable

doubleMatrix residual_var

Vector< double > mean_for_genotype

Matrix< int > switch_table

Matrix< int > switch_table_gmt

doubleMatrix switch_table_prob

Vector< double > Q_freq

Vector< Vector< double > > marker_freq

Vector< Vector< double > > locusFreq

std::string analysis_type

Vector< double > mean_for_pgenotype

Vector< double > P_freq

int P_ndim

Fpmm * F

int founder_allele_counter

Vector< unsigned > connected_groups

Vector< int > allele_vector1

Vector< int > allele_vector2

unsigned connect_counter

Vector< double > RecoVector

unsigned SLlocus

Model * model

int valid_graph

std::vector< std::vector<
int > > founder_allele_neighbors

Matrix< unsigned > re_pdq

Vector< unsigned > temp

double log_lhood

unsigned initial

unsigned numFounders

Matrix< double > recombinationMatrix

SafeSTLVector< GNodeSet * > vectorOfGNsts

GNodeList gNodeList

bool lookOnlyAtLocusToYourLeft

IndexVector haplotypeCoder

unsigned numHaplotypes

bool markerDataIn

Static Public Attributes

double disPenetranceTable [2][3]

Friends

class Model

class GLMM

class Individual

class NuFamily

Constructor & Destructor Documentation

matvec::Population::Population ( void )

Definition at line 34 of file population.cpp.

00035 { 00036 initialize(); 00037 }

matvec::Population::Population ( GeneticDist * D )

Definition at line 39 of file population.cpp.

00040 { 00041 initialize(); 00042 prior = D; 00043 }

matvec::Population::Population ( const int maxsize,

GeneticDist * D

)

Definition at line 51 of file population.cpp.
References resize().

00052 { 00053 resize(maxsize,D); // hashtable needs to resize() afterwards if necessary 00054 }

matvec::Population::Population ( const Population & A )

Definition at line 45 of file population.cpp.
References copyfrom(), and initialize().

00046 { 00047 initialize(); 00048 copyfrom(A); 00049 }

matvec::Population::~Population ( void ) [inline]

Definition at line 116 of file population.h.
References release().

00116 {release();}

Member Function Documentation

void matvec::Population::add_Ga_inv ( SparseMatrix & mme,

const int start_addr,

const double ev_a

)

Definition at line 538 of file pop_mblup.cpp.
References matvec::Matrix< T >::begin(), matvec::Individual::father_id(), matvec::Individual::father_inbcoef(), matvec::getlambda(), matvec::Individual::id(), matvec::SparseMatrix::insert(), matvec::Individual::mother_id(), matvec::Individual::mother_inbcoef(), popmember, and popsize.

00540 { 00541 int i,j,k,ii,jj; 00542 Matrix<double> lambda(3,3); 00543 getlambda(lambda.begin(),3); 00544 for (i=0; i<3; i++) for (j=0; j<3; j++) lambda[i][j] /= ev_a; 00545 Individual *I; 00546 double dval; 00547 unsigned asd[3]; 00548 00549 for (k=0; k<popsize; k++) { 00550 I = popmember[k]; 00551 asd[0] = I->id(); asd[1] = I->father_id(); asd[2] = I->mother_id(); 00552 dval = 4.0/(2.0 - I->father_inbcoef() - I->mother_inbcoef()); 00553 for (i=0; i<3; i++) { 00554 if (asd[i] != 0) { 00555 ii = start_addr + asd[i]; 00556 for (j=0; j<3; j++) { 00557 if (asd[j] != 0) { 00558 jj = start_addr + asd[j]; 00559 if (jj >= ii) mme.insert(ii,jj,lambda[i][j]*dval); 00560 } 00561 } 00562 } 00563 } 00564 } 00565 }

void matvec::Population::add_Gv_inv ( SparseMatrix & mme,

const int start_addr,

const double ev_v,

const double er

)

Definition at line 567 of file pop_mblup.cpp.
References matvec::Locus::allele, matvec::Genome::chromosome, matvec::compute_pdq(), matvec::Individual::father_id(), matvec::Individual::genome0, matvec::Individual::genome1, matvec::Individual::id(), matvec::inva22(), matvec::Chromosome::locus, and matvec::Individual::mother_id().

00574 { 00575 int invert,i,j,k,kk,t,ii,jj,iii,jjj; 00576 double v1,v2; 00577 unsigned sdi[3]; 00578 Matrix<double> pdm(2,4); 00579 Matrix<double> pdq(2,6); 00580 Matrix<double> blk(6,6); 00581 Matrix<double> dinv(2,2); 00582 double inv_v = 1.0/ev_v; // ev_u must not be 0.0 00583 invert = 0; 00584 pdq[0][4]=1.0; pdq[0][5]=0.0; 00585 pdq[1][4]=0.0; pdq[1][5]=1.0; 00586 00587 Individual *I; 00588 for (t=0; t<popsize; t++) { 00589 I = popmember[t]; 00590 if (I->genome0.chromosome[0].locus[0].allele == 0 || 00591 I->genome1.chromosome[0].locus[0].allele == 0) { 00592 throw exception(" Population::add_Gv_inv(...): ind is ungenotyped"); 00593 } 00594 compute_pdm(I,pdm); 00595 compute_pdq(er,pdm,pdq); 00596 00597 dinv[0][0] = 1.0; dinv[0][1] = 0.0; 00598 dinv[1][0] = 0.0; dinv[1][1] = 1.0; 00599 for (i=0; i<2; i++) { 00600 for (j=0; j<2; j++) { 00601 for (k=0; k<4; k++) dinv[i][j] -= pdq[i][k]*pdq[j][k]; 00602 } 00603 } 00604 invert = inva22(dinv); 00605 if (!invert) throw exception(" dinv is not invertable"); 00606 for (i=0; i<6; i++) { 00607 v1 = pdq[0][i]*dinv[0][0] + pdq[1][i]*dinv[1][0]; 00608 v2 = pdq[0][i]*dinv[0][1] + pdq[1][i]*dinv[1][1]; 00609 for (j=0; j<6; j++) { 00610 blk[i][j] = (v1*pdq[0][j] + v2*pdq[1][j])*inv_v; 00611 } 00612 } 00613 sdi[0] = I->father_id(); sdi[1] = I->mother_id(); sdi[2] = I->id(); 00614 kk = start_addr + 1; 00615 for (i=0; i<3; i++) { 00616 if (sdi[i] != 0) { 00617 ii = kk + (sdi[i]-1)*2; 00618 iii = i*2; 00619 for (j=0; j<3; j++) { 00620 if (sdi[j] != 0) { 00621 jj = kk + (sdi[j]-1)*2 ; 00622 if (jj >= ii) { 00623 jjj = j*2; 00624 mme.insert(ii,jj,blk[iii][jjj]); 00625 mme.insert(ii,jj+1,blk[iii][jjj+1]); 00626 mme.insert(ii+1,jj+1,blk[iii+1][jjj+1]); 00627 if (jj > ii) { 00628 mme.insert(ii+1,jj,blk[iii+1][jjj]); 00629 } 00630 } 00631 } 00632 } 00633 } 00634 } 00635 } 00636 }

void matvec::Population::allele_freq ( const char type[],

const unsigned nal,

...

)

Population * matvec::Population::ancestor_of ( const char name[] )

Definition at line 217 of file population1.cpp.
References ancestor_of(), and get_id().

00218 { 00219 unsigned id = get_id(name); 00220 if (id == 0) throw exception("Population::ancestor_of(): cannot found"); 00221 return ancestor_of(id); 00222 }

Population * matvec::Population::ancestor_of ( const unsigned id )

Definition at line 71 of file population1.cpp.
References mark_ancestor_of(), popsize, renum(), stdid, and sub().
Referenced by ancestor_of().

00072 { 00073 ////////////////////////////////////////////////////////////// 00074 // put all ancestors of individual id into a new population 00075 // user is resposible to release this new population 00076 // using either delete or resize(0,0,0); 00077 ////////////////////////////////////////////////////////////// 00078 00079 if (!stdid) renum(); 00080 if (id<1 || id >popsize) throw exception(" Population::mark_ancestor_of(): bad arg"); 00081 unsigned counter = 1; 00082 unsigned nd = mark_ancestor_of(id,counter); 00083 return sub(nd); 00084 }

void matvec::Population::anterior ( Individual * I,

doubleMatrix & wspace

)

Definition at line 777 of file pop_peeling.cpp.
References matvec::Individual::anterior, matvec::Matrix< double >::begin(), matvec::Individual::father(), fullsibs_prob(), matvec::Individual::get_penetrance(), get_posterior(), matvec::Individual::id(), ind_name(), matvec::Individual::mother(), matvec::Individual::nspouse(), tn_genotype, and trans_mat.
Referenced by anterior_posterior(), and partial_iterative_peeling().

00778 { 00779 ///////////////////////////////////////////////////////////////// 00780 // REQUIREMENTS: memory must be allocated for anterior 00781 ///////////////////////////////////////////////////////////////// 00782 if (I->base()) throw exception("Population::anterior(): it shoould initialize first"); 00783 double ai,val,scale,**trme; 00784 unsigned um,uf,ui; 00785 Vector<double> post_vec(tn_genotype); 00786 Vector<double> pen_vec(tn_genotype); 00787 Vector<double> apm(tn_genotype); 00788 Vector<double> apf(tn_genotype); 00789 Individual *father, *mother; 00790 00791 mother = I->mother(); 00792 father = I->father(); 00793 scale = fullsibs_prob(father,mother,I,wspace); //Beware of call changes here 00794 scale += mother->anterior[tn_genotype]; 00795 for (um=0; um<tn_genotype; um++) apm[um] = mother->anterior[um]; 00796 if (mother->nspouse() > 1) { 00797 scale += get_posterior(mother,father,post_vec); 00798 for (um=0; um<tn_genotype; um++) apm[um] *= post_vec[um]; 00799 } 00800 00801 scale += father->anterior[tn_genotype]; 00802 for (uf=0; uf<tn_genotype; uf++) apf[uf] = father->anterior[uf]; 00803 if (father->nspouse() > 1) { 00804 scale += get_posterior(father,mother,post_vec); 00805 for (uf=0; uf<tn_genotype; uf++) apf[uf] *= post_vec[uf]; 00806 } 00807 00808 scale += I->get_penetrance(pen_vec); 00809 for (val=0.0,ui=0; ui<tn_genotype; ui++) { 00810 trme = trans_mat[ui].begin(); 00811 for (ai=0.0,um=0; um<tn_genotype; um++) for(uf=0; uf<tn_genotype; uf++) { 00812 ai += apm[um]*trme[um][uf]*wspace[um][uf]*apf[uf]; 00813 } 00814 val += I->anterior[ui] = ai*pen_vec[ui]; 00815 } 00816 //This is new 00817 if (val == 0.0) throw exception(std::string("Population::anterior(): genotypes conflict. Related animal:") + ind_name(I->id()) ); 00818 for (ui=0; ui<tn_genotype; ui++) I->anterior[ui] /= val; 00819 scale += std::log(val); 00820 I->anterior[tn_genotype] = scale; 00821 }

void matvec::Population::anterior_posterior ( Individual * I,

doubleMatrix & wspace

)

Definition at line 862 of file pop_peeling.cpp.
References anterior().
Referenced by matvec::NuFamily::iterative_peeling().

00863 { 00864 ////////////////////////////////////////////////////////////////////// 00865 // iterative peeling over individual I 00866 // REQUIREMENTS: memories must be allocated for anterior and posterior 00867 /////////////////////////////////////////////////////////////////////// 00868 if (!(I->base() || I->anterior_iw == 2 || I->anterior_iw == 3)) { 00869 anterior(I,wspace); 00870 } 00871 00872 Individual **spouses = I->spouses(); 00873 unsigned nsp = I->nspouse(); 00874 for (unsigned j=0; j<nsp; j++) { 00875 if (!(I->posterior_iw[j] == 2 || I->posterior_iw[j] == 3)) { 00876 posterior(I,spouses[j],wspace,j); 00877 } 00878 } 00879 }

void matvec::Population::break_loop ( void )

Definition at line 349 of file pop_peeling.cpp.
References matvec::compare_nuf_ptr(), matvec::Individual::connect_keeper, matvec::NuFamily::connectors, matvec::NuFamily::cutting(), matvec::Individual::est_GV, matvec::Individual::father(), matvec::NuFamily::in_connector, matvec::Individual::mother(), matvec::Individual::myoffspring, nufamily_vec, num_t_nufamily, matvec::NuFamily::out_connector, matvec::Individual::related_family, matvec::NuFamily::seq_id, stdid, and matvec::NuFamily::workvec.
Referenced by peeling_sequence().

00350 { 00351 ////////////////////////////////////////////////////////////////////// 00352 // search for any circles by visiting neighbors by neighbors 00353 // REQUIREMENT: call detect_loop() first, ie. there must be loops 00354 ///////////////////////////////////////////////////////////////////// 00355 if (num_t_nufamily == num_nufamily) return; 00356 NuFamily *pre_family, *cur_family, *neighbor; 00357 Individual *enter_door, *exit_door; 00358 std::list<Individual *>::iterator connector; 00359 std::list<NuFamily *> family_loop_connector; 00360 std::list<NuFamily *>::iterator Fam; 00361 std::list<NuFamily *>::reverse_iterator rFam; 00362 pre_family = 0; 00363 cur_family = nufamily_vec[num_t_nufamily]; 00364 // cur_family = nufamily_vec[num_nufamily-1]; 00365 cur_family->in_connector = 0; 00366 while (cur_family->seq_id != 0) { // seq_id should not be 0 at this moment 00367 cur_family->seq_id = 0; // 0 indicates the family is visited 00368 family_loop_connector.push_back(cur_family); 00369 neighbor = 0; 00370 connector = cur_family->connectors.begin(); 00371 while ((connector != cur_family->connectors.end()) && ((*connector) == cur_family->in_connector)) { 00372 ++connector; 00373 } 00374 if (connector == cur_family->connectors.end()) throw exception("Population::break_loop(): a bug!"); 00375 while (connector != cur_family->connectors.end()) { 00376 Fam = (*connector)->related_family.begin(); 00377 while (Fam != (*connector)->related_family.end()) { 00378 if (*Fam != cur_family && *Fam != pre_family ) { 00379 neighbor = *Fam; // this may not a good neighbor, so chech 00380 if (neighbor->out_connector != (*connector)) break; 00381 } 00382 ++Fam; 00383 } 00384 if (neighbor) break; 00385 ++connector; 00386 } 00387 if (!neighbor) { // connector=0 too, a short circle has been found 00388 connector= cur_family->connectors.begin(); 00389 while ((connector!= cur_family->connectors.end())&& (*connector)==cur_family->in_connector) { 00390 ++connector; 00391 } 00392 if (connector == cur_family->connectors.end()) throw exception("Population::break_loop(): a bug!"); 00393 cur_family->out_connector = *connector; 00394 pre_family->in_connector = *connector; 00395 cur_family = pre_family; 00396 break; 00397 } 00398 cur_family->out_connector = *connector; 00399 pre_family = cur_family; 00400 cur_family = neighbor; 00401 cur_family->in_connector = *connector; 00402 } 00403 cur_family->seq_id = 1; // cur_family is the first family in the circle 00404 ///////////////////////////////////////////////////////////////////// 00405 // a circle has been found, cur_family is in this circle, 00406 // let's break it at the position of connector 00407 ///////////////////////////////////////////////////////////////////// 00408 /* 00409 Individual *C = connector[0]; 00410 int perfectcut = 0; 00411 double smallest_gp = 1.0e38; 00412 Fam = (NuFamily **)family_loop_connector.bottom(); 00413 while (Fam) { 00414 enter_door = Fam[0]->in_connector; 00415 exit_door = Fam[0]->out_connector; 00416 if (enter_door != Fam[0]->mother() && enter_door != Fam[0]->father()){ 00417 if (enter_door->est_GV < smallest_gp) { 00418 smallest_gp = enter_door->est_GV; 00419 C = enter_door; 00420 cur_family = *Fam; 00421 perfectcut = 1; 00422 } 00423 } 00424 if (exit_door != Fam[0]->mother() && exit_door != Fam[0]->father()){ 00425 if (exit_door->est_GV < smallest_gp) { 00426 smallest_gp = exit_door->est_GV; 00427 C = exit_door; 00428 cur_family = *Fam; 00429 perfectcut = 1; 00430 } 00431 } 00432 if (Fam[0]->seq_id == 1) break; 00433 Fam = (NuFamily **)family_loop_connector.previous(); 00434 } 00435 if (perfectcut == 0) { 00436 Individual *mother,*father; 00437 Fam = (NuFamily **)family_loop_connector.bottom(); 00438 while (Fam) { 00439 enter_door = Fam[0]->in_connector; 00440 exit_door = Fam[0]->out_connector; 00441 mother = Fam[0]->mother(); 00442 father = Fam[0]->father(); 00443 if (!enter_door->connect_keeper) { 00444 if (enter_door == mother) { 00445 if (Fam[0]->workvec[1] < smallest_gp) { 00446 smallest_gp = Fam[0]->workvec[1]; 00447 C = enter_door; 00448 cur_family = *Fam; 00449 } 00450 } 00451 else if (enter_door == father) { 00452 if (Fam[0]->workvec[0] < smallest_gp) { 00453 smallest_gp = Fam[0]->workvec[0]; 00454 C = enter_door; 00455 cur_family = *Fam; 00456 } 00457 } 00458 } 00459 if (!exit_door->connect_keeper) { 00460 if (exit_door == mother) { 00461 if (Fam[0]->workvec[1] < smallest_gp) { 00462 smallest_gp = Fam[0]->workvec[1]; 00463 C = exit_door; 00464 cur_family = *Fam; 00465 } 00466 } 00467 else if (exit_door == father) { 00468 if (Fam[0]->workvec[0] < smallest_gp) { 00469 smallest_gp = Fam[0]->workvec[0]; 00470 C = exit_door; 00471 cur_family = *Fam; 00472 } 00473 } 00474 } 00475 if (Fam[0]->seq_id == 1) break; 00476 Fam = (NuFamily **)family_loop_connector.previous(); 00477 } 00478 } 00479 // wait a minute, let's first update related_family in C 00480 C->related_family.clear(&cur_family,compare_nuf_ptr); 00481 cur_family->cutting(C); 00482 */ 00483 00484 Individual *C = *connector; 00485 int perfectcut = 0; 00486 double smallest_gp = 1.0e38; 00487 rFam = family_loop_connector.rbegin(); 00488 if (stdid >= 0) { 00489 while (rFam != family_loop_connector.rend()) { 00490 enter_door = (*rFam)->in_connector; 00491 exit_door = (*rFam)->out_connector; 00492 if ((*rFam)->numcut == 0) { 00493 if (enter_door != (*rFam)->mother() && enter_door != (*rFam)->father()) { 00494 if (enter_door->est_GV < smallest_gp) { 00495 smallest_gp = enter_door->est_GV; 00496 C = enter_door; 00497 cur_family = *rFam; 00498 perfectcut = 1; 00499 } 00500 } 00501 if (exit_door != (*rFam)->father() && exit_door != (*rFam)->mother()){ 00502 if (exit_door->est_GV < smallest_gp) { 00503 smallest_gp = exit_door->est_GV; 00504 C = exit_door; 00505 cur_family = *rFam; 00506 perfectcut = 1; 00507 } 00508 } 00509 } 00510 if ((*rFam)->seq_id == 1) break; 00511 ++rFam; 00512 } 00513 } 00514 else { 00515 while (rFam != family_loop_connector.rend()) { 00516 if (((*rFam)->mother_id() == 4 && (*rFam)->father_id() == 10) 00517 && (*rFam)->myoffspring[0]->id() == 11) { 00518 C = (*rFam)->myoffspring[0]; 00519 cur_family = *rFam; 00520 perfectcut = 1; 00521 } 00522 if ((*rFam)->seq_id == 1) break; 00523 ++rFam; 00524 } 00525 } 00526 if (perfectcut == 0) { 00527 Individual *father, *mother; 00528 rFam = family_loop_connector.rbegin(); 00529 while (rFam != family_loop_connector.rend()) { 00530 enter_door = (*rFam)->in_connector; 00531 exit_door = (*rFam)->out_connector; 00532 if ((*rFam)->numcut == 0) { 00533 father = (*rFam)->father(); 00534 mother = (*rFam)->mother(); 00535 if (!enter_door->connect_keeper) { 00536 if (enter_door == mother) { 00537 if ((*rFam)->workvec[1] < smallest_gp) { 00538 smallest_gp = (*rFam)->workvec[1]; 00539 C = enter_door; 00540 cur_family = *rFam; 00541 } 00542 } 00543 else if (enter_door == father) { 00544 if ((*rFam)->workvec[0] < smallest_gp) { 00545 smallest_gp = (*rFam)->workvec[0]; 00546 C = enter_door; 00547 cur_family = *rFam; 00548 } 00549 } 00550 } 00551 if (!exit_door->connect_keeper) { 00552 if (exit_door == mother) { 00553 if ((*rFam)->workvec[1] < smallest_gp) { 00554 smallest_gp = (*rFam)->workvec[1]; 00555 C = exit_door; 00556 cur_family = *rFam; 00557 } 00558 } 00559 else if (exit_door == father) { 00560 if ((*rFam)->workvec[0] < smallest_gp) { 00561 smallest_gp = (*rFam)->workvec[0]; 00562 C = exit_door; 00563 cur_family = *rFam; 00564 } 00565 } 00566 } 00567 } 00568 if ((*rFam)->seq_id == 1) break; 00569 ++rFam; 00570 } 00571 } 00572 // wait a minute, let's first update related_family in C 00573 C->related_family.remove(cur_family); 00574 cur_family->cutting(C); 00575 stdid++; 00576 }

int matvec::Population::build_allele_vector ( unsigned lcs,

unsigned connected_group

)

Definition at line 295 of file pop_graph.cpp.
References connected_groups, founder_allele_counter, process_alleles_neighbors(), and valid_graph.
Referenced by calc_prior_descent_graph(), and M_H_sample_ext().

00295 { 00296 // Authors: Mathias Schelling and Rohan L. Fernando 00297 // (1999) 00298 // Contributors: 00299 valid_graph = 1; // graph is valid 00300 for (int i=1; i<=founder_allele_counter; i++) { 00301 if(connected_groups(i) == connected_group) { 00302 process_alleles_neighbors(lcs,i); 00303 return valid_graph; 00304 } 00305 } 00306 return valid_graph; 00307 //return 0; 00308 }

void matvec::Population::build_connected_groups ( unsigned )

Definition at line 45 of file pop_graph.cpp.
References matvec::Locus::allele, matvec::Individual::assigned_founder_allele, matvec::Genome::chromosome, connect_counter, connected_groups, founder_allele_counter, matvec::Individual::genome0, matvec::Chromosome::locus, matvec::Individual::m_founder, matvec::Individual::p_founder, popmember, popsize, and matvec::Vector< unsigned >::resize().
Referenced by calc_prior_descent_graph(), and M_H_sample_ext().

00045 { 00046 // Authors: Mathias Schelling and Rohan L. Fernando 00047 // (1999) 00048 // Contributors: L. Radu Totir 00049 if (!founder_allele_counter) { 00050 throw exception("Individual::set_founder_alleles() must be called first" ); 00051 } 00052 // For some reason, I cannot put a constant as the second argument of resize 00053 // this should be fixed 00054 int zero=0; 00055 connected_groups.resize(founder_allele_counter,zero); 00056 00057 connect_counter = 0; 00058 for (int i=0; i<popsize; i++) { 00059 if(popmember[i]->genome0.chromosome[0].locus[lcs-1].allele) { //genotype known 00060 popmember[i]->assigned_founder_allele=0; 00061 int mfounder = popmember[i]->m_founder(lcs); 00062 int pfounder = popmember[i]->p_founder(lcs); 00063 if (!connected_groups(mfounder) && !connected_groups(pfounder)){ //both zero 00064 connect_counter++; 00065 connected_groups(mfounder) = connect_counter; 00066 connected_groups(pfounder) = connect_counter; 00067 } 00068 else if (connected_groups(mfounder) && connected_groups(pfounder)) {//both non-zero 00069 //set them equal 00070 int ccgpf = connected_groups(pfounder); 00071 int ccgmf = connected_groups(mfounder); 00072 for (int j=1; j<=founder_allele_counter; j++) { 00073 if (connected_groups(j) == ccgpf ) { 00074 connected_groups(j) = ccgmf; 00075 } 00076 } 00077 } 00078 else { //only one non-zero -> let the other inherit it's connected group 00079 if (connected_groups(mfounder)) { 00080 connected_groups(pfounder) = connected_groups(mfounder); 00081 } 00082 else { 00083 connected_groups(mfounder) = connected_groups(pfounder); 00084 } 00085 } 00086 } 00087 } 00088 }

void matvec::Population::build_founder_allele_neighbors ( unsigned lcs )

Definition at line 277 of file pop_graph.cpp.
References matvec::Locus::allele, allele_vector1, allele_vector2, matvec::Genome::chromosome, founder_allele_counter, founder_allele_neighbors, matvec::Individual::genome0, matvec::Chromosome::locus, matvec::Individual::m_founder, matvec::Individual::p_founder, popmember, popsize, and matvec::Vector< int >::resize().
Referenced by calc_prior_descent_graph(), and M_H_sample_ext().

00277 { 00278 // Authors: Mathias Schelling and Rohan L. Fernando 00279 // (1999) 00280 // Contributors: L. Radu Totir 00281 allele_vector1.resize(founder_allele_counter,0); //need to put this somewhere before building 00282 allele_vector2.resize(founder_allele_counter,0); //allele vectors. So just stuck it here 00283 founder_allele_neighbors.resize(0); 00284 founder_allele_neighbors.resize(founder_allele_counter); 00285 for (int i=0; i<popsize; i++){ 00286 if(popmember[i]->genome0.chromosome[0].locus[lcs-1].allele){ // genotype known 00287 int mfounder = popmember[i]->m_founder(lcs); 00288 int pfounder = popmember[i]->p_founder(lcs); 00289 founder_allele_neighbors[mfounder-1].push_back(i); 00290 founder_allele_neighbors[pfounder-1].push_back(i); 00291 } 00292 } 00293 }

void matvec::Population::build_nufamily ( void )

Definition at line 237 of file pop_peeling.cpp.
References build_spouse_info(), matvec::check_ptr(), matvec::NuFamily::father(), matvec::Individual::father(), matvec::NuFamily::father_indx, matvec::Individual::id(), ind_name(), matvec::NuFamily::mother(), matvec::Individual::mother(), matvec::NuFamily::mother_indx, matvec::NuFamily::noffs(), matvec::Individual::noffs_spouse(), matvec::Individual::nspouse(), nuFamiliesDone, nufamily(), nufamily_vec, num_nufamily, num_t_nufamily, matvec::NuFamily::offspring(), matvec::Individual::offspring(), matvec::NuFamily::pop, popmember, popsize, matvec::Individual::put_family(), matvec::Individual::sex(), spouse_info_built, and matvec::Individual::spouses().
Referenced by matvec::Model::multipoint_genot_probs(), peeling_sequence(), and setupRSampler().

00238 { 00239 if (!spouse_info_built){ 00240 build_spouse_info(); 00241 } 00242 Individual *I; 00243 unsigned i,j,nsp; 00244 num_t_nufamily = 0; 00245 num_nufamily = 0; 00246 for (i=0; i<popsize; i++) { 00247 I = popmember[i]; 00248 if (I->father() && !I->mother() || !I->father() && I->mother()) { 00249 j = I->id(); 00250 throw exception(std::string("one of parent is missing: ") + ind_name(j)); 00251 } 00252 nsp = I->nspouse(); 00253 if (I->sex() == 'F' && nsp > 0) num_nufamily += nsp; 00254 } 00255 if (num_nufamily == 0) return; 00256 nufamily_vec = new NuFamily*[num_nufamily]; 00257 check_ptr(nufamily_vec); 00258 for (i=0; i<num_nufamily; i++) nufamily_vec[i] = new NuFamily; 00259 unsigned k,ii,jj,msp,*mi; 00260 Individual *J,**spouses,**offspring; 00261 NuFamily* nufamily; 00262 const unsigned *noffs_sp; 00263 for (ii=0,i=0; i<popsize; i++) { 00264 I = popmember[i]; 00265 nsp = I->nspouse(); 00266 if (I->sex() == 'F' && nsp > 0) { 00267 offspring = I->offspring(); 00268 noffs_sp = I->noffs_spouse(); 00269 for (k=0, j=0; j<nsp; j++) { 00270 J = offspring[k]->father(); 00271 // J = offspring[k]->mother(); //old version 00272 nufamily = nufamily_vec[ii++]; 00273 nufamily->pop = this; 00274 nufamily->mother(I); 00275 nufamily->father(J); 00276 nufamily->offspring(&(offspring[k]),noffs_sp[j]); 00277 nufamily->father_indx = j; 00278 spouses = J->spouses(); 00279 msp = J->nspouse(); 00280 for (jj=0; jj<msp; jj++) if (spouses[jj] == I) { 00281 nufamily->mother_indx = jj; 00282 break; 00283 } 00284 k += noffs_sp[j]; 00285 } 00286 } 00287 } 00288 //This appears only in the old version 00289 for (i=0; i<num_nufamily; i++) { 00290 int noffs = nufamily_vec[i]->noffs(); 00291 offspring = nufamily_vec[i]->offspring(); 00292 for (k=0;k<noffs;k++) offspring[k]->put_family(i); 00293 } 00294 nuFamiliesDone = true; 00295 }

void matvec::Population::build_offs_info ( void )

Definition at line 507 of file population.cpp.
References matvec::Individual::father(), mark_need_remove, matvec::Individual::mother(), matvec::Individual::myoffspring, matvec::Individual::numoffs, offs_info_built, matvec::Individual::offs_tree, popmember, popsize, renum(), and stdid.
Referenced by build_spouse_info(), input_ped(), and sub().

00508 { 00509 if (offs_info_built) return; 00510 if (!stdid) renum(); 00511 Individual *I; 00512 unsigned i,j,noffs; 00513 for (i=0; i<popsize; i++) { 00514 I = popmember[i]; 00515 if (I->mother()) I->mother()->offs_tree.insert(I); 00516 if (I->father()) I->father()->offs_tree.insert(I); 00517 } 00518 for (i=0; i<popsize; i++) { 00519 I = popmember[i]; 00520 if (I->myoffspring) { 00521 delete [] I->myoffspring;I->myoffspring=0; 00522 } 00523 noffs = I->offs_tree.size(); 00524 I->numoffs = noffs; 00525 if (noffs > 0) { 00526 I->myoffspring = new Individual* [noffs]; 00527 std::set<Individual *>::iterator pos; 00528 for (j=0,pos=I->offs_tree.begin(); pos != I->offs_tree.end(); ++pos,++j) { 00529 I->myoffspring[j] = *pos; 00530 } 00531 I->offs_tree.clear(); 00532 } 00533 } 00534 mark_need_remove = 1; 00535 offs_info_built = 1; 00536 }

void matvec::Population::build_pop_gamete ( )

Definition at line 97 of file pop_gibbs.cpp.
References matvec::Locus::allele, matvec::GeneticDist::chrom(), matvec::Genome::chromosome, matvec::Chromosome::freq(), gametebase, gtindex(), matvec::Chromosome::id(), matvec::Chromosome::locus, matvec::ChromStruct::locus, matvec::new_Genome_vec(), matvec::ChromStruct::nloci(), numchrom, pop_gamete, prior, renum(), matvec::Genome::resize(), stdid, and tn_gamete.
Referenced by build_trans_mat(), cond_genotype_config(), genotype_config(), get_genotype_freq(), llhood_phenotype(), and resize_genotype_counter().

00098 { 00099 ///////////////////////////////////////////////////////////// 00100 // It is assumed that population is in equilibrium status 00101 ///////////////////////////////////////////////////////////// 00102 00103 if (pop_gamete) return; 00104 if (!stdid) renum(); 00105 if (pop_gamete) { 00106 delete [] pop_gamete; 00107 pop_gamete=0; 00108 } 00109 pop_gamete = new_Genome_vec(numchrom); 00110 unsigned chr,nl,nc,num,allele_id; 00111 unsigned i,j,k,h; 00112 double gamete_freq; 00113 ChromStruct *Chrom = prior->chrom(); 00114 Genome *G, *GB; 00115 00116 for (tn_gamete=1,chr = 0; chr< numchrom; chr++) { 00117 nl = Chrom[chr].nloci(); 00118 Vector<unsigned> ngtvec(nl); 00119 Vector<unsigned> gtvec(nl); 00120 for (num=1,k=0; k<nl; k++) { 00121 num *= ngtvec[k] = Chrom[chr].locus[k].nallele(); 00122 } 00123 tn_gamete *= num; 00124 G = &(pop_gamete[chr]); 00125 G->resize(num,nl); 00126 for (i=0; i<num; i++) { 00127 gtindex(i,nl,ngtvec,gtvec); 00128 for (gamete_freq=1.0, k=0; k<nl; k++) { 00129 allele_id = gtvec[k]; 00130 G->chromosome[i].locus[k].allele = allele_id; 00131 gamete_freq *= Chrom[chr].locus[k].allele_freq[allele_id]; 00132 } 00133 G->chromosome[i].id(i); 00134 G->chromosome[i].freq(gamete_freq); 00135 } 00136 } 00137 00138 ///////////////////// build gamete base information //////////////////// 00139 // gamete id 0, 1, 2,.... num-1. There are num*(num-1)/2 combinations: 00140 // (1,0) 00141 // (2,0) (2,1) 00142 // (3,0) (3,1) (3,2) 00143 // . . . 00144 /////////////////////////////////////////////////////////////////////// 00145 if(gametebase){ 00146 if (gametebase) delete [] gametebase; 00147 gametebase=0; 00148 } 00149 gametebase = new_Genome_vec(numchrom); 00150 Chromosome gamete0, *gamete1, *gamete2; 00151 double r_k2_kk, freq; 00152 unsigned ngt,k2,kk,gamete_id; 00153 for (chr = 0; chr<numchrom; chr++) { 00154 nl = Chrom[chr].nloci(); 00155 gamete0.resize(nl); 00156 Vector<unsigned> ngtvec(nl); 00157 Vector<unsigned> gtvec(nl); 00158 for (num=1,i=0; i<nl; i++) num *= Chrom[chr].locus[i].nallele(); 00159 GB = &(gametebase[chr]); 00160 GB->resize(num*(num-1)/2, 0); 00161 G = &(pop_gamete[chr]); 00162 00163 for (nc =0, i=1; i<num; i++) { 00164 gamete1 = &(G->chromosome[i]); 00165 for (j=0; j<i; j++) { 00166 gamete2 = &(G->chromosome[j]); 00167 for (ngt=1, k=0; k<nl; k++) { 00168 if (gamete1->locus[k].allele == gamete2->locus[k].allele) 00169 ngtvec[k] = 1; 00170 else { 00171 ngtvec[k] = 2; 00172 ngt *= 2; 00173 } 00174 } 00175 GB->chromosome[nc].resize(ngt); 00176 for (k=0; k<ngt; k++) { 00177 gtindex(k,nl,ngtvec,gtvec); 00178 freq = 0.5; 00179 k2 = 0; // the first 2-different allele case hasn't found 00180 for (kk=0; kk<nl; kk++) { 00181 if (ngtvec[kk] == 2) { 00182 if (k2 > 0 ) { // not the first 2-different allele case 00183 r_k2_kk = prior->recomb_rate(chr+1,k2,kk+1); 00184 if (gtvec[kk] == gtvec[k2-1]) freq *= 1.0 - r_k2_kk; 00185 else freq *= r_k2_kk; 00186 } 00187 k2 = kk+1; 00188 } 00189 } 00190 /////////////////////////////////////////////////////////////// 00191 // GB->chromosome[nc].locus[k].allele = gamete_id in pop_gamete 00192 // GB->chromosome[nc].locus[k].effect = gamete freq 00193 /////////////////////////////////////////////////////////////// 00194 GB->chromosome[nc].locus[k].effect = freq; 00195 for (kk=0; kk<nl; kk++) { 00196 if (gtvec[kk] == 0) { 00197 gamete0.locus[kk].allele = gamete1->locus[kk].allele; 00198 } 00199 else gamete0.locus[kk].allele = gamete2->locus[kk].allele; 00200 } 00201 gamete_id = G->chrom_id(gamete0); 00202 GB->chromosome[nc].locus[k].allele = gamete_id; 00203 } 00204 nc++; // nc = i*(i-1)/2 + j 00205 } 00206 } 00207 } 00208 tn_genotype = tn_gamete*(tn_gamete+1)/2; 00209 00210 ////////////////////////////////////////////////////////// 00211 // build a vector of genotypic value: mean_for_genotype 00212 ////////////////////////////////////////////////////////// 00213 mean_for_genotype.resize(tn_genotype); 00214 unsigned nchrom = prior->nchrom(); 00215 std::cout << "pop nchrom " << nchrom << std::endl; 00216 Vector<unsigned> ngtvec(nchrom); 00217 Vector<unsigned> sire_gtvec(nchrom); 00218 Vector<unsigned> dam_gtvec(nchrom); 00219 for (i=0; i<nchrom; i++) ngtvec[i] = prior->chrom()[i].locus[0].nallele(); 00220 00221 for (kk=0,i=0; i<tn_gamete; i++) { 00222 gtindex(i,nchrom,ngtvec,sire_gtvec); 00223 for (j=0; j<=i; j++) { 00224 gtindex(j,nchrom,ngtvec,dam_gtvec); 00225 for (freq=0.0,k=0; k<nchrom; k++) { 00226 for (h=0; h<nl; h++) { // nl = 1 00227 freq += prior->genotypic_val(k,h)[sire_gtvec[k]][dam_gtvec[k]]; 00228 } 00229 } 00230 mean_for_genotype[kk++] = freq; 00231 } 00232 } 00233 }

void matvec::Population::build_spouse_info ( void )

Definition at line 538 of file population.cpp.
References build_offs_info(), matvec::compare_father_id(), matvec::compare_mother_id(), matvec::Individual::father(), mark_need_remove, matvec::Individual::mother(), matvec::Individual::mysex, matvec::Individual::numoffs, matvec::Individual::numoffs_spouse, matvec::Individual::numspouse, offs_info_built, matvec::Individual::offspring(), popmember, popsize, renum(), spouse_info_built, matvec::Individual::spouselist, and stdid.
Referenced by build_nufamily(), and input_ped().

00539 { 00540 if (spouse_info_built) return; 00541 if (!stdid) renum(); 00542 if (!offs_info_built) build_offs_info(); 00543 00544 Individual *I,*spouse,**offspring; 00545 unsigned i,j,k,s,noffs,nsp; 00546 for (i=0; i<popsize; i++) { 00547 I = popmember[i]; 00548 if (I->numoffs_spouse) { 00549 delete [] I->numoffs_spouse; I->numoffs_spouse = 0; 00550 } 00551 if (I->spouselist) { 00552 delete [] I->spouselist; I->spouselist = 0; 00553 } 00554 noffs = I->numoffs; 00555 offspring = I->offspring(); 00556 if (noffs == 1) { 00557 I->numspouse = 1; 00558 I->spouselist = new Individual*[1]; 00559 if (I->mysex == 'F') { 00560 I->spouselist[0] = offspring[0]->father(); 00561 } 00562 else { 00563 I->spouselist[0] = offspring[0]->mother(); 00564 } 00565 I->numoffs_spouse = new unsigned [1]; 00566 I->numoffs_spouse[0] = 1; 00567 } 00568 else if (noffs > 1) { 00569 j = 1; nsp = 1; 00570 if (I->mysex == 'F') { 00571 qsort(offspring,noffs,sizeof(Individual*),compare_father_id); 00572 spouse = offspring[0]->father(); 00573 while(j<noffs) { 00574 if (spouse != offspring[j]->father()) { 00575 spouse = offspring[j]->father(); 00576 nsp++; 00577 } 00578 j++; 00579 } 00580 I->numspouse = nsp; 00581 if(nsp>0){ 00582 I->numoffs_spouse = new unsigned [nsp]; 00583 } 00584 else { 00585 I->numoffs_spouse = 0; 00586 } 00587 if(nsp>0){ 00588 I->spouselist = new Individual* [nsp]; 00589 } 00590 else { 00591 I->spouselist = 0; 00592 } 00593 spouse = offspring[0]->father(); 00594 I->spouselist[0] = spouse; 00595 j = 1; s = 0; k = 1; 00596 while(j<noffs) { 00597 if (spouse != offspring[j]->father()) { 00598 spouse = offspring[j]->father(); 00599 I->numoffs_spouse[s++] = k; 00600 I->spouselist[s] = spouse; 00601 k = 1; 00602 } 00603 else { 00604 k++; // count # of myoffspring of each spouses 00605 } 00606 j++; 00607 } 00608 I->numoffs_spouse[s] = k; 00609 } 00610 else { 00611 qsort(offspring,noffs,sizeof(Individual*),compare_mother_id); 00612 spouse = offspring[0]->mother(); 00613 while(j<noffs) { 00614 if (spouse != offspring[j]->mother()) { 00615 spouse = offspring[j]->mother(); 00616 nsp++; 00617 } 00618 j++; 00619 } 00620 I->numspouse = nsp; 00621 if(nsp>0){ 00622 I->numoffs_spouse = new unsigned [nsp]; 00623 } 00624 else { 00625 I->numoffs_spouse = 0; 00626 } 00627 if(nsp>0){ 00628 I->spouselist = new Individual* [nsp]; 00629 } 00630 else { 00631 I->spouselist = 0; 00632 } 00633 spouse = offspring[0]->mother(); 00634 I->spouselist[0] = spouse; 00635 j = 1; s = 0; k = 1; 00636 while(j<noffs) { 00637 if (spouse != offspring[j]->mother()) { 00638 spouse = offspring[j]->mother(); 00639 I->numoffs_spouse[s++] = k; 00640 I->spouselist[s] = spouse; 00641 k = 1; 00642 } 00643 else { 00644 k++; // count # of smyoffspring of each spouses 00645 } 00646 j++; 00647 } 00648 I->numoffs_spouse[s] = k; 00649 } 00650 } 00651 } 00652 mark_need_remove = 1; 00653 spouse_info_built = 1; 00654 }

void matvec::Population::build_trans_mat ( void )

Definition at line 168 of file pop_peeling.cpp.
References build_pop_gamete(), matvec::check_ptr(), matvec::GeneticDist::chrom(), matvec::GeneticDist::nchrom(), matvec::ChromStruct::nloci(), pop_gamete, prior, matvec::Matrix< double >::resize(), tn_gamete, tn_genotype, and trans_mat.
Referenced by peeling_sequence().

00169 { 00170 //This if statement appears in the old version 00171 if (prior->chrom()[0].nloci() != 1) { 00172 trans_mat = new doubleMatrix[1]; 00173 return; 00174 } 00175 if (trans_mat) return; 00176 if (!pop_gamete) build_pop_gamete(); 00177 if(tn_genotype>0){ 00178 trans_mat = new doubleMatrix[tn_genotype]; 00179 } 00180 else { 00181 trans_mat = 0; 00182 } 00183 check_ptr(trans_mat); 00184 unsigned s,s1,s2,d,d1,d2,p,p1,p2,ss1,ss2; 00185 d = prior->nchrom(); 00186 for (s=0; s<d; s++) { 00187 if (prior->chrom()[s].nloci() != 1) { 00188 if(trans_mat){ 00189 delete [] trans_mat; 00190 trans_mat=0; 00191 } 00192 throw exception("Population::build_trans_mat(): one locus/chromosome is alowed"); 00193 } 00194 } 00195 for (s=0; s<tn_genotype; s++) trans_mat[s].resize(tn_genotype,tn_genotype,0.0); 00196 for (s=0,s1=0; s1<tn_gamete; s1++) { 00197 ss1 = s1*(s1+1)/2; 00198 for (s2=0; s2<=s1; s2++) { 00199 ss2 = s2*(s2+1)/2; 00200 for (d=0,d1=0; d1<tn_gamete; d1++) { 00201 if (s1 >= d1) { 00202 p1 = ss1 + d1; 00203 } 00204 else { 00205 p1 = d1*(d1+1)/2 + s1; 00206 } 00207 if (s2 >= d1) { 00208 p2 = ss2 + d1; 00209 } 00210 else { 00211 p2 = d1*(d1+1)/2 + s2; 00212 } 00213 for (d2=0; d2<=d1; d2++) { 00214 trans_mat[p1][s][d] += 0.25; trans_mat[p2][s][d] += 0.25; 00215 if (s1 >= d2) { 00216 p = ss1 + d2; 00217 } 00218 else { 00219 p = d2*(d2+1)/2 + s1; 00220 } 00221 trans_mat[p][s][d] += 0.25; 00222 if (s2 >= d2) { 00223 p = ss2 + d2; 00224 } 00225 else { 00226 p = d2*(d2+1)/2 + s2; 00227 } 00228 trans_mat[p][s][d] += 0.25; 00229 d++; 00230 } 00231 } 00232 s++; 00233 } 00234 } 00235 }

void matvec::Population::buildRecombinationMatrix ( void )

Definition at line 1239 of file population.cpp.
References matvec::GeneticDist::get_distance(), matvec::Model::MapF(), model, prior, recombinationMatrix, and matvec::Matrix< double >::resize().
Referenced by setupRSampler().

01239 { 01240 // Authors: L. Radu Totir and Rohan L. Fernando 01241 // (September, 2003) 01242 // Contributors: 01243 unsigned noMarkerLoci = Individual::numLoci-Individual::nQTL; 01244 recombinationMatrix.resize(noMarkerLoci,noMarkerLoci,0.0); 01245 for(unsigned i=1; i<noMarkerLoci; i++) { 01246 for (unsigned j=1; j<noMarkerLoci; j++) { 01247 if(j>=i){ 01248 recombinationMatrix[i-1][j] = model->MapF(prior->get_distance(1,j+1) - prior->get_distance(1,i)); 01249 //recombinationMatrix[i-1][j] = 0.5; 01250 } 01251 } 01252 } 01253 std::cout << "Recombination Matrix:\n" << recombinationMatrix; 01254 }

double matvec::Population::calc_log_q_ratio ( void )

Definition at line 593 of file pop_graph.cpp.
References matvec::Individual::calc_q(), matvec::Individual::m_gamete, matvec::Individual::m_gamete1, member(), matvec::Individual::p_gamete, matvec::Individual::p_gamete1, matvec::Individual::sample_haplotypes(), size(), SLlocus, and matvec::Individual::update_offsprings_founder_alleles().
Referenced by M_H_sample_ext().

00593 { 00594 // Authors: Mathias Schelling and Rohan L. Fernando 00595 // (1999) 00596 // Contributors: Fabiano V. Pita 00597 double accumulator = 0.0; 00598 for (int i=0; i<size(); i++) { 00599 unsigned flag = 0; 00600 if (member(i)->m_gamete(SLlocus) != member(i)->m_gamete1(SLlocus)) { 00601 flag = 1; 00602 member(i)->sample_haplotypes(0); 00603 } 00604 if (member(i)->p_gamete(SLlocus) != member(i)->p_gamete1(SLlocus)) { 00605 flag = 1; 00606 member(i)->sample_haplotypes(1); 00607 } 00608 if (flag != 0) { 00609 accumulator += member(i)->calc_q(); 00610 member(i)->update_offsprings_founder_alleles(); 00611 } 00612 } 00613 return accumulator; 00614 }

double matvec::Population::calc_prior_descent_graph ( unsigned lcs )

Definition at line 343 of file pop_graph.cpp.
References allele_vector1, allele_vector2, matvec::Individual::base(), build_allele_vector(), build_connected_groups(), build_founder_allele_neighbors(), connect_counter, connected_groups, founder_allele_counter, matvec::Individual::get_allele_v1(), initial, locusFreq, member(), and size().
Referenced by descent_graph_log_lhood(), and matvec::Individual::sample_self().

00343 { 00344 // Authors: Mathias Schelling and Rohan L. Fernando 00345 // (1999) 00346 // Contributors: Fabiano V. Pita 00347 build_connected_groups(lcs); 00348 build_founder_allele_neighbors(lcs); 00349 for (int i=1; i<=connect_counter; i++) { 00350 unsigned buildAlleleVectorResult = build_allele_vector(lcs,i); 00351 if (buildAlleleVectorResult == 0){ 00352 return 0.0; 00353 } 00354 } 00355 // FVP 03/12/03 here I check if the order of ordered heterozygous is respected 00356 // if not => DG is rejected 00357 if (initial == 1){ 00358 for (int i=0; i<size(); i++) { 00359 if (member(i) -> ord_heter == lcs){ 00360 int mfounder = member(i) -> m_founder(lcs); 00361 int pfounder = member(i) -> p_founder(lcs); 00362 00363 if (!((member(i) -> initial_order1 == allele_vector1(mfounder)) || 00364 (member(i) -> initial_order2 == allele_vector1(pfounder)))) { 00365 return 0.0; 00366 } 00367 } 00368 } 00369 } 00370 //show_descent_graph_stuff(); 00371 double prior_prob = 1.0; 00372 double prior_prob1; 00373 double prior_prob2; 00374 for (int i=1; i<=connect_counter; i++) { 00375 prior_prob1 = 1.0; 00376 prior_prob2 = 1.0; 00377 for (int j=1; j<=founder_allele_counter; j++){ 00378 if (connected_groups(j) == i) { 00379 prior_prob1 *= locusFreq(lcs)(allele_vector1(j)); 00380 } 00381 } 00382 for (int j=1; j<=founder_allele_counter; j++){ 00383 if (connected_groups(j) == i) { 00384 if (allele_vector2(j) == -1) { 00385 prior_prob2 = 0.0; 00386 break; 00387 } 00388 prior_prob2 *= locusFreq(lcs)(allele_vector2(j)); 00389 } 00390 } 00391 prior_prob *= prior_prob1 + prior_prob2; 00392 } 00393 if (initial == 0){ 00394 for (int i=0; i<size(); i++) { 00395 if (member(i)->base()){ 00396 member(i)->get_allele_v1(lcs, allele_vector1); 00397 } 00398 } 00399 } 00400 return prior_prob; 00401 }

void matvec::Population::CalcFreqHaploFounders ( unsigned numOfSamples )

Definition at line 2130 of file population.cpp.
References matvec::GeneticDist::chrom(), matvec::ChromStruct::locus, matvec::Individual::mHaplotypeFreq, matvec::Individual::mymother, matvec::Individual::pHaplotypeFreq, popmember, popsize, and prior.
Referenced by matvec::Model::RSamplerGibbs(), matvec::Model::RSamplerGibbsMH(), and matvec::Model::RSamplerMH().

02130 { 02131 // Authors: Helene Gilbert and Rohan L. Fernando 02132 // (Spring 2004) 02133 // Contributors: L. Radu Totir 02134 Individual *ind, *mom; 02135 double doubleSamples = numOfSamples ; 02136 // cout<<" nb samples "<<doubleSamples<<endl; 02137 unsigned nLoci = Individual::numLoci; 02138 for(unsigned i=0;i<popsize;i++){ 02139 ind=popmember[i]; 02140 mom=ind->mymother; 02141 if(!mom){ 02142 // cout<<"Individual "<<ind->myid<<endl; 02143 unsigned ihaplo=0; 02144 for (unsigned i1=0; i1<nLoci-1 ; i1++){ // defined until (n-1)th locus 02145 if(prior->chrom()[0].locus[i1].qtl_ml != 'q' ){ 02146 ind->mHaplotypeFreq[ihaplo] /= doubleSamples ; 02147 ind->pHaplotypeFreq[ihaplo] /= doubleSamples ; 02148 ihaplo++; 02149 } 02150 } 02151 // cout<<ind->mHaplotypeFreq[1]<<" "<<ind->pHaplotypeFreq[1]<<" "<<ind->mHaplotypeFreq[2]<<" "<<ind->pHaplotypeFreq[2]<<" "<<endl; 02152 } 02153 } 02154 }

unsigned matvec::Population::calcTotalNumberOfGenotypes ( unsigned locus )

method to do count total number of possible genotypes across the whole pedigree

Definition at line 1216 of file population.cpp.
References matvec::Individual::genotNodeVector, popmember, and popsize.
Referenced by LGGenotypeElimination().

01216 { 01217 // Authors: Joseph Abraham and Rohan L. Fernando 01218 // (September, 2004) 01219 // Contributors: 01220 unsigned count = 0; 01221 Individual *ind; 01222 for(unsigned i=0;i<popsize;i++){ 01223 ind=popmember[i]; 01224 //ind->displayGenotypes(locus); 01225 unsigned numGenotypes = ind->genotNodeVector[locus].genotypeVector.size(); 01226 if (numGenotypes==0){ 01227 throw exception ("Population::calcTotalNumberOfGenotypes: Incompatible genotypes in pedigree"); 01228 } 01229 count += numGenotypes; 01230 } 01231 return count; 01232 }

void matvec::Population::compute_pdm ( const Individual * ind,

Matrix< double > & pdm

)

Definition at line 738 of file pop_mblup.cpp.
References matvec::LocusStruct::allele_freq, matvec::GeneticDist::chrom(), matvec::Individual::father(), matvec::Individual::genotype(), matvec::ChromStruct::locus, matvec::Individual::mother(), nmarker, matvec::pdm_val(), matvec::pr_osd(), and prior.

00739 { 00740 int j,k1,k2; 00741 unsigned gtsire[2],gtdam[2],gtoff[2]; 00742 double pr_missing,p0,pros; 00743 Matrix<double> pdm_i(2,4); 00744 Individual *father, *mother; 00745 LocusStruct *ML = &(prior->chrom()[0].locus[0]); 00746 00747 father = I->father(); 00748 mother = I->mother(); 00749 // int ori = I->p_origin; 00750 // parental origin is unknow: allele1 and allele2 00751 I->genotype(0,0,gtoff); 00752 for (j=0; j<4; j++) { 00753 pdm[0][j] = 0.0; pdm[1][j] = 0.0; 00754 } 00755 if (father && mother) { 00756 father->genotype(0,0,gtsire); 00757 mother->genotype(0,0,gtdam); 00758 if (pdm_val(gtoff,gtsire,gtdam,pdm,0)) 00759 throw exception(" incompatible marker genotypes related to animal"); 00760 } 00761 else if (father && !mother) { 00762 //********************************************************** 00763 // pdm = sum (pdm_d * Pr(d|s,o)) over d 00764 // Pr(d|s,o) = Pr(o|s,d)*Pr(d) / sum(Pr(o|s,d)*Pr(d)) over d 00765 //*********************************************************** 00766 father->genotype(0,0,gtsire); 00767 for (pros=0.0, k1=0; k1<nmarker; k1++) { 00768 gtdam[0] = k1; 00769 p0 = ML->allele_freq[k1]; 00770 for (k2=0; k2<nmarker; k2++) { 00771 gtdam[1] = k2; 00772 pr_missing = pr_osd(gtoff,gtsire,gtdam,0)*p0*ML->allele_freq[k2]; 00773 if (pr_missing != 0.0) { 00774 pros += pr_missing; 00775 pdm_val(gtoff,gtsire,gtdam,pdm_i,0); 00776 for (j=0; j<4; j++) { 00777 pdm[0][j] += pdm_i[0][j] * pr_missing; 00778 pdm[1][j] += pdm_i[1][j] * pr_missing; 00779 } 00780 } 00781 } 00782 } 00783 if (pros > 0.0) { 00784 for (j=0; j<4; j++) {pdm[0][j] /= pros; pdm[1][j] /= pros; } 00785 for (j=2; j<4; j++) { // the last two columns are set to zero 00786 pdm[0][j] = 0; // because mother is unknown 00787 pdm[1][j] = 0; 00788 } 00789 } 00790 else { 00791 throw exception("incompatible marker genotypes related to animal"); 00792 } 00793 } 00794 else if (!father && mother) { 00795 //*********************************************************** 00796 // pdm = sum (pdm_s * Pr(s|d,o)) over s 00797 // Pr(s|d,o) = Pr(o|s,d)*Pr(s) / sum(Pr(o|s,d)*Pr(s)) over s 00798 //************************************************************ 00799 mother->genotype(0,0,gtdam); 00800 for (pros=0.0, k1=0; k1<nmarker; k1++) { 00801 gtsire[0] = k1; 00802 p0 = ML->allele_freq[k1]; 00803 for (k2=0; k2<nmarker; k2++) { 00804 gtsire[1] = k2; 00805 pr_missing = pr_osd(gtoff,gtsire,gtdam,0)*p0*ML->allele_freq[k2]; 00806 if (pr_missing != 0.0) { 00807 pros += pr_missing; 00808 pdm_val(gtoff,gtsire,gtdam,pdm_i,0); 00809 for (j=0; j<4; j++) { 00810 pdm[0][j] += pdm_i[0][j] * pr_missing; 00811 pdm[1][j] += pdm_i[1][j] * pr_missing; 00812 } 00813 } 00814 } 00815 } 00816 if (pros > 0.0) { 00817 for (j=0; j<4; j++) { 00818 pdm[0][j] /= pros; pdm[1][j] /= pros; 00819 } 00820 for (j=0; j<2; j++) { // the first two columns are set to zero 00821 pdm[0][j] = 0; // because father is unknown 00822 pdm[1][j] = 0; 00823 } 00824 } 00825 else { 00826 throw exception("incompatible marker genotypes related to animal"); 00827 } 00828 } 00829 }

void matvec::Population::cond_genotype_config ( )

Definition at line 494 of file pop_gibbs.cpp.
References matvec::Vector< T >::begin(), matvec::Matrix< T >::begin(), build_pop_gamete(), matvec::Genome::chromosome, matvec::Chromosome::freq(), matvec::Individual::maternal_chrom(), matvec::Genome::nchrom(), numchrom, partial_cdist(), matvec::Individual::paternal_chrom(), pop_gamete, popmember, popsize, matvec::sampling(), and size().
Referenced by matvec::Model::genotype_dist_gibbs(), and matvec::Model::genotypic_value_gibbs().

00495 { 00496 //////////////////////////////////////////////////////////////////// 00497 // a conditional genotype configuration G over a given pedigree, 00498 // given y. 00499 /////////////////////////////////////////////////////////////////// 00500 00501 if (!pop_gamete) build_pop_gamete(); 00502 00503 Individual *I; 00504 unsigned i,j,k,k0,k1,size; 00505 Chromosome *C0, *C1; 00506 00507 unsigned maxg = pop_gamete[0].nchrom(); // ng = maxximum number of gametes 00508 for (i=1; i<numchrom; i++) { // among chromosomes 00509 k = pop_gamete[i].nchrom(); 00510 if (k > maxg) maxg = k; 00511 } 00512 ////////////////////////////////////////////////////// 00513 // get working space ready for partial_cdist() 00514 // freq_mat, gid_mat, cdist_value, cdist_prob 00515 //////////////////////////////////////////////////// 00516 00517 Matrix<double> freq_mat(2,maxg); 00518 Matrix<unsigned> gid_mat(2,maxg); 00519 Matrix<unsigned> cdist_value(2,maxg*maxg); 00520 Vector<double> cdist_prob(maxg*maxg); 00521 00522 for (i=0; i<popsize; i++) { 00523 I = popmember[i]; 00524 C0 = I->paternal_chrom(); 00525 C1 = I->maternal_chrom(); 00526 00527 /////////////////////////////////////////////////////////////// 00528 // first randomly assigned genotypes for each chromosome except 00529 // the first chromosome. 00530 // there is other alternatives such as using calculations of 00531 // f(y|Q1), f(y|Q1,Q2) ..., 00532 /////////////////////////////////////////////////////////////// 00533 for (j=1; j<numchrom; j++) { 00534 maxg = pop_gamete[j].nchrom(); 00535 for (k=0; k<maxg; k++) { 00536 freq_mat[0][k] = pop_gamete[j].chromosome[k].freq(); 00537 } 00538 k0 = sampling(freq_mat[0], maxg); 00539 k1 = sampling(freq_mat[0], maxg); 00540 C0[j] = pop_gamete[j].chromosome[k0]; 00541 C1[j] = pop_gamete[j].chromosome[k1]; 00542 } 00543 for (j=0; j<numchrom; j++) { 00544 size = partial_cdist(I,j,cdist_value.begin(),cdist_prob.begin(),gid_mat.begin(),freq_mat.begin()); 00545 k = sampling(cdist_prob.begin(),size); 00546 k0 = cdist_value[0][k]; // paternal gamete id 00547 k1 = cdist_value[1][k]; // maternal gamete id 00548 C0[j] = pop_gamete[j].chromosome[k0]; 00549 C1[j] = pop_gamete[j].chromosome[k1]; 00550 } 00551 } 00552 }

void matvec::Population::confirm_sex ( const int col_sex )

Definition at line 399 of file population.cpp.
References matvec::Individual::father(), ind_name(), matvec::Individual::mother(), matvec::Individual::myfather, matvec::Individual::mymother, matvec::Individual::mysex, popmember, popsize, matvec::Individual::sex(), and temp.
Referenced by input_ped().

00400 { 00401 ////////////////////////////////////////////////////////////////////////// 00402 // if Pedigree has an explicity column for, 00403 // then it overwrite 'mother' and 'father' columns, i.e. mother 00404 // and father just mean parents, not sexuality 00405 // else 'mother' and 'father' columns make their sexuality 00406 ////////////////////////////////////////////////////////////////////////// 00407 00408 Individual *I,*temp,*mother,*father; 00409 unsigned i; 00410 if (col_sex >= 0) { // sex has been already set, now check 00411 for (i=0; i<popsize; i++) { // to see it is confirmable 00412 I = popmember[i]; 00413 mother = I->mother(); 00414 father = I->father(); 00415 if (mother && father) { 00416 if (mother->sex() != '.' && mother->sex() == father->sex()) { 00417 throw exception(std::string("Population::confirm_sex(): ") + ind_name(i+1)); 00418 } 00419 else if (mother->sex() == 'M') { 00420 temp = mother; 00421 I->mymother = father; 00422 I->myfather = temp; 00423 } 00424 else if (father->sex() == 'F') { 00425 temp = father; 00426 I->myfather = mother; 00427 I->mymother = temp; 00428 } 00429 } 00430 else if (father && !mother) { 00431 if (father->sex() == 'F') { 00432 temp = father; 00433 I->myfather = mother; 00434 I->mymother = temp; 00435 } 00436 } 00437 else if (!father && mother) { 00438 if (mother->sex() == 'M') { 00439 temp = mother; 00440 I->mymother = father; 00441 I->myfather = temp; 00442 } 00443 } 00444 } 00445 } 00446 else { 00447 for (i=0; i<popsize; i++) { 00448 I = popmember[i]; 00449 mother = I->mother(); 00450 father = I->father(); 00451 // if (mother && mother->mysex == 'M') throw exception(std::string("its mother has been father earlier: ")+ind_name(i+1)); 00452 // if (father && father->mysex == 'F') throw exception(std::string("its father has been mother earlier") + ind_name(i+1)); 00453 if (mother) { 00454 if (mother->mysex == 'M') { 00455 if(Pedigree::type==Pedigree::raw){ 00456 throw exception(std::string("its mother has been father earlier: ") + ind_name(i+1)); 00457 } 00458 else { 00459 std::cerr << "mother of individual "<< (i+1); 00460 throw exception(" has been used as father earlier"); 00461 } 00462 } 00463 else { 00464 mother->mysex = 'F'; 00465 } 00466 } 00467 if (father) { 00468 if (father->mysex == 'F') { 00469 if(Pedigree::type==Pedigree::raw){ 00470 throw exception(std::string("its father has been mother earlier: ") + ind_name(i+1)); 00471 } 00472 else { 00473 std::cerr << "father of individual "<< (i+1); 00474 throw exception(" has been used as mother earlier"); 00475 } 00476 } 00477 else { 00478 father->mysex = 'M'; 00479 } 00480 } 00481 } 00482 /////////////////////////////////////////////////////////////////// 00483 // if I has never been parent, then set it as a female. 00484 // actually whether it is set to be male or female does not matter 00485 /////////////////////////////////////////////////////////////////// 00486 for (i=0; i<popsize; i++) { 00487 I = popmember[i]; 00488 if (I->sex() == '.') I->mysex = 'F'; 00489 } 00490 } 00491 }

matvec::Population::copyCandidateToAccepted ( string samplerType )

save candidate state as accepted states; allele origins for accepted states are already in (p)m_gameteOld

Definition at line 1558 of file population.cpp.
References matvec::Individual::genotNodeVector, matvec::Individual::malleleStateNodeVector, matvec::Individual::palleleStateNodeVector, popmember, and popsize.
Referenced by matvec::Model::RSamplerGibbsMH(), and matvec::Model::RSamplerMH().

01558 { 01559 // Authors: L. Radu Totir 01560 // (October, 2003) 01561 // Contributors: 01562 Individual *ind; 01563 for(unsigned i=0;i<popsize;i++){ 01564 ind=popmember[i]; 01565 for(unsigned j=0;j<Individual::numLoci;j++){ 01566 if(samplerType=="genotypic"){ 01567 ind->genotNodeVector[j].acceptedGenotypeState = ind->genotNodeVector[j].candidateGenotypeState; 01568 ind->genotNodeVector[j].acceptedGenotypeVector = ind->genotNodeVector[j].genotypeVector; 01569 } 01570 else if(samplerType=="allelic"){ 01571 ind->malleleStateNodeVector[j].acceptedAlleleState = ind->malleleStateNodeVector[j].candidateAlleleState; 01572 ind->malleleStateNodeVector[j].acceptedAlleleStateVector = ind->malleleStateNodeVector[j].alleleStateVector; 01573 ind->palleleStateNodeVector[j].acceptedAlleleState = ind->palleleStateNodeVector[j].candidateAlleleState; 01574 ind->palleleStateNodeVector[j].acceptedAlleleStateVector = ind->palleleStateNodeVector[j].alleleStateVector; 01575 } 01576 else { 01577 cerr << "Sampler type should be either genotypic or allelic;" << endl; 01578 exit(1); 01579 } 01580 } 01581 } 01582 }

void matvec::Population::copyfrom ( const Population & A )

Definition at line 56 of file population.cpp.
References blupsol, evaluate_method, fdone, gametebase, hashtable, mark_need_remove, markersymbol, maxnamelen, maxpopsize, n_markerLoci, nmarker, nuFamiliesDone, numchrom, numtrait, offs_info_built, pm_storage, pop_gamete, popsize, prior, resize(), spouse_info_built, stdid, and trans_mat.
Referenced by operator=(), and Population().

00057 { 00058 if (this == &A) 00059 return; 00060 unsigned i; 00061 this->resize(A.maxpopsize, A.prior); 00062 prior = A.prior; 00063 fdone = A.fdone; 00064 stdid = A.stdid; 00065 popsize = A.popsize; 00066 numchrom = A.numchrom; 00067 numtrait = A.numtrait; 00068 nmarker = A.nmarker; 00069 maxpopsize = A.maxpopsize; 00070 maxnamelen = A.maxnamelen; 00071 //RLF 00072 n_markerLoci = A.n_markerLoci; 00073 nuFamiliesDone = A.nuFamiliesDone; 00074 //RLF 00075 trans_mat = A.trans_mat; 00076 blupsol = A.blupsol; 00077 hashtable = A.hashtable; 00078 evaluate_method = A.evaluate_method; 00079 mark_need_remove = A.mark_need_remove; 00080 00081 spouse_info_built = 0; 00082 offs_info_built = 0; 00083 00084 popsize = A.popsize; 00085 00086 for (i=0; i<popsize; i++) 00087 pm_storage[i] = A.pm_storage[i]; 00088 for (i=0; i<numchrom; i++) { 00089 pop_gamete[i] = A.pop_gamete[i]; 00090 gametebase[i] = A.gametebase[i]; 00091 } 00092 00093 if (nmarker != A.nmarker) { 00094 if(markersymbol){ 00095 delete [] markersymbol; 00096 markersymbol=0; 00097 } 00098 if(A.nmarker>0){ 00099 markersymbol = new std::string [A.nmarker]; 00100 } 00101 else { 00102 markersymbol = 0; 00103 } 00104 nmarker = A.nmarker; 00105 } 00106 for (i=0; i<nmarker; i++) markersymbol[i] = A.markersymbol[i]; 00107 }

matvec::Population::copyGNodeStatesToCandidateStates ( string samplerType )

saves the GNode states (alleles or genotypes) for candidate sample

Definition at line 1587 of file population.cpp.
References matvec::Individual::genotNodeVector, matvec::Individual::malleleStateNodeVector, matvec::Individual::palleleStateNodeVector, popmember, and popsize.
Referenced by matvec::Model::RSamplerGibbsMH(), and matvec::Model::RSamplerMH().

01587 { 01588 // Authors: L. Radu Totir 01589 // (October, 2003) 01590 // Contributors: 01591 Individual *ind; 01592 for(unsigned i=0;i<popsize;i++){ 01593 ind=popmember[i]; 01594 for(unsigned j=0;j<Individual::numLoci;j++){ 01595 if(samplerType=="genotypic"){ 01596 ind->genotNodeVector[j].candidateGenotypeState = ind->genotNodeVector[j].genotypeState; 01597 } 01598 else if(samplerType=="allelic"){ 01599 ind->malleleStateNodeVector[j].candidateAlleleState = ind->malleleStateNodeVector[j].alleleState; 01600 ind->palleleStateNodeVector[j].candidateAlleleState = ind->palleleStateNodeVector[j].alleleState; 01601 } 01602 else { 01603 cerr << "Sampler type should be either genotypic or allelic;" << endl; 01604 exit(1); 01605 } 01606 } 01607 } 01608 }

void matvec::Population::count_genotype ( Individual * I )

Definition at line 53 of file pop_gibbs.cpp.
References matvec::Individual::genotype_counter, matvec::Individual::maternal_chrom_id(), numchrom, and matvec::Individual::paternal_chrom_id().
Referenced by gibbs_iterate().

00054 { 00055 unsigned i,j,k,t; 00056 for (t=0; t<numchrom; t++) { 00057 i = I->paternal_chrom_id(t); 00058 j = I->maternal_chrom_id(t); 00059 if (i >= j) k = i*(i+1)/2+j; 00060 else k = j*(j+1)/2 + i; 00061 I->genotype_counter[t][k] += 1.0; 00062 } 00063 }

void matvec::Population::countFounders ( void )

Definition at line 493 of file population.cpp.
References matvec::Individual::mymother, numFounders, popmember, and popsize.
Referenced by input_ped().

00493 { 00494 // Authors: L. Radu Totir 00495 // (September, 2004) 00496 // Contributors: 00497 Individual *ind; 00498 numFounders=0; 00499 for(int i=0;i<popsize;i++){ 00500 ind=popmember[i]; 00501 if(!(ind->mymother)){ 00502 numFounders++; 00503 } 00504 } 00505 }

void matvec::Population::countGenotypes ( string samplerType )

Definition at line 2327 of file population.cpp.
References matvec::GeneticDist::chromosome, matvec::Individual::genotNodeVector, matvec::ChromStruct::locus, matvec::Individual::malleleStateNodeVector, matvec::Individual::palleleStateNodeVector, popmember, popsize, and prior.
Referenced by matvec::Model::RSamplerMH().

02327 { 02328 // Authors: Rohan L. Fernando 02329 // (October, 2004) 02330 // Contributors: 02331 02332 unsigned nLoci = Individual::numLoci; 02333 unsigned mat, pat, geno; 02334 Individual *ind; 02335 for(unsigned i=0;i<popsize;i++){ 02336 ind=popmember[i]; 02337 for (unsigned j=0; j<nLoci; j++){ 02338 if (samplerType=="genotypic"){ 02339 mat = ind->genotNodeVector[j].getAcceptedMatState(); 02340 pat = ind->genotNodeVector[j].getAcceptedPatState(); 02341 } 02342 else if (samplerType=="allelic"){ 02343 mat = ind->malleleStateNodeVector[j].getAcceptedAlleleState(); 02344 pat = ind->palleleStateNodeVector[j].getAcceptedAlleleState(); 02345 } 02346 geno = mat*prior->chromosome[0].locus[j].nallele() + pat; 02347 ind->genotNodeVector[j].genotypeCount[geno]++; 02348 } 02349 } 02350 }

void matvec::Population::countHaplotypes ( string samplerType )

Definition at line 2261 of file population.cpp.
References matvec::Individual::genotNodeVector, matvec::IndexVector::getIndex(), haplotypeCoder, matvec::Individual::malleleStateNodeVector, matvec::Individual::matHaplotypeCount, matvec::Individual::palleleStateNodeVector, matvec::Individual::patHaplotypeCount, popmember, and popsize.
Referenced by matvec::Model::RSamplerGibbs(), matvec::Model::RSamplerGibbsMH(), and matvec::Model::RSamplerMH().

02261 { 02262 // Authors: L. Radu Totir and Rohan L. Fernando 02263 // (August, 2004) 02264 // Contributors: 02265 std::vector<unsigned> matHaplotype, patHaplotype; 02266 unsigned nLoci = Individual::numLoci; 02267 Individual *ind; 02268 for(unsigned i=0;i<popsize;i++){ 02269 ind=popmember[i]; 02270 matHaplotype.resize(nLoci); 02271 patHaplotype.resize(nLoci); 02272 for (unsigned j=0; j<nLoci; j++){ 02273 if (samplerType=="genotypic"){ 02274 matHaplotype[j] = ind->genotNodeVector[j].getAcceptedMatState(); 02275 patHaplotype[j] = ind->genotNodeVector[j].getAcceptedPatState(); 02276 } 02277 else if (samplerType=="allelic"){ 02278 matHaplotype[j] = ind->malleleStateNodeVector[j].getAcceptedAlleleState(); 02279 patHaplotype[j] = ind->palleleStateNodeVector[j].getAcceptedAlleleState(); 02280 } 02281 } 02282 unsigned matIndex = haplotypeCoder.getIndex(matHaplotype); 02283 unsigned patIndex = haplotypeCoder.getIndex(patHaplotype); 02284 ind->matHaplotypeCount[matIndex]++; 02285 ind->patHaplotypeCount[patIndex]++; 02286 } 02287 }

double matvec::Population::cprob_children ( const Individual * I,

const unsigned jc,

unsigned sg[],

unsigned dg[],

double sgf[],

double dgf[]

)

Definition at line 235 of file pop_gibbs.cpp.
References matvec::Individual::father(), ind_gamete(), matvec::Individual::maternal_chrom_id(), matvec::Individual::mother(), matvec::Individual::myoffspring, matvec::Individual::noffs(), matvec::Individual::nspouse(), matvec::Individual::numoffs_spouse, matvec::Individual::offspring(), matvec::Individual::paternal_chrom_id(), and matvec::Individual::sex().
Referenced by full_cdist().

00237 { 00238 //////////////////////////////////////////////////// 00239 // note that if no children, return 1.0 00240 ////////////////////////////////////////////////// 00241 00242 unsigned i,j,t,k,nc,ns,gid,ssize,dsize,noffs; 00243 Individual *child, **offspring; 00244 ns = I->nspouse(); 00245 double prob,tmp1,tmp2; 00246 noffs = I->noffs(); 00247 offspring = I->offspring(); 00248 00249 prob=1.0; 00250 if (I->sex()== 'F') { // I is a mother 00251 dsize = ind_gamete(I, jc,dg,dgf); 00252 for (k=0; k<noffs; k++) { 00253 child = offspring[k]; 00254 gid = child->maternal_chrom_id(jc); 00255 for (tmp1=0.0, t=0; t < dsize; t++) { 00256 if (gid == dg[t]) { 00257 tmp1 = dgf[t]; 00258 break; 00259 } 00260 } 00261 if (tmp1 == 0.0) return tmp1; 00262 prob *= tmp1; 00263 } 00264 for (k=0, i=0; i<ns; i++) { 00265 ssize = ind_gamete(I->myoffspring[k]->father(), jc, sg,sgf); 00266 nc = I->numoffs_spouse[i]; 00267 for (j=0; j<nc; j++) { 00268 child = offspring[k++]; 00269 gid = child->paternal_chrom_id(jc); 00270 for (tmp2=0.0, t=0; t < ssize; t++) { 00271 if (gid == sg[t]) { 00272 tmp2 = sgf[t]; 00273 break; 00274 } 00275 } 00276 if (tmp2 == 0.0) return tmp2; 00277 prob *= tmp2; 00278 } 00279 } 00280 } 00281 else { // I is a father 00282 ssize = ind_gamete(I, jc,sg,sgf); 00283 for (k=0; k<noffs; k++) { 00284 child = offspring[k]; 00285 gid = child->paternal_chrom_id(jc); 00286 for (tmp1=0.0, t=0; t < ssize; t++) { 00287 if (gid == sg[t]) { 00288 tmp1 = sgf[t]; 00289 break; 00290 } 00291 } 00292 if (tmp1 == 0.0) return tmp1; 00293 prob *= tmp1; 00294 } 00295 00296 for (k=0, i=0; i<ns; i++) { 00297 dsize = ind_gamete(I->myoffspring[k]->mother(), jc, dg,dgf); 00298 nc = I->numoffs_spouse[i]; 00299 for (j=0; j<nc; j++) { 00300 child = offspring[k++]; 00301 gid = child->maternal_chrom_id(jc); 00302 for (tmp2=0.0, t=0; t < dsize; t++) { 00303 if (gid == dg[t]) { 00304 tmp2 = dgf[t]; 00305 break; 00306 } 00307 } 00308 if (tmp2 == 0.0) return tmp2; 00309 prob *= tmp2; 00310 } 00311 } 00312 } 00313 return prob; 00314 }

double matvec::Population::cprob_op ( const unsigned goff[],

const unsigned gparent[],

const int ori

)

Definition at line 870 of file pop_mblup.cpp.
References matvec::pr_osd().

00875 { 00876 LocusStruct *ML = &(prior->chrom()[0].locus[0]); 00877 unsigned gtunknown[2]; 00878 double pros,p0; 00879 int k1,k2; 00880 for (pros=0.0, k1=0; k1<nmarker; k1++) { 00881 gtunknown[0] = k1; 00882 p0 = ML->allele_freq[k1]; 00883 for (k2=0; k2<nmarker; k2++) { 00884 gtunknown[1] = k2; 00885 pros += pr_osd(goff,gparent,gtunknown,ori)*p0*ML->allele_freq[k2]; 00886 } 00887 } 00888 return pros; 00889 }

Population * matvec::Population::descendant_of ( const char name[] )

Definition at line 224 of file population1.cpp.
References descendant_of(), and get_id().

00225 { 00226 unsigned id = get_id(name); 00227 if (id == 0) throw exception(" Population::descendant_of(): cannot found"); 00228 return descendant_of(id); 00229 }

Population * matvec::Population::descendant_of ( const unsigned id )

put all descendants of individual id into a new population. caller is resposible to release this new population using either delete or resize(0,0,0)
Definition at line 123 of file population1.cpp.
References mark_descendant_of(), popsize, renum(), stdid, and sub().
Referenced by descendant_of().

00124 { 00125 if (!stdid) renum(); 00126 if (id < 1 || id > popsize ) throw exception(" Population::mark_descendant_of(): bad arg"); 00127 unsigned counter = 1; 00128 unsigned nd = mark_descendant_of(id,counter); 00129 return sub(nd); 00130 }

void matvec::Population::descent_graph_init_parm ( void )

Definition at line 33 of file pop_graph.cpp.
References matvec::GeneticDist::chrom(), matvec::ChromStruct::locus, locusFreq, prior, and matvec::Vector< Vector< double > >::resize().
Referenced by descent_graph_setup().

00033 { 00034 // Authors: Mathias Schelling and Rohan L. Fernando 00035 // (1999) 00036 // Contributors: L. Radu Totir 00037 unsigned nLoci = Individual::numLoci; 00038 locusFreq.resize(nLoci); 00039 // Reading in frequencies for locus i into vector locusFreq 00040 for (unsigned i=0;i<nLoci; i++){ 00041 locusFreq[i] = prior->chrom()[0].locus[i].allele_freq; 00042 } 00043 }

double matvec::Population::descent_graph_log_lhood ( )

Definition at line 403 of file pop_graph.cpp.
References calc_prior_descent_graph(), matvec::GeneticDist::chrom(), matvec::ChromStruct::locus, matvec::Individual::m_gamete, member(), matvec::Individual::mymother, matvec::Individual::p_gamete, prior, RecoVector, size(), and temp.
Referenced by M_H_sample_ext().

00403 { 00404 // Authors: Mathias Schelling and Rohan L. Fernando 00405 // (1999) 00406 // Contributors: Fabiano V. Pita, L. Radu Totir 00407 //calc prior probability: 00408 double prior_prob = 0.0; 00409 unsigned nLoci = Individual::numLoci; 00410 for (int i=1; i<=nLoci; i++) { 00411 double temp = 1.0; 00412 if(prior->chrom()[0].locus[i-1].qtl_ml!='q'){ 00413 temp = calc_prior_descent_graph(i); 00414 } 00415 // std::cout << "prior marker locus " << i << ": " << temp << std::endl; 00416 if (temp == 0.0) { 00417 //i.e. graph is not valid, return inpossible log lh and return 00418 return 9999.9; 00419 } 00420 prior_prob += std::log(temp); 00421 } 00422 // std::cout << "log prior prob: " << prior_prob<< std::endl; 00423 //calc transmission probability 00424 double trans_prob = 0.0; 00425 for (int i=1; i<size(); i++) { 00426 if (member(i)->mymother) { // if not founder 00427 double temp = 0.25; 00428 for (int j=2; j<=nLoci; j++) { 00429 if (member(i)->m_gamete(j-1) == member(i)->m_gamete(j)) { 00430 temp *= 1 - RecoVector(j-1); 00431 } 00432 else { 00433 temp *= RecoVector(j-1); 00434 } 00435 if (member(i)->p_gamete(j-1) == member(i)->p_gamete(j)) { 00436 temp *= 1 - RecoVector(j-1); 00437 } 00438 else { 00439 temp *= RecoVector(j-1); 00440 } 00441 } 00442 trans_prob += std::log(temp); 00443 } 00444 } 00445 // std::cout << "log transmission prob: " << trans_prob<< std::endl; 00446 //return log lh 00447 return prior_prob + trans_prob; 00448 }

void matvec::Population::descent_graph_setup ( Data * )

Definition at line 634 of file pop_graph.cpp.
References descent_graph_init_parm(), matvec::GeneticDist::get_distance(), input_markerData(), matvec::Model::MapF(), model, prior, RecoVector, and matvec::Vector< double >::resize().
Referenced by matvec::Model::DGSamplerSetup().

00634 { 00635 // Authors: Mathias Schelling and Rohan L. Fernando 00636 // (1999) 00637 // Contributors: Fabiano V. Pita, L. Radu Totir 00638 00639 // set up all the necessary tables for descent graph calculations 00640 input_markerData(data); 00641 descent_graph_init_parm(); 00642 int nLociInt = Individual::numLoci-1; 00643 RecoVector.resize(nLociInt); 00644 for (int i=1; i<=nLociInt; i++) { 00645 RecoVector(i) = model->MapF(prior->get_distance(1,i+1) - prior->get_distance(1,i)); 00646 } 00647 std::cout << "RecoVector:\n" << RecoVector; 00648 }

int matvec::Population::detect_loop ( void )

Definition at line 297 of file pop_peeling.cpp.
References matvec::compare_seq_id(), matvec::NuFamily::connectors, matvec::NuFamily::kernal, matvec::NuFamily::nconnector(), nufamily_vec, num_nufamily, num_t_nufamily, matvec::NuFamily::seq_id, and matvec::NuFamily::terminal.
Referenced by peeling_sequence().

00297 { 00298 ////////////////////////////////////////////////////////////////////////// 00299 // if an individual's group_id > 1, then it is a connector. 00300 // we count # of connectors within each family, each family can have 00301 // more than 1 connectors, however, a terminal family only have 00302 // 1 connector. If a family has no connector at all, it is the 00303 // kernal family. an isolated individial is not a nuclei family 00304 // if we can't clip out terminal to the end of the pedigree, that means 00305 // there is at least one loop in the rest part of the pedigree. 00306 // 00307 // REQUIREMENTS: build_nufamily() has to be called; 00308 ////////////////////////////////////////////////////////////////////////// 00309 00310 unsigned i,nc,keep_going = 2; 00311 std::list<Individual *>::iterator it; 00312 Individual **connector; 00313 NuFamily *family; 00314 unsigned i0 = num_t_nufamily; 00315 while (keep_going) { 00316 keep_going = 2; // it remains 2 until it's changed 00317 for (i=i0; i<num_nufamily; i++) { 00318 family = nufamily_vec[i]; 00319 if ( !family->terminal) { 00320 nc = family->nconnector(); 00321 if (nc == 0) { // found a kernal nuclei family 00322 family->seq_id = num_t_nufamily++; 00323 family->kernal = 1; 00324 family->terminal = 1; 00325 keep_going = 1; // reset to 1 00326 } 00327 else if (nc == 1) { // found a terminal nuclei family 00328 family->seq_id = num_t_nufamily++; 00329 it = family->connectors.begin(); 00330 (*it)->related_family.remove(family); 00331 (*it)->group_id -= 1; 00332 family->terminal = 1; 00333 keep_going = 1; // reset to 1 00334 } 00335 else { 00336 family->seq_id = num_nufamily+nc; // nothing found 00337 } 00338 if (num_t_nufamily == num_nufamily) {keep_going = 0; break;} 00339 } 00340 } 00341 if (keep_going == 2) break; // can't clip farther, loops may be found 00342 } 00343 if (num_t_nufamily > i0 ) { 00344 qsort((char *)nufamily_vec,num_nufamily,sizeof(NuFamily *),compare_seq_id); 00345 } 00346 return (num_nufamily-num_t_nufamily); 00347 }

void matvec::Population::display ( const char key[] )

Definition at line 808 of file population.cpp.
References matvec::Locus::allele, matvec::Vector< T >::begin(), matvec::Genome::chromosome, matvec::Individual::father_id(), matvec::Individual::genotype_counter, matvec::Individual::id(), ind_name(), matvec::Chromosome::locus, matvec::Individual::mother_id(), matvec::Chromosome::nloci(), numchrom, pop_gamete, popmember, popsize, release_genotype_counter(), matvec::Individual::sex(), and matvec::Genome::size().

00809 { 00810 Individual *I; 00811 unsigned i,t,t1,t2,t3,k,nl,num; 00812 if (strcmp(key,"ped") == 0) { 00813 for (i=0; i<popsize; i++) { 00814 I = popmember[i]; 00815 std::cout << ind_name(i+1) << " " << I->id() << " " << I->mother_id(); 00816 std::cout << " " << I->father_id() << " " << I->sex() << "\n"; 00817 } 00818 } 00819 else if(strcmp(key,"genotype") == 0) { 00820 if (!(popmember[0]->genotype_counter)) return; 00821 Chromosome *C1, *C2; 00822 Vector<double> *vv; 00823 for (i=0; i<popsize; i++) { 00824 I = popmember[i]; 00825 std::cout << "individual " << ind_name(i+1) << ":\n"; 00826 for (t=0; t<numchrom; t++) { 00827 vv = &(I->genotype_counter[t]); 00828 std::cout << " chromosome " << t+1 << ":\n"; 00829 num = pop_gamete[t].size(); 00830 nl = pop_gamete[t].chromosome[0].nloci(); 00831 for (k=0,t1=0; t1<num; t1++) { 00832 C1 = &(pop_gamete[t].chromosome[t1]); 00833 for (t2=0; t2<=t1; t2++) { 00834 C2 = &(pop_gamete[t].chromosome[t2]); 00835 std::cout << " ___ "; 00836 for (t3=0; t3<nl; t3++) std::cout << C1->locus[t3].allele << " "; 00837 std::cout << "___\n"; 00838 std::cout << " ~~~ "; 00839 for (t3=0; t3<nl; t3++) std::cout << C2->locus[t3].allele << " "; 00840 std::cout << "~~~: "<< static_cast<unsigned>(vv->begin()[k]) << "\n"; 00841 k++; // k = t1*(t1+1)/2 + t2; 00842 } 00843 } 00844 } 00845 } 00846 } 00847 std::cout << "\n"; 00848 release_genotype_counter(); 00849 }

void matvec::Population::DisplayFreqHaploFounders ( void )

Definition at line 2156 of file population.cpp.
References matvec::GeneticDist::chrom(), ind_name(), matvec::ChromStruct::locus, matvec::Individual::mHaplotypeFreq, model, matvec::Individual::myid, matvec::Individual::mymother, matvec::Model::myRSamplerParms, matvec::Individual::pHaplotypeFreq, popmember, popsize, prior, and matvec::RSamplerParms::resultsFile.
Referenced by matvec::Model::RSamplerGibbs(), matvec::Model::RSamplerGibbsMH(), and matvec::Model::RSamplerMH().

02156 { 02157 // Authors: Helene Gilbert and Rohan L. Fernando 02158 // (Spring 2004) 02159 // Contributors: L. Radu Totir 02160 const char *s = model->myRSamplerParms.resultsFile.c_str(); 02161 ofstream outfile; 02162 outfile.open(s); 02163 // cout<<endl<<" Founders'haplotype Frequencies "<<endl; 02164 Individual *ind, *mom; 02165 unsigned nLoci = Individual::numLoci; 02166 for (unsigned i1=0; i1<nLoci-1 ; i1++){ // defined until (n-1)th locus 02167 if(prior->chrom()[0].locus[i1].qtl_ml != 'q' ){ 02168 outfile << "Number of possible haplotypes at locus " << i1+1 << " = " 02169 << prior->chrom()[0].locus[i1].nHaplotypes << " => "; 02170 for (unsigned ih = 0; ih < prior->chrom()[0].locus[i1].nHaplotypes; ih++){ 02171 outfile << "(" << prior->chrom()[0].locus[i1].al1Haplo[ih] 02172 << "," << prior->chrom()[0].locus[i1].al2Haplo[ih] << ") "; 02173 } 02174 outfile << endl; 02175 } 02176 } 02177 outfile << endl; 02178 for(unsigned i=0;i<popsize;i++){ 02179 ind=popmember[i]; 02180 mom=ind->mymother; 02181 if(!mom){ 02182 if (Pedigree::type==Pedigree::raw){ 02183 outfile << "Individual " << ind_name(i+1) << endl; 02184 } 02185 else { 02186 outfile << "Individual " << ind->myid << endl; 02187 } 02188 unsigned ihaplo=0; 02189 for (unsigned i1=0; i1<nLoci-1 ; i1++){ // defined until (n-1)th locus 02190 outfile << "For interval " << i1+1 << " => " << endl; 02191 if(prior->chrom()[0].locus[i1].qtl_ml != 'q' ){ 02192 for (unsigned ih = 0; ih < prior->chrom()[0].locus[i1].nHaplotypes; ih++){ 02193 if(ind->mHaplotypeFreq[ihaplo][ih]||ind->pHaplotypeFreq[ihaplo][ih]){ 02194 outfile << "(" << prior->chrom()[0].locus[i1].al1Haplo[ih] 02195 << "," << prior->chrom()[0].locus[i1].al2Haplo[ih] << ") => "; 02196 outfile << "maternal = " << ind->mHaplotypeFreq[ihaplo][ih] << " " 02197 << "paternal = " << ind->pHaplotypeFreq[ihaplo][ih] << endl; 02198 } 02199 } 02200 ihaplo++; 02201 } 02202 } 02203 } 02204 } 02205 }

void matvec::Population::displayGenotypeFrequencies ( unsigned numSamples,

std::ostream & outfile = cout

)

Definition at line 2357 of file population.cpp.
References matvec::GeneticDist::chromosome, matvec::Individual::genotNodeVector, ind_name(), matvec::ChromStruct::locus, model, matvec::Individual::myid, matvec::Model::myRSamplerParms, popmember, popsize, prior, and matvec::RSamplerParms::resultsFile.
Referenced by matvec::Model::RSamplerMH().

02357 { 02358 // Authors: Rohan L. Fernando 02359 // (October, 2004) 02360 // Contributors: 02361 const char *s = model->myRSamplerParms.resultsFile.c_str(); 02362 Individual *ind; 02363 unsigned nLoci = Individual::numLoci; 02364 double denominator = 1.0/numSamples; 02365 for (unsigned j=0;j<nLoci;j++){ 02366 unsigned numAlleles = prior->chromosome[0].locus[j].nallele(); 02367 outfile << "Locus: " << j+1 << endl; 02368 for(unsigned i=0;i<popsize;i++){ 02369 ind=popmember[i]; 02370 if(Pedigree::type==Pedigree::raw){ 02371 outfile << "Genotype (maternal/paternal) for: " << setw(5) << ind_name(i+1) << endl; 02372 } 02373 else { 02374 outfile << "Genotype (maternal/paternal) for: " << setw(5) << ind->myid << endl; 02375 } 02376 outfile << " ---------------------------- " << endl; 02377 for (unsigned mat=0;mat<numAlleles;mat++){ 02378 for (unsigned pat=0;pat<numAlleles;pat++){ 02379 unsigned geno = mat*prior->chromosome[0].locus[j].nallele() + pat; 02380 if (ind->genotNodeVector[j].genotypeCount[geno]){ 02381 outfile <<setw(2) << mat+1 << "/" <<setw(2) << pat+1 << " " 02382 <<setw(10) << ind->genotNodeVector[j].genotypeCount[geno]*denominator << endl; 02383 } 02384 } 02385 } 02386 outfile << " ---------------------------- " << endl; 02387 } 02388 } 02389 }

void matvec::Population::displayHaplotypeFrequencies ( unsigned numSamples,

std::ostream & os = cout

)

Definition at line 2292 of file population.cpp.
References matvec::IndexVector::getVector(), haplotypeCoder, ind_name(), matvec::Individual::matHaplotypeCount, matvec::Individual::myid, numHaplotypes, matvec::Individual::patHaplotypeCount, popmember, and popsize.
Referenced by matvec::Model::RSamplerGibbs(), and matvec::Model::RSamplerMH().

02292 { 02293 // Authors: L. Radu Totir and Rohan L. Fernando 02294 // (August, 2004) 02295 // Contributors: 02296 02297 unsigned nLoci = Individual::numLoci; 02298 Individual *ind; 02299 double numerator = 1.0/(1.0*numSamples); 02300 for(unsigned i=0;i<popsize;i++){ 02301 ind=popmember[i]; 02302 if(Pedigree::type==Pedigree::raw){ 02303 os << "Haplotype for:" << setw(5) << ind_name(i+1); 02304 } 02305 else { 02306 os << "Haplotype for:" << setw(5) << ind->myid; 02307 } 02308 os << setw(10) << "maternal" << setw(10) << "paternal" << endl; 02309 os << " ---------------------------- " << endl; 02310 for(unsigned j=0;j<numHaplotypes;j++){ 02311 if (ind->matHaplotypeCount[j]!=0||ind->patHaplotypeCount[j]!=0){ 02312 std::vector<unsigned> outVector = haplotypeCoder.getVector(j); 02313 for (unsigned k=0;k<nLoci;k++){ 02314 os << setw(4) << outVector[k]+1; 02315 } 02316 os << setw(10) << ind->matHaplotypeCount[j]*numerator 02317 << setw(10) << ind->patHaplotypeCount[j]*numerator << endl; 02318 } 02319 } 02320 os << " ---------------------------- " << endl; 02321 } 02322 }

void matvec::Population::displaySegregationIndicators ( std::ostream & outfile = cout )

Definition at line 1714 of file population.cpp.
References ind_name(), matvec::Individual::m_gamete, matvec::Individual::myid, matvec::Individual::mymother, matvec::Individual::p_gamete, popmember, and popsize.
Referenced by matvec::Model::RSamplerMH().

01714 { 01715 // Authors: Rohan L. Fernando 01716 // (November, 2004) 01717 // Contributors: 01718 for(unsigned i=0;i<popsize;i++){ 01719 Individual *ind=popmember[i]; 01720 if(Pedigree::type==Pedigree::raw){ 01721 outfile << setw(5) << ind_name(i+1) << " "; 01722 } 01723 else { 01724 outfile << setw(5) << ind->myid << " "; 01725 } 01726 Individual *mom=ind->mymother; 01727 if (mom){ 01728 for (unsigned locus = 0; locus < Individual::numLoci; locus++){ 01729 outfile << ind->m_gamete[locus] << " " << ind->p_gamete[locus] << " "; 01730 } 01731 } 01732 else { 01733 for (unsigned locus = 0; locus < Individual::numLoci; locus++){ 01734 outfile << " 9 9 "; 01735 } 01736 } 01737 outfile << endl; 01738 } 01739 }

Population * matvec::Population::family_of ( const char name[] )

Definition at line 238 of file population1.cpp.
References get_id(), and relative_of().

00239 { 00240 unsigned id = get_id(name); 00241 if (id == 0) throw exception("Population::family_of(): cannot found"); 00242 return relative_of(id); 00243 }

Population* matvec::Population::family_of ( const unsigned id ) [inline]

Definition at line 189 of file population.h.
References relative_of().

00189 {return relative_of(id);}

unsigned matvec::Population::fetch_families ( std::ostream & stream )

Definition at line 254 of file population1.cpp.
References matvec::DataNode::double_val(), matvec::Individual::father(), matvec::Individual::group_id, matvec::Individual::id(), ind_name(), mark_need_remove, mark_relative_recursive(), maxnamelen, matvec::DataNode::missing, matvec::Individual::mother(), numtrait, popmember, popsize, matvec::Individual::record(), remove_mark(), renum(), and stdid.

00255 { 00256 if (!stdid) renum(); 00257 if (popsize == 0) return 0; 00258 if (mark_need_remove) remove_mark(); 00259 00260 Individual *I; 00261 int k; 00262 unsigned i,familyid = 1; 00263 for (i=0; i<popsize; i++) { 00264 I = popmember[i]; 00265 if (I->group_id == 0) { 00266 mark_relative_recursive(I,familyid++); 00267 } 00268 } 00269 familyid--; 00270 00271 unsigned nf; 00272 unsigned fd = 1; 00273 while (fd <= familyid) { 00274 for (nf=0,i=0; i<popsize; i++) if (popmember[i]->group_id == fd) nf++; 00275 stream << "family " << fd << ": " << nf << "\n"; 00276 for (i=0; i<popsize; i++) { 00277 I = popmember[i]; 00278 if (I->group_id != fd) continue; 00279 if (maxnamelen == 0) { 00280 stream << std::setw(10) << I->id(); 00281 if (I->mother()) { 00282 stream << " " << std::setw(10) << I->mother()->id(); 00283 } 00284 else { 00285 stream << " " << std::setw(10) << "."; // TW 00286 } 00287 if (I->father()) { 00288 stream << " " << std::setw(10) << I->father()->id(); 00289 } 00290 else { 00291 stream << " " << std::setw(10) << "."; // TW 00292 } 00293 } 00294 else { 00295 stream << std::setw(maxnamelen) << ind_name(I->id()); 00296 if (I->mother()) { 00297 stream << " " << std::setw(maxnamelen)<<ind_name(I->mother()->id()); 00298 } 00299 else { 00300 stream << " " << std::setw(maxnamelen) << "."; // TW 00301 } 00302 if (I->father()) { 00303 stream << " " << std::setw(maxnamelen)<<ind_name(I->father()->id()); 00304 } 00305 else { 00306 stream << " " << std::setw(maxnamelen) << "."; // TW 00307 } 00308 } 00309 for (k=0; k<numtrait; k++) { 00310 if (I->record()[k].missing) { 00311 stream << " ."; // TW 00312 } 00313 else { 00314 stream << " " << I->record()[k].double_val(); 00315 } 00316 } 00317 stream << "\n"; 00318 } 00319 fd++; 00320 } 00321 mark_need_remove = 1; 00322 return familyid; 00323 }

unsigned matvec::Population::fetch_families ( const char filename[] )

Definition at line 245 of file population1.cpp.

00246 { 00247 std::ofstream outfile(filename); 00248 if (!outfile) throw exception("Population::fetch_families(): cannot open or already exists"); 00249 unsigned retval = this->fetch_families(outfile); 00250 outfile.close(); 00251 return retval; 00252 }

void matvec::Population::finishUpAlleleStateInit ( unsigned locus )

Definition at line 1256 of file population.cpp.
References matvec::Locus::allele, matvec::Genome::chromosome, matvec::Individual::genome0, matvec::Individual::genome1, matvec::Individual::genotNodeVector, matvec::Chromosome::locus, matvec::Individual::malleleStateNodeVector, matvec::MaternalPaternalAlleles::maternal, matvec::Individual::palleleStateNodeVector, matvec::MaternalPaternalAlleles::paternal, popmember, and popsize.
Referenced by SimpleGenotypeElimination().

01256 { 01257 // Authors: L. Radu Totir and Rohan L. Fernando 01258 // (June, 2003) 01259 // Contributors: 01260 Individual *ind; 01261 MaternalPaternalAlleles matPat; 01262 unsigned j = locus; 01263 for(unsigned i=0;i<popsize;i++) { 01264 ind=popmember[i]; 01265 ind->genotNodeVector[j].genotypeVector.clear(); 01266 unsigned allelePat = ind->genome0.chromosome[0].locus[j].allele; 01267 unsigned alleleMat = ind->genome1.chromosome[0].locus[j].allele; 01268 unsigned n = ind->malleleStateNodeVector[j].alleleStateVector.size(); 01269 for(unsigned k=0;k<n;k++){ 01270 unsigned allelek = ind->malleleStateNodeVector[j].alleleStateVector[k]; 01271 ind->palleleStateNodeVector[j].alleleStateVector.push_back(allelek); 01272 matPat.maternal = allelek-1; 01273 for(unsigned l=0;l<n;l++){ 01274 unsigned allelel = ind->malleleStateNodeVector[j].alleleStateVector[l]; 01275 matPat.paternal = allelel-1; 01276 // both alleles at a marker assumed present or absent 01277 if(allelePat!=0 && allelePat!=alleleMat && matPat.paternal==matPat.maternal){ 01278 continue; 01279 } 01280 ind->genotNodeVector[j].genotypeVector.push_back(matPat); 01281 } 01282 } 01283 } 01284 }

unsigned matvec::Population::full_cdist ( Individual * I,

const unsigned jc,

unsigned ** cdist_value,

double * cdist_prob,

unsigned ** gid_mat,

double ** freq_mat

)

Definition at line 350 of file pop_gibbs.cpp.
References matvec::Matrix< double >::begin(), matvec::Genome::chromosome, cprob_children(), matvec::Individual::father(), matvec::Individual::genome0, matvec::Individual::genome1, ind_gamete(), matvec::Individual::mother(), matvec::Individual::noffs(), pop_gamete, and residual_var.
Referenced by gibbs_iterate().

00353 { 00354 //////////////////////////////////////////////////////////////////// 00355 // working space gid_mat,freq_mat must have at least four rows 00356 //////////////////////////////////////////////////////////////////// 00357 00358 unsigned *s_gamete_id = gid_mat[0]; 00359 unsigned *d_gamete_id = gid_mat[1]; 00360 double *s_gamete_freq = freq_mat[0]; 00361 double *d_gamete_freq = freq_mat[1]; 00362 unsigned ss = ind_gamete(I->father(), jc, s_gamete_id, s_gamete_freq); 00363 unsigned ds = ind_gamete(I->mother(), jc, d_gamete_id, d_gamete_freq); 00364 00365 unsigned i,j,k,sg,dg; 00366 double prob; 00367 for (k=0,i=0; i<ss; i++) for (j=0; j<ds; j++) { 00368 sg = s_gamete_id[i]; 00369 dg = d_gamete_id[j]; 00370 I->genome0.chromosome[jc] = pop_gamete[jc].chromosome[sg]; 00371 I->genome1.chromosome[jc] = pop_gamete[jc].chromosome[dg]; 00372 prob = 1.0; 00373 if (I->noffs() > 0) { 00374 prob = cprob_children(I,jc,gid_mat[2],gid_mat[3],freq_mat[2], 00375 freq_mat[3]); 00376 } 00377 if (prob != 0.0) { 00378 prob *= (*penetrance_f)(I,(const double**)residual_var.begin()); 00379 prob *= s_gamete_freq[i]*d_gamete_freq[j]; 00380 if (prob > 1.0e-10) { 00381 cdist_value[0][k] = sg; 00382 cdist_value[1][k] = dg; 00383 cdist_prob[k++] = prob; 00384 } 00385 } 00386 } 00387 if (k == 0) throw exception(" Population.full_cdist(...): you have probably found a bug!"); 00388 for (prob=0.0,i=0; i<k; i++) prob += cdist_prob[i]; 00389 for (i=0; i<k; i++) cdist_prob[i] /= prob; 00390 return k; 00391 }

double matvec::Population::fullsibs_prob ( Individual * dam,

Individual * sire,

Individual * excludeI,

doubleMatrix & wspace

)

Definition at line 706 of file pop_peeling.cpp.
References matvec::Matrix< double >::assign(), matvec::Individual::get_penetrance(), get_posterior(), matvec::Individual::noffs_spouse(), matvec::Individual::nspouse(), matvec::Individual::offspring(), matvec::Individual::spouses(), tn_genotype, and trans_mat.
Referenced by anterior(), and posterior().

00709 { 00710 wspace.assign(1.0); 00711 unsigned j,k,um,uf,uj,nsp,nfs; 00712 const unsigned *noffs_sp; 00713 double val, scale; 00714 Vector<double> post_vec(tn_genotype); 00715 Vector<double> pen_vec(tn_genotype); 00716 Individual *J,**offspring, **fullsibs,**spouses; 00717 00718 spouses = sire->spouses(); 00719 offspring = sire->offspring(); 00720 nsp = sire->nspouse(); 00721 noffs_sp = sire->noffs_spouse(); 00722 for (k=0, j=0; j<nsp; j++) { 00723 if ( spouses[j] == dam) break; 00724 k += noffs_sp[j]; 00725 } 00726 nfs = noffs_sp[j]; 00727 fullsibs = &offspring[k]; 00728 00729 for (scale = 0.0, j=0; j<nfs; j++) { 00730 J = fullsibs[j]; 00731 if (J != excludeI) { 00732 nsp = J->nspouse(); 00733 if (nsp ) scale += get_posterior(J,0,post_vec); 00734 scale += J->get_penetrance(pen_vec); 00735 for (um=0; um<tn_genotype; um++) { 00736 for (uf=0; uf<tn_genotype; uf++) { 00737 val = 0.0; 00738 if (nsp) { 00739 for (uj=0; uj<tn_genotype; uj++) { 00740 val += trans_mat[uj][um][uf]*pen_vec[uj]*post_vec[uj]; 00741 } 00742 } 00743 else { 00744 for (uj=0; uj<tn_genotype; uj++) { 00745 val += trans_mat[uj][um][uf]*pen_vec[uj]; 00746 } 00747 } 00748 wspace[um][uf] *= val; 00749 } 00750 } 00751 } 00752 } 00753 // std::cout << "scale from fullsibs " << scale << std::endl; 00754 return scale; 00755 }

void matvec::Population::genotype_config ( const char type[] )

Definition at line 440 of file pop_gibbs.cpp.
References matvec::Vector< T >::begin(), build_pop_gamete(), matvec::Genome::chromosome, matvec::Individual::father(), ind_gamete(), matvec::Individual::maternal_chrom(), matvec::Individual::mother(), matvec::Genome::nchrom(), numchrom, matvec::Individual::paternal_chrom(), pop_gamete, popmember, popsize, and matvec::sampling().

00442 { 00443 //////////////////////////////////////////////////////////////// 00444 // a possible genotype configuration G over a given pedigree, 00445 ///////////////////////////////////////////////////////////////// 00446 00447 if (!pop_gamete) build_pop_gamete(); 00448 00449 Individual *I; 00450 unsigned i,j, s,d,nsg,ndg; 00451 Chromosome *C0, *C1; 00452 00453 unsigned maxg = pop_gamete[0].nchrom(); // ng = maxximum number of gametes 00454 for (i=1; i<numchrom; i++) { // among chromosomes 00455 j = pop_gamete[i].nchrom(); 00456 if (j > maxg) maxg = j; 00457 } 00458 ///////////////////////////////////////////////////////////// 00459 // get working space to hold parents gametes information 00460 ///////////////////////////////////////////////////////////// 00461 00462 Vector<unsigned> s_gamete_id(maxg); 00463 Vector<unsigned> d_gamete_id(maxg); 00464 Vector<double> s_gamete_freq(maxg); 00465 Vector<double> d_gamete_freq(maxg); 00466 00467 if (strcmp(type,"qtl") == 0) { 00468 00469 for (i=0; i<popsize; i++) { 00470 I = popmember[i]; 00471 C0 = I->paternal_chrom(); 00472 C1 = I->maternal_chrom(); 00473 for (j=0; j<numchrom; j++) { 00474 nsg = ind_gamete(I->father(), j, s_gamete_id.begin(),s_gamete_freq.begin()); 00475 ndg = ind_gamete(I->mother(), j, d_gamete_id.begin(),d_gamete_freq.begin()); 00476 s = sampling(s_gamete_freq.begin(), nsg); 00477 d = sampling(d_gamete_freq.begin(), ndg); 00478 C0[j] = pop_gamete[j].chromosome[s_gamete_id[s]]; 00479 C1[j] = pop_gamete[j].chromosome[d_gamete_id[d]]; 00480 } 00481 } 00482 } 00483 else if (strcmp(type,"ml") == 0) { 00484 ; 00485 } 00486 else if (strcmp(type,"ml_qtl") == 0) { 00487 ; 00488 } 00489 else { 00490 ; 00491 } 00492 }

void matvec::Population::genotype_dist_peeling ( const int prtlevel = 1,

const int estifreq = 1

)

Definition at line 1097 of file pop_peeling.cpp.
References matvec::Individual::anterior, matvec::Individual::base(), matvec::Vector< T >::begin(), matvec::Individual::est_GV, matvec::Individual::father_id(), matvec::Individual::id(), matvec::Vector< double >::inner_product(), matvec::Individual::mother_id(), matvec::Individual::nspouse(), matvec::Session::output_precision, matvec::Individual::posterior, matvec::Individual::residual_var, matvec::Vector< double >::resize(), matvec::SESSION, and matvec::Individual::sex().
Referenced by matvec::Model::genotype_dist_peeling(), matvec::Model::genotypic_value_peeling(), matvec::Model::graph_log_likelihood_peeling(), matvec::Model::log_likelihood_peeling(), maxant_maxpost(), and maxant_maxpost_old().

01100 { 01101 ////////////////////////////////////////////////////////////////////////// 01102 // this works for arbitrary pedigrees 01103 // based on iterative peeling algorithm: Arendonk(1989), Fernando(1993) 01104 // prtlevel = 0, print nothing 01105 // 1, default, print genotype distribution and genotypic values 01106 // estifreq = 0, not estimate allele frequencies 01107 // 1, default, estimate allele frequencies 01108 ////////////////////////////////////////////////////////////////////////// 01109 build_trans_mat(); 01110 if (!trans_mat) throw exception("Population::genotype_dist_peeling(): fail to build trans matrix"); 01111 Vector<double> post_vec(tn_genotype); 01112 Vector<double> genotype_freq = get_genotype_freq(); 01113 double *gfreq = genotype_freq.begin(); 01114 doubleMatrix wspace(tn_genotype,tn_genotype); 01115 unsigned i,j,k,nallele; 01116 doubleMatrix I_genotype_freq(popsize,tn_genotype); // this is working space 01117 Individual *I; 01118 double *ve,val,diff,maxchange, maxchange_freq; 01119 unsigned inner_niterate,outer_niterate; 01120 unsigned nb = this->nbase(); 01121 double allelefreq; 01122 outer_niterate = 0; 01123 do { 01124 genotype_freq = get_genotype_freq(); 01125 gfreq = genotype_freq.begin(); 01126 01127 for (i=0; i<popsize; i++) { 01128 I = popmember[i]; 01129 I->residual_var = &residual_var; // this is for time-being 01130 //////////////////////////////////////////////////////////////////// 01131 // the last column in anterior and posterior are for scaling factor 01132 //////////////////////////////////////////////////////////////////// 01133 I->initial_anterior(gfreq,tn_genotype); 01134 k = I->nspouse(); 01135 if (!I->posterior) { 01136 if(k>0){ 01137 I->posterior = new Vector<double> [k]; 01138 } 01139 else { 01140 I->posterior = 0; 01141 } 01142 } 01143 for (j=0; j<k; j++) I->posterior[j].resize(tn_genotype+1,1.0); 01144 } 01145 01146 inner_niterate = 0; 01147 do { 01148 inner_niterate++; 01149 maxchange = 0.0; 01150 for (i=0; i<popsize; i++) { 01151 I = popmember[i]; 01152 ve = I_genotype_freq[i]; 01153 anterior_posterior(I, wspace); 01154 get_posterior(I,0,post_vec); 01155 for (val=0.0,j=0; j<tn_genotype; j++) { 01156 val += gfreq[j] = I->anterior[j]*post_vec[j]; 01157 } 01158 if (inner_niterate > 1) { 01159 for (j=0; j<tn_genotype; j++) { 01160 gfreq[j] /= val; 01161 diff = fabs(ve[j] - gfreq[j]); 01162 if (diff > maxchange) maxchange = diff; 01163 } 01164 } 01165 for (j=0; j<tn_genotype; j++) ve[j] = gfreq[j]; 01166 } 01167 } while (inner_niterate == 1 || maxchange >= 0.000001); 01168 01169 if (prtlevel==1) { 01170 std::cout << " the # of inner iterations: " << inner_niterate << "\n"; 01171 } 01172 01173 if (estifreq == 1) { 01174 ////////////////////////////////////////////////////////////////// 01175 // it only works for the case of two-allele and one chromosome 01176 // PLEASE FIX ME 01177 /////////////////////////////////////////////////////////////////// 01178 01179 if (prior->nchrom() != 1 || prior->chrom()[0].nloci() != 1 01180 || prior->chrom()[0].locus[0].nallele() != 2) { 01181 throw exception("genotype_dist_peeling(): works only for the 1-chrom and 2-allele"); 01182 } 01183 outer_niterate++; 01184 if (prtlevel==1) { 01185 std::cout << " outer iteration: " << outer_niterate << "\n"; 01186 } 01187 for (allelefreq=0.0,i=0; i<popsize; i++) { 01188 I = popmember[i]; 01189 if (I->base()) { 01190 ve = I_genotype_freq[i]; 01191 allelefreq += ve[0]; 01192 allelefreq += ve[1] * 0.5; 01193 } 01194 } 01195 allelefreq /= nb; 01196 val = prior->chrom()[0].locus[0].allele_freq[0]; 01197 maxchange_freq = fabs(allelefreq - val); 01198 prior->chrom()[0].locus[0].allele_freq[0] = allelefreq; 01199 prior->chrom()[0].locus[0].allele_freq[1] = 1.0-allelefreq; 01200 if(pop_gamete){ 01201 delete [] pop_gamete; 01202 pop_gamete=0; 01203 } 01204 pop_gamete = 0; 01205 } 01206 else { 01207 maxchange_freq = 0.0; 01208 } 01209 } while (maxchange_freq >= 0.00001); 01210 if (prtlevel == 1) { 01211 std::cout << "allele frequencies\n"; 01212 for (i=0; i<prior->nchrom(); i++) { 01213 for (j=0; j<prior->chrom()[i].nloci(); j++) { 01214 for (k=0; k<prior->chrom()[i].locus[j].nallele(); k++) { 01215 std::cout << prior->chrom()[i].locus[j].allele_freq[k] << " "; 01216 } 01217 std::cout << "\n"; 01218 } 01219 } 01220 } 01221 01222 I = popmember[popsize-1]; 01223 double llhood_val = get_posterior(I,0,post_vec); 01224 val = I->anterior.inner_product(post_vec); 01225 llhood_val += std::log(val) + I->anterior[tn_genotype]; 01226 std::cout << "log_likelihood = " << llhood_val << "\n"; 01227 01228 int W = SESSION.output_precision + 6; // 6 = +.e+00 01229 const double **gv; 01230 01231 unsigned kk,nchrom=prior->nchrom(); 01232 for (kk=0,k=0; k<nchrom; k++) { 01233 gv = prior->genotypic_val(k,0); 01234 nallele = prior->chrom()[k].locus[0].nallele(); 01235 for (i=0; i<nallele; i++) for (j=0; j<=i; j++) { 01236 gfreq[kk++] = gv[i][j]; 01237 } 01238 } 01239 for (i=0; i<popsize; i++) { 01240 I = popmember[i]; 01241 if (prtlevel==1) { 01242 std::cout << ind_name(I->id()) 01243 << " " << I->sex() 01244 << " " << I->id() 01245 << " " << I->mother_id() 01246 << " " << I->father_id() 01247 << "\n"; 01248 } 01249 ve = I_genotype_freq[i]; 01250 for (val=0.0, j=0; j<tn_genotype; j++) { 01251 val += ve[j] * gfreq[j]; 01252 if (prtlevel==1) std::cout << " " << std::setw(W) << ve[j]; 01253 } 01254 I->est_GV = val; 01255 if (prtlevel==1) std::cout << " " << std::setw(W) << val << "\n"; 01256 } 01257 }

Vector< double > matvec::Population::get_genotype_freq ( void )

Definition at line 679 of file pop_peeling.cpp.
References build_pop_gamete(), matvec::Genome::chromosome, gtindex(), matvec::Genome::nchrom(), matvec::GeneticDist::nchrom(), pop_gamete, prior, tn_gamete, and tn_genotype.
Referenced by matvec::NuFamily::cutting(), maxant_maxpost(), and maxant_maxpost_old().

00680 { 00681 if (!pop_gamete) build_pop_gamete(); 00682 unsigned nchrom = prior->nchrom(); 00683 Vector<double> genotype_freq(tn_genotype); 00684 Vector<double> gamete_freq(tn_gamete); 00685 Vector<unsigned> ngtvec(nchrom); 00686 Vector<unsigned> gtvec(nchrom); 00687 double val; 00688 unsigned i,j,k; 00689 for (i=0; i<nchrom; i++) ngtvec[i] = pop_gamete[i].nchrom(); 00690 for (i=0; i<tn_gamete; i++) { 00691 gtindex(i,nchrom,ngtvec,gtvec); 00692 for (val=1.0,j=0; j<nchrom; j++) { 00693 val *= pop_gamete[j].chromosome[gtvec[j]].freq(); 00694 } 00695 gamete_freq[i] = val; 00696 } 00697 for (k=0,i=0; i<tn_gamete; i++) { 00698 for (j=0; j<i; j++) { 00699 genotype_freq[k++] = 2.0*gamete_freq[i]*gamete_freq[j]; 00700 } 00701 genotype_freq[k++] = gamete_freq[i]*gamete_freq[i]; 00702 } 00703 return genotype_freq; 00704 }

unsigned matvec::Population::get_id ( const char name[] ) [inline]

Definition at line 153 of file population.h.
References matvec::HashTable::get_id().
Referenced by ancestor_of(), descendant_of(), family_of(), and relative_of().

00153 {return hashtable.get_id(name);}

double matvec::Population::get_posterior ( Individual * I,

Individual * exclJ,

Vector< double > & vec

)

Definition at line 757 of file pop_peeling.cpp.
References matvec::Vector< double >::begin(), matvec::Individual::nspouse(), matvec::Individual::posterior, matvec::Individual::spouses(), and tn_genotype.
Referenced by anterior(), fullsibs_prob(), maxant_maxpost(), maxant_maxpost_old(), and posterior().

00758 { 00759 unsigned i,j,nsp = I->nspouse(); 00760 Individual **spouses = I->spouses(); 00761 double retval,*ve; 00762 for (i=0; i<tn_genotype; i++) vec[i]=1.0; 00763 for (retval=0.0,i=0; i<nsp; i++) { 00764 if (spouses[i] != exclJ) { 00765 ve = I->posterior[i].begin(); 00766 for (j=0; j<tn_genotype; j++) { 00767 vec[j] *= *ve++; 00768 } 00769 retval += *ve; 00770 // std::cout << i << " retval " << retval << std::endl; 00771 } 00772 } 00773 // std::cout << "retval from get_posterior " << retval << std::endl; 00774 return retval; 00775 }

void matvec::Population::getGNodeListSample ( unsigned maxsize,

unsigned numLoci,

string samplerType

)

Definition at line 2569 of file population.cpp.
References gNodeList, initJointAlleleNodeList(), matvec::GNodeList::peelCutAndSample(), and matvec::GNodeList::reinitGNodeList().

02569 { 02570 // Authors: L. Radu Totir and Rohan L. Fernando 02571 // (September, 2004) 02572 // Contributors: 02573 if(samplerType=="genotypic"){ 02574 cerr << "For joint peeling, the sampler type must be allelic;" << endl; 02575 exit(1); 02576 } 02577 else if(samplerType=="allelic"){ 02578 initJointAlleleNodeList(numLoci); 02579 } 02580 else { 02581 cerr << "Sampler type should be allelic;" << endl; 02582 exit(1); 02583 } 02584 bool notSampled = true; 02585 while(notSampled){ 02586 try { 02587 gNodeList.peelCutAndSample(maxsize); 02588 notSampled = false; 02589 } 02590 catch (matvec::InvalidSample) { 02591 cout << "Inconsistent sample, I will try to sample again" << endl; 02592 gNodeList.reinitGNodeList(); 02593 } 02594 //setSegregationIndex(j,samplerType); 02595 } 02596 }

void matvec::Population::getGNodeListSample ( unsigned maxsize,

unsigned startBlock,

unsigned stopBlock,

unsigned sizeBlock,

string samplerType

)

obtain a new candidate sample using peeling and cutting (if necessary)

Definition at line 1489 of file population.cpp.
References matvec::GeneticDist::chrom(), gNodeList, initAlleleNodeList(), initGenotypeNodeList(), matvec::ChromStruct::nloci(), matvec::GNodeList::peelCutAndSample(), prior, matvec::GNodeList::reinitGNodeList(), and setSegregationIndex().
Referenced by matvec::Model::RSamplerGibbs(), matvec::Model::RSamplerGibbsMH(), and matvec::Model::RSamplerMH().

01489 { 01490 // Authors: L. Radu Totir and Rohan L. Fernando 01491 // (October, 2003) 01492 // Contributors: 01493 for (int j=0;j<prior->chrom()[0].nloci();j++){ 01494 // cout << "Sample locus: " << j+1 << endl; 01495 if(samplerType=="genotypic"){ 01496 initGenotypeNodeList(j); 01497 } 01498 else if(samplerType=="allelic"){ 01499 initAlleleNodeList(j); 01500 } 01501 else { 01502 cerr << "Sampler type should be either genotypic or allelic;" << endl; 01503 exit(1); 01504 } 01505 bool notSampled = true; 01506 while(notSampled){ 01507 try { 01508 gNodeList.peelCutAndSample(maxsize,startBlock,stopBlock,sizeBlock); 01509 notSampled = false; 01510 } 01511 catch (matvec::InvalidSample) { 01512 cout << "Inconsistent sample, I will try to sample again" << endl; 01513 gNodeList.reinitGNodeList(); 01514 } 01515 } 01516 setSegregationIndex(j,samplerType); 01517 } 01518 }

void matvec::Population::getInitialGNodeListSample ( unsigned maxsize,

unsigned numLoci,

string samplerType

)

Definition at line 2535 of file population.cpp.
References matvec::GeneticDist::chrom(), matvec::GNodeList::displayGNodeSets(), gNodeList, initJointAlleleNodeList(), numFounders, matvec::ChromStruct::peelOrder, matvec::GNodeList::peelOrderCutAndSample(), popsize, prior, and matvec::GNodeList::reinitGNodeList().

02535 { 02536 // Authors: L. Radu Totir and Rohan L. Fernando 02537 // (September, 2004) 02538 // Contributors: 02539 if (samplerType=="genotypic"){ 02540 cerr << "For joint peeling, the sampler type must be allelic;" << endl; 02541 exit(1); 02542 } 02543 else if (samplerType=="allelic"){ 02544 prior->chrom()[0].peelOrder.resize(2*numLoci*(popsize + popsize-numFounders)); 02545 initJointAlleleNodeList(numLoci); 02546 gNodeList.displayGNodeSets(); 02547 } 02548 else { 02549 cerr << "Sampler type should be allelic;" << endl; 02550 exit(1); 02551 } 02552 bool notSampled = true; 02553 while(notSampled){ 02554 try { 02555 gNodeList.peelOrderCutAndSample(maxsize); 02556 notSampled = false; 02557 } 02558 catch (matvec::InvalidSample) { 02559 cout << "Inconsistent sample, I will try to sample again" << endl; 02560 gNodeList.reinitGNodeList(); 02561 } 02562 } 02563 }

void matvec::Population::getInitialGNodeListSample ( unsigned maxsize,

unsigned startBlock,

unsigned stopBlock,

unsigned sizeBlock,

string samplerType

)

Definition at line 1451 of file population.cpp.
References matvec::GeneticDist::chrom(), gNodeList, initAlleleNodeList(), initGenotypeNodeList(), matvec::ChromStruct::locus, model, matvec::Model::myRSamplerParms, matvec::ChromStruct::nloci(), matvec::GNodeList::peelOrderCutAndSample(), popsize, matvec::RSamplerParms::printFlag, prior, matvec::GNodeList::reinitGNodeList(), and setSegregationIndex().
Referenced by matvec::Model::RSamplerGibbs(), matvec::Model::RSamplerGibbsMH(), matvec::Model::RSamplerInitialDG(), and matvec::Model::RSamplerMH().

01451 { 01452 // Authors: L. Radu Totir and Rohan L. Fernando 01453 // (February, 2004) 01454 // Contributors: 01455 unsigned printFlag = model->myRSamplerParms.printFlag; 01456 for (int j=0;j<prior->chrom()[0].nloci();j++){ 01457 if (printFlag>0) cout << "Sample locus: " << j+1 << endl; 01458 if (samplerType=="genotypic"){ 01459 prior->chrom()[0].locus[j].peelOrder.resize(popsize); 01460 initGenotypeNodeList(j); 01461 } 01462 else if (samplerType=="allelic"){ 01463 prior->chrom()[0].locus[j].peelOrder.resize(2*popsize); 01464 initAlleleNodeList(j); 01465 //gNodeList.displayGNodeSets(); 01466 } 01467 else { 01468 cerr << "Sampler type should be either genotypic or allelic;" << endl; 01469 exit(1); 01470 } 01471 bool notSampled = true; 01472 while(notSampled){ 01473 try { 01474 gNodeList.peelOrderCutAndSample(maxsize,startBlock,stopBlock,sizeBlock); 01475 notSampled = false; 01476 } 01477 catch (matvec::InvalidSample) { 01478 cout << "Inconsistent sample, I will try to sample again" << endl; 01479 gNodeList.reinitGNodeList(); 01480 } 01481 } 01482 setSegregationIndex(j,samplerType); 01483 } 01484 }

void matvec::Population::getOldGNodeListProbability ( unsigned maxsize,

unsigned startBlock,

unsigned stopBlock,

unsigned sizeBlock,

string samplerType

)

recompute the needed quantities for the existing (old) sample (for the MH step)

Definition at line 1523 of file population.cpp.
References matvec::GeneticDist::chrom(), gNodeList, initAlleleNodeList(), initGenotypeNodeList(), matvec::ChromStruct::nloci(), matvec::GNodeList::peelCutAndCompute(), prior, matvec::GNodeList::reinitGNodeList(), and setSegregationIndex().
Referenced by matvec::Model::RSamplerGibbsMH(), and matvec::Model::RSamplerMH().

01523 { 01524 // Authors: L. Radu Totir and Rohan L. Fernando 01525 // (October, 2003) 01526 // Contributors: 01527 for (int j=0;j<prior->chrom()[0].nloci();j++){ 01528 // cout << "Sample locus: " << j+1 << endl; 01529 if (samplerType=="genotypic"){ 01530 initGenotypeNodeList(j); 01531 } 01532 else if (samplerType=="allelic"){ 01533 initAlleleNodeList(j); 01534 } 01535 else { 01536 cerr << "Sampler type should be either genotypic or allelic;" << endl; 01537 exit(1); 01538 } 01539 bool notSampled = true; 01540 while(notSampled){ 01541 try { 01542 gNodeList.peelCutAndCompute(maxsize,startBlock,stopBlock,sizeBlock); 01543 notSampled = false; 01544 } 01545 catch (matvec::InvalidSample) { 01546 cout << "Inconsistent sample, I will try to sample again" << endl; 01547 gNodeList.reinitGNodeList(); 01548 } 01549 } 01550 setSegregationIndex(j,samplerType); 01551 } 01552 }

void matvec::Population::gibbs_iterate ( const int count_genotype = 0 )

Definition at line 554 of file pop_gibbs.cpp.
References matvec::Vector< T >::begin(), matvec::Matrix< T >::begin(), matvec::Genome::chromosome, count_genotype(), full_cdist(), matvec::Individual::maternal_chrom(), matvec::Genome::nchrom(), numchrom, matvec::Individual::paternal_chrom(), pop_gamete, popmember, popsize, matvec::sampling(), and size().
Referenced by matvec::Model::genotype_dist_gibbs(), and matvec::Model::genotypic_value_gibbs().

00555 { 00556 Individual *I; 00557 Chromosome *C0, *C1; 00558 unsigned i,j,k,k0,k1,size; 00559 00560 unsigned maxg = pop_gamete[0].nchrom(); // maxg = maximum number of gametes 00561 for (i=1; i<numchrom; i++) { // among chromosomes 00562 k = pop_gamete[i].nchrom(); 00563 if (k > maxg) maxg = k; 00564 } 00565 00566 ////////////////////////////////////////////////////// 00567 // get working space ready for full_cdist() 00568 // freq_mat, gid_mat must have at least four rows 00569 // cdist_value, cdist_prob 00570 //////////////////////////////////////////////////// 00571 00572 Matrix<double> freq_mat(4,maxg); 00573 Matrix<unsigned> gid_mat(4,maxg); 00574 Matrix<unsigned> cdist_value(2,maxg*maxg); 00575 Vector<double> cdist_prob(maxg*maxg); 00576 00577 for (i=0; i<popsize; i++) { 00578 I = popmember[i]; 00579 C0 = I->paternal_chrom(); 00580 C1 = I->maternal_chrom(); 00581 for (j=0; j<numchrom; j++) { 00582 size = full_cdist(I,j,cdist_value.begin(),cdist_prob.begin(),gid_mat.begin(),freq_mat.begin()); 00583 k = sampling(cdist_prob.begin(),size); 00584 k0 = cdist_value[0][k]; // paternal gamete id 00585 k1 = cdist_value[1][k]; // maternal gamete id 00586 C0[j] = pop_gamete[j].chromosome[k0]; 00587 C1[j] = pop_gamete[j].chromosome[k1]; 00588 } 00589 if (count_gt) count_genotype(I); 00590 } 00591 }

void matvec::Population::Gibbs_sample ( void )

Definition at line 1090 of file pop_graph.cpp.
References member(), matvec::Individual::sample_self(), and size().

01090 { 01091 // Authors: Mathias Schelling and Rohan L. Fernando 01092 // (1999) 01093 // Contributors: Fabiano V. Pita 01094 for (int i=0; i<size(); i++){ 01095 member(i)->sample_self(0);//maternal 01096 member(i)->sample_self(1);//paternal 01097 } 01098 }

void matvec::Population::gtindex ( const unsigned num,

const unsigned ni,

const Vector< unsigned > & ngt,

Vector< unsigned > & gtvec

) [static]

Definition at line 777 of file population.cpp.
Referenced by build_pop_gamete(), and get_genotype_freq().

00778 { 00779 //////////////////////////////////////////////////////////////////////// 00780 // suppose, 3 individuals, the 1st and 2nd can have 2 genotypes, 00781 // the 3nd can have 3 genotypes. thus ngt = {2,2,3}. 00782 // I want num be transformed into gtvec system according to ngt: 00783 // 00784 // num gtvec[0] gtvec[1] gtvec[2] 00785 // 0 -> 0 0 0 00786 // 1 -> 0 0 1 00787 // 2 -> 0 0 2 00788 // 3 -> 0 1 0 00789 // 4 -> 0 1 1 00790 // 5 -> 0 1 2 00791 // 6 -> 1 0 0 00792 // 7 -> 1 0 1 00793 // 8 -> 1 0 2 00794 // 9 -> 1 1 0 00795 // 10 -> 1 1 1 00796 // 11 -> 1 1 2 00797 ///////////////////////////////////////////////////////////////////// 00798 { 00799 long i; 00800 unsigned ir=num; 00801 for (i=ni-1; i>=0; i--) { 00802 gtvec[i] = static_cast<unsigned>(fmod(static_cast<double>(ir),static_cast<double>(ngt[i]))); 00803 ir /= ngt[i]; 00804 } 00805 } 00806 }

Vector< double > matvec::Population::inbcoef ( void )

Definition at line 455 of file population1.cpp.
References fdone, matvec::Individual::inbcoef(), inbcoef_meuwissen(), popmember, and popsize.

00456 { 00457 if (!fdone) inbcoef_meuwissen(); 00458 Vector<double> fvec(popsize); 00459 for (unsigned i=0; i<popsize; i++) fvec[i] = popmember[i]->inbcoef(); 00460 return fvec; 00461 }

void matvec::Population::inbcoef_meuwissen ( void )

Definition at line 351 of file population1.cpp.
References matvec::Individual::father_id(), fdone, matvec::Individual::inbc, matvec::Individual::mother_id(), popmember, popsize, renum(), matvec::Vector< T >::reserve(), and stdid.
Referenced by inbcoef().

00352 { 00353 if (!stdid) renum(); 00354 Vector<double> Fam; Fam.reserve(popsize+1); // f vector starting from 0 00355 double *L; 00356 if(popsize>0){ 00357 L = new double [popsize]; 00358 } 00359 else { 00360 L = 0; 00361 } 00362 L--; 00363 double *D; 00364 if(popsize>0){ 00365 D = new double [popsize]; 00366 } 00367 else { 00368 D = 0; 00369 } 00370 D--; 00371 unsigned *point; 00372 if(popsize>0){ 00373 point = new unsigned [popsize]; 00374 } 00375 else { 00376 point = 0; 00377 } 00378 memset(point,'\0',sizeof(unsigned)*popsize); 00379 point--; 00380 unsigned i,j,k,s,d,ks,kd,kk,previous_s,previous_d; 00381 double r,fi; 00382 Individual *I; 00383 00384 previous_s = 0; previous_d = 0; 00385 Fam[0] = -1.0; 00386 for (i=1; i<=popsize; i++) { 00387 I = popmember[i-1]; 00388 s = I->father_id(); 00389 d = I->mother_id(); 00390 00391 D[i] = 0.5 - 0.25*(Fam[s] + Fam[d]); 00392 if (s == 0 || d == 0) { 00393 Fam[i] = 0.0; 00394 } 00395 else if ((s == previous_s) && (d == previous_d)) { 00396 Fam[i] = Fam[i-1]; 00397 } 00398 else { 00399 D[i] = 0.5 - 0.25*(Fam[s] + Fam[d]); 00400 fi = -1.0; 00401 L[i] = 1.0; 00402 j = i; 00403 while (j != 0) { 00404 k = j; 00405 r = 0.5*L[k]; 00406 I = popmember[k-1]; 00407 ks = I->father_id(); 00408 kd = I->mother_id(); 00409 if (ks < kd) {kk = ks; ks = kd; kd = kk;} 00410 if (ks > 0) { 00411 while (point[k] > ks) k = point[k]; 00412 L[ks] += r; 00413 if (ks != point[k]) { 00414 point[ks] = point[k]; 00415 point[k] = ks; 00416 } 00417 if (kd > 0) { 00418 while (point[k] > kd) k = point[k]; 00419 L[kd] += r; 00420 if (kd != point[k]) { 00421 point[kd] = point[k]; 00422 point[k] = kd; 00423 } 00424 } 00425 } 00426 fi += L[j]*L[j]*D[j]; 00427 L[j] = 0.0; 00428 k = j; 00429 j = point[j]; 00430 point[k] = 0; 00431 } 00432 Fam[i] = fi; 00433 } 00434 previous_s = s; previous_d = d; 00435 } 00436 L++; 00437 if(L){ 00438 delete [] L; 00439 L=0; 00440 } 00441 D++; 00442 if(D){ 00443 delete [] D; 00444 D=0; 00445 } 00446 point++; 00447 if(point){ 00448 delete [] point; 00449 point=0; 00450 } 00451 for (i=0; i<popsize; i++) popmember[i]->inbc = Fam[i+1]; 00452 fdone = 1; 00453 }

void matvec::Population::inbcoef_quaas ( void )

Definition at line 325 of file population1.cpp.
References matvec::Individual::father_id(), matvec::Individual::father_inbcoef(), fdone, matvec::Individual::inbc, matvec::Individual::mother_id(), matvec::Individual::mother_inbcoef(), popmember, popsize, renum(), matvec::Vector< T >::reserve(), and stdid.

00326 { 00327 if (!stdid) renum(); 00328 unsigned i,j,s,d; 00329 Individual *I; 00330 double t; 00331 Vector<double> fvalue; fvalue.reserve(popsize+1); // f vector starting from 0 00332 fvalue[0] = 0.0; 00333 for (j=0; j<popsize; j++) { 00334 I = popmember[j]; 00335 t = 1.0 - 0.25*(fvalue[I->father_id()] + fvalue[I->mother_id()]); 00336 fvalue[j+1] += t; 00337 I->inbc = std::sqrt(t); 00338 for (i=j+1; i<popsize; i++) { 00339 I = popmember[i]; s = I->father_id(); d = I->mother_id(); 00340 t = 0.0; 00341 if (s > j) t += 0.5*(I->father_inbcoef()); 00342 if (d > j) t += 0.5*(I->mother_inbcoef()); 00343 I->inbc = t; 00344 fvalue[i+1] += t*t; 00345 } 00346 } 00347 for (i=0; i<popsize; i++) popmember[i]->inbc = fvalue[i+1] - 1.0; 00348 fdone = 1; 00349 }

unsigned matvec::Population::ind_gamete ( const Individual * I,

const unsigned jc,

unsigned g_id[],

double fq[]

)

Definition at line 65 of file pop_gibbs.cpp.
References matvec::Locus::allele, matvec::Genome::chromosome, matvec::Locus::effect, matvec::Chromosome::freq(), gametebase, matvec::Chromosome::locus, matvec::Individual::maternal_chrom_id(), matvec::Genome::nchrom(), matvec::Chromosome::nloci(), matvec::Individual::paternal_chrom_id(), and pop_gamete.
Referenced by cprob_children(), full_cdist(), genotype_config(), llhood_phenotype(), and partial_cdist().

00067 { 00068 unsigned i,k,retval,ii,jj; 00069 if (I) { 00070 ii = I->paternal_chrom_id(jc); 00071 jj = I->maternal_chrom_id(jc); 00072 if (ii == jj) { 00073 retval = 1; 00074 g_id[0] = ii; 00075 fq[0] = 1.0; 00076 } 00077 else { 00078 if (ii > jj) k = ii*(ii-1)/2 + jj; 00079 else k = jj*(jj-1)/2 + ii; 00080 Chromosome* C = &(gametebase[jc].chromosome[k]); 00081 retval = C->nloci(); 00082 for (i=0; i<retval; i++) { 00083 g_id[i] = C->locus[i].allele; 00084 fq[i] = C->locus[i].effect; 00085 } 00086 } 00087 } 00088 else { 00089 retval = pop_gamete[jc].nchrom(); 00090 for (i=0; i<retval; i++) { 00091 g_id[i] = i; fq[i] = pop_gamete[jc].chromosome[i].freq(); 00092 } 00093 } 00094 return retval; 00095 }

const char * matvec::Population::ind_name ( const unsigned id ) const [inline]

Definition at line 354 of file population.h.
References matvec::HashTable::find(), hashtable, and matvec::HashTable::size().
Referenced by anterior(), build_nufamily(), confirm_sex(), display(), matvec::NuFamily::display(), DisplayFreqHaploFounders(), displayGenotypeFrequencies(), displayHaplotypeFrequencies(), displaySegregationIndicators(), fetch_families(), multi_geno_dist_peeling(), multi_m_geno_dist_peeling(), posterior(), renum(), matvec::Individual::set_switch(), and sub().

00355 { 00356 if (hashtable.size()==0 ){ 00357 std::ostringstream ostr; 00358 ostr << id; 00359 std::string s = ostr.str(); 00360 return s.c_str(); 00361 } 00362 else { 00363 return (const char *)hashtable.find(id); 00364 } 00365 }

void matvec::Population::initAlleleNodeList ( unsigned whichLocus )

creates the allele node list, and the founder, penetrance, and transmission sets for each allele node

Definition at line 1289 of file population.cpp.
References matvec::GNodeSet::attachMeToMyGnodes(), matvec::GeneticDist::chrom(), matvec::GNodeList::completeSetofGNsts, matvec::GNode::connectFlag, matvec::Individual::currentLocus, gNodeList, matvec::GNodeList::howToSample, matvec::GNode::id, matvec::ChromStruct::locus, matvec::Individual::malleleStateNodeVector, matvec::Individual::myfather, matvec::Individual::mymother, matvec::GNode::numberOfCuts, matvec::TransmissionSet::offspring, matvec::AllelePenetranceSet::owner, matvec::DisAllelePenetranceSet::owner, matvec::Individual::palleleStateNodeVector, matvec::TransmissionSet::paternal, popmember, popsize, prior, matvec::GNodeList::releaseGNsts(), renum(), and vectorOfGNsts.
Referenced by getGNodeListSample(), getInitialGNodeListSample(), and getOldGNodeListProbability().

01289 { 01290 // Authors: L. Radu Totir and Rohan L. Fernando 01291 // (June, 2003) 01292 // Contributors: 01293 Individual::currentLocus=whichLocus; 01294 GNodeSet::currentLocus=whichLocus; 01295 renum(); 01296 Individual *ind, *mom, *dad; 01297 gNodeList.resize(2*popsize); 01298 gNodeList.howToSample = "single"; 01299 GNode::gNodeListPtr = &gNodeList; 01300 vectorOfGNsts.clear(); 01301 gNodeList.releaseGNsts(); 01302 unsigned i,j; 01303 for(i=0;i<popsize;i++){ 01304 j=i*2; 01305 ind=popmember[i]; 01306 mom=ind->mymother; 01307 dad=ind->myfather; 01308 AlleleStateNode* imAllele = &ind->malleleStateNodeVector[ind->currentLocus]; 01309 AlleleStateNode* ipAllele = &ind->palleleStateNodeVector[ind->currentLocus]; 01310 ipAllele->id = j+1; 01311 ipAllele->connectFlag = j+1; 01312 ipAllele->numberOfCuts= 0; 01313 imAllele->id = j+2; 01314 imAllele->connectFlag = j+2; 01315 imAllele->numberOfCuts= 0; 01316 gNodeList[j] = ipAllele; 01317 gNodeList[j+1] = imAllele; 01318 01319 if(prior->chrom()[0].locus[Individual::currentLocus].qtl_ml=='r'){ 01320 DisAllelePenetranceSet *aPenSet = new DisAllelePenetranceSet; 01321 aPenSet->owner = ind; 01322 aPenSet->insert(imAllele); 01323 aPenSet->insert(ipAllele); 01324 vectorOfGNsts.push_back(aPenSet); 01325 gNodeList.completeSetofGNsts.insert(aPenSet); 01326 aPenSet->attachMeToMyGnodes(); 01327 } 01328 else{ 01329 AllelePenetranceSet *aPenSet = new AllelePenetranceSet; 01330 aPenSet->owner = ind; 01331 aPenSet->insert(imAllele); 01332 aPenSet->insert(ipAllele); 01333 vectorOfGNsts.push_back(aPenSet); 01334 gNodeList.completeSetofGNsts.insert(aPenSet); 01335 aPenSet->attachMeToMyGnodes(); 01336 } 01337 01338 if(mom){ 01339 AlleleStateNode *mMAllele = &mom->malleleStateNodeVector[ind->currentLocus]; 01340 AlleleStateNode *mPAllele = &mom->palleleStateNodeVector[ind->currentLocus]; 01341 AlleleStateNode *pMAllele = &dad->malleleStateNodeVector[ind->currentLocus]; 01342 AlleleStateNode *pPAllele = &dad->palleleStateNodeVector[ind->currentLocus]; 01343 01344 TransmissionSet *mTransmSet = new TransmissionSet; 01345 mTransmSet->offspring = ind; 01346 mTransmSet->paternal = false; 01347 mTransmSet->insert(imAllele); 01348 mTransmSet->insert(mMAllele); 01349 mTransmSet->insert(mPAllele); 01350 vectorOfGNsts.push_back(mTransmSet); 01351 gNodeList.completeSetofGNsts.insert(mTransmSet); 01352 mTransmSet->attachMeToMyGnodes(); 01353 01354 TransmissionSet *pTransmSet = new TransmissionSet; 01355 pTransmSet->offspring = ind; 01356 pTransmSet->paternal = true; 01357 pTransmSet->insert(ipAllele); 01358 pTransmSet->insert(pMAllele); 01359 pTransmSet->insert(pPAllele); 01360 vectorOfGNsts.push_back(pTransmSet); 01361 gNodeList.completeSetofGNsts.insert(pTransmSet); 01362 pTransmSet->attachMeToMyGnodes(); 01363 } 01364 else { 01365 AlleleFounderSet *mFoundSet = new AlleleFounderSet; 01366 mFoundSet->insert(imAllele); 01367 vectorOfGNsts.push_back(mFoundSet); 01368 gNodeList.completeSetofGNsts.insert(mFoundSet); 01369 mFoundSet->attachMeToMyGnodes(); 01370 01371 AlleleFounderSet *pFoundSet = new AlleleFounderSet; 01372 pFoundSet->insert(ipAllele); 01373 vectorOfGNsts.push_back(pFoundSet); 01374 gNodeList.completeSetofGNsts.insert(pFoundSet); 01375 pFoundSet->attachMeToMyGnodes(); 01376 } 01377 } 01378 }

void matvec::Population::initGenotypeFreq ( void )

Definition at line 2032 of file population.cpp.
References matvec::GeneticDist::chromosome, matvec::Individual::genotNodeVector, matvec::ChromStruct::locus, popmember, popsize, and prior.
Referenced by matvec::Model::RSamplerMH().

02032 { 02033 Individual *ind; 02034 for(unsigned i=0;i<popsize;i++){ 02035 ind=popmember[i]; 02036 for (unsigned j=0;j<Individual::numLoci;j++){ 02037 unsigned numAlleles = prior->chromosome[0].locus[j].nallele(); 02038 ind->genotNodeVector[j].genotypeCount.resize(numAlleles*numAlleles,0); 02039 } 02040 } 02041 }

void matvec::Population::initGenotypeNodeList ( unsigned whichLocus )

creates the genotype node list, and the founder, penetrance, and transmission sets for each genotype node

Definition at line 1384 of file population.cpp.
References matvec::GNodeSet::attachMeToMyGnodes(), matvec::GeneticDist::chrom(), matvec::GNodeList::completeSetofGNsts, matvec::GNode::connectFlag, matvec::Individual::currentLocus, matvec::Individual::genotNodeVector, gNodeList, matvec::GNodeList::howToSample, matvec::GNode::id, matvec::ChromStruct::locus, matvec::Individual::myfather, matvec::Individual::mymother, matvec::GNode::numberOfCuts, matvec::TransitionSet::offspring, matvec::DisGenoPenetranceSet::owner, matvec::GenoPenetranceSet::owner, popmember, popsize, prior, matvec::GNodeList::releaseGNsts(), renum(), and vectorOfGNsts.
Referenced by getGNodeListSample(), getInitialGNodeListSample(), and getOldGNodeListProbability().

01384 { 01385 // Authors: L. Radu Totir and Rohan L. Fernando 01386 // (June, 2003) 01387 // Contributors: 01388 Individual::currentLocus=whichLocus; 01389 GNodeSet::currentLocus=whichLocus; 01390 renum(); 01391 Individual *ind, *mom, *dad; 01392 gNodeList.resize(popsize); 01393 gNodeList.howToSample = "single"; 01394 GNode::gNodeListPtr = &gNodeList; 01395 vectorOfGNsts.clear(); 01396 gNodeList.releaseGNsts(); 01397 unsigned i,j; 01398 for(i=0;i<popsize;i++){ 01399 ind=popmember[i]; 01400 mom=ind->mymother; 01401 dad=ind->myfather; 01402 GenotypeNode* indGenotype = &ind->genotNodeVector[ind->currentLocus]; 01403 01404 indGenotype->id = i+1; 01405 indGenotype->connectFlag = i+1; 01406 indGenotype->numberOfCuts= 0; 01407 gNodeList[i] = indGenotype; 01408 01409 if(prior->chrom()[0].locus[Individual::currentLocus].qtl_ml=='q'){ 01410 GenoPenetranceSet *gPenSet = new GenoPenetranceSet; 01411 gPenSet->owner = ind; 01412 gPenSet->insert(indGenotype); 01413 vectorOfGNsts.push_back(gPenSet); 01414 gNodeList.completeSetofGNsts.insert(gPenSet); 01415 gPenSet->attachMeToMyGnodes(); 01416 } 01417 else if(prior->chrom()[0].locus[Individual::currentLocus].qtl_ml=='r'){ 01418 DisGenoPenetranceSet *gPenSet = new DisGenoPenetranceSet; 01419 gPenSet->owner = ind; 01420 gPenSet->insert(indGenotype); 01421 vectorOfGNsts.push_back(gPenSet); 01422 gNodeList.completeSetofGNsts.insert(gPenSet); 01423 gPenSet->attachMeToMyGnodes(); 01424 } 01425 if(mom){ 01426 GenotypeNode *momGenotype = &mom->genotNodeVector[ind->currentLocus]; 01427 GenotypeNode *dadGenotype = &dad->genotNodeVector[ind->currentLocus]; 01428 TransitionSet *TransitSet = new TransitionSet; 01429 TransitSet->offspring = ind; 01430 TransitSet->insert(indGenotype); 01431 TransitSet->insert(momGenotype); 01432 TransitSet->insert(dadGenotype); 01433 vectorOfGNsts.push_back(TransitSet); 01434 gNodeList.completeSetofGNsts.insert(TransitSet); 01435 TransitSet->attachMeToMyGnodes(); 01436 } 01437 else { 01438 GenoFounderSet *FoundSet = new GenoFounderSet; 01439 FoundSet->insert(indGenotype); 01440 vectorOfGNsts.push_back(FoundSet); 01441 gNodeList.completeSetofGNsts.insert(FoundSet); 01442 FoundSet->attachMeToMyGnodes(); 01443 } 01444 } 01445 }

void matvec::Population::initialize ( void )

Definition at line 109 of file population.cpp.
References fdone, gametebase, lookOnlyAtLocusToYourLeft, mark_need_remove, markerDataIn, markersymbol, maxnamelen, maxpopsize, n_markerLoci, nmarker, nuFamiliesDone, nufamily_vec, num_nufamily, numchrom, numtrait, offs_info_built, pm_storage, pop_gamete, popmember, popsize, prior, sex_confirmed, spouse_info_built, stdid, and trans_mat.
Referenced by Population().

00110 { 00111 popsize = 0; 00112 maxpopsize = 0; 00113 numchrom = 0; 00114 numtrait = 0; 00115 nmarker = 0; 00116 fdone = 0; 00117 stdid = 0; 00118 num_nufamily = 0; 00119 maxnamelen = 0; 00120 //RLF 00121 n_markerLoci = 0; 00122 nuFamiliesDone = false; 00123 //RLF 00124 sex_confirmed = 0; 00125 mark_need_remove = 0; 00126 spouse_info_built = 0; 00127 offs_info_built = 0; 00128 00129 prior = 0; 00130 trans_mat = 0; 00131 pop_gamete = 0; 00132 gametebase = 0; 00133 pm_storage = 0; 00134 popmember = 0; 00135 nufamily_vec = 0; 00136 markersymbol = 0; 00137 markerDataIn = false; 00138 // LRT 00139 lookOnlyAtLocusToYourLeft = true; 00140 }

void matvec::Population::initJointAlleleNodeList ( unsigned howManyLoci )

Definition at line 2391 of file population.cpp.
References matvec::GNodeSet::attachMeToMyGnodes(), matvec::GeneticDist::chrom(), matvec::GNodeList::completeSetofGNsts, matvec::GNode::connectFlag, matvec::RAlleleFounderSet::forLocus, matvec::RTransmissionSet::forLocus, matvec::RAllelePenetranceSet::forLocus, matvec::RDisAllelePenetranceSet::forLocus, gNodeList, matvec::GNodeList::howToSample, matvec::GNode::id, matvec::ChromStruct::locus, matvec::Individual::malleleOriginNodeVector, matvec::Individual::malleleStateNodeVector, matvec::Individual::myfather, matvec::Individual::mymother, matvec::GNode::numberOfCuts, numFounders, matvec::RTransmissionSet::offspring, matvec::RAllelePenetranceSet::owner, matvec::RDisAllelePenetranceSet::owner, matvec::Individual::palleleOriginNodeVector, matvec::Individual::palleleStateNodeVector, matvec::RTransmissionSet::paternal, popmember, popsize, prior, matvec::RecombinationSet::r, recombinationMatrix, matvec::GNodeList::releaseGNsts(), renum(), and vectorOfGNsts.
Referenced by getGNodeListSample(), and getInitialGNodeListSample().

02391 { 02392 // Authors: L. Radu Totir and Rohan L. Fernando 02393 // (August, 2004) 02394 // Contributors: 02395 renum(); 02396 Individual *ind, *mom, *dad; 02397 gNodeList.resize(2*howManyLoci*(popsize + popsize-numFounders)); 02398 gNodeList.howToSample = "joint"; 02399 GNode::gNodeListPtr = &gNodeList; 02400 vectorOfGNsts.clear(); 02401 gNodeList.releaseGNsts(); 02402 unsigned count = 0; 02403 for (int locus=0;locus<howManyLoci;locus++){ 02404 unsigned i,j,k; 02405 unsigned startPosition = (popsize + popsize-numFounders)*2*locus; 02406 unsigned countOffspring = 0; 02407 for(i=0;i<popsize;i++){ 02408 j=startPosition + i*2; 02409 ind=popmember[i]; 02410 mom=ind->mymother; 02411 dad=ind->myfather; 02412 AlleleStateNode* imAlleleState = &ind->malleleStateNodeVector[locus]; 02413 AlleleStateNode* ipAlleleState = &ind->palleleStateNodeVector[locus]; 02414 ipAlleleState->id = j+1; 02415 ipAlleleState->connectFlag = j+1; 02416 ipAlleleState->numberOfCuts= 0; 02417 imAlleleState->id = j+2; 02418 imAlleleState->connectFlag = j+2; 02419 imAlleleState->numberOfCuts= 0; 02420 gNodeList[j] = ipAlleleState; 02421 gNodeList[j+1] = imAlleleState; 02422 02423 if(prior->chrom()[0].locus[locus].qtl_ml=='r'){ 02424 RDisAllelePenetranceSet *aPenSet = new RDisAllelePenetranceSet; 02425 aPenSet->forLocus = locus; 02426 aPenSet->owner = ind; 02427 vectorOfGNsts.push_back(aPenSet); 02428 aPenSet->insert(imAlleleState); 02429 aPenSet->insert(ipAlleleState); 02430 gNodeList.completeSetofGNsts.insert(aPenSet); 02431 aPenSet->attachMeToMyGnodes(); 02432 } 02433 else{ 02434 RAllelePenetranceSet *aPenSet = new RAllelePenetranceSet; 02435 aPenSet->forLocus = locus; 02436 aPenSet->owner = ind; 02437 vectorOfGNsts.push_back(aPenSet); 02438 aPenSet->insert(imAlleleState); 02439 aPenSet->insert(ipAlleleState); 02440 gNodeList.completeSetofGNsts.insert(aPenSet); 02441 aPenSet->attachMeToMyGnodes(); 02442 } 02443 02444 if(mom){ 02445 k=startPosition + 2*popsize + countOffspring*2; 02446 countOffspring++; 02447 AlleleOriginNode *imAlleleOrigin = &ind->malleleOriginNodeVector[locus]; 02448 AlleleOriginNode *ipAlleleOrigin = &ind->palleleOriginNodeVector[locus]; 02449 ipAlleleOrigin->id = k+1; 02450 ipAlleleOrigin->connectFlag = k+1; 02451 ipAlleleOrigin->numberOfCuts= 0; 02452 imAlleleOrigin->id = k+2; 02453 imAlleleOrigin->connectFlag = k+2; 02454 imAlleleOrigin->numberOfCuts= 0; 02455 02456 gNodeList[k] = ipAlleleOrigin; 02457 gNodeList[k+1] = imAlleleOrigin; 02458 02459 RecombinationSet *mRecomSet = new RecombinationSet; 02460 mRecomSet->r=0.5; 02461 mRecomSet->insert(imAlleleOrigin); 02462 if(locus){ 02463 AlleleOriginNode *imLeftAlleleOrigin = &ind->malleleOriginNodeVector[locus-1]; 02464 mRecomSet->r=recombinationMatrix[locus-1][locus]; 02465 mRecomSet->insert(imLeftAlleleOrigin); 02466 } 02467 vectorOfGNsts.push_back(mRecomSet); 02468 gNodeList.completeSetofGNsts.insert(mRecomSet); 02469 mRecomSet->attachMeToMyGnodes(); 02470 02471 RecombinationSet *pRecomSet = new RecombinationSet; 02472 pRecomSet->r=0.5; 02473 pRecomSet->insert(ipAlleleOrigin); 02474 if(locus){ 02475 AlleleOriginNode *ipLeftAlleleOrigin = &ind->palleleOriginNodeVector[locus-1]; 02476 pRecomSet->r=recombinationMatrix[locus-1][locus]; 02477 pRecomSet->insert(ipLeftAlleleOrigin); 02478 } 02479 vectorOfGNsts.push_back(pRecomSet); 02480 gNodeList.completeSetofGNsts.insert(pRecomSet); 02481 pRecomSet->attachMeToMyGnodes(); 02482 02483 AlleleStateNode *mMAlleleState = &mom->malleleStateNodeVector[locus]; 02484 AlleleStateNode *mPAlleleState = &mom->palleleStateNodeVector[locus]; 02485 AlleleStateNode *pMAlleleState = &dad->malleleStateNodeVector[locus]; 02486 AlleleStateNode *pPAlleleState = &dad->palleleStateNodeVector[locus]; 02487 02488 RTransmissionSet *mTransmSet = new RTransmissionSet; 02489 mTransmSet->forLocus = locus; 02490 mTransmSet->offspring = ind; 02491 mTransmSet->paternal = false; 02492 mTransmSet->insert(imAlleleState); 02493 mTransmSet->insert(imAlleleOrigin); 02494 mTransmSet->insert(mMAlleleState); 02495 mTransmSet->insert(mPAlleleState); 02496 vectorOfGNsts.push_back(mTransmSet); 02497 gNodeList.completeSetofGNsts.insert(mTransmSet); 02498 mTransmSet->attachMeToMyGnodes(); 02499 02500 RTransmissionSet *pTransmSet = new RTransmissionSet; 02501 pTransmSet->forLocus = locus; 02502 pTransmSet->offspring = ind; 02503 pTransmSet->paternal = true; 02504 pTransmSet->insert(ipAlleleState); 02505 pTransmSet->insert(ipAlleleOrigin); 02506 pTransmSet->insert(pMAlleleState); 02507 pTransmSet->insert(pPAlleleState); 02508 vectorOfGNsts.push_back(pTransmSet); 02509 gNodeList.completeSetofGNsts.insert(pTransmSet); 02510 pTransmSet->attachMeToMyGnodes(); 02511 } 02512 else { 02513 RAlleleFounderSet *mFoundSet = new RAlleleFounderSet; 02514 mFoundSet->forLocus = locus; 02515 mFoundSet->insert(imAlleleState); 02516 vectorOfGNsts.push_back(mFoundSet); 02517 gNodeList.completeSetofGNsts.insert(mFoundSet); 02518 mFoundSet->attachMeToMyGnodes(); 02519 02520 RAlleleFounderSet *pFoundSet = new RAlleleFounderSet; 02521 pFoundSet->forLocus = locus; 02522 pFoundSet->insert(ipAlleleState); 02523 vectorOfGNsts.push_back(pFoundSet); 02524 gNodeList.completeSetofGNsts.insert(pFoundSet); 02525 pFoundSet->attachMeToMyGnodes(); 02526 } 02527 } 02528 } 02529 }

unsigned matvec::Population::input_data ( const char fname[],

GeneticDist * G

)

Definition at line 43 of file pop_mblup.cpp.
References matvec::Locus::allele, matvec::GeneticDist::chrom(), matvec::Genome::chromosome, matvec::DataNode::double_val(), matvec::Individual::genome0, matvec::Individual::genome1, matvec::INPUT_LINE_WIDTH, matvec::Chromosome::locus, matvec::ChromStruct::locus, matvec::Individual::myfather, matvec::Individual::myid, matvec::Individual::mymother, matvec::Individual::myrecord, nmarker, matvec::Individual::p_origin, popmember, popsize, resize(), stdid, and matvec::validline().

00044 { 00045 size_t linewidth = INPUT_LINE_WIDTH; 00046 char *line; 00047 if(linewidth>0){ 00048 line = new char [linewidth]; 00049 } 00050 else { 00051 line = 0; 00052 } 00053 std::ifstream in(fname); 00054 if (!in) { 00055 if(line){ 00056 delete [] line; 00057 line=0; 00058 } 00059 throw exception("Population::input_data(): cannot open file"); 00060 } 00061 unsigned numrec = 0; 00062 while (in.getline(line,linewidth)) if (validline(line)) numrec++; 00063 in.close(); 00064 resize(numrec,G); 00065 popsize = numrec; 00066 nmarker = G->chrom()[0].locus[0].nallele(); 00067 stdid = 1; 00068 // for gcc-3.2 replaced 00069 // std::istrstream str_in(line,linewidth); 00070 // with 00071 std::istringstream str_in(line); 00072 00073 00074 /////////////////// user attention begins //////////////////// 00075 // the example of data (pedigree) file 00076 // the requirements: 00077 // (1) animal id's must be standard (starting from 1, parents precede 00078 // their offspring, 0 means missing) 00079 // (2) marker genotype is represented by two columns: marker_allele ids 00080 // allele_id must be standard (starting from 1 and 0 means missing) 00081 // (3) parental origin of marker genotypes is unknown. 00082 // 00083 // 1 0 0 80 1 1 00084 // 2 0 0 120 0 0 00085 // 3 1 2 90 1 2 00086 // 4 0 2 110 1 2 00087 // 5 3 4 115 1 2 00088 // 00089 // note that . is not allowed to represent the missing value 00090 // 00091 unsigned ind,sire,dam; 00092 int a1,a2; 00093 double y; 00094 Individual* I; 00095 numrec = 0; 00096 in.open(fname); 00097 while (in.getline(line,linewidth)) { 00098 if (validline(line) ) { 00099 str_in >> ind >> sire >> dam >> y >> a1 >> a2; 00100 if (a1<0 || a1>nmarker || a2<0 || a2>nmarker) throw exception(" Pop.input_data(): invalid allele_id"); 00101 I = popmember[numrec++]; 00102 I->genome0.chromosome[0].locus[0].allele = a1-1; 00103 I->genome1.chromosome[0].locus[0].allele = a2-1; // locus[0] is ML 00104 I->myid = ind; 00105 if (sire == 0) { I->myfather = 0; } 00106 else { I->myfather = popmember[sire-1]; } 00107 if (dam == 0) { I->mymother = 0; } 00108 else { I->mymother = popmember[dam-1]; } 00109 I->myrecord[0].double_val(y); 00110 I->p_origin = 0; 00111 str_in.seekg(0,std::ios::beg); 00112 } 00113 } 00114 in.close(); 00115 if(line){ 00116 delete [] line; 00117 line=0; 00118 } 00119 return numrec; 00120 /////////////////// user attention end //////////////////// 00121 }

unsigned matvec::Population::input_data ( Data * D )

Definition at line 313 of file population.cpp.
References matvec::Locus::allele, matvec::GeneticDist::chrom(), matvec::Genome::chromosome, matvec::FieldStruct::classi(), matvec::Field::col_struct, matvec::Data::datasheet, matvec::DataNode::double_val(), matvec::HashTable::find(), matvec::Individual::genome0, matvec::HashTable::get_id(), hashtable, matvec::Data::hashtable, matvec::Data::in_memory(), matvec::Data::input_datasheet(), matvec::Chromosome::locus, matvec::ChromStruct::locus, matvec::DataNode::missing, matvec::Individual::myrecord, matvec::FieldStruct::name(), matvec::GeneticDist::name(), matvec::Data::num_cols(), matvec::Data::num_rows(), numtrait, popmember, prior, matvec::Data::rawcol(), matvec::Data::release_datasheet(), and matvec::DataNode::unsigned_val().
Referenced by matvec::Model::genotype_dist_gibbs(), matvec::Model::genotype_dist_peeling(), matvec::Model::genotypic_value_gibbs(), matvec::Model::genotypic_value_peeling(), matvec::Model::graph_log_likelihood_peeling(), and matvec::Model::log_likelihood_peeling().

00314 { 00315 if (!D || D->num_rows() <= 0) throw exception("Population::input_data(): data is empty"); 00316 unsigned numcol = D->num_cols(); 00317 unsigned numrec = D->num_rows(); 00318 unsigned i,j,t,id,dummy_ind; 00319 const char *strid; 00320 00321 if (!D->in_memory()) D->input_datasheet(); 00322 for (t=1; t<numcol; t++) { // 1st column is reserved for intercept 00323 // for (t=0; t<numcol; t++) { Old version 00324 if (D->datasheet[t].col_struct.classi()=='G') break; 00325 } 00326 if (t==numcol) throw exception("Population::input_data(): no pedigree exits"); 00327 DataNode *tcol, *column = D->rawcol(t); 00328 HashTable* dhtable = D->hashtable[t]; 00329 if (strcmp(prior->name(),"GeneticDist") == 0) { 00330 /////////////////////////////////////////////////////////// 00331 // need more work for finding marker-locus in prior 00332 /////////////////////////////////////////////////////////// 00333 Individual *I; 00334 int gtype[2]; 00335 if (prior->chrom()[0].locus[1].qtl_ml == 'm') { 00336 for (t=0; t<numcol; t++) { 00337 if (D->datasheet[t].col_struct.name()== "marker1") break; 00338 } 00339 if (t==numcol) throw exception(" Population::input_data(): no marker1 exits"); 00340 tcol = D->rawcol(t); // marker1 column 00341 for (i=0; i<numrec; i++) { 00342 if (tcol[i].double_val() == 0.0) { 00343 gtype[0]=0; 00344 gtype[1]=0; 00345 } 00346 else if (tcol[i].double_val() == 1.0) { 00347 gtype[0]=0; 00348 gtype[1]=1; 00349 } 00350 else if (tcol[i].double_val() == 2.0) { 00351 gtype[0]=1; 00352 gtype[1]=1; 00353 } 00354 else { 00355 throw exception("Population::input_data(): marker1: invalid genotype"); 00356 } 00357 00358 id = column[i].unsigned_val(); 00359 strid = (const char *)dhtable->find(id); 00360 id = hashtable.get_id(strid); 00361 if (id == 0) continue; 00362 I = popmember[id-1]; 00363 I->genome0.chromosome[0].locus[1].allele = gtype[0]; 00364 I->genome1.chromosome[0].locus[1].allele = gtype[1]; 00365 } 00366 } 00367 } 00368 00369 for (j=0, t=0; t<numcol; t++) { 00370 if (D->datasheet[t].col_struct.classi()=='T') j++; 00371 } 00372 if (j != numtrait) throw exception("Population::input_data(): Pop's ntrait != Data's ntrait"); 00373 // this is where the phenotype from the data file is 00374 // merged with the pedigree informatiojn 00375 for (j=0, t=1; t<numcol; t++) { // 1st column is reserved for intercept 00376 if (D->datasheet[t].col_struct.classi()=='T') { 00377 tcol = D->rawcol(t); 00378 dummy_ind = 0; 00379 for (i=0; i<numrec; i++) { 00380 id = column[i].unsigned_val(); 00381 strid = (const char *)dhtable->find(id); 00382 id = hashtable.get_id(strid); 00383 if (id > 0) { 00384 popmember[id-1]->myrecord[j] = tcol[i]; 00385 } 00386 else { 00387 column[i].missing = 1; 00388 dummy_ind++; 00389 } 00390 } 00391 j++; 00392 } 00393 } 00394 D->release_datasheet(); 00395 if (j==0) throw exception("Population::input_data(): trait(s) not available"); 00396 return (numrec - dummy_ind); 00397 }

unsigned matvec::Population::input_descentGraph ( char * dgfile )

Definition at line 853 of file population.cpp.
References matvec::HashTable::get_id(), hashtable, matvec::GeneticDist::nloci_chrom(), popmember, popsize, prior, matvec::Individual::put_gametes(), and matvec::Vector< T >::resize().

00853 { 00854 // Authors: Rohan L. Fernando and Fabiano V. Pita 00855 // (2003) 00856 // Contributors: L. Radu Totir 00857 std::string strid; 00858 int id, m1, m2; 00859 Individual *I; 00860 Vector <unsigned> m_g , p_g; 00861 unsigned nLoci = prior->nloci_chrom(1); 00862 m_g.resize(nLoci); 00863 p_g.resize(nLoci); 00864 std::ifstream mfile(dgfile); 00865 if(!mfile) { 00866 throw exception("Couldn't open file "); 00867 } 00868 for(int i=0; i<popsize;i++){ 00869 if (!(mfile >> strid)){ 00870 throw exception(" Not enough rows in dgfile "); 00871 } 00872 id = hashtable.get_id(strid.c_str()); 00873 I = popmember[id-1]; 00874 for(int j=1; j<=nLoci;j++){ 00875 if(!(mfile >> m_g(j) )){ 00876 std::cout << " Not enough columns in dgfile " << std::endl; 00877 } 00878 if(!(mfile >> p_g(j))){ 00879 throw exception(" Not enough columns in dgfile "); 00880 } 00881 } 00882 I->put_gametes(m_g, p_g); 00883 } 00884 return 1; 00885 }

unsigned matvec::Population::input_markerData ( Data * D )

Definition at line 908 of file population.cpp.
References matvec::Locus::allele, matvec::GeneticDist::chrom(), matvec::Genome::chromosome, matvec::FieldStruct::classi(), matvec::Field::col_struct, matvec::Data::datasheet, matvec::DataNode::double_val(), matvec::HashTable::find(), matvec::Individual::genome0, matvec::Individual::genome1, matvec::HashTable::get_id(), hashtable, matvec::Data::hashtable, matvec::Data::in_memory(), matvec::Data::input_datasheet(), matvec::Chromosome::locus, matvec::ChromStruct::locus, markerDataIn, matvec::Individual::myrecord, n_markerLoci, matvec::FieldStruct::name(), matvec::GeneticDist::nloci_chrom(), matvec::Data::num_cols(), matvec::Data::num_rows(), matvec::GeneticDist::numMarkerLoci, numtrait, popmember, prior, matvec::Data::rawcol(), matvec::Data::release_datasheet(), matvec::Vector< T >::resize(), matvec::Field::type(), and matvec::DataNode::unsigned_val().
Referenced by descent_graph_setup(), matvec::Model::multipoint_setup(), and setupRSampler().

00909 { 00910 if (markerDataIn) return D->num_rows(); 00911 if (!D || D->num_rows() <= 0) throw exception(" Population::input_data(): data is empty"); 00912 unsigned numcol = D->num_cols(); 00913 unsigned numrec = D->num_rows(); 00914 unsigned i,j,t,id; 00915 const char *strid; 00916 00917 if (!D->in_memory()) D->input_datasheet(); 00918 for (t=0; t<numcol; t++) { 00919 if (D->datasheet[t].col_struct.classi()=='G') break; 00920 } 00921 if (t==numcol) throw exception("Population::input_data(): no pedigree exits"); 00922 char TY = D->datasheet[t].type(); 00923 DataNode *tcol, *column = D->rawcol(t); 00924 HashTable* dhtable = D->hashtable[t]; 00925 Vector<DataNode*> markerCols1, markerCols2; 00926 std::string nm1,nm2; 00927 00928 // number of marker loci is assumed to be = number of loci - 1 00929 // locus(0) is the qtl and the remaining loci are all markers 00930 // n_markerLoci = prior->nloci_chrom(1) - 1; 00931 00932 // I (LRT) have replaced the line above with the following: 00933 n_markerLoci = prior->numMarkerLoci; 00934 markerCols1.resize(prior->nloci_chrom(1)); 00935 markerCols2.resize(prior->nloci_chrom(1)); 00936 00937 // find marker columns in data sheet for each marker locus 00938 // pointers to columns in data sheet are stored in 00939 // markerCols1 and markerCols2 00940 00941 for (i=0; i<prior->nloci_chrom(1); i++){ 00942 if(prior->chrom()[0].locus[i].qtl_ml=='m'){ 00943 nm1 = prior->chrom()[0].locus[i].nameOfcol1; 00944 nm2 = prior->chrom()[0].locus[i].nameOfcol2; 00945 00946 for (t=0; t<numcol; t++) { // column 1 00947 if (nm1 == D->datasheet[t].col_struct.name()) break; 00948 } 00949 if (t==numcol) throw exception("Population::input_data(): no data for nameOfcol1 locus i"); 00950 markerCols1(i+1) = D->rawcol(t); 00951 00952 for (t=0; t<numcol; t++) { // column 2 00953 if (nm2 == D->datasheet[t].col_struct.name()) break; 00954 } 00955 if (t==numcol) throw exception("Population::input_data(): no data for nameOfcol2 locus i"); 00956 markerCols2(i+1) = D->rawcol(t); 00957 } 00958 } 00959 Individual *I; 00960 int gtype[2]; 00961 00962 // Now we transfer the marker data from the data sheet to Individuals 00963 00964 for (i=0; i<numrec; i++) { 00965 if (TY == 'S') { 00966 id = column[i].unsigned_val(); 00967 strid = (const char *)dhtable->find(id); 00968 id = hashtable.get_id(strid); 00969 I = popmember[id-1]; 00970 } 00971 else { 00972 id = (unsigned) column[i].double_val(); 00973 I = popmember[id-1]; 00974 } 00975 for (j=0; j<prior->nloci_chrom(1); j++){ 00976 if(prior->chrom()[0].locus[j].qtl_ml=='m'){ 00977 I->genome0.chromosome[0].locus[j].allele = (int) markerCols1(j+1)[i].double_val(); 00978 I->genome1.chromosome[0].locus[j].allele = (int) markerCols2(j+1)[i].double_val(); 00979 } 00980 } 00981 } 00982 for (j=0, t=0; t<numcol; t++) { 00983 if (D->datasheet[t].col_struct.classi()=='T') j++; 00984 } 00985 if (j != numtrait) throw exception("Population::input_data(): Pop's ntrait != Data's ntrait"); 00986 for (j=0, t=1; t<numcol; t++) { // 1st column is reserved for intercept 00987 if (D->datasheet[t].col_struct.classi()=='T') { 00988 tcol = D->rawcol(t); 00989 for (i=0; i<numrec; i++) { 00990 if (TY == 'S') { 00991 id = column[i].unsigned_val(); 00992 strid = (const char *)dhtable->find(id); 00993 id = hashtable.get_id(strid); 00994 } 00995 else{ 00996 id = (unsigned) column[i].double_val(); 00997 } 00998 popmember[id-1]->myrecord[j] = tcol[i]; 00999 } 01000 j++; 01001 } 01002 } 01003 D->release_datasheet(); 01004 if (j == 0) throw exception("Population::input_data(): trait(s) not available"); 01005 //std::cout << "exit from input marker data \n"; 01006 markerDataIn = true; 01007 return numrec; 01008 }

unsigned matvec::Population::input_ped ( matvec::Pedigree & P,

GeneticDist * D

)

Definition at line 200 of file population.cpp.
References build_offs_info(), build_spouse_info(), confirm_sex(), countFounders(), input_ped0(), matvec::GeneticDist::name(), and matvec::Pedigree::sex_confirmed.

00200 { 00201 input_ped0(P,D); 00202 if (strcmp(D->name(),"UnknownDist") != 0) { 00203 if (P.sex_confirmed) { 00204 confirm_sex(1); 00205 } 00206 else { 00207 confirm_sex(-1); 00208 } 00209 countFounders(); 00210 build_offs_info(); 00211 build_spouse_info(); 00212 } 00213 }

unsigned matvec::Population::input_ped ( const char fname[],

const char recfmt[],

GeneticDist * D

)

Definition at line 301 of file population.cpp.
Referenced by setupRSampler().

00303 { 00304 ////////////////////////////////////////////////// 00305 // return 0 if everything is successful 00306 // 1 something wrong 00307 ///////////////////////////////////////////////// 00308 Pedigree ped; 00309 ped.input(fname,recfmt); 00310 return input_ped(ped,D); 00311 }

unsigned matvec::Population::input_ped0 ( matvec::Pedigree & P,

GeneticDist * D

)

Definition at line 215 of file population.cpp.
References matvec::warning().
Referenced by input_ped().

00216 { 00217 ////////////////////////////////////////////////// 00218 // return 0 if everything is successful 00219 // 1 something wrong 00220 ///////////////////////////////////////////////// 00221 if (P.size() == 0) { 00222 warning("Population::input_ped(pedigree,data): pedigree is empty"); 00223 return 1; // This has been changed from mv12 00224 } 00225 resize(P.size()+P.ngroup(),D); 00226 hashtable.resize(maxpopsize); 00227 popsize = maxpopsize; 00228 maxnamelen = P.maxnamelen; 00229 if (P.in_memory()) P.release_pedsheet(); 00230 unsigned i,s, idssex[4]; 00231 // unsigned i,s, isd[4]; old declaration 00232 00233 size_t recsize = sizeof(unsigned)*4; 00234 unsigned nanim; 00235 nanim=P.size(); 00236 Individual* I; 00237 std::ifstream pedfile(P.diskfname().c_str(),std::ios::in); 00238 //std::cout << P.diskfname().c_str() << " " << recsize << "\n"; 00239 if (!pedfile) throw exception("Population::input_ped(): cannot open file"); 00240 popmember--; 00241 for (i=1; i<=nanim; i++) { 00242 I = popmember[i]; 00243 pedfile.read((char *)idssex,recsize); 00244 I->myid = idssex[0]; 00245 if (idssex[1] > 0) { 00246 I->mymother = popmember[idssex[1]]; 00247 } 00248 else { 00249 I->mymother = 0; 00250 } 00251 if (idssex[2] > 0) { 00252 I->myfather = popmember[idssex[2]]; 00253 } 00254 else { 00255 I->myfather = 0; 00256 } 00257 I->mysex = (char)idssex[3]; 00258 // std::cout << I->myid << " "<< I->mymother << " " << I->myfather << "\n"; 00259 } 00260 for(i=(nanim+1);i<=popsize;i++) { 00261 I = popmember[i]; 00262 I->myid=i; 00263 } 00264 popmember++; 00265 if (P.maxnamelen > 0) { 00266 char *strid = new char [P.maxnamelen+1]; 00267 for (i=0; i<popsize; i++) { 00268 pedfile.read((char *)&s,sizeof(unsigned)); 00269 pedfile.read(strid,s); 00270 hashtable.insert(strid); 00271 } 00272 if(strid){ 00273 delete [] strid; 00274 strid=0; 00275 } 00276 } 00277 stdid = 1; 00278 if(!P.inbcoef_done()) P.

Public Methods
void	mark_descendant_recursive (const Individual *ind, const unsigned counter=1)
void	mark_ancestor_recursive (const Individual *ind, const unsigned counter=1)
void	mark_relative_recursive (Individual *ind, const unsigned counter=1)
void	partial_iterative_peeling (Individual *I, doubleMatrix &wspace)
void	setPenetrance (penetranceType p)
	Population (void)
	Population (GeneticDist *D)
	Population (const int maxsize, GeneticDist *D)
	Population (const Population &A)
	~Population (void)
const Population &	operator= (const Population &A)
void	copyfrom (const Population &A)
void	initialize (void)
void	inbcoef_quaas (void)
void	inbcoef_meuwissen (void)
void	confirm_sex (const int col_sex)
doubleMatrix	rela (void)
doubleMatrix	relv (const double er)
doubleMatrix	relv_nomissing (const double er)
double	cprob_op (const unsigned goff[], const unsigned gparent[], const int ori)
double	llhood_genotype (const unsigned c_id, const unsigned l_id)
double	llhood_phenotype (const unsigned nrep)
doubleMatrix	relv_inv (const double er)
void	allele_freq (const char type[], const unsigned nal,...)
void	penetrance_ml (const char fname[])
void	release (void)
void	resize (const unsigned maxn, GeneticDist *D)
void	build_offs_info (void)
void	build_spouse_info (void)
void	build_trans_mat (void)
void	renum (void)
void	display (const char key[])
unsigned	input_ped (const char fname[], const char recfmt[], GeneticDist *D)
unsigned	input_ped0 (matvec::Pedigree &P, GeneticDist *D)
unsigned	input_ped (matvec::Pedigree &P, GeneticDist *D)
unsigned	input_data (Data *D)
unsigned	input_data (const char fname[], GeneticDist *G)
unsigned	size (void) const
unsigned	nbase (void) const
unsigned	get_id (const char name[])
unsigned	input_markerData (Data *D)
unsigned	input_descentGraph (char *dgfile)
void	output_descentGraph (char *dgfile)
const char *	ind_name (const unsigned id) const
void	compute_pdm (const Individual *ind, Matrix< double > &pdm)
Vector< double >	inbcoef (void)
void	setup_m_ww (SparseMatrix &mme, Vector< double > &rhs, Vector< int > &start_addr, const double ev_e)
void	add_Ga_inv (SparseMatrix &mme, const int start_addr, const double ev_a)
void	add_Gv_inv (SparseMatrix &mme, const int start_addr, const double ev_v, const double er)
SparseMatrix *	setup_m_mme1 (Vector< double > &rhs, const double ev_v, const double ev_u, const double ev_e, const double er)
SparseMatrix *	setup_m_mme (Vector< double > &rhs, const double ev_v, const double ev_u, const double ev_e, const double er)
Individual *	member (unsigned k)
Vector< double > *	mblup1 (const double ev_v, const double ev_u, const double ev_e, const double er)
Vector< double > *	mblup (const double ev_v, const double ev_u, const double ev_e, const double er)
void	remove_mark (void)
unsigned	mark_descendant_of (const unsigned id, const unsigned counter=1)
unsigned	mark_ancestor_of (const unsigned id, const unsigned counter=1)
unsigned	mark_relative_of (const unsigned id, const unsigned counter=1)
unsigned	mark_families (const unsigned start_family_id=1)
Population *	sub (const unsigned subsize)
Population *	ancestor_of (const unsigned id)
Population *	descendant_of (const unsigned id)
Population *	relative_of (const unsigned id)
Population *	family_of (const unsigned id)
Population *	ancestor_of (const char name[])
Population *	descendant_of (const char name[])
Population *	relative_of (const char name[])
Population *	family_of (const char name[])
unsigned	fetch_families (const char filename[])
unsigned	fetch_families (std::ostream &stream)
unsigned	n_nufamily (void) const
unsigned	ind_gamete (const Individual *I, const unsigned jc, unsigned g_id[], double fq[])
double	cprob_children (const Individual *I, const unsigned jc, unsigned sg[], unsigned dg[], double sgf[], double dgf[])
unsigned	full_cdist (Individual I, const unsigned jc, unsigned cdist_value, double cdist_prob, unsigned gid_mat, double freq_mat)
unsigned	partial_cdist (Individual I, const unsigned jc, unsigned cdist_value, double cdist_prob, unsigned gid_mat, double freq_mat)
int	peeling_sequence (void)
int	detect_loop (void)
void	maxant_maxpost (void)
void	break_loop (void)
void	cond_genotype_config ()
void	genotype_config (const char type[])
void	build_pop_gamete ()
void	gibbs_iterate (const int count_genotype=0)
void	count_genotype (Individual *I)
void	release_genotype_counter (void)
void	resize_genotype_counter (void)
void	build_nufamily (void)
void	anterior (Individual *I, doubleMatrix &wspace)
void	posterior (Individual I, Individual J, doubleMatrix &wspace, const unsigned pj)
void	anterior_posterior (Individual *I, doubleMatrix &wspace)
void	genotype_dist_peeling (const int prtlevel=1, const int estifreq=1)
double	log_likelihood_peeling (const unsigned maxit=3)
double	get_posterior (Individual I, Individual exclJ, Vector< double > &vec)
double	fullsibs_prob (Individual dam, Individual sire, Individual *excludeI, doubleMatrix &wspace)
Vector< double >	get_genotype_freq (void)
void	maxant_maxpost_old (void)
NuFamily **	nufamily (void)
const DataNode *	record_ind (const unsigned i) const
void	set_switches (void)
double	multipoint_likelihood (int maxit)
double	multipoint_init_parm (Fpmm &Farg)
void	multi_geno_dist_peeling (const unsigned prtlevel)
double	multi_m_log_likelihood_peeling (int maxit)
void	multi_m_geno_dist_peeling (const unsigned prtlevel)
void	build_connected_groups (unsigned)
int	update_allele_vectors (unsigned, Individual *ind)
void	build_founder_allele_neighbors (unsigned lcs)
int	build_allele_vector (unsigned lcs, unsigned connected_group)
void	process_alleles_neighbors (unsigned lcs, unsigned founder_allele)
void	show_descent_graph_stuff (void)
double	calc_prior_descent_graph (unsigned lcs)
void	descent_graph_init_parm (void)
double	descent_graph_log_lhood ()
void	descent_graph_setup (Data *)
double	M_H_sample_ext (double l_hood, unsigned &option)
double	MH_ibd_sample (unsigned l, double log_lhood, unsigned &option)
double	calc_log_q_ratio (void)
void	show_change (const int, const int)
void	sum_descentState (unsigned)
void	sum_genotype_freq1 (unsigned)
void	sum_genotype_freq2 (unsigned)
void	sum_genotype_freq3 (unsigned)
void	output_pdq (int n_samples, char *fname)
void	Gibbs_sample (void)
void	update_ibdMatrix (unsigned locus)
void	output_ibdMatrix (unsigned n_samples, char *fname)
void	countFounders (void)
void	initJointAlleleNodeList (unsigned howManyLoci)
void	initAlleleNodeList (unsigned whichLocus)
	creates the allele node list, and the founder, penetrance, and transmission sets for each allele node
void	initGenotypeNodeList (unsigned whichLocus)
	creates the genotype node list, and the founder, penetrance, and transmission sets for each genotype node
void	getInitialGNodeListSample (unsigned maxsize, unsigned startBlock, unsigned stopBlock, unsigned sizeBlock, string samplerType)
void	getInitialGNodeListSample (unsigned maxsize, unsigned numLoci, string samplerType)
void	getOldGNodeListProbability (unsigned maxsize, unsigned startBlock, unsigned stopBlock, unsigned sizeBlock, string samplerType)
	recompute the needed quantities for the existing (old) sample (for the MH step)
void	getGNodeListSample (unsigned maxsize, unsigned startBlock, unsigned stopBlock, unsigned sizeBlock, string samplerType)
	obtain a new candidate sample using peeling and cutting (if necessary)
void	getGNodeListSample (unsigned maxsize, unsigned numLoci, string samplerType)
void	setSegregationIndex (unsigned atLocus, string samplerType)
	set the segregation indicators when these can be infered from the ordered genotypes that have been sampled
void	setupRSampler (Pedigree &P, Data D, GeneticDist G)
	method to setup the structures needed by the R-sampler
void	SimpleGenotypeElimination (unsigned locus)
	eliminates genotypes using the individuals phenotypes
void	LGGenotypeElimination (unsigned locus)
	method to do Lange and Goradia genotype elimination
void	setAlleleStateVectors (unsigned locus)
void	buildRecombinationMatrix (void)
void	finishUpAlleleStateInit (unsigned locus)
void	switchStateBack (string samplerType)
void	copyGNodeStatesToCandidateStates (string samplerType)
	saves the GNode states (alleles or genotypes) for candidate sample
void	copyCandidateToAccepted (string samplerType)
	save candidate state as accepted states; allele origins for accepted states are already in (p)m_gameteOld
void	storeSampledGametes (void)
	store sampled gametes to be able to retreive them if we retreive the gametes sampled by Gibbs if MH rejects the sample
void	retreiveSampledGametes (void)
void	initGenotypeFreq (void)
void	countHaplotypes (string samplerType)
void	countGenotypes (string samplerType)
void	displayHaplotypeFrequencies (unsigned numSamples, std::ostream &os=cout)
void	displayGenotypeFrequencies (unsigned numSamples, std::ostream &outfile=cout)
void	displaySegregationIndicators (std::ostream &outfile=cout)
unsigned	calcTotalNumberOfGenotypes (unsigned locus)
	method to do count total number of possible genotypes across the whole pedigree
void	sampleSegregationIndicators (void)
	sample the unknown segregation indicators of non-founders based on the ordered parental genotypes
void	PrintSegregationIndicators (void)
	print the segregation indicators of the non-founders
void	resizeSegregationIndicators (void)
double	MH_ibd_sample_map (double log_lhood, unsigned option)
void	sum_descentState_map (void)
void	setHaploSegregationIndex (unsigned atLocus, unsigned atMarker)
void	ListAlleleFounders (void)
void	SetPossibleHaplotypes (void)
void	SetFreqHaploFounders (void)
void	UpdateFreqHaploFounders (void)
void	CalcFreqHaploFounders (unsigned numOfSamples)
void	DisplayFreqHaploFounders (void)
Static Public Methods
void	gtindex (const unsigned num, const unsigned ni, const Vector< unsigned > &ngt, Vector< unsigned > &gtvec)
Public Attributes
Matrix< float >	ibdMatrix
int	stdid
int	mark_need_remove
int	sex_confirmed
int	offs_info_built
int	spouse_info_built
int	fdone
bool	nuFamiliesDone
double	scaling_val
NuFamily **	nufamily_vec
unsigned	num_nufamily
unsigned	num_t_nufamily
unsigned	nonloop_t_nufamily
unsigned	tn_genotype
unsigned	tn_gamete
unsigned	maxnamelen
unsigned	tn_qtl
unsigned	n_markerLoci
doubleMatrix *	trans_mat
unsigned	popsize
unsigned	maxpopsize
unsigned	numchrom
unsigned	numtrait
Individual *	pm_storage
Individual **	popmember
Vector< double >	blupsol
std::string	evaluate_method
Genome *	pop_gamete
Genome *	gametebase
GeneticDist *	prior
int	nmarker
std::string *	markersymbol
penetranceType	penetrance_f
SparseMatrix	hmmec
HashTable	hashtable
doubleMatrix	residual_var
Vector< double >	mean_for_genotype
Matrix< int >	switch_table
Matrix< int >	switch_table_gmt
doubleMatrix	switch_table_prob
Vector< double >	Q_freq
Vector< Vector< double > >	marker_freq
Vector< Vector< double > >	locusFreq
std::string	analysis_type
Vector< double >	mean_for_pgenotype
Vector< double >	P_freq
int	P_ndim
Fpmm *	F
int	founder_allele_counter
Vector< unsigned >	connected_groups
Vector< int >	allele_vector1
Vector< int >	allele_vector2
unsigned	connect_counter
Vector< double >	RecoVector
unsigned	SLlocus
Model *	model
int	valid_graph
std::vector< std::vector< int > >	founder_allele_neighbors
Matrix< unsigned >	re_pdq
Vector< unsigned >	temp
double	log_lhood
unsigned	initial
unsigned	numFounders
Matrix< double >	recombinationMatrix
SafeSTLVector< GNodeSet * >	vectorOfGNsts
GNodeList	gNodeList
bool	lookOnlyAtLocusToYourLeft
IndexVector	haplotypeCoder
unsigned	numHaplotypes
bool	markerDataIn
Static Public Attributes
double	disPenetranceTable [2][3]
Friends
class	Model
class	GLMM
class	Individual
class	NuFamily