From 11b308e7363e937966b035b4891db32b4eece3bf Mon Sep 17 00:00:00 2001 From: Tor Aamodt Date: Fri, 1 Oct 2010 08:55:28 -0800 Subject: integrating recent changes from fermi-test into fermi (i'll use "fermi" for more disruptive changes to the pipeline model such as updating the MSHRs and getting rid of the warp tracker, ripping out DWF, etc...) [git-p4: depot-paths = "//depot/gpgpu_sim_research/fermi/distribution/": change = 7805] --- benchmarks/CUDA/MUM/suffix-tree.cpp | 1784 ----------------------------------- 1 file changed, 1784 deletions(-) delete mode 100644 benchmarks/CUDA/MUM/suffix-tree.cpp (limited to 'benchmarks/CUDA/MUM/suffix-tree.cpp') diff --git a/benchmarks/CUDA/MUM/suffix-tree.cpp b/benchmarks/CUDA/MUM/suffix-tree.cpp deleted file mode 100644 index 872fcc3..0000000 --- a/benchmarks/CUDA/MUM/suffix-tree.cpp +++ /dev/null @@ -1,1784 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include - -#define ulong4 uint32_t -#include "mummergpu.h" - -#define MPOOL 1 - -#include "PoolMalloc.hh" - - -#include - -using namespace std; - -// Enable verification/debug options -#define VERIFY 0 -#define VERBOSE 0 -const bool DEBUG = 0; - -// Setting for linear time alg -bool FORCEROOT = false; -bool DOJUMP = true; -bool DOINTERNALSKIP = true; -bool DOPHASETRICK = true; - -// Statistics -int skippedbases = 0; -int skippedextensions = 0; - -char substrbuffer[1024]; -const char * substr(const char * str, int start, int len) -{ - if (len > 1024) { len = 1024; } - strncpy(substrbuffer, str+start, len); - substrbuffer[len] = '\0'; - - return substrbuffer; -} - - -// Helper to convert from ascii to single byte -unsigned char b2i(char base) -{ - switch (base) - { - case 'A' : return 0; - case 'C' : return 1; - case 'G' : return 2; - case 'T' : return 3; - case '$' : return 4; - - default: - cerr << "Unknown base: " << base << endl; - return b2i('A'); - }; -} - -#include -#include - -class EventTime_t -{ -public: - /// Constructor, starts the stopwatch - EventTime_t() - { - start(); - memset(&m_end, 0, sizeof(struct timeval)); - } - - - /// Explicitly restart the stopwatch - void start() - { - gettimeofday(&m_start, NULL); - } - - - /// Explicitly stop the stopwatch - void stop() - { - gettimeofday(&m_end, NULL); - } - - - /// Return the duration in seconds - double duration() - { - if ((m_end.tv_sec == 0) && (m_end.tv_usec == 0)) { stop(); } - return ((m_end.tv_sec - m_start.tv_sec)*1000000.0 + (m_end.tv_usec - m_start.tv_usec)) / 1e6; - } - - - /** \brief Pretty-print the duration in seconds. - ** If stop() has not already been called, uses the current time as the end - ** time. - ** \param format Controls if time should be enclosed in [ ] - ** \param precision Controls number of digits past decimal pt - **/ - std::string str(bool format = true, - int precision=2) - { - double r = duration(); - - char buffer[1024]; - sprintf(buffer, "%0.*f", precision, r); - - if (format) - { - string s("["); - s += buffer; - s += "s]"; - return s; - } - - return buffer; - } - - -private: - /// Start time - struct timeval m_start; - - /// End time - struct timeval m_end; -}; - - -// A node in the suffix tree -class SuffixNode -{ -public: - static int s_nodecount; - -#ifdef MPOOL - void *operator new( size_t num_bytes, PoolMalloc_t *mem) - { - return mem->pmalloc(num_bytes); - } -#endif - - SuffixNode(int s, int e, int leafid, - SuffixNode * p, SuffixNode * x) - : m_start(s), m_end(e), - m_nodeid(++s_nodecount), - m_leafid(leafid), - m_parent(p), - m_suffix(x) - { - for (int i = 0; i < basecount; i++) - { m_children[i] = NULL; } - - m_depth = len(); - if (p) - m_depth += p->m_depth; - } - - ~SuffixNode() - { - for (int i = 0; i < basecount; i++) - { - if (m_children[i]) { delete m_children[i]; } - } - } - - int id() - { - if (this) { return m_nodeid; } - return 0; - } - - bool isLeaf() - { - for (int i = 0; i < basecount; i++) - { - if (m_children[i]) { return false; } - } - - return true; - } - - const char * str(const char * refstr) - { - return substr(refstr, m_start, m_end-m_start+1); - } - - int len(int i=-1) - { - if (i != -1) - { - if (i < m_end) - { - return i - m_start + 1; - } - } - - return m_end - m_start + 1; - } - - int depth() - { - return m_depth; - } - - ostream & printLabel(ostream & os, const char * refstr) - { - if (m_start == m_end && m_start == 0) - { - os << "\"ROOT\""; - } - else - { - os << "\"" << str(refstr) << "\""; - - // << " [" << m_start - // << "," << m_end - // << "(" << m_nodeid << ")\""; - } - - return os; - } - - - ostream & printNodeLabel(ostream & os) - { - os << m_nodeid; - return os; - } - - ostream & printEdgeLabel(ostream & os, const char * refstr) - { - string seq = substr(refstr, m_start, m_end-m_start+1); - os << "\"" << seq << "\""; - //os << "\"" << seq << " [" << m_start << "," << m_end << "]\""; - return os; - } - - int m_start; // start pos in string - int m_end; // end pos in string - int m_nodeid; // the id for this node - int m_leafid; // For leafs, the start position of the suffix in the string - SuffixNode * m_children [basecount]; // children nodes - SuffixNode * m_parent; // parent node - SuffixNode * m_suffix; // suffixlink - int m_depth; -#if VERIFY - string m_pathstring; // string of path to node -#endif -}; - -int SuffixNode::s_nodecount(0); - -ostream & operator<< (ostream & os, SuffixNode * n) -{ - return n->printNodeLabel(os); -} - - -// Encapsulate the tree with some helper functions -class SuffixTree -{ -public: - SuffixTree(const char * s) : m_string(s) - { - m_strlen = strlen(s); -#ifdef MPOOL - m_root = new (&m_pool) SuffixNode(0,0,0,NULL,NULL); // whole tree -#else - m_root = new SuffixNode(0,0,0,NULL,NULL); // whole tree -#endif - m_root->m_suffix = m_root; - } - - ~SuffixTree() - { -#ifdef MPOOL -#else - delete m_root; -#endif - } - - SuffixNode * m_root; - const char * m_string; - int m_strlen; - -#ifdef MPOOL - PoolMalloc_t m_pool; -#endif - - // Print a node for dot - void printNodeDot(SuffixNode * node, ostream & dfile) - { - int children = 0; - for (int i = 0; i < basecount; i++) - { - SuffixNode * child = node->m_children[i]; - if (child) - { - children++; - - dfile << " " << node << "->" << child; - - //node->printNodeLabel(dfile, m_string) << " -> "; - //child->printNodeLabel(dfile, m_string); - - //dfile << " [minlen=" << child->len() << ", label="; - dfile << " [minlen=1, label="; - child->printEdgeLabel(dfile, m_string) << "]" << endl; - - printNodeDot(child, dfile); - } - } - - if (node->m_suffix) - { - dfile << " " << node << " -> " << node->m_suffix - << " [style=dotted, constraint=false]" << endl; - - //node->printLabel(dfile, m_string) << " -> "; - //node->m_suffix->printLabel(dfile, m_string) << " [style=dotted, constraint=false]" << endl; - } - - if (children == 0) - { - //dfile << " " << node << " [shape=box, label="; - //node->printLabel(dfile, m_string) << "]" << endl; - - dfile << " " << node << " [shape=box,width=.2,height=.2,label=\"" << node->id() << ":" << node->m_leafid << "\"]" << endl; - } - else - { - //dfile << " " << node << " [label="; - //node->printLabel(dfile, m_string) << "]" << endl; - dfile << " " << node << " [width=.2,height=.2,label=\"" << node->id() << "\"]" << endl; - } - } - - // Print the whole tree for dot - void printDot(const char * dotfilename) - { - ofstream dfile; - dfile.open(dotfilename, ofstream::out | ofstream::trunc); - - cerr << "Printing dot tree to " << dotfilename << endl; - - dfile << "digraph G {" << endl; - dfile << " size=\"7.5,10\"" << endl; - dfile << " center=true" << endl; - dfile << " label=\"Suffix tree of \'" << m_string << "\' len:" - << m_strlen-1 << " nc:" - << SuffixNode::s_nodecount << "\"" << endl; - - printNodeDot(m_root, dfile); - dfile << "}" << endl; - } - - // Print a node in text format - void printNodeText(ostream & out, SuffixNode * n, int depth) - { - for (int b = 0; b < basecount; b++) - { - if (n->m_children[b]) - { - for (int i = 0; i < depth; i++) - { - out << " "; - } - out << " "; - out << n->m_children[b]->str(m_string) << endl; - printNodeText(out, n->m_children[b], depth+1); - } - } - } - - // Print the tree in Text - void printText(ostream & out) - { - out << "Suffix Tree len=" << m_strlen-1 << endl; - out << "String: \"" << m_string << "\"" << endl; - out << "+" << endl; - printNodeText(out, m_root, 0); - } - - // Print the tree as list of sorted suffixes - void printTreeSorted(ostream & out, SuffixNode * node, const string & pathstring) - { - bool isLeaf = true; - - string ps(pathstring); - if (node != m_root) { ps.append(node->str(m_string)); } - - for (int i = 0; i < basecount; i++) - { - if (node->m_children[i]) - { - isLeaf = false; - printTreeSorted(out, node->m_children[i], ps); - } - } - - if (isLeaf) { out << ps << endl; } - } - - void printTreeFlat(ostream & out) - { - cerr << "nodeid\tparent\tSL\tstart\tend\t$\tA\tC\tG\tT\tnodestring" << endl; - cout << "0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0" << endl; - printNodeFlat(out, m_root); - } - - void printNodeFlat(ostream & out, SuffixNode * node) - { - out << node->id() << "\t" - << node->m_parent->id() << "\t" - << node->m_suffix->id() << "\t" - << node->m_start << "\t" - << node->m_end << "\t"; - - for (int i = 0; i < basecount; i++) - { - out << node->m_children[i]->id() << "\t"; - } - - out << node->m_start << "\t" << node->m_end << endl; - - if (node == m_root) { out << "ROOT" << endl; } - else { out << node->str(m_string) << endl; } - - for (int i = 0; i < basecount; i++) - { - if (node->m_children[i]) { printNodeFlat(out, node->m_children[i]); } - } - } - -#if VERIFY - void setNodePath(SuffixNode * node, const string & parentString) - { - node->m_pathstring = parentString; - - if (node != m_root) - { - node->m_pathstring.append(m_string, node->m_start, node->m_end - node->m_start + 1); - } - - for (int b = 0; b < basecount; b++) - { - if (node->m_children[b]) - { - setNodePath(node->m_children[b], node->m_pathstring); - } - } - } - - int verifyNodeSuffixLinks(SuffixNode * node, int & linkcount) - { - int errs = 0; - if (node != m_root && node->m_suffix) - { - const string & np = node->m_pathstring; - const string & sp = node->m_suffix->m_pathstring; - if (np.substr(1, np.length() -1) != sp) - { - cerr << "Suffix Link Mismatch!!" << endl; - node->printLabel(cerr, m_string) << ": " << np << endl; - node->m_suffix->printLabel(cerr, m_string) << ": " << sp << endl; - errs++; - } - - linkcount++; - } - - if (node == m_root && node->m_suffix != m_root) - { - cerr << "Error m_root suffix != m_root !!!" << endl; - errs++; - } - - int childcount = 0; - for (int b = 0; b < basecount; b++) - { - if (node->m_children[b]) - { - childcount++; - errs += verifyNodeSuffixLinks(node->m_children[b], linkcount); - } - } - - if (childcount && !node->m_suffix) - { - errs++; - node->printLabel(cerr, m_string) << " has no suffix link!!!" << endl; - } - - return errs; - } - - void verifySuffixLinks() - { - cerr << endl; - cerr << "Verifing links" << endl; - setNodePath(m_root, ""); - int linkcount = 0; - int err = verifyNodeSuffixLinks(m_root, linkcount); - cerr << err << " suffix link errors detected" << endl; - cerr << linkcount << " suffix links checked" << endl; - - if (err) { exit(1); } - } -#endif - - - void buildUkkonen() - { - int len = m_strlen - 1; // length of the string, not of the buffer (remove s) - char base = m_string[1]; - - if (DEBUG) - { - cerr << "Building Ukkonen Tree for " << m_string << endl - << "Len: " << len << endl; - } - - // Construct T1 -#ifdef MPOOL - SuffixNode * node = new (&m_pool) SuffixNode(1, len, 1, m_root, NULL); // leaf: 1 -#else - SuffixNode * node = new SuffixNode(1, len, 1, m_root, NULL); // leaf: 1 -#endif - m_root->m_children[b2i(base)] = node; - SuffixNode * firstleaf = node; - SuffixNode * lastleaf = node; - - if (DEBUG) - { cerr << "Phase 1 Child: "; node->printLabel(cerr, m_string) << endl; } - - int startj = 2; - - // phase i+1 - for (int i = 2; i <= len; i++) - { - // Start at the last leaf created which will allow easy - // access to the node for startj - node = lastleaf; - int nodewalk = 0; - - // Keep track of last internal nodes created in split so we can add suffix links - SuffixNode * splitnode = NULL; - - if (!DOPHASETRICK) - { - startj = 2; - node = firstleaf; - } - - if (DEBUG) - { - char next = m_string[i]; - cerr << endl; - cerr << i << ".0 " << "Phase " << i << " adding " << next << " starting with " << startj << endl; - - string beta = substr(m_string, 1, i); - cerr << i << ".1" << " Extension 1: \"" << beta << "\" [implicit]" << endl; - } - - for (int j = startj; j <= i; j++) - { - // Goal: Ensure S[j .. i] (beta) is in the suffix tree - // Precondition: S[j-1 .. i] (alpha) is in the suffix tree "near" node - // All Internal nodes have a suffix link - - // Idea: 1) Remember where alpha is in the tree relative to node - // 2) Walk up the tree w bases until we get to a node with a suffix link. - // 3) Follow suffix link which shifts the path from S[j-1..i] to S[j..i] - // 4) Walk down tree in new location ensuring S[i-w .. i] is in tree - - // Notes: 1) All internal nodes have a suffix link by next extension - // 2) Any time we walk up to root, have to check S[j..i] - // 3) Suffix [1..i] is always present so start extension j with 2 - - int betapos = i; // The first position in string we need to check in tree - - if (DEBUG) - { - cerr << endl; - string beta = substr(m_string, j, i-j+1); - cerr << i << "." << j << " Phase " << i << " Extension " << j << ": \"" << beta << "\" bp:" << betapos << endl; - - cerr << i << "." << j << " Walking up from n:"; - node->printLabel(cerr, m_string) << " nw: " << nodewalk << endl; - } - - if (node == m_root) - { - // If we are at root, we have to check the full string s[j..i] anyways - } - else - { - if (nodewalk) - { - // partially walked down node->child, but didn't switch to child - // Match at i=6 on left... nodewalk=2, at 5 after suffix link - // 5 = i-2+1 - // o ----- o - // 5 A A 5 <- - // -> 6 T T 6 - - betapos -= nodewalk-1; - - if (DEBUG) - { - cerr << i << "." << j << " Adjusted nw: " << nodewalk << endl; - } - } - else - { - // Exactly at a node or leaf. - // Walk up to parent, subtracting length of that edge - int len = node->len(i); - betapos -= len-1; - node = node->m_parent; - - if (DEBUG) - { - cerr << i << "." << j << " Adjusted len: " << len << endl; - } - } - - if (DEBUG) - { - cerr << i << "." << j << " parent bp: " << betapos << " n:"; - node->printLabel(cerr, m_string) << endl; - } - - if (node->m_suffix == NULL) - { - // Subtract entire edge length - betapos -= node->len(i); - node = node->m_parent; - - if (DEBUG) - { - cerr << i << "." << j << " grandparent bp: " << betapos << " n:"; - node->printLabel(cerr, m_string) << endl; - } - - #if VERIFY - if (node->m_suffix == NULL) - { - cerr << "Missing suffix link!!! "; - exit(1); - } - #endif - } - } - - // jump across suffix link - node = node->m_suffix; - if (node == m_root) { betapos = j; } // have to check full string - - if (DEBUG) - { - cerr << i << "." << j << " Starting to walk down from bp: " << betapos << " to " << i << " n:"; - node->printLabel(cerr, m_string) << endl; - } - - if (FORCEROOT && node != m_root) - { - node = m_root; - betapos = j; - - if (DEBUG) - { - cerr << i << "." << j << " AtRoot bp: " << betapos << endl; - } - } - - bool done = false; - startj = j+1; // assume this extension should be skipped in the next phase - - while ((betapos <= i) && !done) - { - char base = m_string[betapos]; - unsigned char b = b2i(base); - SuffixNode * child = node->m_children[b]; - - if (DEBUG) - { - cerr << i << "." << j << " node betapos: " << betapos << "[" << base << "] n:"; - node->printLabel(cerr, m_string) << " "; - if (child) { cerr << "c: "; child->printLabel(cerr, m_string); } - cerr << endl; - } - - if (!child) - { - if (splitnode && betapos == splitnode->m_start) - { - if (DEBUG) - { - cerr << i << "." << j << " Add SL1: "; - splitnode->m_parent->printLabel(cerr, m_string) << " sl-> "; - node->printLabel(cerr, m_string) << endl; - } - - splitnode->m_parent->m_suffix = node; - splitnode = NULL; - } - -#ifdef MPOOL - SuffixNode * newnode = new (&m_pool) SuffixNode(betapos, len, j, node, NULL); // leaf: j -#else - SuffixNode * newnode = new SuffixNode(betapos, len, j, node, NULL); // leaf: j -#endif - node->m_children[b] = newnode; - lastleaf = newnode; - - if (DEBUG) - { - cerr << i << "." << j << " New Node: "; - newnode->printLabel(cerr, m_string) << endl; - } - - node = newnode; - - // This is the first base that differs, but the edgelength to - // i may be longer. Therefore set nodewalk to 0, so the entire - // edge is subtracted. - nodewalk = 0; - done = true; - break; - } - else - { - int nodepos = child->m_start; - nodewalk = 0; - - char nodebase = m_string[nodepos]; - - #if VERIFY - if (nodebase != base) - { - char nb = m_string[nodepos]; - cerr << "ERROR: first base on edge doesn't match edge label" << endl; - cerr << " nb: " << nb << " base: " << base << endl; - exit(1); - } - #endif - - // By construction, the string from j-1 to betapos to i-1 - // must already by present in the suffix tree - // Therefore, we can skip checking every character, and zoom - // to exactly the right character, possibly skipping the entire edge - - if (DOJUMP) - { - int mustmatch = i-1 - betapos + 1; - int childlen = child->len(i); - - if (mustmatch >= childlen) - { - betapos += childlen; - nodepos += childlen; - - skippedbases += childlen; - - if (DEBUG) - { - cerr << i << "." << j << " Edge Jump by: " << childlen << " new bp: " << betapos << " np: " << nodepos << endl; - } - - #if VERIFY - if (nodepos != child->m_end+1) - { - cerr << "ERROR: jump should have skipped entire edge, but didn't!" << endl; - exit(1); - } - #endif - } - else if (mustmatch) - { - betapos += mustmatch; - nodepos += mustmatch; - nodewalk += mustmatch; - - skippedbases += mustmatch; - - if (DEBUG) - { - cerr << i << "." << j << " Partial Jump by: " << mustmatch << " new bp: " << betapos << " np: " << nodepos << endl; - } - - #if VERIFY - if (VERIFY) - { - if (m_string[betapos-1] != m_string[nodepos-1]) - { - cerr << "ERROR: jump should have matched at least the mustmatch-1 characters" << endl; - cerr << "s[bp-1]: " << m_string[betapos-1] << " s[np-1]: " << m_string[nodepos-1] << endl; - exit(1); - } - } - #endif - } - } - - while (nodepos <= child->m_end && betapos <= i) - { - nodebase = m_string[nodepos]; - - #if VERBOSE - cerr << i << "." << j << " child bp: " << betapos << "[" << m_string[betapos] - << "] nb [" << nodebase << "]" << endl; - #endif - - if (m_string[betapos] == nodebase) - { - if (splitnode && betapos == splitnode->m_start) - { - if (DEBUG) - { - cerr << i << "." << j << " Add SL2: "; - splitnode->m_parent->printLabel(cerr, m_string) << " sl-> "; - node->printLabel(cerr, m_string) << endl; - } - - splitnode->m_parent->m_suffix = node; - splitnode = NULL; - } - - nodepos++; betapos++; nodewalk++; - - if (betapos == i+1) - { - if (DEBUG) - { - cerr << i << "." << j << " Internal edge match nw: " << nodewalk << endl; - } - - if ((nodewalk == child->len(i)) && (child->m_end == len)) - { - // we walked the whole edge to leaf, implicit rule I extension - if (DEBUG) - { - cerr << i << "." << j << " Leaf Node, Implicit Rule I Extension" << endl; - } - } - else - { - // "Real" rule III implicit extension - - // The j-1 extension was the last explicit extension in this round - // Start the next round at the last explicit extension - if (DOPHASETRICK) - { - startj = j; - - int skip = startj - 2; - - if (DEBUG) - { - cerr << i << "." << j << " Implicit Extension... start next phase at " << startj << ", saved " << skip << endl; - } - - skippedextensions += skip; - } - - if (DOINTERNALSKIP) - { - // Since we hit an internal match on a non-leaf, we know every other - // extension in this phase will also hit an internal match. - - // Have to be careful since leafs get the full string immediately, but - // they really have a Rule 1 extension - - int skip = i-j; - - if (DEBUG) - { - cerr << i << "." << j << " Implicit Extension... skipping rest of phase, saved " << skip << endl; - } - - skippedextensions += skip; - j = i+1; - } - } - - done = true; - } - } - else - { - if (DEBUG) { cerr << i << "." << j << " Spliting "; child->printLabel(cerr, m_string); } - - // Split is a copy of the child with the end shifted - // Then adjust start of child -#ifdef MPOOL - SuffixNode * split = new (&m_pool) SuffixNode(child->m_start, nodepos-1, 0, node, NULL); // internal -#else - SuffixNode * split = new SuffixNode(child->m_start, nodepos-1, 0, node, NULL); // internal -#endif - - split->m_children[b2i(nodebase)] = child; - child->m_start = nodepos; - child->m_parent = split; - - if (DEBUG) - { - cerr << " => "; - split->printLabel(cerr, m_string) << " + "; - child->printLabel(cerr, m_string) << endl; - } - - node->m_children[b] = split; - node = split; - - if (splitnode && betapos == splitnode->m_start) - { - if (DEBUG) - { - cerr << i << "." << j << " Add SL3: "; - splitnode->m_parent->printLabel(cerr, m_string) << " sl-> "; - node->printLabel(cerr, m_string) << endl; - } - - splitnode->m_parent->m_suffix = split; - splitnode = NULL; - } - - // Now create the new node -#ifdef MPOOL - SuffixNode * newnode = new (&m_pool) SuffixNode(betapos, len, j, split, NULL); // leaf j -#else - SuffixNode * newnode = new SuffixNode(betapos, len, j, split, NULL); // leaf j -#endif - lastleaf = newnode; - - split->m_children[b2i(m_string[betapos])] = newnode; - splitnode = newnode; - - node = newnode; - - if (DEBUG) - { - cerr << i << "." << j << " Split New Node: "; - newnode->printLabel(cerr, m_string) << endl; - } - - // This is the first base that differs, but the edgelength to - // i may be longer. Therefore set nodewalk to 0, so the entire - // edge is subtracted. - nodewalk = 0; - done = true; - break; - } - } - } - - if (!done) { node = child; } - } - } - - #if VERIFY - if (VERIFY) { verifySuffixLinks(); } - #endif - } - } -}; - - -SuffixTree * gtree = NULL; - -void buildUkkonenSuffixTree(const char * str) -{ - gtree = new SuffixTree(str); - gtree->buildUkkonen(); -} - -static const int MAX_TEXTURE_DIMENSION = 4096; -static const int BLOCKSIZE = 32; - -inline TextureAddress id2addr(int id) -{ - TextureAddress retval; - - int bigx = id & 0x1FFFF; - int bigy = id >> 17; - retval.y = (bigy << 5) + (bigx & 0x1F); - retval.x = bigx >> 5; - - return retval; -} - -void buildNodeTexture(SuffixNode * node, - PixelOfNode * nodeTexture, - PixelOfChildren * childrenTexture, - AuxiliaryNodeData aux_data[], - const char * refstr) -{ - int id = node->id(); - - aux_data[id].length = node->len(); - aux_data[id].depth = node->depth(); - aux_data[id].leafid = node->m_leafid; - aux_data[id].parent = id2addr(node->m_parent->id()); - - if (aux_data[id].leafid != 0) - { - aux_data[id].leafchar = refstr[aux_data[id].leafid-1]; - } - else - { - aux_data[id].leafchar = 0; - } - - TextureAddress myaddress(id2addr(id)); - id = myaddress.x + myaddress.y*MAX_TEXTURE_DIMENSION; - - nodeTexture[id].start = node->m_start; - nodeTexture[id].end = node->m_end; - nodeTexture[id].suffix = id2addr(node->m_suffix->id()); - - - for (int i = 0; i < basecount; i++) - { - if (node->m_children[i]) - { - TextureAddress childaddr = id2addr(node->m_children[i]->id()); - - // Unfortunately, the $ link doesn't fit into PixelOfChildren - if (i == b2i('$')) - { - nodeTexture[id].childD = childaddr; - } - else - { - childrenTexture[id].children[i] = childaddr; - } - - buildNodeTexture(node->m_children[i], nodeTexture, childrenTexture, aux_data, refstr); - } - } -} - -void buildSuffixTreeTexture(PixelOfNode** nodeTexture, - PixelOfChildren **childrenTexture, - unsigned int* width, unsigned int* height, - AuxiliaryNodeData **aux_data, - const char * refstr) -{ - assert(SuffixNode::s_nodecount < MAX_TEXTURE_DIMENSION*MAX_TEXTURE_DIMENSION); - assert(sizeof(PixelOfNode) == 16); - assert(sizeof(PixelOfChildren) == 16); - - // Leave space for NULL node - int allnodes = SuffixNode::s_nodecount+1; - - *width = MAX_TEXTURE_DIMENSION; - *height = (int)ceil((allnodes+0.0) / MAX_TEXTURE_DIMENSION)+BLOCKSIZE; - - // allocate space for the node and children textures - *nodeTexture = (PixelOfNode*) calloc((*width)*(*height), sizeof(PixelOfNode)); - *childrenTexture = (PixelOfChildren*) calloc((*width)*(*height), sizeof(PixelOfChildren)); - - *aux_data = (AuxiliaryNodeData*)calloc(SuffixNode::s_nodecount + 1, sizeof(AuxiliaryNodeData)); - - if (!*nodeTexture || !*childrenTexture || !*aux_data) - { - printf("arg. texture allocation failed.\n"); - exit(-1); - } - - buildNodeTexture(gtree->m_root, *nodeTexture, *childrenTexture, *aux_data, refstr); -}; - - -void printTreeTexture(const char * texfilename, - PixelOfNode * nodeTexture, - PixelOfChildren * childrenTexture, - int nodecount) -{ - cerr << "Printing tree texture to " << texfilename << endl; - - ofstream texfile; - texfile.open(texfilename, ofstream::out | ofstream::trunc); - - texfile << "id\tx\ty\tstart\tend\ta.x\ta.y\tc.x\tc.y\tg.x\tg.y\tt.x\tt.y\t$.x\t$.y" << endl; - for (int i = 0; i < nodecount; i++) - { - TextureAddress myaddress(id2addr(i)); - - texfile << i << "\t" - << myaddress.x << "\t" - << myaddress.y << "\t" - << nodeTexture[i].start << "\t" - << nodeTexture[i].end << "\t"; - - for (int j = 0; j < 4; j++) - { - texfile << childrenTexture[i].children[j].x << "\t"; - texfile << childrenTexture[i].children[j].y << "\t"; - } - - texfile << nodeTexture[i].childD.x << "\t"; - texfile << nodeTexture[i].childD.y << endl; - } - - texfile.close(); -} - -void renumberTree() -{ - queue > nodequeue; - - nodequeue.push(make_pair(gtree->m_root,0)); - int nodecount = 0; - - while(!nodequeue.empty()) - { - pair npair = nodequeue.front(); nodequeue.pop(); - - SuffixNode * node = npair.first; - int depth = npair.second; - - node->m_nodeid = ++nodecount; - - if (depth < 16) - { - for (int i = 0; i < basecount; i++) - { - SuffixNode * child = node->m_children[i]; - if (child) { nodequeue.push(make_pair(child,depth+1)); } - } - } - else - { - for (int i = 0; i < basecount; i++) - { - SuffixNode * child = node->m_children[i]; - if (child) - { - child->m_nodeid = ++nodecount; - - for(int j = 0; j < basecount; j++) - { - SuffixNode * gchild = child->m_children[j]; - - if (gchild) - { - gchild->m_nodeid = ++nodecount; - for (int k = 0; k < basecount; k++) - { - SuffixNode * ggchild = gchild->m_children[k]; - - if (ggchild) - { - ggchild->m_nodeid = ++nodecount; - - for (int l = 0; l < basecount; l++) - { - SuffixNode * gggchild = ggchild->m_children[l]; - - if (gggchild) - { - gggchild->m_nodeid = ++nodecount; - - for (int m = 0; m < basecount; m++) - { - SuffixNode * ggggchild = gggchild->m_children[m]; - if (ggggchild){ nodequeue.push(make_pair(ggggchild, depth+5)); } - } - } - } - } - } - } - } - } - } - } - } -} - - -extern "C" -void createTreeTexture(const char * refstr, - PixelOfNode** nodeTexture, - PixelOfChildren** childrenTexture, - unsigned int* width, unsigned int* height, - AuxiliaryNodeData** aux_data, - int* num_nodes, - const char * dotfilename, - const char * texfilename) -{ - cerr << " Creating Suffix Tree... "; - EventTime_t btimer; - SuffixNode::s_nodecount = 0; - buildUkkonenSuffixTree(refstr); - cerr << SuffixNode::s_nodecount << " nodes " - << btimer.str(true, 5) << endl; - - cerr << " Renumbering tree... "; - EventTime_t rtimer; - renumberTree(); - cerr << rtimer.str(true, 5) << endl; - - - EventTime_t ftimer; - cerr << " Flattening Tree... "; - // indexTree(); - buildSuffixTreeTexture(nodeTexture, - childrenTexture, - width, height, - aux_data, - gtree->m_string); - - *num_nodes = SuffixNode::s_nodecount; - cerr << ftimer.str(true, 5) << endl; - - if (dotfilename) - { - gtree->printDot(dotfilename); - } - - if (texfilename) - { - printTreeTexture(texfilename, *nodeTexture, *childrenTexture, SuffixNode::s_nodecount+1); - } - - delete gtree; - gtree = NULL; -} - - -extern "C" -void getReferenceString(const char * filename, char** refstr, size_t* reflen) -{ - EventTime_t timer; - cerr << "Loading ref: " << filename << "... "; - - string S="s"; - - ifstream file; - file.open(filename); - - if (!file) - { - cerr << "Can't open " << filename << endl; - exit (1); - } - - // Skip over the reference header line - char refline[2048]; - file.getline(refline, sizeof(refline)); - - if (refline[0] != '>') - { - cerr << endl - << "ERROR: Reference file is not in FASTA format" - << endl; - } - - // Now read the reference string - string buffer; - while (file >> buffer) - { - if (buffer[0] == '>') - { - cerr << endl - << "ERROR: Only a single reference sequence is supported!" - << endl; - - exit (1); - } - else - { - for (unsigned int i = 0; i < buffer.length(); i++) - { - char b = toupper(buffer[i]); - if (b == 'A' || b == 'C' || b == 'G' || b=='T') - { - S += b; - } - else - { - S += 'A'; - } - } - } - } - - S += "$"; - - *refstr = strdup(S.c_str()); - *reflen = strlen(*refstr) + 1; - - cerr << *reflen-3 << " bp. " << timer.str(true, 5) << endl; -} - -inline void addChar(char **buf, int * size, int * pos, char c) -{ - if (*pos == *size) - { - (*size) *= 2; // double the size of the buffer - *buf = (char *) realloc(*buf, *size); - if (!*buf) - { - cerr << "ERROR: Realloc failed, requested: " << *size << endl; - } - } - - (*buf)[*pos] = c; - (*pos)++; -} - -inline size_t bytesNeededOnGPU(unsigned int querylen, int min_match_len) -{ - if (min_match_len == -1) - return sizeof(MatchCoord) + (querylen + 10); - else - return sizeof(MatchCoord) * (querylen - min_match_len + 1) + (querylen + 10); -} - -//Gets up to set_size queries. -extern "C" -void getQueriesTexture(int qfile, - char** queryTexture, - size_t* queryTextureSize, - int** queryAddrs, - char*** queryNames, - int** queryLengths, - unsigned int* numQueries, - size_t memory_avail, - int min_match_length, - bool rc) -{ - EventTime_t timer; -//fprintf(stderr,"1"); - int qstringpos = 0; - int qstringsize = 1024*1024; - char * qstring = (char *) malloc(qstringsize); - - bool resetAmbiguity = true; - - // offset of query i in qstring - int offsetspos = 0; - int offsetssize = 1024; - int * offsets = (int *) malloc(offsetssize * sizeof(int)); - int * lengths = (int *) malloc(offsetssize * sizeof(int)); - - int qrylen = 0; - int this_qrylen = 0; - - int bytes_read; - unsigned char buf[32*1024]; - - vector names; - string header; - bool inheader = false; - int total_read = 0; - - unsigned char dnachar [256]; - - bool set_full = false; -//fprintf(stderr,"2"); - // tracks the GPU memory needed by the queries read so far. - size_t curr_mem_usage = 0; - - for (int i = 0; i < 256; i++) - { - dnachar[i] = 0; - } - - dnachar[(unsigned char) 'A'] = 1; - dnachar[(unsigned char) 'C'] = 1; - dnachar[(unsigned char) 'G'] = 1; - dnachar[(unsigned char) 'T'] = 1; -//fprintf(stderr,"3"); - while ((bytes_read = read(qfile, buf, sizeof(buf))) != 0) - { - // cerr << "bytes_read: " << bytes_read << endl; - - if (bytes_read == -1) - { - cerr << "ERROR: Error reading file: " << errno << endl; - exit(1); - } - - int i = 0; - - if (inheader) - { - // Handle case where last read was inside a header - for (; i < bytes_read; i++) - { - if (buf[i] == '\n') - { - inheader = false; - i++; - char* name = strdup(header.c_str()); - names.push_back(name); - header.clear(); - break; - } - else - { - header.insert(header.end(), buf[i]); - } - } - } -// fprintf(stderr,"4"); - for (; i < bytes_read; i++) - { - unsigned char b = toupper(buf[i]); - - if (b == '>') - { - - if (curr_mem_usage >= memory_avail) - { - set_full = true; - off_t seek = lseek(qfile, -(bytes_read - i), SEEK_CUR); - if (seek == (off_t)-1) - { - cerr<< "lseek failed: "<< errno< %s\n", names.back()); - if (rc) - printf("> %s Reverse\n", names.back()); - names.pop_back(); - --offsetspos; - qstringpos -= this_qrylen + 1; - } - else - { - addChar(&qstring, &qstringsize, &qstringpos, '\0'); - lengths[offsetspos - 1] = this_qrylen; - curr_mem_usage += bytesNeededOnGPU(this_qrylen, min_match_length); - } - } -// fprintf(stderr,"6"); - if (offsetspos == offsetssize) - { - offsetssize *= 2; - offsets = (int *) realloc(offsets, sizeof(int)*offsetssize); - lengths = (int *) realloc(lengths, sizeof(int)*offsetssize); - if (!offsets || !lengths) - { - cerr << endl - << "ERROR: Realloc failed: requested " - << sizeof(int) * offsetssize << endl; - exit(1); - } - } - - offsets[offsetspos++] = qstringpos; - - inheader = true; - - // Try to walk out of header - for (i++; i < bytes_read; i++) - { - if (buf[i] == '\n') - { - inheader = false; - char* name = strdup(header.c_str()); - names.push_back(name); - header.clear(); - break; - } - else - { - header.insert(header.end(), buf[i]); - } - } - - addChar(&qstring, &qstringsize, &qstringpos, 'q'); - this_qrylen = 0; - } - else if (dnachar[b]) - { - addChar(&qstring, &qstringsize, &qstringpos, b); - qrylen++; - this_qrylen++; - } - else if (isspace(b)) - { - - } - else if (resetAmbiguity) - { - addChar(&qstring, &qstringsize, &qstringpos, 'x'); - this_qrylen++; - } - else - { - cerr << endl - << "ERROR: Unexpected character: " << buf[i] - << " in query file at byte: " << total_read+i << endl; - exit(1); - } - } -// fprintf(stderr,"7"); - if (set_full) - break; - - total_read += bytes_read; - } - - if (qstringpos) - { - if (this_qrylen < min_match_length) - { - printf("> %s\n", names.back()); - if (rc) - printf("> %s Reverse\n", names.back()); - names.pop_back(); - --offsetspos; - qstringpos -= this_qrylen + 1; - } - else - { - addChar(&qstring, &qstringsize, &qstringpos, '\0'); - lengths[offsetspos - 1] = this_qrylen; - curr_mem_usage += bytesNeededOnGPU(this_qrylen, min_match_length); - } - } - - *numQueries = offsetspos; - - if (offsetspos == 0) - { - free(offsets); - free(lengths); - free(qstring); - *queryAddrs = NULL; - *queryTexture = NULL; - *queryTextureSize = 0; - *queryNames = NULL; - - return; - } - - - *queryAddrs = offsets; - - *queryTexture = qstring; - *queryTextureSize = qstringpos; - *queryNames = (char**)malloc(names.size() * sizeof(char*)); - *queryLengths = lengths; -//fprintf(stderr,"8"); - for (unsigned int i = 0; i < *numQueries; ++i) - { - *(*queryNames + i) = names[i]; - } - cerr << offsetspos << " queries (" - << qrylen << " bp), need " - << curr_mem_usage << " bytes on the GPU " - << timer.str(true, 5) << endl; -} - -struct pathblock -{ - TextureAddress node_addr; - int string_depth; -}; - -#define __USE_BUFFERED_IO__ - -static const size_t output_buf_limit = 32*1024; -char output_buf[output_buf_limit]; - -//FIXME: needs to be reinitialized to zero at the beginning of each round of printing. -size_t bytes_written = 0; - -int addToBuffer(char* string) -{ - size_t buf_length = strlen(string); - - if (buf_length + bytes_written>= output_buf_limit) - { - size_t chunk = (output_buf_limit - bytes_written - 1); - strncpy(output_buf + bytes_written, string, chunk); - output_buf[bytes_written + chunk] = 0; - printf("%s", output_buf); - strncpy(output_buf, string + chunk, buf_length - chunk); - bytes_written = buf_length - chunk; - } - else - { - strncpy(output_buf + bytes_written, string, buf_length); - bytes_written += buf_length; - } - return 0; -} - - -inline int addr2id(TextureAddress addr) -{ - int blocky = addr.y & 0x1F; - int bigy = addr.y >> 5; - int bigx = (addr.x << 5) + blocky; - return bigx + (bigy << 17); -} - -#define CHILDREN(node_addr) ((((PixelOfChildren*)(page->ref.h_children_tex_array)) + (node_addr.x) + ((node_addr.y) * MAX_TEXTURE_DIMENSION))->children) -#define DOLLAR(node_addr) ((((PixelOfNode*) (page->ref.h_node_tex_array)) + (node_addr.x) + ((node_addr.y) * MAX_TEXTURE_DIMENSION))->childD) -#define LEAFID(x) (page->ref.aux_data[x].leafid) -#define LEAFCHAR(x) (page->ref.aux_data[x].leafchar) - -char buf[256]; - -void printNodeAlignments(const char* ref, - const ReferencePage* page, - const char queryflankingbase, - const TextureAddress node, - const int qrypos, - int qrylen, - const pathblock path[], - int path_idx, - bool rc) -{ - int nodeid = addr2id(node); - char isLeaf = LEAFCHAR(nodeid); - - if (path[path_idx].node_addr.data == node.data) - { - qrylen = path[path_idx].string_depth; - path_idx--; - } - - if (isLeaf) - { - if (isLeaf != queryflankingbase) - { - int leafid = LEAFID(nodeid); - int left_in_ref = (leafid - 1) + page->begin; - int right_in_ref = left_in_ref + qrylen; - - if ((left_in_ref != page->begin || page->shadow_left == -1) && - (right_in_ref != page->end || page->shadow_right == -1)) - { - if (!(left_in_ref > page->begin && right_in_ref < page->shadow_left)) - { - //sprintf(buf, "\t%d\t%d\t%d\n", node->m_leafid, qrypos, qrylen); - sprintf(buf, "%8d%10d%10d\n", left_in_ref, qrypos, qrylen); - addToBuffer(buf); - } - } - } - } - else - { - TextureAddress* children = CHILDREN(node); - for (int i = 0; i < basecount - 1; ++i) - { - if ((children + i)->data) - { - printNodeAlignments(ref, page, queryflankingbase, *(children+i), - qrypos, qrylen, path, path_idx, rc); - } - } - - TextureAddress dollar = DOLLAR(node); - if (dollar.data) - { - printNodeAlignments(ref, page, queryflankingbase, dollar, - qrypos, qrylen, path, path_idx, rc); - } - } -} - -void flushOutput() -{ - if (bytes_written) - { - output_buf[bytes_written] = 0; - printf("%s", output_buf); - bytes_written = 0; - } -} - - -//FIXME: hardcoded path buffer, needs to be as long as the longest query in the query set. -pathblock path[8192]; - -#define NODE_SDEPTH(x) (page->ref.aux_data[x].depth) -#define NODE_LENGTH(x) (page->ref.aux_data[x].length) -#define NODE_PARENT(x) (page->ref.aux_data[x].parent) - -char RC(char c) -{ - switch(c) - { - case 'A': return 'T'; - case 'C': return 'G'; - case 'G': return 'C'; - case 'T': return 'A'; - case 'q': return '\0'; - default: return c; - }; -} - -void printAlignments(char* ref, - ReferencePage* page, - char* query, - int qrylen, - int nodeid, - int qrypos, - int edge_match, - int min_match, - bool rc, - bool forwardcoordinates) -{ - TextureAddress node_addr = id2addr(nodeid); - TextureAddress prev; - prev.data = 0; - - int path_idx = 0; - int string_depth = NODE_SDEPTH(nodeid) - 1; - - if (edge_match > 0) - { - string_depth = NODE_SDEPTH(nodeid) - (NODE_LENGTH(nodeid) - edge_match) - 1; - } - else - { - edge_match = NODE_LENGTH(nodeid); - } - - path[path_idx].node_addr = node_addr; - path[path_idx].string_depth = string_depth; - path_idx++; - string_depth -= edge_match; - prev = node_addr; - - node_addr = NODE_PARENT(nodeid); - - while ((node_addr.data) && string_depth >= min_match) - { - nodeid = addr2id(node_addr); - path[path_idx].node_addr = node_addr; - path[path_idx].string_depth = string_depth; - path_idx++; - string_depth -= NODE_LENGTH(nodeid); - - prev = node_addr; - node_addr = NODE_PARENT(nodeid); - } - - char flankingbase = query[qrypos]; - - if (rc) - { - flankingbase = RC(query[strlen(query)-qrypos]); - if (forwardcoordinates) { qrypos = qrylen - 1 - qrypos; } - } - - printNodeAlignments(ref, page, flankingbase, prev, qrypos + 1, - NODE_SDEPTH(addr2id(prev)), path, path_idx - 1, rc); -} -- cgit v1.3