C++ Genome类代码示例

double Genome::GetGeneticalDistanceFrom(const Genome& other) const {
	double totalWeightDifference = 0.0;
	size_t numberOfOverlapingGenes = 0;

	size_t sizeOfSmallerGenome = min(this->GetGeneCount(), other.GetGeneCount());
	auto IsHistoricalMarkingSameAt = [&](size_t i) {
		return this->GetGeneAt(i).historicalMarking == other[i].historicalMarking;
	};

	for (size_t i = 0; i < sizeOfSmallerGenome && IsHistoricalMarkingSameAt(i); ++i) {
		totalWeightDifference += (double)abs(this->GetGeneAt(i).weight - other[i].weight);
		++numberOfOverlapingGenes;
	}

	auto numberOfDisjointGenes = this->GetGeneCount() + other.GetGeneCount() - (size_t)2 * numberOfOverlapingGenes;
	auto sizeOfBiggerGenome = max(this->GetGeneCount(), other.GetGeneCount());
	// TODO jnf: Think how we'll handle the next line
	auto disjointGenesInfluence = (double)numberOfDisjointGenes /* / (double)sizeOfBiggerGenome*/;

	auto averageWeightDifference = totalWeightDifference / (double)numberOfOverlapingGenes;

	disjointGenesInfluence *= (double)parameters.advanced.speciation.importanceOfDisjointGenes;
	averageWeightDifference *= (double)parameters.advanced.speciation.importanceOfAverageWeightDifference;

	return disjointGenesInfluence + averageWeightDifference;
}

Genome*
AverageCombinator::combine(Genome *x, Genome *y) {
	srand ( time(NULL) );

	int size = x->get_size( );
	Genome *z = new Genome();

	double chanceFactor = 0.5;
	for(int i = 0; i < size; i++) {
		double chance = (((double) rand()) / RAND_MAX+1);
		if(chance < chanceFactor) {
			double totalTime = x->get_gene(i)->getTime().getTimeInMinutes() + y->get_gene(i)->getTime().getTimeInMinutes();

			Time t;
			t.addMinute(floor((totalTime / 2) + 0.5));
			z->add_gene(x->get_gene(i)->getPlane(), t);
		} else if(chance < 0.75) {
			z->add_gene( x->get_gene(i)->getPlane(), x->get_gene(i)->getTime());
		} else {
			z->add_gene( y->get_gene(i)->getPlane(), y->get_gene(i)->getTime());
		}
	}

	return z;
}

void DownPropagator::propagateFitnesses(Individual ** population, int populationSize) {
	int tempFitness;
	Genome * tempGenome;
	GeneNode ** tempPools;
	int tempGenomeLength;
	int * tempIndexes;

	for (int i = 0; i < populationSize; i++) {
		tempFitness = population[i]->getFitness();
		
		tempGenome = population[i]->getGenome();
		tempPools = tempGenome->getGeneNodes();
		tempGenomeLength = tempGenome->getGenomeLength();
		tempIndexes = tempGenome->getGenome();		

		for (int k = 0; k < tempGenomeLength; k++) {
			if (tempPools[k]->getFitnessAtIndex(tempIndexes[k]) < tempFitness) {
				tempPools[k]->setFitnessAtIndex(tempIndexes[k], tempFitness);
			}

			//TODO: Find a more efficient (read: not n^2) way to do
			//this
			tempPools[k]->propagateFitnesses();
		}
	}
}

Genome Breeder::replicateGenome(const DNA& parent)
{
    Genome newGenome;

    const real32* parentWeights = parent.genome.readWeights();

    real32* weights = newGenome.editWeights();

    int32 mutationCount = 0;

    for (uint32 i = 0; i < newGenome.getLength(); i++) {

        const int32 dice = RandomGen::randomInt(0, parent.traits.mutationRate);
        if (dice >= 670) {
            weights[i] = parentWeights[i];
        }
        else if (dice >= 335) {
            weights[i] = parentWeights[i] + RandomGen::randomFloat(-0.5f, 0.5f);
            mutationCount++;
        }
        else {
            weights[i] = Genome::randomGenomeWeight();
            mutationCount++;
        }
    }

    std::stringstream sb;
    sb << "mutation count: " << mutationCount;
    Console::write(sb.str());

    return newGenome;
}

void LodExtract::writeSequences(const Genome* inParent,
                                const vector<const Genome*>& inChildren)
{
  vector<const Genome*> inGenomes = inChildren;
  inGenomes.push_back(inParent);
  const Genome* outParent = _outAlignment->openGenome(inParent->getName());
  (void)outParent;
  assert(outParent != NULL && outParent->getNumBottomSegments() > 0);
  string buffer;

  for (hal_size_t i = 0; i < inGenomes.size(); ++i)
  {
    const Genome* inGenome = inGenomes[i];
    Genome* outGenome = _outAlignment->openGenome(inGenome->getName());
    if (inGenome == inParent || outGenome->getNumChildren() == 0)
    {   
      SequenceIteratorConstPtr inSeqIt = inGenome->getSequenceIterator();
      SequenceIteratorConstPtr end = inGenome->getSequenceEndIterator();
      for (; inSeqIt != end; inSeqIt->toNext())
      {
        const Sequence* inSequence = inSeqIt->getSequence();
        if (inSequence->getSequenceLength() > 0)
        {
          Sequence* outSequence = outGenome->getSequence(inSequence->getName());
          assert(outSequence != NULL);
          inSequence->getString(buffer);
          outSequence->setString(buffer);
        }
      }
    }
  }
}

void Genome::crossover(const Genome &with,Genome &baby1,Genome &baby2)
{
	baby1.setFitness(0.0);
	baby2.setFitness(0.0);

	int crossoverPoint = NeuralNetRandom::RandomInt(0,numberOfGenes()-1);

	std::vector<double> newChromosome1;
	std::vector<double> newChromosome2;

	std::vector<double>::iterator myGenes = _chromosome.begin();
	std::vector<double> otherChromosome = with.chromosome();
	std::vector<double>::iterator otherGenes = otherChromosome.begin();

	if (crossoverPoint > 0)
	{
		newChromosome1.assign(myGenes,myGenes+crossoverPoint+1);
		newChromosome2.assign(otherGenes,otherGenes+crossoverPoint+1);
	}
	newChromosome1.insert(newChromosome1.end(),otherGenes+crossoverPoint,otherChromosome.end());
	newChromosome2.insert(newChromosome2.end(),myGenes+crossoverPoint,_chromosome.end());

	baby1.setChromosome(newChromosome1);
	baby2.setChromosome(newChromosome2);
}

void FullRun::displayResult(void) {
    unsigned long numRecombinant = 0;
    unsigned long numLongRec = 0;


    for (unsigned long i = 0; i < childDataset->getSize(); i++) {
        Genome* child = childDataset->getGenome(i);

        if (child->getRecombinantType() != Genome::Na_REC) {
            numRecombinant++;
        }

        if (child->getRecombinantType() == Genome::LONG_REC) {
            numLongRec++;
            if (fileLongRecombinants != NULL) {
                fileLongRecombinants->writeLine(child->getAccession());
            }
        }
    }


    Interface::instance()
            << "Number of triples tested :              " << getNumTotalTriplets() << "\n"
            << "Number of p-values computed exactly :   " << numComputedExactly << "\n"
            << "Number of p-values approximated (HS) :  " << numApproximated << "\n"
            << "Number of p-values not computed :       " << numSkipped << "\n"
            << endl
            << "Number of recombinant triplets :                               \t" << numRecombinantTriplets << "\n"
            << "Number of distinct recombinant sequences :                     \t" << numRecombinant << "\n";

    if (cloStopAtFirstRecombinant) {
        Interface::instance() << "[[ search stopped when first one found ]]\n" << endl;
    }

    if (!cloNoBpCalculated) {
        Interface::instance()
                << "Number of distinct recombinant sequences at least "
                << cloMinLongRecombinantLength << "nt long : \t" << numLongRec << "\n"
                << "Longest of short recombinant segments :                        \t"
                << longestRecombinantSegment << "nt\n";
    }
    Interface::instance().showOutput(true);

    Interface::instance() << Interface::SEPARATOR << endl;
    Interface::instance().showLog(true);

    char formatedPVal[20];
    sprintf(formatedPVal,
            "%1.3e",
            static_cast<double> (stats::correction::dunnSidak(minPVal, getNumTripletsForCorrection()))
            );
    Interface::instance()
            << "Rejection of the null hypothesis of clonal evolution at p = "
            << stats::correction::dunnSidak(minPVal, getNumTripletsForCorrection()) << "\n"
            << "                                                        p = " << formatedPVal << "\n"
            << "                                            Uncorrected p = " << minPVal << "\n"
            << "                                            Bonferroni  p = "
            << stats::correction::bonferroni(minPVal, getNumTripletsForCorrection()) << "\n";
    Interface::instance().showOutput(true);
}

void GenomeMetaTest::createCallBack(Alignment *alignment) {
    hal_size_t alignmentSize = alignment->getNumGenomes();
    CuAssertTrue(_testCase, alignmentSize == 0);

    Genome *ancGenome = alignment->addRootGenome("AncGenome", 0);

    MetaData *ancMeta = ancGenome->getMetaData();
    ancMeta->set("Young", "Jeezy");
}

SolverEvolver::Genome SolverEvolver::random_genome()
{
	Genome ret;
	for (int i = 0; i < ret.size(); ++i)
	{
		ret[i] = random_double(0.0f);
	}
	return std::move(ret);
}

void Trainer::init()
{
    for (unsigned int i = 0; i < N_AGENTS; i++)
    {
        Genome genome;
        genome.initRandomData();
        Agent *agent = new Agent(genome);
        agents.push_back(agent);
    }
}

void Genome::mutate(Genome &theGenome)
{
	std::vector<double> chromosome = theGenome.chromosome();
	for (int i=0;i<(int)chromosome.size();++i)
	{
		if (NeuralNetRandom::RandomFloat() < MutationRate)
			chromosome[i] += NeuralNetRandom::RandomClamped()*MaxMutationPerturbation;
	}
	theGenome.setChromosome(chromosome);
}

void TopSegmentSequenceTest::createCallBack(Alignment *alignment) {
    Genome *ancGenome = alignment->addRootGenome("Anc0", 0);
    vector<Sequence::Info> seqVec(1);
    seqVec[0] = Sequence::Info("Sequence", 1000000, 5000, 700000);
    ancGenome->setDimensions(seqVec);

    ancGenome->setSubString("CACACATTC", 500, 9);
    TopSegmentIteratorPtr tsIt = ancGenome->getTopSegmentIterator(100);
    tsIt->getTopSegment()->setCoordinates(500, 9);
}

int main() {
    const unsigned mutations=2;
    Genome::set_mutations_at_cloning(mutations);
    std::cout<<"[Cloning with "<<mutations<<" mutations.]"<<std::endl;
    
    Genome A;
    std::cout<<"HEALTHY GENOME: "<<std::endl<<" "<<A<<std::endl;
    std::cout<<" Bad mutations within first 10 years = "<<A.bad_mutations(10)<<std::endl;
    Genome B=A.clone();
    std::cout<<"AFTER 1st CLONING: "<<std::endl<<" "<<B<<std::endl;
    std::cout<<" Bad mutations within the whole Genome = "<<B.bad_mutations(Genome::max_age-1)<<std::endl;
    for (unsigned i=0; i<5; ++i) {
        A=B.clone();
        B=A.clone();
    }
    std::cout<<"AFTER 10 MORE CLONINGS: "<<std::endl<<" "<<A<<std::endl;
    std::cout<<" Bad mutations within first 5 years = "<<A.bad_mutations(5)<<std::endl;
    std::cout<<"AFTER 1 MORE CLONING: "<<std::endl<<" "<<B<<std::endl;
    std::cout<<" Bad mutations within first 10 years = "<<A.bad_mutations(10)<<std::endl;
    for (unsigned i=0; i<50; ++i) {
        A=B.clone();
        B=A.clone();
    }
    std::cout<<"AFTER 100 MORE CLONINGS: "<<std::endl<<" "<<A<<std::endl;
    std::cout<<" Bad mutations within first 2 years = "<<A.bad_mutations(2)<<std::endl;
    std::cout<<"AFTER 1 MORE CLONING: "<<std::endl<<" "<<B<<std::endl;
    std::cout<<" Bad mutations within first 2 years = "<<A.bad_mutations(2)<<std::endl;
    
    return 0;
}

Agent * add_agent(char * filename, Agent *mother) {
	Genome * myGenome;
	myGenome = new Genome(filename);
	Agent * newAgent;
	char portNum[256]; sprintf(portNum, "%i", curPort); 
	int pid = myGenome->buildOrganismFromGenome(portNum);
	newAgent = new Agent(myGenome, pid, curPort++);
	agents.push_back(newAgent);
	printf("loaded %s on port %i (%p)\n", filename, newAgent->get_port(), newAgent);
	return newAgent;
}

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