Merge branch 'master' of https://github.com/jasperlinthorst/reveal

Conflicts: reveallib/reveal.c
jasperlinthorst · Aug 20, 2021 · d5f5b9a · d5f5b9a
2 parents 8beaf90 + 6a0d70d
commit d5f5b9a
Show file tree

Hide file tree

Showing 4 changed files with 41 additions and 12 deletions.
diff --git a/reveal/rem.py b/reveal/rem.py
@@ -323,16 +323,24 @@ def graphalign(index,mum):
         isize=index.n
         mns=[]
         topop=[]
-        logging.debug("Nodes in subgraph:")
-        for node in nodes:
-            logging.debug("%s"%str(node))
+        #logging.debug("Nodes in subgraph:")
+        #for node in nodes:
+        #    logging.debug("%s"%str(node))
 
         sp=[sp for gid,sp in spd]
         matching=set()
         for pos in sp:
             matching.add((pos,pos+l))
-            # logging.debug("Lookup node for sp=%d"%pos)
+            logging.debug("Lookup node for sp=%d"%pos)
             old=t[pos].pop()
+            logging.debug("Node=%s"%str(old))
+
+            if old.end-pos<l:
+                import pickle
+                open("T.txt",'w').write(index.T)
+                pickle.dump(index.LCP,open("LCP.pickle",'w'))
+                pickle.dump(t,open("tree.pickle",'w'))
+
             assert(old.end-old.begin>=l)
             mn,other=breaknode(old,pos,l)
             mns.append(mn)
@@ -566,7 +574,6 @@ def align_genomes(args):
 
     logging.info("Constructing index...")
     idx.construct()
-
     logging.info("Done.")
 
     if len(args.inputfiles)==2 and not graph:

diff --git a/reveal/schemes.py b/reveal/schemes.py
@@ -137,6 +137,20 @@ def lookup(mum):
     qrpoint=dict()
     for pos in sp:
         t=ts[pos]
+        if len(t)!=1:
+            import pickle
+            logging.debug("Dumping tree to file...")
+            pickle.dump(ts, open('tree.pickle','w'))
+            logging.debug("mum: %s"%str(mum))
+            print(type(pos))
+            print(ts[pos+1])
+            print(ts[pos-1])
+            logging.debug("l: %s"%str(l))
+            logging.debug("pos: %s"%str(pos))
+            logging.debug("sp: %s"%str(sp))
+            logging.debug("t: %s"%str(t))
+            #logging.debug("ts: %s"%str(ts))
+
         assert(len(t)==1)
         node=iter(t).next()
         ndata=G.node[node]
@@ -230,7 +244,7 @@ def graphmumpicker(mums,idx,precomputed=False,minlength=0):
                 logging.debug("No MUMS that span all input genomes, segment genomes.")
                 mmums=segment(mums)
                 logging.debug("Segmented genomes/graphs into %s, now %d MUMS for chaining."%(mmums[0][2],len(mmums)))
-
+            
             if args.trim:
                 logging.debug("Trimming overlap between mums.")
                 mmums=trim_overlap(mmums)

diff --git a/reveallib/interface.c b/reveallib/interface.c
@@ -269,9 +269,18 @@ static PyObject *construct(RevealIndex *self, PyObject *args, PyObject *keywds)
             return NULL;
          };
     }
-
-    //TODO: sanity check!
-
+
+#ifdef REVEALDEBUG
+//    int j=0;
+//    for (i=1; i<self->n; i++){
+//        for (j=0; i<self->LCP[i]; j++){
+//            assert(self->T[self->SA[i]+j]!='$'); //Make sure lcp intervals dont span sentinels
+//            assert(self->T[self->SA[i-1]+j]!='$');
+//        }
+//    }
+//
+//    //TODO: sanity check!
+#endif
     //if caching is specified write sa and lcp to disk
     if (self->cache==1){
         fprintf(stderr,"Writing LCP and SA to disk...");

diff --git a/reveallib/reveal.c b/reveallib/reveal.c
@@ -309,8 +309,8 @@ PyObject * getmultimems(RevealIndex *index, PyObject *args, PyObject *keywds) {
     int maxdepth=1000;
     int *flag_so=calloc(mainidx->nsamples,sizeof *flag_so);
     lcp_t *stack_lcp=malloc(maxdepth * sizeof *stack_lcp);
-    lcp_t *stack_lb=malloc(maxdepth * sizeof *stack_lb);
-    lcp_t *stack_ub=malloc(maxdepth * sizeof *stack_ub);
+    saidx_t *stack_lb=malloc(maxdepth * sizeof *stack_lb);
+    saidx_t *stack_ub=malloc(maxdepth * sizeof *stack_ub);
     lcp_t i_lcp;
     int depth=0;
     saidx_t i,lb,i_lb,i_ub;
@@ -456,7 +456,6 @@ PyObject * getmultimums(RevealIndex *index, PyObject *args, PyObject *keywds) {
     lcp_t *stack_lcp=malloc(maxdepth * sizeof *stack_lcp);
     saidx_t *stack_lb=malloc(maxdepth * sizeof *stack_lb);
     saidx_t *stack_ub=malloc(maxdepth * sizeof *stack_ub);
-
     lcp_t i_lcp;
     int depth=0;
     saidx_t i,lb,i_lb,i_ub;