individual.cpp Source File

00001 //****************************************************
00002 //  April, 1993, University of Illinois
00003 // Copyright (C) 1993, 1994 Tianlin Wang
00004 /* Copyright (C) 1994-2003 Matvec Development Team. 
00005 
00006 This program is free software; you can redistribute it and/or
00007 modify it under the terms of the GNU Library General Public
00008 License as published by the Free Software Foundation; either
00009 version 2 of the License, or (at your option) any later version.
00010 
00011 This program is distributed in the hope that it will be useful,
00012 but WITHOUT ANY WARRANTY; without even the implied warranty of
00013 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00014 Library General Public License for more details.
00015 
00016 You should have received a copy of the GNU Library General Public
00017 License along with this library; if not, write to the Free
00018 Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
00019 MA 02111-1307, USA 
00020 */
00021 
00022 #include <iomanip>
00023 #include "population.h"
00024 #include "chromosome.h"
00025 #include "individual.h"
00026 #include "model.h"
00027 namespace matvec {
00028         //BRS: Declaration of static variables
00029         Population *Individual::population;
00030         doubleMatrix Individual::g_weight;
00031         std::vector <std::vector<double> > Individual::vec_prob;
00032         std::vector <std::vector<double> > Individual::vec_cutsetval;
00033         SafeSTLVector<unsigned> Individual::QTLPosVector;
00034         unsigned Individual::numLoci;
00035         unsigned Individual::currentLocus;
00036         unsigned Individual::nQTL;
00037         double Individual::isqrt2pi = 1.0/std::sqrt(6.283185308);
00038         extern double ranf(void);
00039         
00040         int compare_ind_pt(const void* x, const void* y)
00041         {
00042                 Individual **U1,**U2;
00043                 U1 = (Individual **)x;
00044                 U2 = (Individual **)y;
00045                 int retval = 1;
00046                 if (*U1 < *U2)
00047                         retval = -1;
00048                 else if (*U1 == *U2) {
00049                         retval = 0;
00050                 }
00051                 return retval;
00052         }
00053         
00054         int compare_ind_id(const void* a, const void* b)
00055                 //  sort the popmember table which contains just address of each individal
00056                 // in popmember table, it isn't time consuming
00057         {
00058                 Individual **x,**y;
00059                 x = (Individual **)a; y = (Individual **)b;
00060                 return ((*x)->id() - (*y)->id());
00061         }
00062         
00063         int compare_ind_gid(const void* a, const void* b)
00064                 //  sort the popmember table which contains just address of each individal
00065                 // in popmember table, it isn't time consuming
00066         {
00067                 Individual **x,**y;
00068                 x = (Individual **)a; y = (Individual **)b;
00069                 return ((*x)->gid() - (*y)->gid());
00070         }
00071         
00072         int compare_mother_id(const void* x, const void* y)
00073         {
00074                 Individual **U1,**U2;
00075                 U1 = (Individual **)x;
00076                 U2 = (Individual **)y;
00077                 return ((*U1)->mother_id() - (*U2)->mother_id());
00078         }
00079         
00080         int compare_father_id(const void* x, const void* y)
00081         {
00082                 Individual **U1,**U2;
00083                 U1 = (Individual **)x;
00084                 U2 = (Individual **)y;
00085                 return ((*U1)->father_id() - (*U2)->father_id());
00086         }
00087         
00088         Individual::Individual(void)
00089         {
00090                 population       = 0;
00091                 genotype_id      = -1;
00092                 myid             = 0;
00093                 mysex            = '.';  // Changed from zero
00094                 numtrait         = 0;
00095                 numchrom         = 0;
00096                 myrecord         = 0;
00097                 posterior        = 0;
00098                 p_origin         = 1;
00099                 RBV              = 0.0;
00100                 MBV              = 0.0;
00101                 inbc             = 0.0;
00102                 est_GV           = 0.0;
00103                 true_GV          = 0.0;
00104                 xbzu_val         = 0.0;
00105                 group_id         = 0;
00106                 spouselist       = 0;
00107                 myoffspring      = 0;
00108                 residual_var     = 0;
00109                 loop_connector   = 0;
00110                 connect_keeper   = 0;
00111                 numoffs_spouse   = 0;
00112                 genotype_counter = 0;
00113                 marker_index     = 0;
00114                 assigned_founder_allele = 0;
00115         }
00116         
00117         Individual::Individual(Population *P)
00118         {
00119                 if (population == P) return;
00120                 population   = P;
00121                 genotype_id  = -1;
00122                 myid         = 0;
00123                 mysex        = '.';
00124                 myrecord     = 0;
00125                 posterior    = 0;
00126                 p_origin     = 1;
00127                 RBV          = 0.0;
00128                 MBV          = 0.0;
00129                 inbc         = 0.0;
00130                 est_GV       = 0.0;
00131                 true_GV      = 0.0;
00132                 xbzu_val     = 0.0;
00133                 group_id     = 0;
00134                 
00135                 spouselist   = 0;
00136                 myoffspring  = 0;
00137                 residual_var = 0;
00138                 loop_connector = 0;
00139                 connect_keeper = 0;
00140                 numoffs_spouse = 0;
00141                 genotype_counter = 0;
00142                 marker_index     = 0;
00143                 assigned_founder_allele = 0;
00144                 remodel(P);
00145         }
00146         
00147         Individual::Individual(const Individual& A)
00148         {
00149                 numoffs        = 0;
00150                 myrecord       = 0;
00151                 posterior      = 0;
00152                 myoffspring    = 0;
00153                 spouselist     = 0;
00154                 numoffs_spouse = 0;
00155                 copyfrom(A);
00156                 assigned_founder_allele = 0;
00157         }
00158         
00159         void Individual::copyfrom(const Individual& A)
00160         {
00161                 if (this == &A) return;
00162                 population  = A.population;
00163                 genotype_id = A.genotype_id;
00164                 myid      = A.myid;
00165                 mysex       = A.mysex;                inbc      = A.inbc;
00166                 myfather    = A.myfather;             mymother  = A.mymother;
00167                 MBV         = A.MBV;                  RBV       = A.RBV;
00168                 true_GV     = A.true_GV;              est_GV    = A.est_GV;
00169                 xbzu_val    = A.xbzu_val;             numtrait  = A.numtrait;
00170                 numspouse   = A.numspouse;            numoffs   = A.numoffs;
00171                 group_id    = A.group_id;             offs_tree = A.offs_tree;
00172                 p_origin    = A.p_origin;             numchrom  = A.numchrom;
00173                 genome0     = A.genome0;              genome1   = A.genome1;
00174                 loop_connector = A.loop_connector;  anterior  = A.anterior;
00175                 residual_var = A.residual_var;      related_family = A.related_family;
00176                 offs_tree = A.offs_tree;            connect_keeper = A.connect_keeper;
00177                 marker_index = A.marker_index;
00178                 assigned_founder_allele = A.assigned_founder_allele;
00179                 //BRS
00180                 index_sw=A.index_sw;
00181                 bet_sw=A.bet_sw;
00182                 eps_sw=A.eps_sw;
00183                 family=A.family;
00184                 gprobs=A.gprobs;
00185                 m_anterior=A.m_anterior;
00186                 m_posterior=A.m_posterior;
00187                 m_anterior_scale=A.m_anterior_scale;
00188                 m_posterior_scale=A.m_posterior_scale;
00189                 //BRS
00190                 
00191                 unsigned i;
00192                 if (posterior) { delete [] posterior; posterior=0;}
00193                 if (A.posterior) {
00194                         if (numspouse>0){
00195                                 posterior = new Vector<double> [numspouse];
00196                         }
00197                         else {
00198                                 posterior = 0;
00199                         }
00200                         for (i=0; i<numspouse; i++) posterior[i] = A.posterior[i];
00201                 }
00202                 if (myrecord) {delete [] myrecord; myrecord = 0;}
00203                 if (A.myrecord) {
00204                         if(numtrait>0){
00205                                 myrecord = new DataNode [numtrait];
00206                         }
00207                         else {
00208                                 myrecord = 0;
00209                         }
00210                         for (i=0; i<numtrait; i++) myrecord[i] = A.myrecord[i];
00211                 }
00212                 if (spouselist) {delete [] spouselist; spouselist= 0;}
00213                 if (A.spouselist) {
00214                         if (numspouse>0){
00215                                 spouselist = new Individual* [numspouse];
00216                         }
00217                         else {
00218                                 spouselist = 0;
00219                         }
00220                         for (i=0; i<numspouse; i++) spouselist[i] = A.spouselist[i];
00221                 }
00222                 if (myoffspring) {delete [] myoffspring; myoffspring=0;}
00223                 if (A.myoffspring) {
00224                         if(numoffs>0){
00225                                 myoffspring = new Individual* [numoffs];
00226                         }
00227                         else {
00228                                 myoffspring = 0;
00229                         }
00230                         for (i=0; i<numoffs; i++) myoffspring[i] = A.myoffspring[i];
00231                 }
00232                 if (numoffs_spouse) {
00233                         delete [] numoffs_spouse; numoffs_spouse = 0;
00234                 }
00235                 if (A.numoffs_spouse) {
00236                         if(numspouse>0){
00237                                 numoffs_spouse = new unsigned [numspouse];
00238                         }
00239                         else {
00240                                 numoffs_spouse = 0;
00241                         }
00242                         for (i=0; i<numspouse; i++) numoffs_spouse[i] = A.numoffs_spouse[i];
00243                 }
00244                 if (A.genotype_counter) {
00245                         if(numchrom>0){
00246                                 genotype_counter = new Vector<double> [numchrom];
00247                         }
00248                         else {
00249                                 genotype_counter = 0;
00250                         }
00251                         for (i=0; i<numchrom; i++) genotype_counter[i] = A.genotype_counter[i];
00252                 }
00253                 else {
00254                         genotype_counter = 0;
00255                 }
00256         }
00257         
00258         void Individual::remodel(Population *P)
00259         {
00260                 population   = P;
00261                 genotype_id  = -1;
00262                 numtrait = population->prior->ntrait();
00263                 numchrom  = population->prior->nchrom();
00264                 
00265                 genome0.remodel(population->prior);
00266                 genome1.remodel(population->prior);
00267                 
00268                 if (myrecord) {
00269                         delete [] myrecord;
00270                         myrecord=0;
00271                 }
00272                 if(numtrait>0){
00273                         myrecord = new DataNode [numtrait];
00274                 }
00275                 else {
00276                         myrecord = 0;
00277                 }
00278                 
00279                 myfather = 0;
00280                 mymother = 0;
00281                 if (numoffs_spouse) {
00282                         delete [] numoffs_spouse; numoffs_spouse = 0;
00283                 }
00284                 numspouse = 0;
00285                 if (myoffspring) { delete [] myoffspring; myoffspring = 0;}
00286                 numoffs  = 0;
00287                 RBV          = 0.0;
00288                 MBV          = 0.0;
00289                 inbc         = 0.0;
00290                 mysex        = '.';
00291                 est_GV       = 0.0;
00292                 true_GV      = 0.0;
00293                 group_id     = 0;
00294                 xbzu_val     = 0.0;
00295                 genotype_counter = 0;
00296                 if (posterior) {delete [] posterior; posterior=0;}
00297                 anterior_iw = 0;
00298                 anterior.resize(0);
00299         }
00300         
00301         const Individual& Individual::operator=(const Individual& A)
00302         {
00303                 copyfrom(A);
00304                 return *this;
00305         }
00306         
00307         unsigned Individual::father_id(void) const
00308         {
00309                 if (myfather) return myfather->myid;
00310                 else  return 0;
00311         }
00312         
00313         unsigned Individual::mother_id(void) const
00314         {
00315                 if (mymother) return mymother->myid;
00316                 else return 0;
00317         }
00318         
00319         unsigned Individual::father_gid(void) const
00320         {
00321                 if (myfather) return myfather->group_id;
00322                 else  return 0;
00323         }
00324         
00325         unsigned Individual::mother_gid(void) const
00326         {
00327                 if (mymother) return mymother->group_id;
00328                 else  return 0;
00329         }
00330         
00331         double Individual::father_inbcoef(void) const
00332         {
00333                 if (myfather) return myfather->inbcoef();
00334                 else  return -1.0;
00335         }
00336         
00337         double Individual::mother_inbcoef(void) const
00338         {
00339                 if (mymother) return mymother->inbcoef();
00340                 else  return -1.0;
00341         }
00342         
00343         void Individual::initial_anterior(const double *gfreq,const unsigned tng,
00344                                                                           const int cond)
00345         {
00346                 /////////////////////////////////////////////////////////////////////////
00347                 // Conditional initialization for anterior
00348                 // cond = -1  is the default, means initialize anterior un-conditionally
00349                 // cond = 1, associated with extra data resulting from iterating
00350                 /////////////////////////////////////////////////////////////////////
00351                 
00352                 if (cond == anterior_iw) return;      // anterior_iw is the condition
00353                 unsigned j;
00354                 if (anterior.size() != tng+1) anterior.resize(tng+1);
00355                 if (myrecord[0].missing) {
00356                         for (j=0; j<tng; j++) anterior[j] = gfreq[j];
00357                         anterior[tng] = 0.0;
00358                         return;
00359                 }
00360                 double s;
00361                 Vector<double> pen_vec(tng);
00362                 anterior[tng] = get_penetrance(pen_vec);
00363                 for (s=0.0,j=0; j<tng; j++) s += (anterior[j] = pen_vec[j]*gfreq[j]);
00364                 for (j=0; j<tng; j++) anterior[j] /= s;
00365                 anterior[tng] += std::log(s);
00366         }
00367         
00368         void Individual::initial_posterior(const unsigned tng, const int cond)
00369         {
00370                 /////////////////////////////////////////////////////////////////////////
00371                 // Conditional initialization for posterior
00372                 // cond = -1  is the default, means initialize anterior un-conditionally
00373                 // cond = 1, associated with extra data resulting from iterating
00374                 /////////////////////////////////////////////////////////////////////
00375                 int i;
00376                 if (posterior) {
00377                         for (i=0; i<numspouse; i++) {
00378                                 if (posterior_iw[i] != cond) {
00379                                         posterior[i].resize(tng+1,1.0);
00380                                         posterior[i][tng] = 0.0;
00381                                 }
00382                         }
00383                 }
00384                 else {
00385                         if(numspouse>0){
00386                                 posterior = new Vector<double> [numspouse];
00387                                 posterior_iw.resize(numspouse,0.0); //BRS
00388                         }
00389                         else {
00390                                 posterior = 0;
00391                         }
00392                         for (i=0; i<numspouse; i++) {
00393                                 posterior[i].resize(tng+1,1.0);
00394                                 posterior[i][tng] = 0.0;
00395                         }
00396                 }
00397         }
00398         
00399         void Individual::pretend_missing(int on,const double *gfreq,const unsigned tng)
00400         {
00401                 if (on) {
00402                         if (p_origin) {
00403                                 record()[0].missing += 1;
00404                                 initial_anterior(gfreq,tng);
00405                                 if (connect_keeper == 1) {
00406                                         initial_posterior(tng,1);
00407                                 }
00408                                 else {
00409                                         initial_posterior(tng);
00410                                 }
00411                         }
00412                         else {
00413                                 initial_posterior(tng,1);
00414                         }
00415                 }
00416                 else {
00417                         if (p_origin) record()[0].missing -= 1;
00418                 }
00419         }
00420         
00421         //RLF modified to work with pop->mean_for_genotype[i]
00422         
00423         double Individual::get_penetrance(Vector<double> &pen)
00424         {
00425                 ////////////////////////////////////////////////////////////////
00426                 // pen should have enough space to hold all penetrance values
00427                 // this works for multi-chromosomes, but not for multi-loci.
00428                 // For two chromosomes with two alleles each, there are 10 genotypes
00429                 //     AB Ab aB ab
00430                 // AB  0
00431                 // Ab  1  2
00432                 // aB  3  4  5
00433                 // ab  6  7  8  9
00434                 ////////////////////////////////////////////////////////////////
00435                 unsigned i,j,t;
00436                 unsigned tn_genotype = pen.size();
00437                 if (myrecord[0].missing) {
00438                         for (i=0; i<tn_genotype; i++) pen[i] = 1.0;
00439                         return 0.0;
00440                 }
00441                 
00442                 //*
00443                 
00444                  unsigned tn_gamete = static_cast<unsigned>(0.5*std::sqrt(static_cast<double>(8*tn_genotype + 1)) - 0.5);
00445                  for (genotype_id=0,i=0; i<tn_gamete; i++) {
00446                          for (j=0; j<=i; j++) {
00447                                  pen[genotype_id] = (*penetrance_f)(this,(const double**)residual_var->begin());
00448                                  genotype_id++;
00449                          }
00450                  }
00451                  
00452                  //RLF
00453                  double inverse_sqrt2pi = 1.0/std::sqrt(4.0*std::asin(1.0));   // pi =2.0*std::asin(1.0)
00454                  double mu, v = *residual_var[0][0];
00455                  double y = record()[0].double_val();
00456                  
00457                  for (i=0; i < tn_genotype; i++) {
00458                          mu = population->mean_for_genotype[i];
00459                          pen[i] = inverse_sqrt2pi/std::sqrt(v) * std::exp(-(y-mu)*(y-mu)/(2*v));
00460                          
00461                  }
00462                  //RLF
00463                  
00464                  double scale = pen.sum();
00465                  if (fabs(scale) < 1.0e-8) {
00466                          warning("this record is not conformable:");
00467                          this->display();
00468                          scale = 3.0;
00469                          for (i=0; i<tn_genotype; i++) pen[i] = 1.0;
00470                  }
00471                  else {
00472                          for (i=0; i<tn_genotype; i++) pen[i] /= scale;
00473                  }
00474                  return std::log(scale);
00475         }
00476                  //RLF
00477                  
00478                  void Individual::gamete(Chromosome *C,const int n,const double r) const
00479                  {
00480                          /************************************************
00481                          * C is non marked genome,  D is marked chromosome
00482                          ****************************************************/
00483                          if (numchrom != n) throw exception("Individual::gamete(): bad arg");
00484                          const Chromosome *A1, *A2;
00485                          A1 = genome0.chromosome;
00486                          A2 = genome1.chromosome;
00487                          for (unsigned i=0; i<numchrom; i++) {
00488                                  /*
00489                                   if (ranf() <= 0.5) {
00490                                           C[i].qtl = A1[i].qtl;
00491                                           if (ranf() <= r) {
00492                                                   C[i].marker = A2[i].marker;
00493                                           }
00494                                           else {
00495                                                   C[i].marker = A1[i].marker;
00496                                           }
00497                                   }
00498                                   else {
00499                                           C[i].qtl = A2[i].qtl;
00500                                           if (ranf() <= r) {
00501                                                   C[i].marker = A1[i].marker;
00502                                           }
00503                                           else {
00504                                                   C[i].marker = A2[i].marker;
00505                                           }
00506                                   }
00507                                   */
00508                          }
00509                  }
00510                  
00511                  void Individual::release(void)
00512                  {
00513                          if (myrecord)    {delete [] myrecord; myrecord = 0;}
00514                          if (spouselist)  {delete [] spouselist; spouselist = 0;}
00515                          if (posterior)   {delete [] posterior; posterior = 0; }
00516                          if (myoffspring) {delete [] myoffspring; myoffspring = 0;}
00517                          if (numoffs_spouse) {
00518                                  delete [] numoffs_spouse; numoffs_spouse = 0;
00519                          }
00520                          if (genotype_counter) {
00521                                  delete [] genotype_counter; genotype_counter = 0;
00522                          }
00523                  }
00524                  
00525                  double Individual::genotypic_val(void) const
00526                  {
00527                          /////////////////////////////////////////////////////////////////////
00528                          //  assuming there are two alleles [0, 1],
00529                          //  then there are a total of 4 genotypes
00530                          //    ----------------------
00531                          //      |    0        1
00532                          //   -----------------------
00533                          //    0 |  (0 0)   (0 1)
00534                          //    1 |  (1 0)   (1 1)
00535                          //  ------------------------
00536                          //    allele id must be an integer starting from 0
00537                          ///////////////////////////////////////////////////////////////////////
00538                          double retval = 0;
00539                          if (genotype_id >= 0) {
00540                                  retval = population->mean_for_genotype[genotype_id];
00541                          }
00542                          else {
00543                                  warning("Individual::genotypic_val(): no genotype available");
00544                          }
00545                          return retval;
00546                          /*
00547                           unsigned i,j,nc, nl,sd[2];
00548                           nc = prior->nchrom();
00549                           ChromStruct *Chrom = prior->chrom();
00550                           const double** gv;
00551                           double retval = 0;
00552                           for (i=0; i<nc; i++) {
00553                                   nl = Chrom[i].nloci();
00554                                   for (j=0; j<nl; j++) {
00555                                           genotype(i,j,sd);
00556                                           gv = prior->genotypic_val(i,j);
00557                                           retval +=  gv[sd[0]][sd[1]];
00558                                   }
00559                           }
00560                           return retval;
00561                           */
00562                  }
00563                  
00564                  void Individual::set_genotype(const unsigned c,const unsigned l,
00565                                                                            const unsigned a0,const unsigned a1)
00566                  {
00567                          genome0.chromosome[c].locus[l].allele = a0;
00568                          genome1.chromosome[c].locus[l].allele = a1;
00569                  }
00570                  
00571                  void Individual::genotype(const unsigned c, const unsigned l,unsigned gtype[])
00572                  const
00573                  {
00574                          gtype[0] = genome0.chromosome[c].locus[l].allele;
00575                          gtype[1] = genome1.chromosome[c].locus[l].allele;
00576                  }
00577                  
00578                  void Individual::save(std::ostream& out)
00579                  {
00580                          out << "  I " << myid << "   F " << father_id() << "   M " << mother_id();
00581                          out << " Marker Info: ";
00582                          out << genome0.chromosome[0].locus[1].allele << " " ;
00583                          out << genome1.chromosome[0].locus[1].allele << " " ;
00584                          out << genome0.chromosome[0].locus[2].allele << " " ;
00585                          out << genome1.chromosome[0].locus[2].allele;
00586                          //out << "   Switch " << index_sw << "   Beta " << bet_sw << "   Epl " << eps_sw;
00587                          out << std::endl;
00588                  }
00589                  
00590                  //BRS
00591                  void Individual::set_switch(int Nloci)
00592                  
00593                  {
00594                          /* This assigns the switches and gametes to each individual in the population.
00595                          By definition, the maternal gamete gets the lower allele number unless it is known
00596                          which source the allele came from.
00597                          
00598                          An locus is switchable only in two circumstances:
00599                          1)  founder animal is heterozgyous at that locus.
00600 2)  both parents are heterozygous for the same alleles at that locus
00601 Under all other cases, the source of the allele can be tracked.
00602 Thus, each locus needs to be checked to determine the source of the allele.
00603 
00604 */
00605                          int i,j,k,l,allele1,allele2,sire,dam,sireF,sireM,damF,damM,levelS=1,levelEB=1,Ixswitch=0,epl=0,beta=0,temp,tallele;
00606                          if (myfather && (mymother == 0)) {
00607                                  std::cout << "error father known, mother not! for " << population->ind_name(id()) << std::endl;
00608                                  std::cout << "This program is not designed to handle this yet!\n Please supply a mother with records." << std::endl;
00609                                  exit(1);
00610                          }
00611                          if (mymother && (myfather == 0)) {
00612                                  std::cout << "error mother known, father not! for " << population->ind_name(id()) << std::endl;
00613                                  std::cout << "This program is not designed to handle this yet! Please supply a father with records." << std::endl;
00614                                  exit(1);
00615                          }
00616                          
00617                          if (myfather == 0)  { // animal is a founder
00618                                  for (j=(Nloci); j > 0; j--)  {
00619                                          if (genome0.chromosome[0].locus[j].allele > genome1.chromosome[0].locus[j].allele) {  // swap alleles
00620                                                  temp=genome0.chromosome[0].locus[j].allele;
00621                                                  genome0.chromosome[0].locus[j].allele=genome1.chromosome[0].locus[j].allele;
00622                                                  genome1.chromosome[0].locus[j].allele =temp;
00623                                          }
00624                                          if (genome0.chromosome[0].locus[j].allele != genome1.chromosome[0].locus[j].allele) {   // switchable
00625                                                  Ixswitch += levelS;
00626                                          }
00627                                          levelS *=2;
00628                                  }
00629                                  index_sw=Ixswitch;
00630                          }
00631                          else { // not a founder so find all the alleles from individual and parents
00632                                  for (j=(Nloci); j>0; j--)  {
00633                                          allele1=genome0.chromosome[0].locus[j].allele;
00634                                          allele2=genome1.chromosome[0].locus[j].allele;
00635                                          sireM=myfather->genome0.chromosome[0].locus[j].allele;
00636                                          sireF=myfather->genome1.chromosome[0].locus[j].allele;
00637                                          damM=mymother->genome0.chromosome[0].locus[j].allele;
00638                                          damF=mymother->genome1.chromosome[0].locus[j].allele;
00639                                          if ((allele1 != allele2) && (((sireF == damF) && (sireM == damM)) || ((sireF == damM) && (sireM == damF)))) {
00640                                                  Ixswitch += levelS; // Locus is switchable !!!
00641                                                                                          // Check if alleles in correct order if not swap them
00642                                                  if (allele1 > allele2) {
00643                                                          genome1.chromosome[0].locus[j].allele=allele1;
00644                                                          genome0.chromosome[0].locus[j].allele=allele2;
00645                                                  }
00646                                          }
00647                                          // Not switchable so determine source of each allele
00648                                          else if ((allele1 == sireM) && (allele2 == damM)) {
00649                                                  genome0.chromosome[0].locus[j].allele=allele2;
00650                                                  genome1.chromosome[0].locus[j].allele=allele1;
00651                                          }
00652                                          else if ((allele1 == sireM) && (allele2 == damF))  {
00653                                                  genome0.chromosome[0].locus[j].allele=allele2;
00654                                                  genome1.chromosome[0].locus[j].allele=allele1;
00655                                          }
00656                                          else if ((allele1 == sireF) && (allele2 == damM)) {
00657                                                  genome0.chromosome[0].locus[j].allele=allele2;
00658                                                  genome1.chromosome[0].locus[j].allele=allele1;
00659                                          }
00660                                          else if ((allele1 == sireF) && (allele2 == damF))  {
00661                                                  genome0.chromosome[0].locus[j].allele=allele2;
00662                                                  genome1.chromosome[0].locus[j].allele=allele1;
00663                                          }
00664                                          else if ((allele2 == sireM) && (allele1 == damM))  {
00665                                                  genome0.chromosome[0].locus[j].allele=allele1;
00666                                                  genome1.chromosome[0].locus[j].allele=allele2;
00667                                          }
00668                                          else if ((allele2 == sireM) && (allele1 == damF))  {
00669                                                  genome0.chromosome[0].locus[j].allele=allele1;
00670                                                  genome1.chromosome[0].locus[j].allele=allele2;
00671                                          }
00672                                          else if ((allele2 == sireF) && (allele1 == damM))  {
00673                                                  genome0.chromosome[0].locus[j].allele=allele1;
00674                                                  genome1.chromosome[0].locus[j].allele=allele2;
00675                                          }
00676                                          else if ((allele2 == sireF) && (allele1 == damF))  {
00677                                                  genome0.chromosome[0].locus[j].allele=allele1;
00678                                                  genome1.chromosome[0].locus[j].allele=allele2;
00679                                          }
00680                                          else    {
00681                                                  std::cerr << "ERROR ident alleles different from parents - perhaps parent not in pedigree?\n"
00682                                                  << "Id " <<  allele1 << " " << allele2 << " "  << sireM  << " " << sireF  << " " << damM << " "  << damF << std::endl;
00683                                                  throw exception("ERROR ident alleles different from parents - perhaps parent not in pedigree?");
00684                                          }
00685                                          // Now determine epsilon
00686                                          if (sireM == sireF)
00687                                                  epl += (levelEB*2);
00688                                          else if ( genome1.chromosome[0].locus[j].allele == sireF)
00689                                                  epl += levelEB; // Paternal allele from father's father
00690                                                                                  // Now determine  beta
00691                                          if (damM == damF)
00692                                                  beta += (levelEB*2);
00693                                          else if (genome0.chromosome[0].locus[j].allele == damF)
00694                                                  beta += levelEB; // maternal allele from dam's father
00695                                          
00696                                          levelS *=2;
00697                                          levelEB *=3;
00698                                  }
00699                          }
00700                          index_sw=Ixswitch;
00701                          eps_sw=epl;
00702                          bet_sw=beta;
00703                          //std::cout << "Switch " << Ixswitch << " Epli " << epl << " beta " << beta << std::endl;
00704 }
00705 
00706 unsigned Individual::n_switches(void){
00707         return (population->switch_table[index_sw][0]);
00708 }
00709 
00710 
00711 void Individual::initial_multi_anterior(doubleMatrix& penetrance)
00712 {
00713         unsigned i,j,k,a1,a2,ndim;
00714         double freq, scale, sum=0.0, freq_mark=1.0;
00715         unsigned n_switch = population->switch_table[index_sw][0];
00716         ndim=population->P_ndim;
00717         if (m_anterior.get_nrow() != n_switch) m_anterior.resize(ndim,n_switch,4,0.0);
00718         m_anterior_iw = 0;
00719         
00720         for (k=1;k<=population->n_markerLoci; k++){
00721                 a1 = genome0.chromosome[0].locus[k].allele;
00722                 a2 = genome1.chromosome[0].locus[k].allele;
00723                 freq_mark *= population->marker_freq(k)(a1)*population->marker_freq(k)(a2);
00724         }
00725         if (myrecord[0].missing) {
00726                 m_anterior_scale =0.0;
00727                 for (k=0;k<ndim; k++){
00728                         for (j=0; j<4; j++) {
00729                                 freq = (population->Q_freq[j])*freq_mark*(population->P_freq[k]);
00730                                 for (i=0; i<n_switch; i++) {
00731                                         m_anterior[k][i][j] = freq;
00732                                         sum += freq; // accumulating scaling factor
00733                                 }
00734                         }
00735                 }
00736         }
00737         else {
00738                 //    std::cout << " before call Pen " <<  std::endl;
00739                 m_anterior_scale = get_m_penetrance(penetrance);
00740                 //   std::cout << "scale in intial_anterior "<< m_anterior_scale << std::endl;
00741                 //    std::cout << " Pen " << penetrance << std::endl;
00742                 for (k=0;k<ndim; k++){
00743                         for (j=0; j<4; j++) {
00744                                 freq = (population->Q_freq[j])*freq_mark*(population->P_freq[k])*penetrance[k][j];
00745                                 //  std::cout << "k " << k << " j " << j << " Q_f " <<  population->Q_freq[j];
00746                                 //   std::cout << " M_f " << freq_mark << " Pen " << penetrance[k][j] << " freq " << freq << std::endl;
00747                                 for (i=0; i<n_switch; i++) {
00748                                         m_anterior[k][i][j] = freq;
00749                                         sum += freq; // accumulating scaling factor
00750                                 }
00751                         }
00752                 }
00753         }
00754         //  std::cout << "sum " << sum << std::endl;
00755         
00756         // Need to rescale due to probability of observing the marker genotypes
00757         for (k=0;k<ndim; k++){
00758                 for (i=0; i<n_switch; i++) {
00759                         for (j=0; j<4; j++) {
00760                                 m_anterior[k][i][j] /= sum;
00761                         }
00762                 }
00763         }
00764         m_anterior_scale += std::log(sum);
00765         //  std::cout << "scale " << m_anterior_scale << std::endl;
00766 }
00767 
00768 void Individual::initial_multi_posterior(const int cond) {
00769         /////////////////////////////////////////////////////////////////////////
00770         // Conditional initialization for posterior
00771         // cond = -1  is the default, means initialize anterior un-conditionally
00772         // cond = 1, associated with extra data resulting from iterating
00773         /////////////////////////////////////////////////////////////////////
00774         int i,ndim;
00775         unsigned n_switch = population->switch_table[index_sw][0];
00776         ndim=population->P_ndim;
00777         if (numspouse) {
00778                 if (m_posterior.size() != numspouse ) {
00779                         m_posterior.resize(numspouse);
00780                         m_posterior_scale.resize(numspouse,0.0);
00781                         posterior_iw.resize(numspouse,0.0); //BRS
00782                 }
00783                 if (posterior){
00784                         m_posterior_scale.resize(numspouse,0.0);
00785                         posterior_iw.resize(numspouse,0.0); //BRS
00786                         for (i=0; i<numspouse; i++) {
00787                                 if (posterior_iw[i] != cond) {
00788                                         m_posterior[i].resize(ndim,n_switch,4,1.0);
00789                                 }
00790                         }
00791                 }
00792                 else {
00793                         for (i=0; i<numspouse; i++) {
00794                                 m_posterior[i].resize(ndim,n_switch,4,1.0);
00795                         }
00796                         m_posterior_scale.resize(numspouse,0.0);
00797                         posterior_iw.resize(numspouse,0.0); //BRS
00798                 }
00799         }
00800 }
00801 
00802 double Individual::get_m_penetrance(doubleMatrix& pen) {
00803         
00804         unsigned QTL,PGN;
00805         
00806         if (myrecord[0].missing) {
00807                 for (PGN=0; PGN<population->P_ndim; PGN++) {
00808                         for (QTL=0; QTL<4; QTL++) {
00809                                 pen[PGN][QTL] = 1.0;
00810                         }
00811                 }
00812                 return 0.0;
00813         }
00814         
00815         double inverse_sqrt2pi = 1.0/std::sqrt(4.0*std::asin(1.0));  // pi = 2.0*std::asin(1.0)
00816         double mu, v = *residual_var[0][0],scale;
00817         double y = record()[0].double_val();
00818         
00819         for (PGN=0; PGN<population->P_ndim; PGN++) {
00820                 for (QTL=0; QTL<4; QTL++) {
00821                         mu = population->mean_for_genotype[QTL]+population->mean_for_pgenotype[PGN];
00822                         pen[PGN][QTL] = inverse_sqrt2pi/std::sqrt(v) * std::exp(-(y-mu)*(y-mu)/(2*v));
00823                         //  std::cout << QTL << " Q_mu " << population->mean_for_genotype[QTL] << " P_mu " << population->mean_for_pgenotype[PGN] <<" mu " << mu << " y " << y << " v " << v << " pen is " << setprecision(12) << pen[PGN][QTL] << std::endl;
00824                 }
00825         }
00826         //   std::cout << "pen in get_m_pen " << pen;
00827         scale = (pen.sum()).sum();
00828         //  std::cout << "Pen scale " << scale << std::endl;
00829         if (fabs(scale) < 1.0e-8) {
00830                 std::cout << id() << " has record " << y << std::endl;
00831                 warning("this record is not conformable");
00832                 scale = 3.0;
00833                 for (PGN=0; PGN<population->P_ndim; PGN++) {
00834                         for (QTL=0; QTL<4; QTL++)
00835                                 pen[PGN][QTL] = 1.0;
00836                 }
00837         }
00838         else {
00839                 for (PGN=0; PGN<population->P_ndim; PGN++) {
00840                         for (QTL=0; QTL<4; QTL++) {
00841                                 pen[PGN][QTL] /= scale;
00842                         }
00843                 }
00844                 //    std::cout << "pen in get pen " << pen << std::endl;
00845         }
00846         return std::log(scale);
00847 }
00848 
00849 void Individual::pretend_multi_missing(int on, doubleMatrix& pen)
00850 {
00851         if (on) {
00852                 if (p_origin) {
00853                         record()[0].missing += 1;
00854                         initial_multi_anterior(pen);
00855                         //     initial_anterior(gfreq,tng);
00856                         if (connect_keeper == 1) {
00857                                 initial_multi_posterior(1);
00858                                 // initial_posterior(tng,1);
00859                         }
00860                         else {
00861                                 initial_multi_posterior(-1);
00862                                 // initial_posterior(tng);
00863                         }
00864                 }
00865                 else {
00866                         initial_multi_posterior(1);
00867                         //   initial_posterior(tng,1);
00868                 }
00869         }
00870         else {
00871                 if (p_origin) record()[0].missing -= 1;
00872         }
00873 }
00874 //RLF had  modified to work with pop->mean_for_genotype[i]
00875 
00876 double Individual::get_m_posterior(int excJ , Dblock& post_mat) {
00877         
00878         /////////////////////////////////////////////////////////
00879         //  if excJ < 0, say -1, means taking all posteriors
00880         /////////////////////////////////////////////////////////
00881         
00882         // make sure post_mat is initialized properly !!!!!!!!!!!!!!!!
00883         
00884         unsigned nswitches = n_switches();
00885         int spouse,i,j,k,ndim;
00886         double retval=0.0;
00887         ndim=post_mat.get_ndim();
00888         
00889         for (spouse=0; spouse<numspouse; spouse++) {
00890                 if (spouse != excJ) {
00891                         for (k=0; k < ndim; k++){
00892                                 for (i=0;i<nswitches; i++){
00893                                         for (j=0;j<4;j++){
00894                                                 post_mat[k][i][j] *= (m_posterior[spouse])[k][i][j];
00895                                         }
00896                                 }
00897                         }
00898                         retval += m_posterior_scale[spouse];
00899                 }
00900         }
00901         return retval;
00902 }
00903 
00904 
00905 void Individual::cal_gprobs(Dblock& post_mat_f) {
00906         double val;
00907         int ii,jj,kk;
00908         
00909         get_m_posterior(-1,post_mat_f);
00910         gprobs.resize(4);
00911         val = 0;
00912         for (jj=0;jj<4; jj++) {
00913                 for (kk=0; kk < post_mat_f.get_ndim() ; kk++) {
00914                         for (ii=0;ii< n_switches(); ii++) {
00915                                 gprobs[jj] += post_mat_f[kk][ii][jj]*m_anterior[kk][ii][jj];
00916                         }
00917                 }
00918                 val += gprobs[jj];
00919         }
00920         for (jj=0; jj<4; jj++) {
00921                 gprobs[jj] /= val;
00922         }
00923 }
00924 
00925 void Individual::collapse_antpost() {
00926         int i,j, k, nswitch,ndim;
00927         double sum=0.0;
00928         anterior.resize(4,0.0);
00929         initial_posterior(3,-1);
00930         nswitch=n_switches();
00931         ndim=m_anterior.get_ndim();
00932         // std::cout << nswitch << m_anterior << anterior << std::endl;
00933         for (k=0; k<ndim; k++) {
00934                 for (i=0;i<nswitch;i++){
00935                         //  std::cout << i << " " << nswitch << " " << (m_anterior[i][0]) << " " << anterior[i] << std::endl;
00936                         anterior[0] += m_anterior[k][i][0];
00937                         anterior[1] += m_anterior[k][i][1] + m_anterior[k][i][2];
00938                         anterior[2] += m_anterior[k][i][3];
00939                         for (j=0;j<numspouse;j++){
00940                                 (posterior[j])[0] += (m_posterior[j])[k][i][0];
00941                                 (posterior[j])[1] += (m_posterior[j])[k][i][1] + (m_posterior[j])[k][i][2];
00942                                 (posterior[j])[2] += (m_posterior[j])[k][i][3];
00943                         }
00944                 }
00945         }
00946 }
00947 
00948 void Individual::pretend_multi_m_missing(int on, int tng) {
00949         if (on) {
00950                 if (p_origin) {
00951                         record()[0].missing += 1;
00952                         initial_multi_m_anterior(tng);
00953                         //     initial_anterior(gfreq,tng);
00954                         if (connect_keeper == 1) {
00955                                 initial_multi_m_posterior(tng,1);
00956                                 // initial_posterior(tng,1);
00957                         }
00958                         else {
00959                                 initial_multi_m_posterior(tng,-1);
00960                                 // initial_posterior(tng);
00961                         }
00962                 }
00963                 else {
00964                         initial_multi_m_posterior(tng,1);
00965                         //   initial_posterior(tng,1);
00966                 }
00967         }
00968         else {
00969                 if (p_origin) record()[0].missing -= 1;
00970         }
00971 }
00972 
00973 void Individual::initial_multi_m_anterior(const unsigned tn_qtl) {
00974         int i_q, i_sw, i_switches,k, a1, a2;
00975         double freq_mark=1.0;
00976         double P_var = population->F->var; // polygenic varaince
00977         i_switches=n_switches();
00978         
00979         // compute frequency of marker information
00980         for (k=1;k<=population->n_markerLoci; k++){
00981                 a1 = genome0.chromosome[0].locus[k].allele;
00982                 a2 = genome1.chromosome[0].locus[k].allele;
00983                 freq_mark *= population->marker_freq(k)(a1)*population->marker_freq(k)(a2);
00984         }
00985         
00986         if (tn_qtl != mix_anterior.size())
00987                 mix_anterior.resize(tn_qtl);
00988         for (i_q=0; i_q < tn_qtl; i_q++){
00989                 mix_anterior[i_q].resize(i_switches);
00990                 for (i_sw=0; i_sw < i_switches; i_sw++) {
00991                         mix_anterior[i_q][i_sw].nu  = 0.0;
00992                         mix_anterior[i_q][i_sw].tsq = 1.0/P_var;
00993                         mix_anterior[i_q][i_sw].k   = std::log(1.0/std::sqrt(4.0*std::asin(1.0)*P_var)*population->Q_freq[i_q]*freq_mark);
00994                         //       std::cout << id() << " " << i_q << " " << i_sw << " nu " << mix_anterior[i_q][i_sw].nu << " tsq " << mix_anterior[i_q][i_sw].tsq << " k " << mix_anterior[i_q][i_sw].k << std::endl;
00995                 }
00996         }
00997         //  std::cout << "anterior initialized for " << id() << std::endl;
00998 }
00999 
01000 void Individual::initial_multi_m_posterior(const unsigned tn_qtl, const int cond) {
01001         
01002         int i_q, i_sw,spouse,i_switches, i_s;
01003         i_switches=n_switches();
01004         if (numspouse) { 
01005                 if (mix_posterior.size() != numspouse) { 
01006                         mix_posterior.resize(numspouse);
01007                         posterior_iw.resize(numspouse,0.0);
01008                         for (spouse=0; spouse<numspouse; spouse++) {
01009                             mix_posterior[spouse].done = 0;
01010                             mix_posterior[spouse].postvec.resize(tn_qtl);
01011                             for (i_q=0; i_q < tn_qtl; i_q++) {
01012                                         mix_posterior[spouse].postvec[i_q].resize(i_switches); 
01013                             } 
01014                         } 
01015                 } 
01016                 else { 
01017                         if (posterior){ 
01018                                 mix_posterior.resize(numspouse); 
01019                                 posterior_iw.resize(numspouse,0.0); 
01020                                 for (spouse=0; spouse<numspouse; spouse++) { 
01021                                         if (posterior_iw[spouse] != cond) { 
01022                                                 mix_posterior[spouse].done = 0; 
01023                                                 mix_posterior[spouse].postvec.resize(tn_qtl); 
01024                                                 for (i_q=0; i_q < tn_qtl; i_q++) { 
01025                                                         mix_posterior[spouse].postvec[i_q].resize(i_switches); 
01026                                                 } 
01027                                         } 
01028                                 } 
01029                         } 
01030                 } 
01031         } 
01032         else { 
01033                 mix_posterior.resize(numspouse); 
01034                 for (spouse=0; spouse < numspouse; spouse++){ 
01035                         mix_posterior[spouse].done = 0; 
01036                         mix_posterior[spouse].postvec.resize(tn_qtl); 
01037                         for (i_q=0; i_q < tn_qtl; i_q++) { 
01038                                 mix_posterior[spouse].postvec[i_q].resize(i_switches); 
01039                         } 
01040                 } 
01041         } 
01042 }
01043 
01044 
01045 
01046 
01047 void Individual::cal_m_gprobs(int tng) {
01048         double val,y,v_y,k_y,nu_y, a11, a12, a22, k, v, nu, inverse_sqrt2pi, sum1, sum2;
01049         int i_q, i_sw, i_sp, i_switches,j;
01050         i_switches=n_switches();
01051         inverse_sqrt2pi = 1.0/std::sqrt(4.0*std::asin(1.0));  // pi = 2.0*std::asin(1.0)
01052         gprobs.resize(tng);
01053         val = 0.0;
01054         sum1=0.0;
01055         for (i_q=0;i_q<tng;i_q++){
01056                 // contributions from penetrance function
01057                 // skip if phenotype if missing
01058                 
01059                 if (!(record()[0].missing)) {
01060                         y    = record()[0].double_val();
01061                         v_y    = 1.0/(*residual_var[0][0]);
01062                         k_y    = std::log(inverse_sqrt2pi*std::sqrt(v_y));
01063                         nu_y   = population->mean_for_genotype[i_q] - y;
01064                 }
01065                 else{
01066                         nu_y = 0.0;
01067                         v_y = 0.0;
01068                         k_y = 0.0;
01069                 }
01070                 
01071                 for (i_sw=0; i_sw < i_switches; i_sw++){   // loop for genotypes of J
01072                         
01073                         a11 = nu_y*nu_y*v_y;
01074                         a12 = nu_y*v_y;
01075                         a22 = v_y;
01076                         k   = k_y;
01077                         // contributions from anterior
01078                         
01079                         v   = mix_anterior[i_q][i_sw].tsq;
01080                         nu  = mix_anterior[i_q][i_sw].nu;
01081                         k   += mix_anterior[i_q][i_sw].k;
01082                         a11 += nu*nu*v;
01083                         a12 += -nu*v;
01084                         a22 += v;
01085                         
01086                         // Contributions from posteriors.
01087                         // Skip over any posteriors not yet available
01088                         unsigned nspouses =nspouse();
01089                         for (j=0; j<nspouses; j++){
01090                                 if (mix_posterior[j].done)  {
01091                                         v     = mix_posterior[j].postvec[i_q][i_sw].tsq;
01092                                         nu    = mix_posterior[j].postvec[i_q][i_sw].nu;
01093                                         k    += mix_posterior[j].postvec[i_q][i_sw].k;
01094                                         a11  += nu*nu*v;
01095                                         a12  += -nu*v;
01096                                         a22  += v;
01097                                 }
01098                         }
01099                         
01100                         // compute "log likelihood for gi"
01101                         
01102                         v          = 1.0/a22;
01103                         g_weight[i_q][i_sw]   = -0.5*(a11 - a12*a12*v)
01104                                 + k + std::log(std::sqrt(4.0*std::asin(1.0)*v));  // pi = 2.0*std::asin(1.0)
01105                         sum1 += g_weight[i_q][i_sw];
01106                 }
01107         }
01108         
01109         // compute sum of the "likelihoods"
01110         
01111         sum1 /= (tng*i_switches);
01112         for (i_q=0;i_q<tng;i_q++){
01113                 sum2=0.0;
01114                 for (i_sw=0; i_sw < i_switches; i_sw++){   // loop for genotypes of I
01115                         gprobs[i_q] += (std::exp(g_weight[i_q][i_sw] - sum1));
01116                 }
01117                 sum2 += gprobs[i_q];
01118         }
01119         for (i_q=0;i_q<tng;i_q++){
01120                 gprobs[i_q] /= sum2;
01121         }
01122 }
01123 
01124 void Individual::collapse_mix_antpost(){
01125         int tng=4;
01126         anterior.resize(tng,0.0);
01127         initial_posterior(tng,-1);
01128         double val,y,v_y,k_y,nu_y, a11, a12, a22, k, v, nu, inverse_sqrt2pi, sum1, sum2;
01129         int i_q, i_sw, i_sp, i_switches,j;
01130         i_switches=n_switches();
01131         inverse_sqrt2pi = 1.0/std::sqrt(4.0*std::asin(1.0));   // pi = 2.0*std::asin(1.0)
01132         gprobs.resize(3);
01133         val = 0.0;
01134         sum1=0.0;
01135         for (i_q=0;i_q<tng;i_q++){
01136                 // contributions from penetrance function, skip if phenotype is missing
01137                 if (!(record()[0].missing)) {
01138                         y    = record()[0].double_val();
01139                         v_y    = 1.0/(*residual_var[0][0]);
01140                         k_y    = std::log(inverse_sqrt2pi*std::sqrt(v_y));
01141                         nu_y   = population->mean_for_genotype[i_q] - y;
01142                 }
01143                 else{
01144                         nu_y = 0.0;
01145                         v_y = 0.0;
01146                         k_y = 0.0;
01147                 }
01148                 
01149                 for (i_sw=0; i_sw < i_switches; i_sw++){   // loop for genotypes of J
01150                         
01151                         a11 = nu_y*nu_y*v_y;
01152                         a12 = nu_y*v_y;
01153                         a22 = v_y;
01154                         k   = k_y;
01155                         // contributions from Father's anterior
01156                         
01157                         v   = mix_anterior[i_q][i_sw].tsq;
01158                         nu  = mix_anterior[i_q][i_sw].nu;
01159                         k   += mix_anterior[i_q][i_sw].k;
01160                         a11 += nu*nu*v;
01161                         a12 += -nu*v;
01162                         a22 += v;
01163                         
01164                         // compute "log likelihood for gi"
01165                         v    = 1.0/a22;
01166                         g_weight[i_q][i_sw]   = -0.5*(a11 - a12*a12*v)
01167                                 + k + std::log(std::sqrt(4.0*std::asin(1.0)*v));    // pi = 2.0*std::asin(1.0)
01168                         sum1 += g_weight[i_q][i_sw];
01169                 }
01170         }
01171         
01172         // compute sum of the "likelihoods"
01173         
01174         sum1 /= (tng*i_switches);
01175         for (i_q=0;i_q<tng;i_q++){
01176                 sum2=0.0;
01177                 for (i_sw=0; i_sw < i_switches; i_sw++){   // loop for genotypes of I
01178                         anterior[i_q] += (std::exp(g_weight[i_q][i_sw] - sum1));
01179                 }
01180                 sum2 += anterior[i_q];
01181         }
01182     anterior[0] /= sum2;
01183     anterior[1] /= sum2;
01184     anterior[2] /= sum2;
01185     anterior[3] /= sum2;
01186     anterior[1] += anterior[2];
01187     anterior[2]= anterior[3];
01188         
01189         // Now do the posteriors:
01190         unsigned nspouses =nspouse();
01191         for (j=0; j<nspouses; j++){
01192                 sum1=0.0;
01193                 for (i_q=0;i_q<tng;i_q++){
01194                         for (i_sw=0; i_sw < i_switches; i_sw++){   // loop for genotypes of J
01195                                                                                                            // Contributions from posteriors
01196                                 
01197                                 if (mix_posterior[j].done) {
01198                                         v     = mix_posterior[j].postvec[i_q][i_sw].tsq;
01199                                         nu    = mix_posterior[j].postvec[i_q][i_sw].nu;
01200                                         k    += mix_posterior[j].postvec[i_q][i_sw].k;
01201                                         a11  += nu*nu*v;
01202                                         a12  += -nu*v;
01203                                         a22  += v;
01204                                 }
01205                                 // compute "log likelihood for gi"
01206                                 v    = 1.0/a22;
01207                                 g_weight[i_q][i_sw] = -0.5*(a11 - a12*a12*v)
01208                                         + k + std::log(std::sqrt(4.0*std::asin(1.0)*v));    // pi = 2.0*std::asin(1.0)
01209                                 sum1 += g_weight[i_q][i_sw];
01210                         }
01211                 }
01212                 
01213                 // compute sum of the "likelihoods"
01214                 
01215                 sum1 /= (tng*i_switches);
01216                 for (i_q=0;i_q<tng;i_q++){
01217                         sum2=0.0;
01218                         for (i_sw=0; i_sw < i_switches; i_sw++){   // loop for genotypes of I
01219                                 posterior[j][i_q] += (std::exp(g_weight[i_q][i_sw] - sum1));
01220                         }
01221                         sum2 += posterior[j][i_q];
01222                 }
01223                 posterior[j][0] /= sum2;
01224                 posterior[j][1] /= sum2;
01225                 posterior[j][2] /= sum2;
01226                 posterior[j][3] /= sum2;
01227                 posterior[j][1]  +=posterior[j][2];
01228                 posterior[j][2]=posterior[j][3];
01229         }
01230 }
01231 
01232 
01233 double Individual::get_mix_penetrance(Vector<double> &pen) {
01234         
01235         unsigned QTL,PGN;
01236         
01237         if (myrecord[0].missing) {
01238                 for (QTL=0; QTL<3; QTL++) {
01239                         pen[QTL] = 1.0;
01240                 }
01241                 return 0.0;
01242         }
01243         
01244         double inverse_sqrt2pi = 1.0/std::sqrt(4.0*std::asin(1.0));       // pi = 2.0*std::asin(1.0)
01245         double mu, v = *residual_var[0][0],scale;
01246         double y = record()[0].double_val();
01247         
01248     mu = population->mean_for_genotype[0];
01249     pen[0] = inverse_sqrt2pi/std::sqrt(v) * std::exp(-(y-mu)*(y-mu)/(2*v));
01250     mu = population->mean_for_genotype[1];
01251     pen[1] = inverse_sqrt2pi/std::sqrt(v) * std::exp(-(y-mu)*(y-mu)/(2*v));
01252     mu = population->mean_for_genotype[2];
01253     pen[1] += inverse_sqrt2pi/std::sqrt(v) * std::exp(-(y-mu)*(y-mu)/(2*v));
01254     mu = population->mean_for_genotype[3];
01255     pen[2] = inverse_sqrt2pi/std::sqrt(v) * std::exp(-(y-mu)*(y-mu)/(2*v));
01256         
01257         scale = pen[0]+ pen[1] + pen[2];
01258         
01259         if (fabs(scale) < 1.0e-8) {
01260                 warning("this record is not conformable");
01261                 scale = 3.0;
01262                 for (QTL=0; QTL<3; QTL++)
01263                         pen[QTL] = 1.0;
01264         }
01265         else {
01266                 for (QTL=0; QTL<3; QTL++) {
01267                         pen[QTL] /= scale;
01268                 }
01269                 //    std::cout << "pen in get pen " << pen << std::endl;
01270         }
01271         return std::log(scale);
01272 }
01273 
01274 //BRS
01275 
01276 //MJS
01277 // moved to matvec 1.0.2d RLF
01278 //descent graph stuff
01279 void Individual::put_gametes(Vector <unsigned> m_g, Vector <unsigned> p_g) {
01280         // Authors: Matthias Schelling and Rohan L. Fernando 
01281         // (1999) 
01282         // Contributors:
01283         m_gamete = m_g;
01284         p_gamete = p_g;
01285         if(m_gamete.size()!= p_gamete.size()){
01286                 throw exception("Individual::put_gametes: m_gamete.size()!= p_gamete.size()");
01287         }
01288         //init q
01289         q_prob = 1.0;
01290         calc_q();
01291 }
01292 
01293 double Individual::calc_q(void) {
01294         // Authors: Matthias Schelling and Rohan L. Fernando 
01295         // (1999) 
01296         // Contributors: L. Radu Totir 
01297         if ((m_gamete.size() == 0) || (p_gamete.size() == 0)) {
01298                 throw exception ("Individual::calc_q: no graph present.");
01299         }
01300         double old_q_prob = q_prob; //store old q
01301         if (mymother) { // if not founder
01302                 q_prob = 0.25;
01303                 for (int j=0; j<numLoci; j++) {
01304                         if (m_gamete[j] == m_gamete[j+1]) {
01305                                 q_prob *= 1 - population->RecoVector[j];
01306                         }
01307                         else {
01308                                 q_prob *= population->RecoVector[j];
01309                         }
01310                         if (p_gamete[j] == p_gamete[j+1]) {
01311                                 q_prob *= 1 - population->RecoVector[j];
01312                         }
01313                         else {
01314                                 q_prob *= population->RecoVector[j];
01315                         }
01316                 }
01317         }
01318         else {
01319                 q_prob = 1.0;
01320         }
01321         return std::log(old_q_prob) - std::log(q_prob);
01322 }
01323 /*! \fn double Individual::calc_q(void)
01324 *  \brief method to calc q prob for individual; q is stored in 1 variable
01325 */
01326 
01327 void Individual::sample_haplotypes(const unsigned origin) {
01328         //sample haplotypes for an SLlocus,
01329         //origin 0 = maternal, 1 = paternal
01330         Vector <unsigned>* gamete = &m_gamete;
01331         if (origin) {
01332                 gamete = &p_gamete;
01333         }
01334         //first, go from SLlocus down to locus 1
01335         for (unsigned i = population->SLlocus - 1; i > 0; i--) {
01336                 if (ranf() < population->RecoVector(i)) { //do crossover
01337                         if ((*gamete)(i+1) == 0) { (*gamete)(i) = 1; }
01338                         else { (*gamete)(i) = 0; }
01339                 }
01340                 else { //no crossover
01341                         if ((*gamete)(i+1) == 0) { (*gamete)(i) = 0; }
01342                         else { (*gamete)(i) = 1; }
01343                 }
01344         }
01345         //then, go from SLlocus up to last locus
01346         for (unsigned i = population->SLlocus; i < (*gamete).size(); i++) {
01347                 if (ranf() < population->RecoVector(i)) { //do crossover
01348                         if ((*gamete)(i) == 0) { (*gamete)(i+1) = 1; }
01349                         else { (*gamete)(i+1) = 0; }
01350                 }
01351                 else { //no crossover
01352                         if ((*gamete)(i) == 0) { (*gamete)(i+1) = 0; }
01353                         else { (*gamete)(i+1) = 1; }
01354                 }
01355         }
01356 }
01357 
01358 void Individual::set_founder_alleles(bool update) {
01359         // Authors: Matthias Schelling and Rohan L. Fernando 
01360         // (1999) 
01361         // Contributors: L. Radu Totir 
01362         
01363         //  update == 0  ->  first pass
01364         //  update == 1  ->  already done
01365         //                   Done
01366         //                T        F
01367         //             ________________
01368         //           T | yes      yes |
01369         // Update    |              |
01370         //           F | no       yes |
01371         //           |______________|
01372         //                            
01373         // This means that we set the founder_alleles from
01374         // the founder_allele_counter only if not already done
01375         // When we are updating, we change only offsprings
01376         
01377         // get space if needed
01378         if (m_founder.size() == 0) { 
01379                 int zero=0;
01380                 m_founder.resize(numLoci, zero);
01381                 p_founder.resize(numLoci, zero);
01382         }
01383         
01384         // if not already done
01385         //if (!m_founder(1)) {
01386         if (!(!update && m_founder(1))) {
01387                 if (!m_founder(1) && update) {
01388                         std::cout << "should not happen: update = 1, m_founder(1) set" << std::endl;
01389                 }
01390                 // if founder set m/p_founder to founder_allele_counter
01391                 if (mymother == NULL) {
01392                         if (!update) {
01393                                 std::cout << "founder_allele_counter++" << std::endl;
01394                                 population->founder_allele_counter++;
01395                                 for (int i=1; i<= numLoci; i++){
01396                                         m_founder(i) = population->founder_allele_counter;
01397                                         p_founder(i) = population->founder_allele_counter+1;
01398                                 }
01399                                 population->founder_allele_counter++;
01400                         }
01401                 }
01402                 // if not founder 
01403                 else {
01404                         Vector <unsigned> &mm_founder = mymother->get_founder("mother",update);
01405                         Vector <unsigned> &mp_founder = mymother->get_founder("father",update);
01406                         Vector <unsigned> &pm_founder = myfather->get_founder("mother",update);
01407                         Vector <unsigned> &pp_founder = myfather->get_founder("father",update);  
01408                         for (int i=1; i <= numLoci; i++){
01409                                 if (m_gamete(i) == 0) {
01410                                         m_founder(i) = mm_founder(i);
01411                                 }
01412                                 else {
01413                                         m_founder(i) = mp_founder(i);
01414                                 }
01415                                 if (p_gamete(i) == 0) {
01416                                         p_founder(i) = pm_founder(i);
01417                                 }
01418                                 else {
01419                                         p_founder(i) = pp_founder(i);
01420                                 }
01421                         }
01422                 } 
01423         }
01424         }
01425 
01426 Vector <unsigned> &Individual::get_founder(std::string parent, bool update) {
01427         set_founder_alleles(update);
01428         
01429         if (parent == "mother") {
01430                 return m_founder;
01431         }
01432         else {
01433                 return p_founder;
01434         }
01435 }
01436 
01437 void Individual::update_offsprings_founder_alleles(void) {
01438         set_founder_alleles(1); //1 invokes updating
01439                                                         //do all ofsprings
01440         if (numoffs) {
01441                 for (int i=0; i<numoffs; i++) {
01442                         myoffspring[i]->update_offsprings_founder_alleles();
01443                 }
01444         }
01445 }
01446 
01447 double Individual::graph_trans_prob(unsigned uj,
01448                                                                         unsigned um,
01449                                                                         unsigned uf){
01450     double m_prob=0, f_prob=0;
01451     unsigned FQTable[4][4]={2,3,2,3,0,1,0,1,1,0,1,0,3,2,3,2};
01452     unsigned MQTable[4][4]={2,2,3,3,0,0,1,1,1,1,0,0,3,3,2,2};
01453     unsigned FQindex = FQTable[uf][uj];
01454     unsigned MQindex = MQTable[um][uj];
01455     unsigned MMindex = 0, FMindex = 0, size = m_gamete.size();
01456     
01457     //calc marker indexes (bin -> dec)
01458     for(int k = size; k > 0; k--){
01459                 /*std::cout << "m_marker" << k << '\t' << m_gamete(k)
01460                 << "\tp_marker" << k << '\t' << p_gamete(k)
01461                 << std::endl;
01462                 */
01463                 MMindex+=(unsigned)((1 << (size-k)))*m_gamete(k);
01464                 FMindex+=(unsigned)((1 << (size-k)))*p_gamete(k);
01465     }
01466         
01467     if(FQindex==3 || MQindex==3){return 0.0;}
01468     if(FQindex==2){f_prob=1;}
01469     else{
01470                 f_prob=population->model->gamete_prob_table[FMindex][FQindex];
01471                 //std::cout << "FMindex: " << FMindex << "\tFQindex: " << FQindex << std::endl;
01472     }
01473     if(MQindex==2){m_prob=1;}
01474     else{
01475                 m_prob=population->model->gamete_prob_table[MMindex][MQindex];
01476                 //std::cout << "MMindex: " << MMindex << "\tMQindex: " << MQindex << std::endl;
01477     }
01478     //return transition (transmission*transmission) probability
01479     //std::cerr << "return transition " << f_prob*m_prob << std::endl;
01480     return f_prob*m_prob;
01481 }
01482 
01483 
01484 void Individual::t0(unsigned SLlocus, unsigned gender_parent) {
01485         //Sobel & Lange T0 rule
01486         if (gender_parent) {
01487                 if (p_gamete(SLlocus)){p_gamete(SLlocus)=0;}
01488                 else {p_gamete(SLlocus)=1;}
01489         }
01490         else {
01491                 if (m_gamete(SLlocus)){m_gamete(SLlocus)=0;}
01492                 else {m_gamete(SLlocus)=1;}
01493         }
01494         update_offsprings_founder_alleles();
01495         return;
01496 }
01497 
01498 void Individual::t1(unsigned SLlocus) {
01499         //Sobel & Lange T1 rule
01500         //get offsprings
01501         for (int i=0; i<numoffs; i++) {
01502                 if (myoffspring[i]->mymother == this) {
01503                         myoffspring[i]->t0(SLlocus, 0);
01504                 }
01505                 else {
01506                         myoffspring[i]->t0(SLlocus, 1);
01507                 }
01508         }
01509         update_offsprings_founder_alleles();
01510         return;
01511 }
01512 
01513 void Individual::t2a(unsigned SLlocus) {
01514         //Sobel & Lange T2a rule
01515         //choose spouse
01516         Individual *my_spouse = spouselist[unsigned(ranf() * numspouse)];
01517         //get offsprings
01518         for (int i=0; i<numoffs; i++) {
01519                 //check if related with myspouse
01520                 if (myoffspring[i]->myfather == my_spouse
01521                         || myoffspring[i]->mymother == my_spouse) {
01522                         
01523                         //do t0 for myoffspring for both alleles if they originate in the opposite gender
01524                         if (myoffspring[i]->m_gamete(SLlocus) != myoffspring[i]->p_gamete(SLlocus)) {
01525                                 myoffspring[i]->t0(SLlocus,0);
01526                                 myoffspring[i]->t0(SLlocus,1);
01527                         }
01528                         
01529                         //do t1 for myoffspring
01530                         myoffspring[i]->t1(SLlocus);
01531                 }
01532         }
01533         return;
01534 }
01535 
01536 void Individual::t2b(unsigned SLlocus) {
01537         //Sobel & Lange T2b rule
01538         //choose spouse
01539         Individual *my_spouse = spouselist[unsigned(ranf() * numspouse)];
01540         //get offsprings
01541         for (int i=0; i<numoffs; i++) {
01542                 //check if related with myspouse
01543                 if (myoffspring[i]->myfather == my_spouse
01544                         || myoffspring[i]->mymother == my_spouse) {
01545                         
01546                         //do t0 for myoffspring for both alleles if they originate in the same gender
01547                         if (myoffspring[i]->m_gamete(SLlocus) == myoffspring[i]->p_gamete(SLlocus)) {
01548                                 myoffspring[i]->t0(SLlocus,0);
01549                                 myoffspring[i]->t0(SLlocus,1);
01550                         }
01551                         
01552                         //do t1 for myoffspring
01553                         myoffspring[i]->t1(SLlocus);
01554                 }
01555         }
01556         return;
01557 }   
01558 
01559 void Individual::apply_SL_transition(const unsigned rule, const unsigned locus) {
01560         if (rule > 3) {
01561                 throw exception("Population::apply_SL_transition: nonexistent rule!, abort");
01562         }
01563         switch (rule) {
01564                 case 0: {       //T0 rule
01565                         t0(locus,unsigned(ranf()*2)); //choose gender at random
01566                         break;
01567                 }
01568                 case 1: {       //T1 rule
01569                         t1(locus);
01570                         break;
01571                 }
01572                 case 2: {       //T2a rule
01573                         t2a(locus);
01574                         break;
01575                 }
01576                 case 3: {       //T2b rule
01577                         t2b(locus);
01578                         break;
01579                 }
01580         }
01581         return;
01582 }
01583 
01584 void Individual::apply_SL_cascade(const unsigned rule, const unsigned SLlocus) {
01585         //going from SL locus, we apply SL to the other loci depending on distance
01586         //first, go from SLlocus down to locus 1
01587         unsigned rule_applied = 1;
01588         for (unsigned locus = SLlocus - 1; locus > 0; locus--) {
01589                 double u = ranf();
01590                 if ((u > population->RecoVector(locus) && rule_applied) ||
01591                         (u < population->RecoVector(locus) && !rule_applied)) { 
01592                         //apply T rule
01593                         apply_SL_transition(rule, locus);
01594                         rule_applied = 1;
01595                 }
01596                 else {
01597                         rule_applied = 0;
01598                 }
01599         }
01600         
01601         //then go from SLlocus up to last locus
01602         rule_applied = 1;
01603         for (unsigned locus = SLlocus; locus < m_gamete.size(); locus++) {
01604                 double u = ranf();
01605                 if ((u > population->RecoVector(locus) && rule_applied) ||
01606                         (u < population->RecoVector(locus) && !rule_applied)) { 
01607                         //apply T rule
01608                         apply_SL_transition(rule, locus + 1);
01609                         rule_applied = 1;
01610                 }
01611                 else {
01612                         rule_applied = 0;
01613                 }
01614         }
01615         return;
01616 }
01617 
01618 //E. Thompson meiosis sampler
01619 
01620 void Individual::sample_self(int parent_indicator){
01621         // Authors: Matthias Schelling and Rohan L. Fernando 
01622         // (1999) 
01623         // Contributors: L. Radu Totir 
01624         
01625         // parent_indicator = 0 -> maternal
01626         // parent_indicator = 1 -> paternal
01627         Vector <unsigned> *gamete_p;
01628         
01629         if(parent_indicator==0){
01630                 gamete_p=&m_gamete;
01631         }
01632         else {
01633                 gamete_p=&p_gamete;
01634         }
01635         
01636         // Computing the priors
01637         
01638         for ( int j=1; j<=numLoci; j++){
01639                 (*gamete_p)(j)=0;
01640                 update_offsprings_founder_alleles();
01641                 vec_prob[j-1][0]=population->calc_prior_descent_graph(j);
01642                 (*gamete_p)(j)=1;
01643                 update_offsprings_founder_alleles();
01644                 vec_prob[j-1][1]=population->calc_prior_descent_graph(j); 
01645         }
01646         // Computing the "cutset" values
01647         int j=0;
01648         vec_cutsetval[j][0]=vec_prob[j][0];
01649         vec_cutsetval[j][1]=vec_prob[j][1];
01650         for ( int j=1; j<numLoci; j++){ 
01651                 vec_cutsetval[j][0]= vec_prob[j][0]*
01652                 (vec_cutsetval[j-1][0]*(1-population-> RecoVector(j))
01653                  + vec_cutsetval[j-1][1]*(population-> RecoVector(j)));
01654                 vec_cutsetval[j][1]= vec_prob[j][1]*
01655                         (vec_cutsetval[j-1][1]*(1-population-> RecoVector(j))
01656                          + vec_cutsetval[j-1][0]*(population-> RecoVector(j)));
01657         }
01658         
01659         // Sample the last locus from the marginal
01660         
01661         j=numLoci-1;
01662         double u = ranf();
01663         double p1 = vec_cutsetval[j][0];
01664         double p2 = vec_cutsetval[j][1];
01665         if (u<(p1/(p1+p2))){
01666                 (*gamete_p)(j+1)=0;
01667         }
01668         else{
01669                 (*gamete_p)(j+1)=1;
01670         }
01671         update_offsprings_founder_alleles();
01672         
01673         // sample the other loci
01674         
01675         for ( int j=numLoci-2; j>=0; j--){ 
01676                 double u = ranf();
01677                 double recomb=0.0;
01678                 if((*gamete_p)(j+2)==0){
01679                         recomb=(1-population-> RecoVector(j+1));
01680                 }
01681                 else{
01682                         recomb=population-> RecoVector(j+1);
01683                 }
01684                 double p1= recomb*vec_cutsetval[j][0];
01685                 double p2=(1-recomb)*vec_cutsetval[j][1];
01686                 if (u<(p1/(p1+p2))){
01687                         (*gamete_p)(j+1)=0;
01688                 }
01689                 else{
01690                         (*gamete_p)(j+1)=1;
01691                 }
01692                 update_offsprings_founder_alleles();
01693         }
01694 }
01695 
01696 //MJS
01697 
01698 //FVP 02/20/03
01699 Vector<unsigned>  Individual::get_id_pdq(void){
01700         // Authors: Fabiano V. Pita and Rohan L. Fernando 
01701         // (2003) 
01702         // Contributors:
01703         return id_pdq;
01704 }
01705 
01706 void Individual::search_heteroz(unsigned qtl){
01707         // Authors: Fabiano V. Pita and Rohan L. Fernando 
01708         // (2003) 
01709         // Contributors: L. Radu Totir 
01710         
01711         int mult = -1;           //  __1___2___3___qtl___5____6__
01712         unsigned step = 1;
01713         unsigned endl = 0;
01714         unsigned endr = 0;
01715         unsigned lcs = qtl;
01716         
01717         ord_heter = 0;  // 0 means individual is not an ordered heterozygous
01718         
01719         if (mymother != NULL){
01720                 return;
01721         }
01722         
01723         // first goes to the closest left marker
01724         lcs += step * mult;
01725         mult *= -1;
01726         
01727         for (;;) {
01728                 if (genome0.chromosome[0].locus[lcs-1].allele != genome1.chromosome[0].locus[lcs-1].allele){
01729                         ord_heter = lcs;
01730                         break;
01731                 }
01732                 
01733                 if (lcs == 1){  //last locus of left side
01734                         endl = 1;
01735                 }
01736                 
01737                 if (lcs == numLoci){  //last locus of right side
01738                         endr = 1;
01739                 }
01740                 
01741                 step += 1;  
01742                 lcs += step*mult;
01743                 mult *= -1;
01744                 
01745                 if (lcs > numLoci){
01746                         break;
01747                 }
01748                 
01749                 //working just with the right side
01750                 if (endl == 1){
01751                         for (int i=lcs; i<= numLoci; i++){
01752                                 if (genome0.chromosome[0].locus[lcs-1].allele != genome1.chromosome[0].locus[lcs-1].allele){
01753                                         ord_heter = lcs;
01754                                         break;
01755                                 }
01756                         }
01757                         break;
01758                 }
01759                 
01760                 //working just with the left side
01761                 if (endr == 1){
01762                         for (int i=lcs; i>=1; i--){
01763                                 if (genome0.chromosome[0].locus[lcs-1].allele != genome1.chromosome[0].locus[lcs-1].allele){
01764                                         ord_heter = lcs;
01765                                         break;
01766                                 }
01767                         }
01768                         break;
01769                 }
01770         }
01771         
01772 }
01773 
01774 void Individual::get_allele_v1(unsigned lcs, Vector<int> allele_vector1){
01775         if (ord_heter == lcs){
01776                 initial_order1 = allele_vector1(m_founder(lcs));
01777                 initial_order2 = allele_vector1(p_founder(lcs));
01778                 std::cout<< id() <<" "<<lcs<<" "<<initial_order1<<" "<<initial_order2<< std::endl;
01779         }
01780 }
01781 
01782 // FVP 06/02/03 ->accumulates the haplotype type for each individual 
01783 void Individual::count_haplotype(void){
01784         // Authors: Fabiano V. Pita and Rohan L. Fernando 
01785         // (2003) 
01786         // Contributors: L. Radu Totir  
01787         if (m_haplotype.empty()){
01788                 int zero = 0;
01789                 m_haplotype.resize(4,numLoci-1,zero);
01790                 p_haplotype.resize(4,numLoci-1,zero);
01791         }
01792         
01793         if (mymother == NULL) { // if founder
01794                 for (int j=1; j<=numLoci-1;j++){
01795                         m_haplotype(1,j) = 0;
01796                         m_haplotype(2,j) = 0;
01797                         m_haplotype(3,j) = 0;
01798                         m_haplotype(4,j) = 0;
01799                         p_haplotype(1,j) = 0;
01800                         p_haplotype(2,j) = 0;
01801                         p_haplotype(3,j) = 0;
01802                         p_haplotype(4,j) = 0;
01803                 }
01804         }
01805         else {
01806                 for (int j=2; j<=numLoci; j++) {
01807                         if ((m_gamete(j-1) == 0) && (m_gamete(j) == 0)) {
01808                                 m_haplotype(1,j-1)+=1;
01809                         }      
01810                         if ((m_gamete(j-1) == 0) && (m_gamete(j) == 1)) {
01811                                 m_haplotype(2,j-1)+=1;
01812                         }
01813                         if ((m_gamete(j-1) == 1) && (m_gamete(j) == 0)) {
01814                                 m_haplotype(3,j-1)+=1;
01815                         }
01816                         if ((m_gamete(j-1) == 1) && (m_gamete(j) == 1)) {
01817                                 m_haplotype(4,j-1)+=1;
01818                         }
01819                         if ((p_gamete(j-1) == 0) && (p_gamete(j) == 0)) {  
01820                                 p_haplotype(1,j-1)+=1;
01821                         }      
01822                         if ((p_gamete(j-1) == 0) && (p_gamete(j) == 1)) {
01823                                 p_haplotype(2,j-1)+=1;
01824                         }
01825                         if ((p_gamete(j-1) == 1) && (p_gamete(j) == 0)) {
01826                                 p_haplotype(3,j-1)+=1;
01827                         }
01828                         if ((p_gamete(j-1) == 1) && (p_gamete(j) == 1)) {
01829                                 p_haplotype(4,j-1)+=1;
01830                         }
01831                 }
01832         }
01833 }  
01834 
01835 // HG : to print the haplotype frequencies used in additive.cpp
01836 void Individual::display_freq_haplotype(unsigned n)
01837 { 
01838         // Authors: Helene Gilbert and Rohan L. Fernando 
01839         // (2004) 
01840         // Contributors: L. Radu Totir
01841         double doublen = n;
01842         cout << myid << std::setw(8);
01843         for (int j=1; j<numLoci; j++) {
01844                 double m_haplo1 = m_haplotype(1,j);
01845                 double m_haplo2 = m_haplotype(2,j);
01846                 double m_haplo3 = m_haplotype(3,j);
01847                 double m_haplo4 = m_haplotype(4,j);
01848                 double p_haplo1 = p_haplotype(1,j);
01849                 double p_haplo2 = p_haplotype(2,j);
01850                 double p_haplo3 = p_haplotype(3,j);
01851                 double p_haplo4 = p_haplotype(4,j);
01852                 cout <<m_haplo1/doublen<<" "<<m_haplo2/doublen<<" "<<m_haplo3/doublen<<" "<<m_haplo4/doublen<<" "
01853                         <<p_haplo1/doublen<<" "<<p_haplo2/doublen<<" "<<p_haplo3/doublen<<" "<<p_haplo4/doublen<<"       ";
01854         }
01855         cout<<endl;
01856 }
01857 
01858 // Fabiano
01859 // method for calculating the pdq from individual's haplotype probabilities (2 markers interval)
01860 // mm_pdq = considering the individual maternal haplotype, pdq from the mother's maternal haplotype
01861 // mp_pdq = considering the individual maternal haplotype, pdq from the mother's paternal haplotype 
01862 
01863 void Individual::map_pdq(int samples, int interval, BG r_aq, BG r_bq){
01864         double doubleSample = samples;
01865         BG r_ab = (1-r_aq)*r_bq + r_aq*(1-r_bq);
01866         if (mymother == NULL) { // if founder
01867                 bgDMPDQ = 0.0;
01868                 bgSMPDQ = 0.0;    
01869         }
01870         else {  
01871                 //     bgDMPDQ = (1-r_aq)*(1-r_bq)*m_haplotype(1,interval)/samples + (1-r_aq)*(r_bq)*m_haplotype(2,interval)/samples
01872                 //             + (r_aq)*(1-r_bq)*m_haplotype(3,interval)/samples + (r_aq)*(r_bq)*m_haplotype(4,interval)/samples; 
01873                 
01874                 //     bgSMPDQ = (1-r_aq)*(1-r_bq)*p_haplotype(1,interval)/samples + (1-r_aq)*(r_bq)*p_haplotype(2,interval)/samples
01875                 //             + (r_aq)*(1-r_bq)*p_haplotype(3,interval)/samples + (r_aq)*(r_bq)*p_haplotype(4,interval)/samples; 
01876                 //                        !!!
01877                 //     bgDMPDQ = m_haplotype(1,interval)/doubleSample + r_bq/(r_aq+r_bq)*m_haplotype(2,interval)/doubleSample
01878                 //       + (r_aq/(r_aq+r_bq))*m_haplotype(3,interval)/doubleSample ;
01879                 
01880                 //     bgSMPDQ = p_haplotype(1,interval)/doubleSample + r_bq/(r_aq+r_bq)*p_haplotype(2,interval)/doubleSample
01881                 //             + ((r_aq)/(r_aq+r_bq))*p_haplotype(3,interval)/doubleSample; 
01882                 bgDMPDQ = (1-r_aq)*(1-r_bq) / (1-r_ab)     * m_haplotype(1,interval)/doubleSample 
01883                 + (1-r_aq)*(r_bq)   / r_ab         * m_haplotype(2,interval)/doubleSample
01884                 + (r_aq)  *(1-r_bq) / r_ab         * m_haplotype(3,interval)/doubleSample 
01885                 + (r_aq)*(r_bq)     / (1-r_ab)     * m_haplotype(4,interval)/doubleSample; 
01886                 
01887                 bgSMPDQ = (1-r_aq)*(1-r_bq) / (1-r_ab)     * p_haplotype(1,interval)/doubleSample 
01888             + (1-r_aq)*(r_bq)   / r_ab         * p_haplotype(2,interval)/doubleSample
01889             + (r_aq)  *(1-r_bq) / r_ab         * p_haplotype(3,interval)/doubleSample 
01890             + (r_aq)*(r_bq)     / (1-r_ab)     * p_haplotype(4,interval)/doubleSample; 
01891                 //                        !!!
01892                 //        cout<<interval<<" Maternal haplotype pty "<< m_haplotype(1,interval)<<" "<<m_haplotype(2,interval)<<" "<<m_haplotype(3,interval)<<" "<<bgDMPDQ.f<<endl;
01893                 //        cout<<interval<<" Paternal haplotype pty "<< p_haplotype(1,interval)<<" "<<p_haplotype(2,interval)<<" "<<p_haplotype(3,interval)<<" "<<bgSMPDQ.f<<endl;
01894         }
01895 }
01896 
01897 void Individual::setFounderHaplotype(void)
01898 {
01899         // Authors: L. Radu Totir 
01900         // (August, 2004) 
01901         // Contributors:
01902         if(mymother){
01903                 throw exception("Individual::setFounderHaplotype(void) called for a non-founder");
01904         }
01905         else {
01906                 for (unsigned j=0;j<numLoci;j++){
01907                         unsigned patAllele = genome0.chromosome[0].locus[j].allele;
01908                         unsigned matAllele = genome1.chromosome[0].locus[j].allele;
01909                         if (matAllele > patAllele){
01910                                 genome0.chromosome[0].locus[j].allele=matAllele;
01911                                 genome1.chromosome[0].locus[j].allele=patAllele;
01912                         }
01913                 }
01914         }
01915 }
01916 /*! \fn double Individual::setFounderHaplotype(void)
01917 *  \brief makes sure that in founders the maternal allele 
01918 *  is asigned the lower allele state 
01919 */
01920 
01921 double Individual::getAllelePenetrance(void)
01922 {
01923         // Authors: L. Radu Totir and Rohan L. Fernando 
01924         // (September, 2003) 
01925         // Contributors: 
01926         unsigned mState, pState;
01927         if(population->prior->chrom()[0].locus[currentLocus].qtl_ml=='q'){
01928                 double mu=0.0,d,pen;
01929                 double v = *residual_var[0][0];
01930                 unsigned genoState, QLocus;
01931                 if (myrecord[0].missing) {
01932                         return 1.0;
01933                 }
01934                 for (unsigned i=0;i<nQTL;i++) {
01935                         QLocus = QTLPosVector[i];
01936                         mState    = malleleStateNodeVector[QLocus].getState();
01937                         pState    = palleleStateNodeVector[QLocus].getState();
01938                         mu+=population->prior->chrom()[0].locus[QLocus].genotypic_val_mat(pState+1,mState+1);
01939                 }
01940                 double y = record()[0].double_val();
01941                 pen = isqrt2pi/std::sqrt(v)*std::exp(-(y-mu)*(y-mu));
01942                 return pen;
01943         }
01944         else{
01945                 unsigned allelePat = genome0.chromosome[0].locus[currentLocus].allele;
01946                 unsigned alleleMat = genome1.chromosome[0].locus[currentLocus].allele;
01947                 if(allelePat!=0 && allelePat!=alleleMat){
01948                         mState    = malleleStateNodeVector[currentLocus].getMyAlleleState();
01949                         pState    = palleleStateNodeVector[currentLocus].getMyAlleleState();
01950                         if(mState != pState){
01951                                 return 1.0;
01952                         }
01953                         else {
01954                                 return 0.0;
01955                         }
01956                 }
01957                 else {
01958                         return 1.0;
01959                 }
01960         }
01961 }
01962 /*! \fn double Individual::getAllelePenetrance(void)
01963 *  \brief returns the allele penetrance for the RSampler
01964 */
01965 
01966 double Individual::getAllelePenetrance(unsigned forLocus)
01967 {
01968         // Authors: L. Radu Totir and Rohan L. Fernando 
01969         // (September, 2004) 
01970         // Contributors: 
01971         unsigned mState, pState;
01972         if(population->prior->chrom()[0].locus[forLocus].qtl_ml=='q'){
01973                 double mu=0.0,d,pen;
01974                 double v = *residual_var[0][0];
01975                 unsigned genoState, QLocus;
01976                 if (myrecord[0].missing) {
01977                         return 1.0;
01978                 }
01979                 for (unsigned i=0;i<nQTL;i++) {
01980                         QLocus = QTLPosVector[i];
01981                         mState    = malleleStateNodeVector[QLocus].getState();
01982                         pState    = palleleStateNodeVector[QLocus].getState();
01983                         mu+=population->prior->chrom()[0].locus[QLocus].genotypic_val_mat(pState+1,mState+1);
01984                 }
01985                 double y = record()[0].double_val();
01986                 pen = isqrt2pi/std::sqrt(v)*std::exp(-(y-mu)*(y-mu));
01987                 return pen;
01988         }
01989         else{
01990                 unsigned allelePat = genome0.chromosome[0].locus[forLocus].allele;
01991                 unsigned alleleMat = genome1.chromosome[0].locus[forLocus].allele;
01992                 if(allelePat!=0 && allelePat!=alleleMat){
01993                         mState    = malleleStateNodeVector[forLocus].getMyAlleleState();
01994                         pState    = palleleStateNodeVector[forLocus].getMyAlleleState();
01995                         if(mState != pState){
01996                                 return 1.0;
01997                         }
01998                         else {
01999                                 return 0.0;
02000                         }
02001                 }
02002                 else {
02003                         return 1.0;
02004                 }
02005         }
02006 }
02007 /*! \fn double Individual::getAllelePenetrance(void)
02008 *  \brief returns the allele penetrance for the RSampler
02009 */
02010 
02011 double Individual::getDisAllelePenetrance(void)
02012 {
02013         // Authors: L. Radu Totir and Rohan L. Fernando 
02014         // (September, 2004) 
02015         // Contributors: 
02016         double mu=0.0,d,pen;
02017         unsigned mState, pState;
02018         if (myrecord[0].missing) {
02019                 return 1.0;
02020         }
02021         mState    = malleleStateNodeVector[currentLocus].getState();
02022         pState    = palleleStateNodeVector[currentLocus].getState();
02023         int phenotype = int(record()[0].double_val());
02024         // cout << phenotype << "   " << mState+pState << "  " 
02025         // << population->disPenetranceTable[phenotype][mState+pState] 
02026         // << endl;
02027         return population->disPenetranceTable[phenotype][mState+pState];
02028 }
02029 
02030 double Individual::getDisAllelePenetrance(unsigned forLocus)
02031 {
02032         // Authors: L. Radu Totir and Rohan L. Fernando 
02033         // (September, 2004) 
02034         // Contributors: 
02035         double mu=0.0,d,pen;
02036         unsigned mState, pState;
02037         if (myrecord[0].missing) {
02038                 return 1.0;
02039         }
02040         mState    = malleleStateNodeVector[forLocus].getState();
02041         pState    = palleleStateNodeVector[forLocus].getState();
02042         int phenotype = int(record()[0].double_val());
02043         // cout << phenotype << "   " << mState+pState << "  " 
02044         // << population->disPenetranceTable[phenotype][mState+pState] 
02045         // << endl;
02046         return population->disPenetranceTable[phenotype][mState+pState];
02047 }
02048 
02049 double Individual::getGenoPenetrance(void)
02050 {
02051         // Authors: L. Radu Totir and Rohan L. Fernando 
02052         // (June, 2003) 
02053         // Contributors: 
02054         double mu=0.0,d,pen;
02055         double v = *residual_var[0][0];
02056         unsigned mState, pState, QLocus;
02057         if (myrecord[0].missing) {
02058                 return 1.0;
02059         }
02060         for (unsigned i=0;i<nQTL;i++) {
02061                 QLocus = QTLPosVector[i];
02062                 mState    = genotNodeVector[QLocus].getmState();
02063                 pState    = genotNodeVector[QLocus].getpState();
02064                 mu+=population->prior->chrom()[0].locus[QLocus].genotypic_val_mat(pState+1,mState+1);
02065         }
02066         double y = record()[0].double_val();
02067         pen = isqrt2pi/std::sqrt(v)*std::exp(-(y-mu)*(y-mu));
02068         return pen;
02069 }
02070 /*! \fn double Individual::getGenoPenetrance(void)
02071 *  \brief returns the genotype penetrance for the RSampler
02072 */
02073 
02074 double Individual::getDisGenoPenetrance(void)
02075 {
02076         // Authors: L. Radu Totir and Rohan L. Fernando 
02077         // (August, 2004) 
02078         // Contributors: 
02079         double mu=0.0,d,pen;
02080         unsigned mState, pState;
02081         if (myrecord[0].missing) {
02082                 return 1.0;
02083         }
02084         mState    = genotNodeVector[currentLocus].getmState();
02085         pState    = genotNodeVector[currentLocus].getpState();
02086         int phenotype = int(record()[0].double_val());
02087         // cout << phenotype << "   " << mState+pState << "  " 
02088         // << population->disPenetranceTable[phenotype][mState+pState] 
02089         // << endl;
02090         return population->disPenetranceTable[phenotype][mState+pState];
02091 }
02092 
02093 double Individual::getAllelePTransmissionProb(void){
02094         // Authors: L. Radu Totir and Rohan L. Fernando 
02095         // (June, 2003) 
02096         // Contributors:
02097         // cout << currentLocus << endl;
02098         unsigned ip = palleleStateNodeVector[currentLocus].getMyAlleleState();
02099         unsigned fm = myfather->malleleStateNodeVector[currentLocus].getMyAlleleState();
02100         unsigned fp = myfather->palleleStateNodeVector[currentLocus].getMyAlleleState();
02101         // cout << ip << " " << fm << " " << fp <<  endl;
02102         return getTransmissionProb(ip,fm,fp,p_gamete);
02103 }
02104 /*! \fn double Individual::getAllelePTransmissionProb()
02105 *  \brief returns the transmission probability of a paternal allele
02106 in the RSampler
02107 */
02108 
02109 double Individual::getAlleleMTransmissionProb(void){
02110         // Authors: L. Radu Totir and Rohan L. Fernando 
02111         // (June, 2003) 
02112         // Contributors:
02113         unsigned im = malleleStateNodeVector[currentLocus].getMyAlleleState();
02114         // cout << "Maternal allele for individual: " << myid << " is " << im << endl;
02115         unsigned mm = mymother->malleleStateNodeVector[currentLocus].getMyAlleleState();
02116         unsigned mp = mymother->palleleStateNodeVector[currentLocus].getMyAlleleState();
02117         return getTransmissionProb(im,mm,mp,m_gamete);
02118 }
02119 /*! \fn double Individual::getAlleleMTransmissionProb()
02120 *  \brief returns the transmission probability of a maternal allele
02121 in the RSampler
02122 */
02123 
02124 double Individual::getAllelePRTransmissionProb(unsigned forLocus){
02125         // Authors: L. Radu Totir and Rohan L. Fernando 
02126         // (September, 2004) 
02127         // Contributors:
02128         unsigned ipState = palleleStateNodeVector[forLocus].getMyAlleleState();
02129         unsigned fmState = myfather->malleleStateNodeVector[forLocus].getMyAlleleState();
02130         unsigned fpState = myfather->palleleStateNodeVector[forLocus].getMyAlleleState();
02131         unsigned ipOrigin = palleleOriginNodeVector[forLocus].getMyAlleleOrigin();
02132         // cout << ipState << " " << fmState << " " << fpState 
02133         // << " " << ipOrigin << endl;
02134         if(ipOrigin==0){
02135                 if(ipState==fmState){
02136                         return 1.0;
02137                 }
02138                 else {
02139                         return 0.0;
02140                 }
02141         }
02142         if(ipOrigin==1){
02143                 if(ipState==fpState){
02144                         return 1.0;
02145                 }
02146                 else {
02147                         return 0.0;
02148                 }
02149         }
02150 }
02151 /*! \fn double Individual::getAllelePTransmissionProb()
02152 *  \brief returns the transmission probability of a paternal allele 
02153 when the allele origin is assumed known. Used when we peel across 
02154 pedigree and loci jointly.
02155 */
02156 
02157 double Individual::getAlleleMRTransmissionProb(unsigned forLocus){
02158         // Authors: L. Radu Totir and Rohan L. Fernando 
02159         // (September, 2004) 
02160         // Contributors:
02161         unsigned imState = malleleStateNodeVector[forLocus].getMyAlleleState();
02162         unsigned mmState = mymother->malleleStateNodeVector[forLocus].getMyAlleleState();
02163         unsigned mpState = mymother->palleleStateNodeVector[forLocus].getMyAlleleState();
02164         unsigned imOrigin = malleleOriginNodeVector[forLocus].getMyAlleleOrigin();
02165         //cout << imState << " " << mmState << " " << mpState 
02166         //     << " " << imOrigin << endl;
02167         if(imOrigin==0){
02168                 if(imState==mmState){
02169                         return 1.0;
02170                 }
02171                 else {
02172                         return 0.0;
02173                 }
02174         }
02175         if(imOrigin==1){
02176                 if(imState==mpState){
02177                         return 1.0;
02178                 }
02179                 else {
02180                         return 0.0;
02181                 }
02182         }
02183 }
02184 /*! \fn double Individual::getAlleleMTransmissionProb()
02185 *  \brief returns the transmission probability of a maternal allele
02186 when the allele origin is assumed known. Used when we peel across 
02187 pedigree and loci jointly.
02188 */
02189 
02190 double Individual::getTransitionProb(void){
02191         // Authors: L. Radu Totir and Rohan L. Fernando 
02192         // (June, 2003) 
02193         // Contributors:
02194         double result;
02195         result = getMTransmissionProb()*getPTransmissionProb();
02196         return result;
02197 }
02198 /*! \fn double Individual::getTransitionProb()
02199 *  \brief returns the genotype transition probability in the
02200 RSampler
02201 */
02202 
02203 double Individual::getPTransmissionProb(void){
02204         // Authors: L. Radu Totir and Rohan L. Fernando 
02205         // (June, 2003) 
02206         // Contributors:
02207         unsigned ip = genotNodeVector[currentLocus].getpState();
02208         unsigned fm = myfather->genotNodeVector[currentLocus].getmState();
02209         unsigned fp = myfather->genotNodeVector[currentLocus].getpState();
02210         return getTransmissionProb(ip,fm,fp,p_gamete);
02211 }
02212 /*! \fn double Individual::getPTransmissionProb(void)
02213 *  \brief returns the paternal allele transmission probability for
02214 the RSampler for genotype peeling
02215 */
02216 
02217 double Individual::getMTransmissionProb(void){
02218         // Authors: L. Radu Totir and Rohan L. Fernando 
02219         // (June, 2003) 
02220         // Contributors:
02221         unsigned im = genotNodeVector[currentLocus].getmState();
02222         unsigned mm = mymother->genotNodeVector[currentLocus].getmState();
02223         unsigned mp = mymother->genotNodeVector[currentLocus].getpState();
02224         return getTransmissionProb(im,mm,mp,m_gamete);
02225 }
02226 /*! \fn double Individual::getMTransmissionProb(void)
02227 *  \brief returns the maternal allele transmission probability for
02228 the RSampler for genotype peeling
02229 */
02230 
02231 //double Individual::getTransmissionProb(unsigned iA, unsigned mA, unsigned pA,
02232 //                                                                         Vector <unsigned> &gamete){
02233 //      // Authors: L. Radu Totir and Rohan L. Fernando 
02234 //      // (June, 2003) 
02235 //      // Contributors:
02236 //      double temp = 1.0;
02237 //      if(mA==pA){
02238 //              if(iA==mA){
02239 //                      return 1.0;
02240 //              }
02241 //              else {
02242 //                      return 0.0;
02243 //              }
02244 //      }
02245 //      int segregationIndex = gamete[currentLocus]; //LRT 6/28/2004
02246 //      int leftLocus,rightLocus=-1;
02247 //      if(gamete.size()>1){
02248 //              leftLocus  = findLeftLocus(gamete);
02249 //              rightLocus = findRightLocus(gamete);
02250 //      }
02251 //      else{
02252 //              leftLocus = -1;
02253 //      }
02254 //      if(leftLocus==-1){
02255 //              if(iA==mA || iA==pA){
02256 //                      return 0.5;
02257 //              }
02258 //              else {
02259 //                      return 0.0;
02260 //              }
02261 //      }
02262 //      else {
02263 //              gamete[currentLocus] = 0;
02264 //              if(currentLocus>0){
02265 //                      if(gamete[leftLocus]==gamete[currentLocus]){
02266 //                              temp*=1 - population->recombinationMatrix[leftLocus][currentLocus];
02267 //                      }
02268 //                      else {
02269 //                              temp*=population->recombinationMatrix[leftLocus][currentLocus];
02270 //                      }
02271 //              }
02272 //              if(currentLocus<gamete.size()){
02273 //                      if(!(population->lookOnlyAtLocusToYourLeft) && rightLocus!=-1){
02274 //                              if(gamete[currentLocus]==gamete[rightLocus]){
02275 //                                      temp*=1 - population->recombinationMatrix[currentLocus][rightLocus];
02276 //                              }
02277 //                              else {
02278 //                                      temp*=population->recombinationMatrix[currentLocus][rightLocus];
02279 //                              }
02280 //                      }
02281 //              }
02282 //              double temp0 = temp;
02283 //              temp = 1.0;
02284 //              gamete[currentLocus] = 1;
02285 //              if(currentLocus>0){
02286 //                      if(gamete[leftLocus]==gamete[currentLocus]){
02287 //                              temp*=1 - population->recombinationMatrix[leftLocus][currentLocus];
02288 //                      }
02289 //                      else {
02290 //                              temp*=population->recombinationMatrix[leftLocus][currentLocus];
02291 //                      }
02292 //              }
02293 //              if(currentLocus<gamete.size()){
02294 //                      if(!(population->lookOnlyAtLocusToYourLeft) && rightLocus!=-1){
02295 //                              if(gamete[currentLocus]==gamete[rightLocus]){
02296 //                                      temp*=1 - population->recombinationMatrix[currentLocus][rightLocus];
02297 //                              }
02298 //                              else {
02299 //                                      temp*=population->recombinationMatrix[currentLocus][rightLocus];
02300 //                              }
02301 //                      }
02302 //              }
02303 //              gamete[currentLocus] = segregationIndex; //LRT 6/28/2004
02304 //              double temp1 = temp;
02305 //              double iSum  = 1.0/(temp0+temp1);
02306 //              if(iA==mA){
02307 //                      return temp0*iSum;
02308 //              }
02309 //              else if(iA==pA){
02310 //                      return temp1*iSum;
02311 //              }
02312 //              else {
02313 //                      return 0.0;
02314 //              }
02315 //      }
02316 //      return temp;
02317 //}
02318 ///*! \fn double Individual::getTransmissionProb(unsigned im, unsigned
02319 //mm, unsigned mp, Vector
02320 //<unsigned> &gamete)
02321 //*  \brief computes and returns the transmission probability for a
02322 //given trio of alleles taking in account flanking loci in the
02323 //RSampler
02324 //*/
02325 
02326 double Individual::getTransmissionProb(unsigned iA, unsigned mA, unsigned pA,
02327                                                                            Vector <unsigned> &gamete){
02328         // Authors: L. Radu Totir and Rohan L. Fernando 
02329         // (June, 2003) 
02330         // new version,  November, 2004 RLF
02331         // Contributors:
02332         double temp = 1.0;
02333         if(mA==pA){
02334                 if(iA==mA){
02335                         return 1.0;
02336                 }
02337                 else {
02338                         return 0.0;
02339                 }
02340         }
02341         int segregationIndex = gamete[currentLocus]; //LRT 6/28/04
02342         int     leftLocus  = findLeftLocus(gamete);
02343         int     rightLocus = findRightLocus(gamete);
02344         if(leftLocus==-1 && rightLocus==-1){
02345                 if(iA==mA || iA==pA){
02346                         return 0.5;
02347                 }
02348                 else {
02349                         return 0.0;
02350                 }
02351         }
02352         else {
02353                 gamete[currentLocus] = 0;
02354                 if (leftLocus!=-1) {
02355                         if(gamete[leftLocus]==gamete[currentLocus]){
02356                                 temp*=1 - population->recombinationMatrix[leftLocus][currentLocus];
02357                         }
02358                         else {
02359                                 temp*=    population->recombinationMatrix[leftLocus][currentLocus];
02360                         }
02361                 }
02362                 if(rightLocus!=-1){
02363                         if(gamete[currentLocus]==gamete[rightLocus]){
02364                                 temp*=1 - population->recombinationMatrix[currentLocus][rightLocus];
02365                         }
02366                         else {
02367                                 temp*=    population->recombinationMatrix[currentLocus][rightLocus];
02368                         }
02369                 }
02370                 double temp0 = temp;
02371                 temp = 1.0;
02372                 
02373                 gamete[currentLocus] = 1;
02374                 if (leftLocus!=-1) {
02375                         if(gamete[leftLocus]==gamete[currentLocus]){
02376                                 temp*=1 - population->recombinationMatrix[leftLocus][currentLocus];
02377                         }
02378                         else {
02379                                 temp*=    population->recombinationMatrix[leftLocus][currentLocus];
02380                         }
02381                 }
02382                 if(rightLocus!=-1){
02383                         if(gamete[currentLocus]==gamete[rightLocus]){
02384                                 temp*=1 - population->recombinationMatrix[currentLocus][rightLocus];
02385                         }
02386                         else {
02387                                 temp*=    population->recombinationMatrix[currentLocus][rightLocus];
02388                         }
02389                 }
02390                 double temp1 = temp;
02391                 gamete[currentLocus] = segregationIndex; //LRT 6/28/2004
02392                 double iSum  = 1.0/(temp0+temp1);
02393                 if(iA==mA){
02394                         return temp0*iSum;
02395                 }
02396                 else if(iA==pA){
02397                         return temp1*iSum;
02398                 }
02399                 else {
02400                         return 0.0;
02401                 }
02402         }
02403 }
02404 
02405         /*! \fn double Individual::getTransmissionProb(unsigned im, unsigned mm, unsigned mp, Vector<unsigned> &gamete)
02406 *   \brief computes and returns the transmission probability for a given trio of alleles taking in account flanking loci in the
02407 RSampler
02408 */
02409 
02410 
02411 double Individual::getMatthiasTransitionProb(void){
02412         // Authors: L. Radu Totir and Rohan L. Fernando 
02413         // (August, 2004) 
02414         // Contributors:
02415         double result;
02416         result = getMatthiasMTransmissionProb()*getMatthiasPTransmissionProb();
02417         return result;
02418 }
02419 /*! \fn double Individual::getMatthiasTransitionProb()
02420 *  \brief returns the genotype transition probability for the proposal
02421 *  in the RSampler. The cascading origin trick of Matthias Schelling.
02422 */
02423 
02424 double Individual::getMatthiasPTransmissionProb(void){
02425         // Authors: L. Radu Totir and Rohan L. Fernando 
02426         // (August, 2004) 
02427         // Contributors:
02428         unsigned ip = genotNodeVector[currentLocus].getpState();
02429         unsigned fm = myfather->genotNodeVector[currentLocus].getmState();
02430         unsigned fp = myfather->genotNodeVector[currentLocus].getpState();
02431         return getMatthiasTransmissionProb(ip,fm,fp,p_gamete,p_gameteOld);
02432 }
02433 /*! \fn double Individual::getMatthiasPTransmissionProb(void)
02434 *  \brief returns the paternal allele transmission probability for
02435 *   the proposal in the RSampler for genotype peeling. 
02436 */
02437 
02438 double Individual::getMatthiasMTransmissionProb(void){
02439         // Authors: L. Radu Totir and Rohan L. Fernando 
02440         // (August, 2004) 
02441         // Contributors:
02442         unsigned im = genotNodeVector[currentLocus].getmState();
02443         unsigned mm = mymother->genotNodeVector[currentLocus].getmState();
02444         unsigned mp = mymother->genotNodeVector[currentLocus].getpState();
02445         return getMatthiasTransmissionProb(im,mm,mp,m_gamete,m_gameteOld);
02446 }
02447 /*! \fn double Individual::getMatthiasMTransmissionProb(void)
02448 *  \brief returns the maternal allele transmission probability for
02449 *   the proposal in the RSampler for genotype peeling.
02450 */
02451 
02452 double Individual::getMatthiasTransmissionProb(unsigned iA, unsigned mA, unsigned pA,
02453                                                                                            Vector <unsigned> &gamete, 
02454                                                                                            Vector <unsigned> &gameteOld){
02455         // Authors: L. Radu Totir and Rohan L. Fernando 
02456         // (August, 2004) 
02457         // Contributors:
02458         double temp = 1.0;
02459         double rCOLeft,rCORight;
02460         if(mA==pA){
02461                 if(iA==mA){
02462                         return 1.0;
02463                 }
02464                 else {
02465                         return 0.0;
02466                 }
02467         }
02468         int segregationIndex = gamete[currentLocus]; //LRT 6/28/2004
02469         int leftLocus,rightLocus=-1;
02470         if(gamete.size()>1){
02471                 leftLocus  = findLeftLocus(gamete);
02472                 rightLocus = findRightLocus(gamete);
02473         }
02474         else{
02475                 leftLocus = -1;
02476         }
02477         if(leftLocus==-1){
02478                 if(iA==mA || iA==pA){
02479                         return 0.5;
02480                 }
02481                 else {
02482                         return 0.0;
02483                 }
02484         }
02485         else {
02486                 if(gameteOld[leftLocus]==gameteOld[currentLocus]){
02487                         rCOLeft = population->recombinationMatrix[leftLocus][currentLocus];
02488                 }
02489                 else {
02490                         rCOLeft = 1 - population->recombinationMatrix[leftLocus][currentLocus];
02491                 }
02492                 gamete[currentLocus] = 0;
02493                 if(currentLocus>0){
02494                         if(gamete[leftLocus]==gamete[currentLocus]){
02495                                 temp *= 1 - rCOLeft;
02496                         }
02497                         else {
02498                                 temp *= rCOLeft;
02499                         }
02500                 }
02501                 if(currentLocus<gamete.size()){
02502                         if(!(population->lookOnlyAtLocusToYourLeft) && rightLocus!=-1){
02503                                 if(gamete[currentLocus]==gamete[rightLocus]){
02504                                         temp*=1 - population->recombinationMatrix[currentLocus][rightLocus];
02505                                 }
02506                                 else {
02507                                         temp*=population->recombinationMatrix[currentLocus][rightLocus];
02508                                 }
02509                         }
02510                 }
02511                 double temp0 = temp;
02512                 temp = 1.0;
02513                 gamete[currentLocus] = 1;
02514                 if(currentLocus>0){
02515                         if(gamete[leftLocus]==gamete[currentLocus]){
02516                                 temp *= 1 - rCOLeft;
02517                         }
02518                         else {
02519                                 temp *= rCOLeft;
02520                         }
02521                 }
02522                 if(currentLocus<gamete.size()){
02523                         if(!(population->lookOnlyAtLocusToYourLeft) && rightLocus!=-1){
02524                                 if(gamete[currentLocus]==gamete[rightLocus]){
02525                                         temp*=1 - population->recombinationMatrix[currentLocus][rightLocus];
02526                                 }
02527                                 else {
02528                                         temp*=population->recombinationMatrix[currentLocus][rightLocus];
02529                                 }
02530                         }
02531                 }
02532                 gamete[currentLocus] = segregationIndex; //LRT 6/28/2004
02533                 double temp1 = temp;
02534                 double iSum  = 1.0/(temp0+temp1);
02535                 if(iA==mA){
02536                         return temp0*iSum;
02537                 }
02538                 else if(iA==pA){
02539                         return temp1*iSum;
02540                 }
02541                 else {
02542                         return 0.0;
02543                 }
02544         }
02545         return temp;
02546 }
02547 /*! \fn double Individual::getMatthiasTransmissionProb(unsigned im, 
02548 unsigned mm, unsigned mp, Vector <unsigned> &gamete, 
02549 Vector <unsigned> &gameteOld)
02550 *  \brief computes and returns the transmission probability for a
02551 given trio of alleles taking in account flanking loci in the
02552 RSampler
02553 */
02554 
02555 int Individual::findLeftLocus(Vector <unsigned> &gamete){ 
02556         for (int i = currentLocus-1; i>=0; i--){
02557                 if(population->prior->chrom()[0].locus[i].qtl_ml!='q'){
02558                         if(gamete[i]!=9999){
02559                                 return i;
02560                         }
02561                 }
02562         }
02563         return -1;
02564 }
02565 /*! \fn inline int Individual::findLeftLocus(Vector <unsigned> &gamete)
02566 *  \brief finds the first heterozygous marker locus to the left
02567 of the current locus 
02568 */
02569 
02570 int Individual::findRightLocus(Vector <unsigned> &gamete){ 
02571         if (population->lookOnlyAtLocusToYourLeft) return -1;
02572         for (unsigned i = currentLocus+1; i < numLoci ; i++){
02573                 if(population->prior->chrom()[0].locus[i].qtl_ml!='q'){
02574                         if(gamete[i]!=9999){
02575                                 return i;
02576                         }
02577                 }
02578         }
02579         return -1;
02580 }
02581 /*! \fn inline int Individual::findRightLocus(Vector <unsigned> &gamete)
02582 *  \brief finds the first heterozygous marker locus to the right 
02583 of the current locus
02584 */
02585 
02586 ostream& operator<<(ostream& stream, Individual& i) {
02587         std::vector<unsigned> matHaplotype, patHaplotype;
02588         matHaplotype.resize(i.numLoci);
02589         patHaplotype.resize(i.numLoci);
02590         for (unsigned j=0; j<i.numLoci; j++){  
02591                 //     if (i.population->samplerType=="genotypic"){
02592                 matHaplotype[j] = i.genotNodeVector[j].getmState()+1;
02593                 patHaplotype[j] = i.genotNodeVector[j].getpState()+1;
02594                 //     }
02595                 //     else if (i.population->samplerType=="allelic"){
02596                 //       matHaplotype[j] = i.malleleStateNodeVector[j].getMyAlleleState();
02597                 //       patHaplotype[j] = i.palleleStateNodeVector[j].getMyAlleleState();
02598                 //     }
02599         }
02600         for (unsigned j=0; j<i.numLoci; j++){  
02601                 stream << matHaplotype[j] << "  ";
02602         }
02603         stream << endl;
02604         for (unsigned j=0; j<i.numLoci; j++){  
02605                 stream << "-" << "  ";
02606         }
02607         stream << endl;
02608         for (unsigned j=0; j<i.numLoci; j++){  
02609                 stream << patHaplotype[j] << "  ";
02610         }
02611         stream << endl;
02612         return stream;
02613 }
02614 
02615 void Individual::displayGenotypes(unsigned locus){
02616         unsigned numGenos = genotNodeVector[locus].genotypeVector.size();
02617         cout << "Individual: " << myid << endl;
02618         for (unsigned i=0;i<numGenos;i++){
02619                 cout << genotNodeVector[locus].genotypeVector[i].maternal+1;
02620                 cout << "/";
02621                 cout << genotNodeVector[locus].genotypeVector[i].paternal+1 << endl;
02622         }
02623 }
02624 
02625 void Individual::displayAlleleVectors(unsigned locus){
02626         cout << "Individual: " << myid << " Locus: " << locus << endl;  
02627         cout << "maternal allele states" << endl;
02628         unsigned sizeOfVector = malleleStateNodeVector[locus].alleleStateVector.size();
02629         for (unsigned i=0; i<sizeOfVector; i++){
02630                 cout << malleleStateNodeVector[locus].alleleStateVector[i] << endl;
02631         }
02632         cout << "paternal allele states" << endl;
02633         sizeOfVector = palleleStateNodeVector[locus].alleleStateVector.size();
02634         for (unsigned i=0; i<sizeOfVector; i++){
02635                 cout << palleleStateNodeVector[locus].alleleStateVector[i] << endl;
02636         }
02637         
02638 }
02639 
02640 void Individual::setAlleleStateVectors(unsigned locus){
02641         set<unsigned> matAlleles, patAlleles;
02642         malleleStateNodeVector[locus].alleleStateVector.clear();
02643         palleleStateNodeVector[locus].alleleStateVector.clear();        
02644         unsigned numGenos = genotNodeVector[locus].genotypeVector.size();
02645         for (unsigned i=0;i<numGenos;i++){
02646                 matAlleles.insert(genotNodeVector[locus].genotypeVector[i].maternal+1);
02647                 patAlleles.insert(genotNodeVector[locus].genotypeVector[i].paternal+1);
02648         }
02649         set<unsigned>::iterator it;
02650         for (it=matAlleles.begin();it!=matAlleles.end();it++){
02651                 malleleStateNodeVector[locus].alleleStateVector.push_back(*it);
02652         }
02653         for (it=patAlleles.begin();it!=patAlleles.end();it++){
02654                 palleleStateNodeVector[locus].alleleStateVector.push_back(*it);
02655         }               
02656 }
02657 
02658 void Individual::setAlleleOriginVectors(unsigned locus){
02659         SafeSTLVector<unsigned> intersectionVector;
02660         unsigned msize, psize;
02661         
02662         set_intersection(myfather->malleleStateNodeVector[locus].alleleStateVector.begin(),
02663                                          myfather->malleleStateNodeVector[locus].alleleStateVector.end(),
02664                                          palleleStateNodeVector[locus].alleleStateVector.begin(),
02665                                          palleleStateNodeVector[locus].alleleStateVector.end(),
02666                                          inserter(intersectionVector,intersectionVector.begin())  );
02667         msize = intersectionVector.size();
02668         intersectionVector.clear();
02669         set_intersection(myfather->palleleStateNodeVector[locus].alleleStateVector.begin(),
02670                                          myfather->palleleStateNodeVector[locus].alleleStateVector.end(),
02671                                          palleleStateNodeVector[locus].alleleStateVector.begin(),
02672                                          palleleStateNodeVector[locus].alleleStateVector.end(),
02673                                          inserter(intersectionVector,intersectionVector.begin())  );
02674         psize = intersectionVector.size();
02675         
02676         if (msize==0 && psize>0) {
02677                 palleleOriginNodeVector[locus].alleleOriginVector.clear();
02678                 palleleOriginNodeVector[locus].alleleOriginVector.push_back(1);
02679         }
02680         else if(msize>0 && msize==0) {
02681                 palleleOriginNodeVector[locus].alleleOriginVector.clear();
02682                 palleleOriginNodeVector[locus].alleleOriginVector.push_back(0);
02683         }
02684         else {
02685                 throw exception("Individual::setAlleleOriginVectors(locus) possible bug in program");
02686         }
02687         
02688         intersectionVector.clear();
02689         set_intersection(mymother->malleleStateNodeVector[locus].alleleStateVector.begin(),
02690                                          mymother->malleleStateNodeVector[locus].alleleStateVector.end(),
02691                                          malleleStateNodeVector[locus].alleleStateVector.begin(),
02692                                          malleleStateNodeVector[locus].alleleStateVector.end(),
02693                                          inserter(intersectionVector,intersectionVector.begin())  );
02694         msize = intersectionVector.size();
02695         intersectionVector.clear();
02696         set_intersection(mymother->palleleStateNodeVector[locus].alleleStateVector.begin(),
02697                                          mymother->palleleStateNodeVector[locus].alleleStateVector.end(),
02698                                          malleleStateNodeVector[locus].alleleStateVector.begin(),
02699                                          malleleStateNodeVector[locus].alleleStateVector.end(),
02700                                          inserter(intersectionVector,intersectionVector.begin())  );
02701         psize = intersectionVector.size();
02702         if (msize==0 && psize>0) {
02703                 malleleOriginNodeVector[locus].alleleOriginVector.clear();
02704                 malleleOriginNodeVector[locus].alleleOriginVector.push_back(1);
02705         }
02706         else if(msize>0 && msize==0) {
02707                 malleleOriginNodeVector[locus].alleleOriginVector.clear();
02708                 malleleOriginNodeVector[locus].alleleOriginVector.push_back(0);
02709         }
02710         else {
02711                 throw exception("Individual::setAlleleOriginVectors(locus) possible bug in program");
02712         }
02713 }
02714 
02715 
02716 
02717 } /// end of matvec namespace