C++ Patron::getIndividual方法代码示例

void Restaurant::initializeAdmixtureProportions(void) {

	// set up vectors for admixture
	if (assumingAdmixture == true && assumingDpp == false)
		{
		// what type of prior are we placing on the mixing proportions?
		if (settingsPtr->getAdmixtureDist() == "fixed")
			admixtureVarianceParm = settingsPtr->getAdmixtureParm1();
		else if (settingsPtr->getAdmixtureDist() == "exponential")
			admixtureVarianceParm = ranPtr->exponentialRv( settingsPtr->getAdmixtureParm1() );
		else if (settingsPtr->getAdmixtureDist() == "uniform")
			admixtureVarianceParm = ranPtr->uniformRv( settingsPtr->getAdmixtureParm1(), settingsPtr->getAdmixtureParm2() );
		else if (settingsPtr->getAdmixtureDist() == "gamma")
			admixtureVarianceParm = ranPtr->gammaRv( settingsPtr->getAdmixtureParm1(), settingsPtr->getAdmixtureParm2() );
			
		// get the upper and lower limits for the admixture variance parameter
		if (settingsPtr->getAdmixtureDist() == "uniform")
			{
			admixtureLowerLimit = settingsPtr->getAdmixtureParm1();
			admixtureUpperLimit = settingsPtr->getAdmixtureParm2();
			}
		else if (settingsPtr->getAdmixtureDist() == "exponential")
			{
			admixtureLowerLimit = 0.001;
			admixtureUpperLimit = ( 1.0 / settingsPtr->getAdmixtureParm1() ) * 100.0;
			}
		else if (settingsPtr->getAdmixtureDist() == "gamma")
			{
			admixtureLowerLimit = 0.001;
			admixtureUpperLimit = ( settingsPtr->getAdmixtureParm1() / settingsPtr->getAdmixtureParm2() ) * 100.0;
			}
			
		// initialize the counts and probabilities for the mixing proportions
		popCounts = new MbVector<int>[observationsPtr->getNumIndividuals()];
		popProbs  = new MbVector<double>[observationsPtr->getNumIndividuals()];
		for (int i=0; i<observationsPtr->getNumIndividuals(); i++)
			{
			popCounts[i] = MbVector<int>(modelPtr->getNumPopulations());
			popProbs[i] = MbVector<double>(modelPtr->getNumPopulations());
			MbVector<double> alp( modelPtr->getNumPopulations() );
			for (int j=0; j<modelPtr->getNumPopulations(); j++)
				alp[j] = admixtureVarianceParm / modelPtr->getNumPopulations();
			ranPtr->dirichletRv(alp, popProbs[i]);
			for (int j=0; j<modelPtr->getNumPopulations(); j++)
				popCounts[i][j] = 0;
			}
		
		for (int i=0; i<getNumTables(); i++)
			{
			Table *tp = getTable(i);
			for (int j=0; j<tp->getNumPatrons(); j++)
				{
				Patron *p = tp->getPatron(j);
				popCounts[ p->getIndividual() ][ i ]++;
				}
			}
			
		}
}

void Restaurant::updateSeatingForFixedTables(void) {

	// allocate memory for reassignment probabilities
	vector<double> probs;
	
	// loop over all patrons, attempting to reseat each
	for (int idx=0; idx<numPatrons; idx++)
		{
		// get a pointer to the element and its current table
		Patron *p = patrons[idx];
		Table *tablePtr = p->getTable();

		// remove the element from the table
		tablePtr->removePatronFromTable(p);
		if (assumingAdmixture == true)
			popCounts[ p->getIndividual() ][ getIndexForTable(tablePtr) ]--;
		
		// calculate the likelihood when the element is seated at all possible tables
		for (int i=0; i<getNumTables(); i++)
			{
			Table *tp = getTable(i);
			double lnC = 0.0;
			if (assumingAdmixture == true)
				lnC = log( popProbs[ p->getIndividual() ][ i ] );
			probs.push_back( tp->lnLikelihood(p) + lnC );
			}
		normalizeVector( probs );

		// pick a table for reassignment
		tablePtr = getTable( multinomialRv(probs) );
			
		// add the element to that table
		tablePtr->addPatronToTable(p);
		if (assumingAdmixture == true)
			popCounts[ p->getIndividual() ][ getIndexForTable(tablePtr) ]++;
			
		// clear the probs vector
		probs.clear();
		}
}

double Restaurant::lnPriorProbability(void) {

	if (assumingAdmixture == false && assumingDpp == false)
		{
		// fixed number of populations with no admixture
		return -numPatrons * log(getNumTables());
		}
	else if (assumingAdmixture == true && assumingDpp == false)
		{
		// fixed number of poulations with admixture
		double lnP = 0.0;
		for (int idx=0; idx<numPatrons; idx++)
			{
			Patron *p = patrons[idx];
			Table *tablePtr = p->getTable();
			lnP += log(popProbs[ p->getIndividual() ][ getIndexForTable(tablePtr) ]);
			}
		return lnP;
		}
	else if (assumingAdmixture == false && assumingDpp == true)
		{
		// number of populations is a random variable, but no admixture
		/* f({\mathbf z}, k | \alpha, c) = \alpha^k {\prod_{i=1}^k(\eta_i - 1)! \over \prod_{i=1}^c(\alpha + i - 1)} */
		double lnSum1 = 0.0;
		for (vector<Table *>::iterator p=tables.begin(); p != tables.end(); p++)
			{
			lnSum1 += lnFactorial[ (*p)->getNumPatrons() - 1 ];
			}
		double lnSum2 = 0.0;
		for (int i=1; i<=numPatrons; i++)
			lnSum2 += log(concentrationParm + i - 1.0);
		double lnP = getNumTables() * log(concentrationParm) + lnSum1 - lnSum2;
		return lnP;
		}
	return 0.0;
}

void Restaurant::initializePatrons(void) {

	// how many patrons?
	if (assumingAdmixture == false)
		numPatrons = observationsPtr->getNumIndividuals();
	else if (assumingAdmixture == true && individualsInRestaurant == -1)
		numPatrons = observationsPtr->getNumIndividuals() * observationsPtr->getTotalNumAlleleCopies();
	else if (assumingAdmixture == true && individualsInRestaurant >= 0)
		numPatrons = observationsPtr->getTotalNumAlleleCopies();
		
	// allocate the patrons
	patrons = new Patron*[numPatrons];

	// set the index information for each patron
	if (assumingAdmixture == false)
		{
		for (int i=0; i<numPatrons; i++)
			{
			patrons[i] = new PatronIndividual( observationsPtr, i, observationsPtr->getNumLoci(), observationsPtr->getNumUniqueAlleles() );
			patrons[i]->setIndex(i);
			patrons[i]->setIndividual(i);
			}
		}
	else if (assumingAdmixture == true && individualsInRestaurant == -1)
		{
		for (int i=0, m=0; i<observationsPtr->getNumIndividuals(); i++)
			{
			for (int j=0; j<observationsPtr->getNumLoci(); j++)
				{
				int ploidy = 1;
				if (observationsPtr->getIsDiploid(j) == true)
					ploidy = 2;
				for (int k=0; k<ploidy; k++)
					{
					patrons[m] = new PatronAllele( observationsPtr, m, i, j, k, observationsPtr->getNumUniqueAllelesAtLocus(j) );
					m++;
					}
				}
			}
		}
	else if (assumingAdmixture == true && individualsInRestaurant >= 0)
		{
		for (int j=0, m=0; j<observationsPtr->getNumLoci(); j++)
			{
			int ploidy = 1;
			if (observationsPtr->getIsDiploid(j) == true)
				ploidy = 2;
			for (int k=0; k<ploidy; k++)
				{
				patrons[m] = new PatronAllele( observationsPtr, m, individualsInRestaurant, j, k, observationsPtr->getNumUniqueAllelesAtLocus(j) );
				m++;
				}
			}
		}
		
	// set the data for each patron
	for (int i=0; i<numPatrons; i++)
		{
		Patron *p = patrons[i];
		int elemIdx = p->getIndividual();
		p->setData( observationsPtr->getPtrToIndividual(elemIdx) );
		}

	// set the ln likelihood for each data element when placed in a population by itself
	for (int i=0; i<numPatrons; i++)
		{
		Patron *p = patrons[i];
		double x = lnLikelihoodNewTable(p);
		p->setLnLikelihoodAlone(x);
		}		
}