C++ Sequence::getSubsequence方法代码示例

  // jointProb() precodition: neither nmer is degenerate!
  double jointProb(Sequence &parentNmer,Sequence &childNmer,
		   NthOrdSubstMatrix &noPt) {
    double logP=0;
    Sequence parentContext, childContext;
    for(int i=0 ; i<numCols ; ++i) {
      Symbol parentSymbol=parentNmer[i], childSymbol=childNmer[i];
      parentNmer.getSubsequence(0,i,parentContext);
      childNmer.getSubsequence(0,i,childContext);
      SubstitutionMatrix &Pt=*noPt.lookup(parentContext);
      logP+=Pt(parentSymbol,childSymbol);
    }
    return logP;
  }

  RCO_Likelihooder(int numNodes,int numAlpha,Alphabet &alpha,
		   MultSeqAlignment &A,int column,
		   RootNode *root,int order,AlphabetMap &alphabetMap,
		   const BitSet &gapSymbols)
    : A(A), numNodes(numNodes), numAlpha(numAlpha),
      alpha(alpha), L(numNodes,numAlpha), column(column),
      alphabetMap(alphabetMap), gapSymbols(gapSymbols)
  {
    // ctor
    
    // Get the root's left-context for this column (for RCO)
    int rootID=root->getID();
    AlignmentSeq &rootTrack=A.getIthTrack(rootID);
    int contextBegin=column-order;
    if(contextBegin<0) contextBegin=0;
    int contextLen=column-contextBegin;
    rootTrack.getSeq().getSubsequence(contextBegin,contextLen,
				      rootContext);
    int pos=MultSeqAlignment::rightmostGapPos(rootContext,gapSymbols);
    if(pos>=0) {
      Sequence temp;
      rootContext.getSubsequence(pos+1,contextLen-pos-1,temp);
      rootContext=temp;
    }
    rootContextCode=rootContext.asInt(alphabetMap);
    rootContextLength=contextLen;
  }

  double processInternalChild(Symbol parentSymbol,int parentID,int child,
			      NthOrdSubstMatrix &noPt) {
    AlignmentSeq &track=A.getIthTrack(parentID);
    AlignmentSeq &childTrack=A.getIthTrack(child);
    int contextBegin=column-order;
    if(contextBegin<0) contextBegin=0;
    int contextLen=column-contextBegin;
    Sequence context;
    track.getSeq().getSubsequence(contextBegin,contextLen,context);//ACO
    //A.getIthTrack(root->getID()).getSeq().getSubsequence(contextBegin,contextLen,context);//TRCO
    int pos=MultSeqAlignment::rightmostGapPos(context,gapSymbols);
    if(pos>=0) {
      throw "this should not happen (FitchFelsenstein)::processInternalChild";
      Sequence temp;
      context.getSubsequence(pos+1,contextLen-pos-1,temp);
      context=temp;
    }
    SubstitutionMatrix &Pt=*noPt.lookup(context,0,-1);
    Array2D<double>::RowIn2DArray<double> row=L[child];
    Symbol childSym=childTrack.getSeq()[column];
    double ll;
    if(gapSymbols.isMember(childSym)){
      V.setAllTo(NEGATIVE_INFINITY);
      for(Symbol a=0 ; a<numAlpha ; ++a) {
	if(!gapSymbols.isMember(a)) {
	  V[(int)a]=safeAdd(row[a],Pt(parentSymbol,a));
	}
      }
      ll=sumLogProbs(V);
    }
    else {
      ll=safeAdd(row[childSym],Pt(parentSymbol,childSym));
    }
    return ll;
  }

  void getContext(PhylogenyNode *node,int &contextCode,int &contextLength)
  {
    int ID=node->getID();
    AlignmentSeq &track=A.getIthTrack(ID);
    int contextBegin=column-order;
    if(contextBegin<0) contextBegin=0;
    contextLength=column-contextBegin;
    if(contextLength==0) {contextCode=0;return;}
    Sequence context;
    track.getSeq().getSubsequence(contextBegin,contextLength,context);
    int pos=MultSeqAlignment::rightmostGapPos(context,gapSymbols);
    if(pos>=0) {
      Sequence temp;
      int begin=pos+1, len=contextLength-pos-1;
      if(begin<order && len>0)
	context.getSubsequence(pos+1,contextLength-pos-1,temp);
      context=temp;
      contextLength=context.getLength();
    }
    contextCode=context.asInt(*alphabetMap.getRange());
  }

void GenomeFeature::initSeq(const Sequence &from)
{
  from.getSubsequence(begin,end-begin,seq);
}