C++ Gene类代码示例

    std::pair<Chromosome, Breeder::BreedStats> Breeder::breed(const Chromosome& mum, const Chromosome& dad, uint32_t crossOverPerMillion)
    {
        BreedStats stats;
        stats.crossOvers = 0;

        std::random_device rd;
        std::uniform_int_distribution<> randomMillion(0, 999999);
        std::uniform_int_distribution<> mumOrDad(0, 1);
        bool startWithMum = mumOrDad(rd) == 0;

        const Gene* currentGene = startWithMum ? &mum.gene() : &dad.gene();
        const Gene* otherGene = startWithMum ? &dad.gene() : &mum.gene();

        Gene resultingGene;
        resultingGene.resize(currentGene->size());

        for(uint32_t i = 0; i < resultingGene.size(); i++)
        {
            resultingGene[i] = (*currentGene)[i];

            if(crossOverPerMillion > randomMillion(rd))
            {
                std::swap(currentGene, otherGene);
                ++stats.crossOvers;
            }
        }

        return {Chromosome(mum.maxCistronIds(), mum.informationDensity(), mum.startSequence(), mum.stopSequence(), resultingGene), stats};
    }

    void Chromosome::replaceCistronValues(const std::vector<Cistron>& cistrons)
    {
        CistronMap cistronMap;

        for(const auto& cistron : cistrons)
        {
            cistronMap.emplace(cistron.id(), cistron);
        }

        size_t start = 0;

        while((start = mGene.find(mStartSequence, start)) != std::string::npos)
        {
            size_t stop = mGene.find(mStopSequence, start);

            if(stop != std::string::npos)
            {
                size_t cistronStart = start + mStartSequence.size();
                Cistron cistron(mGene.substr(cistronStart, stop - cistronStart), mMaxCistronIds, mInformationDensity);

                if(cistronMap.count(cistron.id()))
                {
                    Gene toInsert = mStartSequence;
                    toInsert.append(cistronMap.at(cistron.id()).gene());
                    toInsert.append(mStopSequence);

                    mGene.replace(start, toInsert.size(), toInsert);
                }
            }

            start = stop;
        }
    }

void Host::RestoreHost(const string& sline)
{
	stringstream ssline(sline);

	ssline >> lociKIR;
	ssline >> lociMHC;
	ssline >>age;
	ssline >> tuning;
	ssline >> dead;
	ssline >> mutationRateHost;
	ssline >> inhibitoryKIRs;
	ssline >> activatingKIRs;
	int totalInfections;
	ssline >> totalInfections;

	/*use the streamstring from the position after "totalInfections"
	 *for that I copy the rest of the line into a string, and then reassign the streamstring with the content of
	 *for "restLine. For it to work, I have to clear the string first!
	 */
	string restLine;
	getline(ssline,restLine);
	ssline.clear();
	ssline.str(restLine);
	string infectionString;
	int type; double inf_time; double immune_time; double clearance_time;
	for(int i=0; i<totalInfections; i++)
	{
		Infection dummyInfection;
		infectionString = dummyInfection.RestoreInfection(ssline);
		infections.push_back(dummyInfection);
		ssline.clear();
		ssline.str(infectionString);
	}

	string geneString;
	//cout <<"mhc string>>>>>"<<ssline.str() <<endl;
	for(int i=0; i<lociMHC*TWO; i++)
	{
		Gene mhc;
		geneString = mhc.RestoreGenes(ssline);
		mhcGenes.push_back(mhc);
		ssline.clear();
		ssline.str(geneString);
	}

	string kirString;
	for(int i=0; i<lociKIR*TWO; i++)
	{
		//cout <<"kir string>>>>>"<<ssline.str() <<endl;
		KIRGene kir;
		kirString = kir.RestoreGenes(ssline);
		kirGenes.push_back(kir);
		ssline.clear();
		ssline.str(kirString);

	}
}

void SmallWorld::initialize() {
	
	for (int i = 0; i < bugNum; i++) {
		Gene gene;
		gene.initialize(minSpeed, maxSpeed, nodeNum * turn);
		genes[i] = gene;
	}
	
	generations++;
}

void FONSEModel::calculateLogLikelihoodRatioPerGroupingPerCategory(std::string grouping, Genome& genome, std::vector<double> &logAcceptanceRatioForAllMixtures)
{
	int numGenes = genome.getGenomeSize();
	//int numCodons = SequenceSummary::GetNumCodonsForAA(grouping);
	double likelihood = 0.0;
	double likelihood_proposed = 0.0;

	double mutation[5];
	double selection[5];
	double mutation_proposed[5];
	double selection_proposed[5];

	std::string curAA;

	Gene *gene;
	SequenceSummary *sequenceSummary;
	unsigned aaIndex = SequenceSummary::AAToAAIndex(grouping);

#ifdef _OPENMP
//#ifndef __APPLE__
	#pragma omp parallel for private(mutation, selection, mutation_proposed, selection_proposed, curAA, gene, sequenceSummary) reduction(+:likelihood,likelihood_proposed)
#endif
	for (unsigned i = 0u; i < numGenes; i++)
	{
		gene = &genome.getGene(i);
		sequenceSummary = gene->getSequenceSummary();
		if (sequenceSummary->getAACountForAA(aaIndex) == 0) continue;

		// which mixture element does this gene belong to
		unsigned mixtureElement = parameter->getMixtureAssignment(i);
		// how is the mixture element defined. Which categories make it up
		unsigned mutationCategory = parameter->getMutationCategory(mixtureElement);
		unsigned selectionCategory = parameter->getSelectionCategory(mixtureElement);
		unsigned expressionCategory = parameter->getSynthesisRateCategory(mixtureElement);
		// get phi value, calculate likelihood conditional on phi
		double phiValue = parameter->getSynthesisRate(i, expressionCategory, false);


		// get current mutation and selection parameter
		parameter->getParameterForCategory(mutationCategory, FONSEParameter::dM, grouping, false, mutation);
		parameter->getParameterForCategory(selectionCategory, FONSEParameter::dOmega, grouping, false, selection);

		// get proposed mutation and selection parameter
		parameter->getParameterForCategory(mutationCategory, FONSEParameter::dM, grouping, true, mutation_proposed);
		parameter->getParameterForCategory(selectionCategory, FONSEParameter::dOmega, grouping, true, selection_proposed);
		
		likelihood += calculateLogLikelihoodRatioPerAA(*gene, grouping, mutation, selection, phiValue);
		likelihood_proposed += calculateLogLikelihoodRatioPerAA(*gene, grouping, mutation_proposed, selection_proposed, phiValue);
	}
	//likelihood_proposed = likelihood_proposed + calculateMutationPrior(grouping, true);
	//likelihood = likelihood + calculateMutationPrior(grouping, false);

	logAcceptanceRatioForAllMixtures[0] = (likelihood_proposed - likelihood);
}

/* -----------------------------------------------------------------------
 * makeGeneVector():
 *
 * Returns a vector of n randomly generated genes of length CHROMOSOME_LENGTH. Gives each
 * object in the vector a fitness level of 0 to start it out.
 * -----------------------------------------------------------------------
 */
vector<Gene> makeGeneVector(int n)
{
    vector<Gene> geneVector;
    
    /* Make INITIAL_POPULATION amount of random genes and push it to firstNGenes */
    for (int i = 0; i < INITIAL_POPULATION; i++)
    {
        Gene g;
        g.setValues(makeRandomGene(), 0);
        geneVector.push_back(g);
    }
    
    return geneVector;
}

genotypes::genotypes(int genomSize) {

	Gene n;
	
	for(int i = 0; i < genomSize; i++){
		for(int j = 0; j < SEGMENT_SIZE; j++){
			if(rand()%2){
				n.flip(j);
			}
		}
		genom.push_back(n);
		n.reset();
	}

}

unsigned int Neuron::loadWeights(Gene const &gene, unsigned int index)
{
  unsigned int i(0);
  for (auto &s : _inputs)
    {
      s.second = gene.getWeight(index + i);
      i++;
    }
  return index + i;
}

	bool operator<(const Gene& g) const {
		if( (compareChar(chrom.c_str(), g.getChr().c_str()) < 0) ||
				(hs == g.getHS() && entrez_ID < g.getEID()) ||
				(hs == g.getHS() && entrez_ID == g.getEID() && start_p < g.getStart()) ||
				(hs == g.getHS() && entrez_ID == g.getEID() && start_p == g.getStart() && stop_p < g.getStop())
		)
			return true;
		else
			return false;
	}

Genome* SelfAdaptiveMutation::mutateProper(Genome* target) {
	std::vector<Gene*> genes = target->getGenome();
	std::vector<Gene*> results;
	for (unsigned int i = 0; i < this->initialGenomeLength; i++)
		results.push_back(genes[i]->copy());

	double rhoZero = HierRNG::gaussian(0, 1);

	for (unsigned int index: this->stdDevIndices) {
		Gene* target = genes[index];
		double stdDev = target->getIndex();
		double tauResult = this->tauPrime * rhoZero
			+ this->tau * HierRNG::gaussian(0, 1);
		stdDev *= this->useTauPlusOneCorrection ?
			pow(stdDev, tauResult) : exp(tauResult);
		results.push_back(target->copy(stdDev));
	}

	for (unsigned int i = 0; i < this->initialGenomeLength; i++) {
		Gene* original = results[i];
		double newIndex = original->getIndex();
		newIndex = newIndex
			+ results[this->initialGenomeLength + i]->getIndex()
			* HierRNG::gaussian(0, 1);
		results[i] = original->copy(newIndex);
		delete(original);
	}

	return new Genome(results, target->getSpeciesNode());
}

void main()
{
	CDiophantine dp(1, 2, 3, 4, 30);

	int ans;
	ans = dp.Solve();
	if (ans == -1)
	{
		cout << "No solution found." << endl;
	}
	else
	{
		Gene gn = dp.GetGene(ans);

		cout << "The solution set to a+2b+3c+4d=30 is:\n";
		cout << "a = " << gn.getAlleles()[0] << "." << endl;
		cout << "b = " << gn.getAlleles()[1] << "." << endl;
		cout << "c = " << gn.getAlleles()[2] << "." << endl;
		cout << "d = " << gn.getAlleles()[3] << "." << endl;
	}
	getchar();
}

void genetic::Generation::Generate(const genetic::Generation& ancestor)
{
  int& nCells = gainput.max_cells;
  m_cells = vector<Cell>(nCells, Cell());

  int iCell = 0;

  // Chose Elites
  for(; iCell < gainput.max_elite; ++iCell) m_cells[iCell] = ancestor.m_cells[iCell];

  // Selection to Next Generation
  while(iCell < nCells)
  {
    Gene child;

    if(drand(1.0) > gainput.thre_cloning)
    {
      child = ancestor.Select().Gen();
    }
    else
    {
      Cell father(ancestor.Select());
      Cell mother(ancestor.Select());
//    if(father.Fitness() > mother.Fitness())
        child = father.Gen() * mother.Gen();
//    else
//      child = mother.Gen() * father.Gen();
    }
    // Mutation
    if(drand(1.0) < gainput.thre_mutation) child.pMutate();
    if(drand(1.0) < gainput.thre_mutation) child.gMutate();

    m_cells[iCell++] = Cell(child);
  }

  ComputeProbability();
}

void Gene::singleCross(Gene& g1,Gene& g2,unsigned position){
	Gene g1temp = g1;
	unsigned word = position / bits_per_word;
	unsigned bit = position % bits_per_word;
	unsigned byte_before = word * sizeof(g1.data[0]);
	memcpy(g1temp.data, g2.data, byte_before);
	memcpy(g2.data, g1.data, byte_before);
	for(int i = 0 ; i < bit ; i++){
		if( (g1.data[word] & (BIT_MASK << i)) != (g2.data[word] & (BIT_MASK << i)) ){
			unsigned pos = i + word * bits_per_word;
			g1temp.toggle(pos);
			g2.toggle(pos);
		}
	}
	g1 = g1temp;
};

double genetic::Evaluate(const double& scale, const double& np, const Gene& gene, const Matrix& K)
{
  Matrix Kp;
  Permute(K, gene.Sequence(), Kp);

  double weight = 0.0;
  for(int i = 0; i < Kp.Nrows(); ++i)
    for(int j = i + 1; j < Kp.Ncols(); ++j)
    {
      double d = j - i;
//    weight += Kp.element(i, j) * pow(d, np);
      weight += Kp.element(i, j) * d * d;
    }

  return scale * weight;
}

void SmallWorld::breed() {
	
	if (bugs.size() > 2) {
		
		vector<Gene> tmp;
		tmp.resize(bugs.size());
		
		tmp[0] = bugs[0].getGene();
		tmp[1] = bugs[1].getGene();
		
		for (int i = 2; i < bugs.size() - 1; i++) {
			Gene gene;
			gene.crossover(tmp[0].dna, tmp[1].dna, nodeNum * turn);
			
			for (int j = 0; j < int(nodeNum / 4); j++) {
				gene.mutation(minSpeed, maxSpeed, probability);
			}
			
			tmp[i] = gene;
		}
		
		for (int i = bugs.size() - 2; i < bugs.size() - 1; i++) {
			Gene gene;
			gene.crossover(tmp[0].dna, tmp[1].dna, nodeNum * turn);
			
			for (int j = 0; j < int(nodeNum / 4); j++) {
				gene.mutation(minSpeed, maxSpeed, probability * 10);
			}
			
			tmp[i] = gene;
		}
		
		for (int i = bugs.size() - 1; i < bugs.size(); i++) {
			Gene gene;
			gene.initialize(minSpeed, maxSpeed, nodeNum * turn);
			tmp[i] = gene;
		}
		
		genes.swap(tmp);
	}
	
	generations++;
}