Add containment

include/oxli/hashgraph.hh

@@ -294,6 +294,7 @@ public:
 
     void update_from(const Nodegraph &other);
     double similarity(const Nodegraph &other);
+    double containment(const Nodegraph &other);
 };
 
 }

include/oxli/storage.hh

@@ -227,6 +227,7 @@ public:
 
     void update_from(const BitStorage&);
     double similarity(const BitStorage&);
+    double containment(const BitStorage&);
 };
 
 

khmer/_oxli/graphs.pxd

@@ -205,6 +205,7 @@ cdef extern from "oxli/hashgraph.hh" namespace "oxli" nogil:
 
         void update_from(const CpNodegraph &) except +oxli_raise_py_error
         double similarity(const CpNodegraph &) except +oxli_raise_py_error
+        double containment(const CpNodegraph &) except +oxli_raise_py_error
 
 
 cdef extern from "oxli/labelhash.hh" namespace "oxli":

khmer/_oxli/graphs.pyx

@@ -869,3 +869,6 @@ cdef class Nodegraph(Hashgraph):
 
     def similarity(self, Nodegraph other):
         return deref(self._ng_this).similarity(deref(other._ng_this))
+
+    def containment(self, Nodegraph other):
+        return deref(self._ng_this).containment(deref(other._ng_this))

src/oxli/hashgraph.cc

@@ -923,6 +923,23 @@ double Nodegraph::similarity(const Nodegraph &otherBASE)
     }
 }
 
+double Nodegraph::containment(const Nodegraph &otherBASE)
+{
+    if (_ksize != otherBASE._ksize) {
+        throw oxli_exception("both nodegraphs must have same k size");
+    }
+    BitStorage * myself = dynamic_cast<BitStorage *>(this->store);
+    const BitStorage * other;
+    other = dynamic_cast<const BitStorage*>(otherBASE.store);
+
+    // if dynamic_cast worked, then the pointers will be not null.
+    if (myself && other) {
+        return myself->containment(*other);
+    } else {
+        throw oxli_exception("similarity failed with incompatible objects");
+    }
+}
+
 template void Hashgraph::consume_seqfile_and_tag<read_parsers::FastxReader>(
     std::string const &filename,
     unsigned int &total_reads,

src/oxli/storage.cc

@@ -123,6 +123,34 @@ double BitStorage::similarity(const BitStorage& other)
     return double(intersection) / double(union_size);
 }
 
+double BitStorage::containment(const BitStorage& other)
+{
+    if (_tablesizes != other._tablesizes) {
+        throw oxli_exception("both nodegraphs must have same table sizes");
+    }
+
+    uint64_t intersection = 0;
+    uint64_t self_size = 0;
+    for (unsigned int table_num = 0; table_num < _n_tables; table_num++) {
+        Byte * me = _counts[table_num];
+        Byte * ot = other._counts[table_num];
+        uint64_t tablesize = _tablesizes[table_num];
+        uint64_t tablebytes = tablesize / 8 + 1;
+
+        for (uint64_t index = 0; index < tablebytes; index++) {
+            // First, get how many values in common we have
+            intersection += __builtin_popcountll(me[index] & ot[index]);
+            self_size += __builtin_popcountll(me[index]);
+        }
+    }
+
+    if (self_size == 0) {
+        self_size = 1;
+    }
+
+    return double(intersection) / double(self_size);
+}
+
 void BitStorage::save(std::string outfilename, WordLength ksize)
 {
     if (!_counts[0]) {