Merge pull request #21 from kb1dds/kernel_proj

Kernel projection for Sheaf.minimalExtend()

pysheaf/pysheaf.py

@@ -872,8 +872,45 @@ def minimalExtend(self,assignment, activeCells=None, testSupport=None, method='n
             ac=[idx for idx in range(len(self.cells)) if idx not in support]
         else:
             ac=activeCells
+
+        if method == 'KernelProj':
+            if not self.isLinear():
+                raise NotImplementedError('KernelProj only works for sheaves of vector spaces')
+
+            if ord != 2:
+                warn('Kernel projection requires order 2 in minimalExtend')
+
+            # Compile dictionary of rows
+            rowstarts=dict()
+            rowidx=0
+            for i in support:
+                rowstarts[i]=rowidx
+                rowidx+=self.cells[i].stalkDim
+
+            newassignment = Section([sc for sc in assignment.sectionCells])
+
+            # Optimize each active cell independently
+            for i in ac:
+                if self.cells[i].stalkDim > 0:
+                    # Matrix of all restrictions out of this cell into the support
+                    mat=np.zeros((sum([self.cells[j].stalkDim for j in support]),
+                                  self.cells[i].stalkDim))
+
+                    for cf in self.cofaces(i): # Iterate over all cofaces of this activeCell
+                        try:
+                            supportidx=support.index(cf.index)
+                            mat[rowstarts[supportidx]:rowstarts[supportidx]+self.cells[supportidx].stalkDim,:]=cf.restriction.matrix
+                        except ValueError:
+                            pass
 
-        if self.isNumeric():
+                # Use least squares to solve for assignment rooted at this cell given the existing assignment
+                asg,bnds=self.serializeAssignment(assignment,activeCells=support) # Confusingly, activeSupport here refers *only* to the support of the assignment
+                result=np.linalg.lstsq(mat,asg)
+
+                newassignment.sectionCells.append(SectionCell(i,result[0]))
+
+            return newassignment
+        elif self.isNumeric():
             initial_guess, bounds = self.serializeAssignment(assignment,ac)
             res=scipy.optimize.minimize( fun = lambda sec: self.consistencyRadius(self.deserializeAssignment(sec,ac,assignment), testSupport=testSupport, ord=ord),
                                          x0 = initial_guess,