C++ Pattern::at方法代码示例

void ParallelInfo::sendZi(const Pattern &toSend)
{
   int P_Tag = 0;

#if defined(TIMING_MODE)
   long long elapsed, start;
   elapsed = 0;
#endif

   // 10-2-2002, Replaced MPI_Send with MPI_Rsend
   //    this assumes that the comp nodes have already
   //    called MPI_Recv, which they should have.
   for (unsigned int i = 1; i < P_NumNodes; i++) {
#if defined(TIMING_MODE)
      start = rdtsc();
#endif
      int * tmpToSend = new int[totalNumNrns];
      for (unsigned int nrn = 0; nrn < totalNumNrns; nrn++)
         tmpToSend[nrn] = toSend.at(nrn);
      MPI_Rsend(&tmpToSend[0], totalNumNrns, MPI_INT, i, P_Tag, MPI_COMM_WORLD);
#if defined(TIMING_MODE)
      elapsed += rdtsc() - start;
#endif
   }

#if defined(TIMING_MODE)
   Output::Out() << "Elapsed root_snd time = " << elapsed * 1.0 / TICKS_PER_SEC
                 << " seconds" << endl;
#endif   
}

void RateMeyerHaeseler::computeFuncDerv(double value, double &df, double &ddf) {
	int nseq = phylo_tree->leafNum;
	int nstate = phylo_tree->getModel()->num_states;
	int i, j, state1, state2;
//	double lh = 0.0;
	double trans, derv1, derv2;
	ModelSubst *model = phylo_tree->getModel();
	Pattern *pat = & phylo_tree->aln->at(optimizing_pattern);
	df = ddf = 0.0;
	for (i = 0; i < nseq-1; i++) if ((state1 = pat->at(i)) < nstate) 
		for (j = i+1; j < nseq; j++) if ((state2 = pat->at(j)) < nstate) {
			double dist = dist_mat[i*nseq + j];
			trans = model->computeTrans(value * dist, state1, state2, derv1, derv2);
//			lh -= log(trans);
			double t1 = derv1 / trans;
			double t2 = derv2 / trans;
			df -= t1 * dist;
			ddf -= dist * dist * (t2 - t1*t1);
		}
//	return lh;
}

double RateMeyerHaeseler::computeFunction(double value) {
	if (!rate_mh) {
		if (value != cur_scale) {
			ptn_tree->scaleLength(value/cur_scale);
			cur_scale = value;
			ptn_tree->clearAllPartialLH();
		}
		return -ptn_tree->computeLikelihood();
	}
	int nseq = phylo_tree->leafNum;
	int nstate = phylo_tree->getModel()->num_states;
	int i, j, state1, state2;
	double lh = 0.0;
	ModelSubst *model = phylo_tree->getModel();
	Pattern *pat = & phylo_tree->aln->at(optimizing_pattern);
	
	for (i = 0; i < nseq-1; i++) if ((state1 = pat->at(i)) < nstate) 
		for (j = i+1; j < nseq; j++) if ((state2 = pat->at(j)) < nstate) 
			lh -= log(model->computeTrans(value * dist_mat[i*nseq + j], state1, state2));
	return lh;
}

double RateMeyerDiscrete::computeFunction(double value) {
	if (!is_categorized) return RateMeyerHaeseler::computeFunction(value);
	if (!rate_mh) {
		if (value != cur_scale) {
			ptn_tree->scaleLength(value/cur_scale);
			cur_scale = value;
			ptn_tree->clearAllPartialLH();
		}
		return -ptn_tree->computeLikelihood();
	}

	double lh = 0.0;
	int nseq = phylo_tree->leafNum;
	int nstate = phylo_tree->getModel()->num_states;
	int i, j, k, state1, state2;
	ModelSubst *model = phylo_tree->getModel();
    int trans_size = nstate * nstate;
	double *trans_mat = new double[trans_size];
	int *pair_freq = new int[trans_size];

	for (i = 0; i < nseq-1; i++) 
		for (j = i+1; j < nseq; j++) {
			memset(pair_freq, 0, trans_size * sizeof(int));
			for (k = 0; k < size(); k++) {
				if (ptn_cat[k] != optimizing_cat) continue;
				Pattern *pat = & phylo_tree->aln->at(k);
				if ((state1 = pat->at(i)) < nstate && (state2 = pat->at(j)) < nstate)
					pair_freq[state1*nstate + state2] += pat->frequency;
			}
			model->computeTransMatrix(value * dist_mat[i*nseq + j], trans_mat);
			for (k = 0; k < trans_size; k++) if (pair_freq[k])
				lh -= pair_freq[k] * log(trans_mat[k]);
		}
	delete [] pair_freq;
	delete [] trans_mat;
	return lh;
}

void ActiveConnect(ArgListType &arg) //AT_FUN
{
   // process the function arguments
   static string FunctionName = "ActiveConnect";
   static IntArg StartNeuron("-Nstart", "start neuron", 1);
   static IntArg EndNeuron("-Nend", "end neuron{-1 gives last neuron}", -1);
   static IntArg StartPat("-Pstart", "start pattern(or vector)", 1);
   static IntArg EndPat("-Pend", "end pattern(or vector) {-1 gives end time}",
                        -1);
   static StrArg SeqName("-name", "sequence name");
   static StrArg ResetSeq("-reset", "reset sequence", "{zeros}");
   static StrArg External("-ext", "external sequence", "{zeros}");
   static FlagArg DoSummary("-sum", "-nosum", "print a summary", 1);
   static int argunset = true;
   static CommandLine ComL(FunctionName);
   if (argunset) {
      ComL.HelpSet("@ActiveConnect() \n"
                   "\tGenerates the data about conditional probabilities of firing\n"
                   "\tand saves the results into the ActiveConBuffer.\n"
                   "\tThe reset sequence is used for the 0 time step and\n"
                   "\tif externals are given then they are ignored.\n"
                   "\tActiveConBuffer is a matrix with 6 Vectors:\n"
                   "\t\t1: TimeStep(Pattern Number)\n"
                   "\t\t2: Absolute Activity Last Time\n"
                   "\t\t3: Absolute Activity This Time - Externals\n"
                   "\t\t4: # of Conns between CoActive Neurons\n"
                   "\t\t5: 4/3 = Ave # of Active Conns for Active Neurons\n"
                   "\t\t6: 4/(2*3) = Prob cij=1 given zi=1,zj=1\n"
                   "\tActiveConnect also updates the variable Pczz which\n"
                   "\tholds the Prob(cij=1 | zi=1,zj=1) for the whole sequence\n");
      ComL.IntSet(4, &StartNeuron, &EndNeuron, &StartPat, &EndPat);
      ComL.StrSet(3, &SeqName, &ResetSeq, &External);
      ComL.FlagSet(1, &DoSummary);
      argunset = false;
   }
   ComL.Process(arg, Output::Err());

   if (!program::Main().getNetworkCreated()) {
      CALL_ERROR << "Error: You must call @CreateNetwork() before you can "
                    "call @ActiveConnect." << ERR_WHERE;
      exit(EXIT_FAILURE);
   }
   UIPtnSequence Seq = SystemVar::getSequence(SeqName, FunctionName, ComL);
   int SeqSize = Seq.size();
   UIPtnSequence ReSeq;
   if (ResetSeq.getValue() != "{zeros}") {
      ReSeq = SystemVar::getSequence(ResetSeq, FunctionName, ComL);
   } else {
      ReSeq.resize(SeqSize, UIVector(0));
   }
   int lastNeuron = 0;
   for (UIPtnSequenceCIt it = Seq.begin(); it != Seq.end(); it++) {
      updateMax(lastNeuron, it->back());
   }
   if (EndNeuron.getValue() == -1) {
      EndNeuron.setValue(lastNeuron);
   }
   if (EndPat.getValue() == -1) {
      EndPat.setValue(SeqSize);
   }
   if (StartNeuron.getValue() < 1 ||
         StartNeuron.getValue() > EndNeuron.getValue() ||
         EndNeuron.getValue() > static_cast<int>(ni)) {
      CALL_ERROR << "Error in " << FunctionName << ": Neuron range of "
          << StartNeuron.getValue() << "..." << EndNeuron.getValue()
          << " is invalid." << ERR_WHERE;
      ComL.DisplayHelp(Output::Err());
   }
   if (EndPat.getValue() < 1 || EndPat.getValue() > SeqSize ||
          StartPat.getValue() < 1 || StartPat.getValue() > EndPat.getValue()) {
      CALL_ERROR << "Error in " << FunctionName << ": Pattern range of "
         << StartPat.getValue() << "..." << EndPat.getValue() << " is invalid."
         << ERR_WHERE;
      ComL.DisplayHelp(Output::Err());
   }

   // get externals
   UIMatrix ExtSeq;
   if (External.getValue() != "{zeros}") {
      ExtSeq = SystemVar::getSequence(External, FunctionName, ComL);
   } else {
      ExtSeq.resize(SeqSize, UIVector(0));
   }

   // Do the loops
   int NumPats = EndPat.getValue() - StartPat.getValue() + 1;
   DataList Ave(NumPats);
   DataList Sum(NumPats);
   DataList Pats(NumPats);
   DataList jNumFired(NumPats);
   DataList iNumFired(NumPats);
   DataList PczzV(NumPats);
   int index = 0;
   UIPtnSequenceCIt SeqConstIterator = Seq.begin();
   UIPtnSequenceCIt ExtConstIterator = ExtSeq.begin();   
   UIVector pat_b4;
   if (StartPat.getValue() == 1) {
      pat_b4 = ReSeq.front();
   } else {
      for (int inc = 0; inc < StartPat.getValue() - 2; ++inc) ++SeqConstIterator;
//.........这里部分代码省略.........

void RateMeyerDiscrete::computeFuncDerv(double value, double &df, double &ddf) {
	if (!is_categorized) {
		RateMeyerHaeseler::computeFuncDerv(value, df, ddf);
		return;
	}

//	double lh = 0.0;
	int nseq = phylo_tree->leafNum;
	int nstate = phylo_tree->getModel()->num_states;
	int i, j, k, state1, state2;
	ModelSubst *model = phylo_tree->getModel();
    int trans_size = nstate * nstate;
	double *trans_mat = new double[trans_size];
	double *trans_derv1 = new double[trans_size];
	double *trans_derv2 = new double[trans_size];
	df = ddf = 0.0;

	int *pair_freq = new int[trans_size];

	for (i = 0; i < nseq-1; i++) 
		for (j = i+1; j < nseq; j++) {
			memset(pair_freq, 0, trans_size * sizeof(int));
			for (k = 0; k < size(); k++) {
				if (ptn_cat[k] != optimizing_cat) continue;
				Pattern *pat = & phylo_tree->aln->at(k);
				if ((state1 = pat->at(i)) < nstate && (state2 = pat->at(j)) < nstate)
					pair_freq[state1*nstate + state2] += pat->frequency;
			}
			double dist = dist_mat[i*nseq + j];
			double derv1 = 0.0, derv2 = 0.0;
			model->computeTransDerv(value * dist, trans_mat, trans_derv1, trans_derv2);
			for (k = 0; k < trans_size; k++) if (pair_freq[k]) {
				double t1 = trans_derv1[k] / trans_mat[k];
				double t2 = trans_derv2[k] / trans_mat[k];
				trans_derv1[k] = t1;
				trans_derv2[k] = (t2 - t1*t1);
//				lh -= log(trans_mat[k]) * pair_freq[k];
				derv1 += trans_derv1[k] * pair_freq[k];
				derv2 += trans_derv2[k] * pair_freq[k];
			}
			df -= derv1 * dist;
			ddf -= derv2 * dist * dist;
		}
	delete [] pair_freq;
	delete [] trans_derv2;
	delete [] trans_derv1;
	delete [] trans_mat;
//	return lh;

/*	double lh = 0.0, derv1, derv2;
	df = 0.0; ddf = 0.0;	
	for (int i = 0; i < size(); i++)
		if (ptn_cat[i] == optimizing_cat) {
			optimizing_pattern = i;
			int freq =  phylo_tree->aln->at(i).frequency;
			lh += RateMeyerHaeseler::computeFuncDerv(value, derv1, derv2) * freq;
			df += derv1 * freq;
			ddf += derv2 * freq;
		}
	return lh;*/
}

void Preprocessor::correctSlanting ()
{
	boost::timer timer;
	timer.restart();

	const double PI = 3.1415926535;

	for ( std::vector<Pattern>::iterator i = patterns_.begin(); i != patterns_.end(); ++i )
	{
	    const unsigned int rotationLimit = 20;
		double shearingFactor;;

	    // This vector stores the greatest number of ink pixels counted in any column of a rotated pattern.
		std::vector<unsigned int> columnPixelMaximalCounts(rotationLimit);

		// Rotate the pattern and count the number of ink pixels per column
		for ( double angle = 0; angle < rotationLimit; ++angle )
		{
			std::map<unsigned int, std::set<unsigned int> > pixelsPerColumn;
			shearingFactor = tan(-angle * PI / 180.0);

			// Store the columns of rotated pixels
			for ( unsigned int j = 0; j < i->height(); ++j )
			{
				for ( unsigned int k = 0; k < i->width(); ++k )
				{
					if ( i->at(j,k) == 1 )  // Ink pixel
					{
						double jp = j;
						double kp = k;
						kp = round(kp - jp * shearingFactor);

						if ( kp < i->width() && kp > 0 )
							pixelsPerColumn[kp].insert(j);
					}
				}
			}

			unsigned int columnPixelMaximalCount = 0;
			for ( std::map<unsigned int, std::set<unsigned int> >::iterator row = pixelsPerColumn.begin(); row != pixelsPerColumn.end(); ++row )
			{
				if ( row->second.size() > columnPixelMaximalCount )
					columnPixelMaximalCount = row->second.size();
			}

			columnPixelMaximalCounts.at(angle) = columnPixelMaximalCount;
		}

		// Get the angle to rotate selecting the rotation that generates the maximal number of pixels in a column
		double targetAngle = distance( columnPixelMaximalCounts.begin(), std::max_element(columnPixelMaximalCounts.begin(), columnPixelMaximalCounts.end()) );
		shearingFactor = tan(-targetAngle * PI / 180.0);

		// Correct slanting of the actual pattern
		if ( targetAngle != 0 )
		{
			Pattern rotatedPattern;
			rotatedPattern.clean();

			for ( unsigned int j = 0; j < i->height(); ++j )
			{
				for ( unsigned int k = 0; k < i->width(); ++k )
				{
					if ( i->at(j,k) == 1 )
					{
						double jp = j;
						double kp = k;
						kp = round(kp - jp * shearingFactor);

						if ( kp < i->width() && kp > 0 )
							rotatedPattern.at(jp, kp) = 1;
					}
				}
			}
			*i = rotatedPattern;
		}
	}

	statistics_.slantingCorrectionTime(timer.elapsed());;
}

void Preprocessor::buildPatterns ()
{
	boost::timer timer;
	timer.restart();

	patterns_.clear();
	patterns_.reserve(regions_.size());
	Pattern pattern;

	// Traverse the list of iterators to regions creating a Pattern object for each region
	for( RegionLines::iterator k = inlineRegions_.begin(); k != inlineRegions_.end(); ++k )
	{
		std::list<RegionIterator> line(k->second);
		for( std::list<RegionIterator>::iterator i = line.begin(); i != line.end(); ++i )
		{
			(*i)->normalizeCoordinates();

			// Normalize region using Magick++ facilities
			Magick::Image image( Magick::Geometry((*i)->width(), (*i)->height()), Magick::ColorGray(1.0) );
			image.type( Magick::BilevelType );

			Magick::Pixels view(image);
			Magick::PixelPacket *originPixel = view.get(0, 0, (*i)->width(), (*i)->height());
			Magick::PixelPacket *pixel;

			for ( unsigned int j = 0; j < (*i)->size(); ++j )
			{
				pixel	= originPixel + ((*i)->at(j).first * view.columns() + (*i)->at(j).second);
				*pixel	= Magick::ColorGray (0.0);
			}
			
			view.sync();
			image.syncPixels();
			image.scale( Magick::Geometry(Pattern::planeSize(), Pattern::planeSize()) );

			// Preprocess the normalized region
			Preprocessor temporalPreprocessor (image, 0, 0, image.rows(), image.columns());
			temporalPreprocessor.applyGlobalThresholding();
			temporalPreprocessor.isolateRegions();

			Region normalizedRegion;
			if ( ! temporalPreprocessor.regions_.empty() )
			{
				// Merge subregions if preprocessing split the original region
				if ( temporalPreprocessor.regions_.size() > 1 )
				{
					for ( RegionIterator j = temporalPreprocessor.regions_.begin(); j != temporalPreprocessor.regions_.end(); ++j )
						normalizedRegion = normalizedRegion + *j;

					temporalPreprocessor.regions_.clear();
					temporalPreprocessor.regions_.push_back(normalizedRegion);
				}
				else
					normalizedRegion = temporalPreprocessor.regions_.front();
			}
			
			// Build the pattern
			pattern.clean();
			for ( unsigned int i = 0; i < normalizedRegion.size(); ++i )
				pattern.at(normalizedRegion.at(i).first, normalizedRegion.at(i).second) = 1;

			// Correct shifting
			if ( image.rows() < Pattern::planeSize() )	// Shift rows from top to the center
			{
				unsigned int offset = (Pattern::planeSize() - image.rows()) / 2;

				while ( offset != 0 )
				{
					for ( int i = Pattern::planeSize()-2; i >= 0; --i )
					{
						for ( unsigned int j = 0; j < Pattern::planeSize(); ++j )
							pattern.at(i+1, j) = pattern.at(i, j);
					}
					for ( unsigned int j = 0; j < Pattern::planeSize(); ++j )
						pattern.at(0, j) = 0;

					--offset;
				}
			}

			if ( image.columns() < Pattern::planeSize() )	// Shift columns from left to center
			{
				unsigned int offset = (Pattern::planeSize() - image.columns()) / 2;

				while ( offset != 0 )
				{
					for ( unsigned int i = 0; i < Pattern::planeSize(); ++i )
					{
						for ( int j = Pattern::planeSize()-2; j >= 0; --j )
							pattern.at(i, j+1) = pattern.at(i, j);
					}
					for ( unsigned int i = 0; i < Pattern::planeSize(); ++i )
						pattern.at(i, 0) = 0;

					--offset;
				}
			}

			patterns_.push_back( pattern );
		}
//.........这里部分代码省略.........