Very preliminary general simplicial complex.

2 files changed, 196 insertions, 9 deletions

diff --git a/phstuff/diphawrapper.py b/phstuff/diphawrapper.py
index 2cead35..b2c9444 100644
--- a/phstuff/diphawrapper.py
+++ b/phstuff/diphawrapper.py
@@ -8,6 +8,7 @@ import tempfile
 import subprocess
 
 import phstuff.barcode as bc
+import phstuff.simplicial as simpl
 
 """Handle marshalling data into and out of DIPHA's formats, and
 optionally manage running the DIPHA executable.
@@ -156,7 +157,7 @@ def save_edge_list(fname, edge_list):
             for j in range(0, len(edge_list[i])):
                 f.write(struct.pack('<d', factor*edge_list[i][j][1]))
 
-def save_cubical(fname, complex):
+def save_cubical(fname, array):
     """Save a (any-dimensional) array that will be interpreted as a
     cubical complex with top-dimensional cells entering the filtration
     at the scale given by the array values.
@@ -168,7 +169,7 @@ def save_cubical(fname, complex):
     
     fname: Name of file to write.
 
-    complex: Any-dimensional array with floating-point values giving
+    array: Any-dimensional array with floating-point values giving
     the filtration values of the top-dimensional cells.
 
     """
@@ -176,14 +177,52 @@ def save_cubical(fname, complex):
     with open(fname, "wb") as f:
         f.write(struct.pack("<q", DIPHA_MAGIC))
         f.write(struct.pack("<q", DIPHA_IMAGE_DATA))
-        f.write(struct.pack("<q", np.product(complex.shape)))
-        f.write(struct.pack("<q", len(complex.shape)))
+        f.write(struct.pack("<q", np.product(array.shape)))
+        f.write(struct.pack("<q", len(array.shape)))
         
-        for n in complex.shape:
+        for n in array.shape:
             f.write(struct.pack("<q", n))
 
-        for w in complex.flat:
-            f.write(struct.pack("<d", w))
+        for w in array.flat:
+            f.write(struct.pack("<d", float(w)))
+
+def save_simplicial(fname, complex):
+    complex.order()
+    assert(complex.is_ordered())
+    
+    with open(fname, "wb") as f:
+        f.write(struct.pack("<q", DIPHA_MAGIC))
+        f.write(struct.pack("<q", DIPHA_WEIGHTED_BOUNDARY_MATRIX))
+        f.write(struct.pack("<q", 0))
+        f.write(struct.pack("<q", complex.size()))
+        f.write(struct.pack("<q", complex.top_dim()))
+
+        # This can be made a whole lot more efficient. We're iterating
+        # in a stupid way.
+        
+        offsets = []
+        total_nonzero = 0
+        for node in complex:
+            d = node.dimension()
+            f.write(struct.pack("<q", d))
+            offsets.append(total_nonzero)
+            if d != 0:
+                total_nonzero += d+1                
+
+        for node in complex:
+            f.write(struct.pack("<d", node.w))            
+
+        for off in offsets:
+            f.write(struct.pack("<q", off))
+
+        f.write(struct.pack("<q", total_nonzero))
+        for node in complex:
+            for bnode in node.boundary_nodes():
+                f.write(struct.pack("<q", bnode.i))
+            
+            
+    
+    
             
 
 def load_barcode(fname, top_dim = None):
@@ -368,10 +407,14 @@ class DiphaRunner:
         save_edge_list(self.ifile, edge_list)
         self.ready = True
 
-    def cubical(self, complex):
+    def cubical(self, array):
         """Initialize to run with a cubical complex. See `save_cubical`."""
 
-        save_cubical(self.ifile, complex)
+        save_cubical(self.ifile, array)
+        self.ready = True
+
+    def simplicial(self, complex):
+        save_simplicial(self.ifile, complex)
         self.ready = True
 
     def run(self):
diff --git a/phstuff/simplicial.py b/phstuff/simplicial.py
new file mode 100644
index 0000000..884a543
--- /dev/null
+++ b/phstuff/simplicial.py
@@ -0,0 +1,144 @@
+class Node:
+    def __init__(self, x, i, w, parent):
+        self.x = x
+        self.i = i
+        self.w = w
+        self.parent = parent
+        self.children = dict()
+
+    def is_root(self):
+        return self.parent is None
+
+    def simplex(self):
+        s = []
+        node = self
+        if node.is_root():
+            return None
+
+        while True:
+            s.append(node.x)
+            node = node.parent
+            if node.is_root():
+                break
+            
+        s.reverse()
+        return s
+
+    def boundary_nodes(self):
+        if self.is_root():
+            return None
+        if self.parent.is_root():
+            return []
+
+        ret = []
+        up = []
+        skip = self
+
+        while not skip.is_root():
+            down = list(up)
+            node = skip.parent
+            while len(down) != 0:
+                node = node.children[down.pop()]
+            ret.append(node)
+            up.append(skip.x)
+            skip = skip.parent
+
+        assert(len(ret) == self.dimension() + 1)
+        return ret
+
+    def boundary_simplices(self):
+        ret = []
+        for node in self.boundary_nodes():
+            ret.append(node.simplex())
+        return ret
+
+    def dimension(self):
+        ret = 0
+        if self.is_root():
+            return None
+
+        node = self
+        while not node.parent.is_root():
+            ret += 1
+            node = node.parent
+        return ret
+
+class ComplexIterator:
+    def __init__(self, start):
+        self.__to_go = start.children.values()
+
+    def next(self):
+        if len(self.__to_go) == 0:
+            raise StopIteration
+        else:
+            ret = self.__to_go.pop(0)
+            self.__to_go.extend(ret.children.values())
+            return ret
+
+class Complex:
+    def __init__(self):
+        self.__root = Node(None, None, None, None)
+        self.__count = 0
+        self.__top_dim = -1
+        self.__ordered = False
+
+    def size(self):
+        return self.__count
+
+    def top_dim(self):
+        return self.__top_dim
+
+    def is_ordered(self):
+        return self.__ordered
+
+    def order(self):
+        i = 0
+        for node in self:
+            node.i = i
+            i += 1
+        self.__ordered = True
+
+    def find(self, simplex): # The simplex must be sorted!
+        p = len(simplex) - 1
+
+        if p < -1:
+            return None
+
+        node = self.__root
+
+        for v in simplex:
+            if v in node.children:
+                node = node.children[v]
+            else:
+                return None
+        return node
+
+    def add(self, simplex, weight):
+        self.__ordered = False
+        simplex = sorted(simplex)
+
+        p = len(simplex) - 1
+        if p < 0:
+            return None
+
+        last = simplex[-1]
+        pre = simplex[0:p]
+
+        parent = self.find(pre)
+
+        node = Node(last, None, weight, parent)
+        self.__count += 1
+        self.__top_dim = max(self.__top_dim, p)
+        
+        return parent.children.setdefault(last, node)
+
+
+    def __contains__(self, simplex):
+        return self.find(simplex) is not None
+
+    def __iter__(self):
+        return ComplexIterator(self.__root)
+
+
+        
+