parseDotBracket.cpp

#include <string>
#include <cassert>
#include <vector>
#include <map>
#include <iostream>

//unsigned long basePairCount; will be computed
//unsigned long maxDepth; will be computed
bool isDotBracket(const std::string &string, unsigned long &basePairCount,
            unsigned long &maxDepth) {
  unsigned long leftBrackets = 0;

  std::string::const_iterator it;
  for (it = string.begin(); it != string.end(); it++) {
    switch (*it) {
    case '.':
      break;
    case '(':
      leftBrackets++;
      basePairCount++;
      if(leftBrackets>0)
        maxDepth = std::max(maxDepth,leftBrackets);
      break;
    case ')':
      if (leftBrackets) {
        leftBrackets--;
        break;
      }
    default:
      return false;
    }
  }
  // if no brackets remain open, string is valid dotBracket
  return !leftBrackets;
}

// basepair map is built traversing the dotBracket string, with aid of a stack of 
// leftBrackets

void buildBPIndex(const std::string &base, const std::string &dotBracket) {
    unsigned long basePairCount = 0, maxDepth = 0, stackPtr = 0;
    unsigned int structureLength, node;

    // base string can contain anything!
    isDotBracket(dotBracket, basePairCount, maxDepth);
    assert(base.length() == dotBracket.length());
    int length = dotBracket.length();

    std::vector<unsigned int> bracketStack(maxDepth + 1);

    int * baseIndex = new int[length];
    stackPtr = 0;

    // read dotBracket string char by char
    for (int i = 0; i < length; i++) {
        switch (dotBracket[i]) {
            case '.':
                // dot pairs with itself
                //std::cout << "dot " << i << std::endl;
                baseIndex[i] = i;
                break;

            case '(':
                //descend

                // We will see the pairing partner later, for now put -1
                baseIndex[i] = -1;

                // descend - push node stack
                stackPtr++;
                bracketStack[stackPtr] = i;
                break;

            case ')':
                // ascend
                // we know the pairing partner now
                int leftBracket = bracketStack[stackPtr];
                //std::cout << "we remember the partner, it was " << i << " and " << leftBracket << std::endl;
                if (baseIndex[leftBracket] != -1) {
                  std::cout << "Something went wrong, unbalanced structure?" << std::endl;
                  exit(0);
                }
                baseIndex[i] = leftBracket;
                baseIndex[leftBracket] = i;

                // ascend
                stackPtr--;
                break;
        }
    }
    for (int i = 0; i < length; i++) {
      std::cout << baseIndex[i] << " ";
    }
    std::cout << std::endl;

}

int main(int argc, char** argv) {
  //std::cout << "Calling buildBPIndex(" << argv[1] << "," << argv[2] << ");" << std::endl;
  buildBPIndex(argv[1],argv[2]);
}