C++ TimeSeries::push_back方法代码示例

// Iterate from an initial state
void DynSysModel::Iterate( NTuple& initial_state, size_t length, TimeSeries& result )
{
	// Make sure the result is empty
	result.clear();

	// Save the initial state as the first time series entry
	result.push_back( initial_state );

	mState = initial_state;
	NTuple next_state;
	PolyModelIter iter;

	size_t i,k;
	for( i = length-1; i > 0; --i )
	{
		k = 1;
		// For each polynomial;
		next_state.Reset();
		iter = mModel.begin();
		while( iter != mModel.end() )
		{
			// Evaluate the k'th polynomial at the current state
			//   to produce a new value for the k'th variable
			next_state.Assign( k, mModel[k-1].Evaluate( mState ) );

			// next polynomial
			++iter; ++k;
		}
		// Update the current state to be the newly compute state
		mState = next_state;
		result.push_back( next_state );
	}
}

Datum blocker(const RealMatrix& Ua, const RealMatrix sa, const vGroup& ensemble, const uint blocksize, uint repeats, ExtractPolicy& policy) {


  
  TimeSeries<double> coverlaps;

  for (uint i=0; i<repeats; ++i) {
    vector<uint> picks = pickFrames(ensemble.size(), blocksize);
    
    if (debug) {
      cerr << "***Block " << blocksize << ", replica " << i << ", picks " << picks.size() << endl;
      dumpPicks(picks);
    }
    
    vGroup subset = subgroup(ensemble, picks);
    boost::tuple<RealMatrix, RealMatrix> pca_result = pca(subset, policy);
    RealMatrix s = boost::get<0>(pca_result);
    RealMatrix U = boost::get<1>(pca_result);

    if (length_normalize)
      for (uint j=0; j<s.rows(); ++j)
        s[j] /= blocksize;

    coverlaps.push_back(covarianceOverlap(sa, Ua, s, U));
  }

  return( Datum(coverlaps.average(), coverlaps.variance(), coverlaps.size()) );

}

// Iterate from an initial state with the k'th function knocked out
void DynSysModel::KoIterate( NTuple& initial_state, size_t length, TimeSeries& result, size_t kov )
{
	// Make sure the result is empty
	result.clear();

	// Force k'th entry to zero in the initial state -  no longer needed, corrected when file is read in
	// NTuple state1 = initial_state;
	// state1.Reset( kov );

	// Save the initial state as the first time series entry
	result.push_back( initial_state );

	mState = initial_state;
	NTuple next_state;
	PolyModelIter iter;

	size_t i,k;
	for( i = length-1; i > 0; --i )
	{
		k = 1;
		// For each polynomial
		next_state.Reset();
		iter = mModel.begin();
		while( iter != mModel.end() )
		{
			// Force the knockout function result to zero
			if( k == kov )
			{
				next_state.Assign( k, 0 );
			}
			else
			{
				// Evaluate the k'th polynomial at the current state
				//   to produce a new value for the k'th variable
				next_state.Assign( k, mModel[k-1].Evaluate( mState ) );
			}

			// next polynomial
			++iter; ++k;
		}
		// Update the current state to be the newly compute state
		mState = next_state;
		result.push_back( next_state );
	}
}

void BFFindMotif::FindMotifSub(std::deque<Point> &window)
{
	size_t i = 0;
	for(i = 0; i < window.size() - m_MotifLength; ++i)
	{
		double distance = 0.0;
		bool newMotif = true;
		TimeSeries ts;
		ts.reserve(m_MotifLength);
		if(m_SlideWindow.size() >= m_MotifLength)		//	Only process slide window larger than motif length
		{
			//	Get time series
			for(size_t j = i; j < i + m_MotifLength; ++j)
			{
				ts.push_back(window[j].second);
			}

			//	Compare with candidate motif
			for(size_t j = 0; j < m_CandidateMotif.size(); ++j)
			{
				distance = EuclideanDistance(m_CandidateMotif[j].second, ts);

				if((2 * m_Radius > distance) && (m_Radius < distance))	//	Neither new motif nor similar motif
				{
					newMotif = false;
				}
				else if(m_Radius > distance)	//	Similar motif
				{
					m_CandidateMotif[j].first++;
					i += m_Step;
					newMotif = false;
					break;		//	Impossible to be similar with other candidates
				}
			}

			if(true == newMotif)	//	New motif
			{
				m_CandidateMotif.push_back(make_pair<long long, TimeSeries>(1, ts));
				i += m_Step;
			}
			
		}
		else
		{
			cerr << "Window size:" << m_SlideWindow.size() << endl;
		}
	
	}
	
	for(size_t k = 0; k < i; ++k)
	{
		window.pop_front();
	}
	
}

void SIMMA:: FindMotifSub(long &bufferCount)
{
	long lLastSize = m_K;			//	last size of candidate motif buffer
	long step = 32;

	size_t i = 0;
	for(i = 0; i < m_SlideWindow.size() - m_MotifLength; ++i)
	{
		double distance = 0.0;
		bool newMotif = true;
		TimeSeries ts;
		ts.reserve(m_MotifLength);

		++bufferCount;

		//	Check buffer when size big enough
		if(true == m_bBufferCheck)
		{
			if((2 * m_K < bufferCount) && (5 * m_K < (int)m_CandidateMotif.size()))
			{
				this->BufferCheck();
				lLastSize = m_SlideWindow.size();
				bufferCount = 0;
			}
		}

		if(m_SlideWindow.size() >= m_MotifLength)	//	Only process slide window larger than motif length
		{
			//	Get time series
			for(size_t j = i; j < i + m_MotifLength; ++j)
			{
				ts.push_back(m_SlideWindow[j].second);
			}

			//	Compare with candidate motif
			for(size_t j = 0; j < m_CandidateMotif.size(); ++j)
			{
				distance = EuclideanDistance(m_CandidateMotif[j].second, ts);

				if((2 * m_Radius > distance) && (m_Radius < distance))	//	Neither new motif nor similar motif
				{
					newMotif = false;
				}
				else if(m_Radius > distance)	//	Similar motif
				{
					m_CandidateMotif[j].first++;
					long jump = this->NonTrivialStep(i);

					if(0 != jump)
						i += jump;
					else
						i = m_SlideWindow.size() - m_MotifLength;

					newMotif = false;
					break;		//	Impossible to be similar with other candidates
				}
			}
			
			//	Check whether current time series is new motif
			if(true == newMotif)
			{
				m_CandidateMotif.push_back(make_pair<long long, TimeSeries>(1, ts));
				long jump = this->NonTrivialStep(i);

				if(0 != jump)
					i += jump;
				else
					i = m_SlideWindow.size() - m_MotifLength;
			}
		}
	}
	
	//	Remove used elements
	for(size_t k = 0; k < i; ++k)
	{
		m_SlideWindow.pop_front();
	}
}

int main(int argc, char *argv[]) {
  
  if (argc < 4 || argc > 6) {
    cerr << "Usage- block_avgconv model traj sel [range [1 = do not align trajectory]]\n";
    cerr << fullHelpMessage();
    exit(-1);
  }

  string hdr = invocationHeader(argc, argv);

  int k = 1;
  AtomicGroup model = createSystem(argv[k++]);
  pTraj traj = createTrajectory(argv[k++], model);
  AtomicGroup subset = selectAtoms(model, argv[k++]);
  
  vector<uint> sizes;
  bool do_align = true;

  if (argc == k) {
    uint step = traj->nframes() / default_starting_number_of_blocks;
    for (uint i=step; i<traj->nframes() * default_fraction_of_trajectory; i += step)
      sizes.push_back(i);
  } else {
    sizes = parseRangeList<uint>(argv[k++]);
    if (argc == k+1)
      do_align = (argv[k][0] != '1');
  }


  cout << "# " << hdr << endl;
  cout << "# n\tavg\tvar\tblocks\tstderr\n";

  vector<AtomicGroup> ensemble;
  cerr << "Reading trajectory...\n";
  readTrajectory(ensemble, subset, traj);

  if (do_align) {
    cerr << "Aligning trajectory...\n";
    boost::tuple<vector<XForm>, greal, int> result = iterativeAlignment(ensemble);
  } else
    cerr << "Trajectory is already aligned!\n";

  cerr << "Processing- ";
  for (uint block = 0; block < sizes.size(); ++block) {
    if (block % 50)
      cerr << ".";

    uint blocksize = sizes[block];

    vector<AtomicGroup> averages;
    for (uint i=0; i<ensemble.size() - blocksize; i += blocksize) {

      vector<uint> indices(blocksize);
      for (uint j=0; j<blocksize; ++j)
        indices[j] = i+j;
      
      averages.push_back(averageSelectedSubset(ensemble, indices));
    }
    
    TimeSeries<double> rmsds;
    for (uint j=0; j<averages.size() - 1; ++j)
      for (uint i=j+1; i<averages.size(); ++i) {
        AtomicGroup left = averages[j];
        AtomicGroup right = averages[i];
        left.alignOnto(right);
        rmsds.push_back(left.rmsd(right));
      }
        

    double v = rmsds.variance();
    uint n = averages.size();
    cout << boost::format("%d\t%f\t%f\t%d\t%f\n") % blocksize % rmsds.average() % v % n % sqrt(v/n);
  }
  
  cerr << "\nDone!\n";
}