From 250bf31ab28ff2a1f33a3cba6b3b9bc09854ec3e Mon Sep 17 00:00:00 2001 From: Michael Robinson Date: Tue, 19 Jun 2018 10:47:20 -0400 Subject: [PATCH 1/4] Outline of tasks for faster linear algorithm for extending partial assignments --- pysheaf/pysheaf.py | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/pysheaf/pysheaf.py b/pysheaf/pysheaf.py index 1a444e6..a9b8e27 100644 --- a/pysheaf/pysheaf.py +++ b/pysheaf/pysheaf.py @@ -872,8 +872,21 @@ 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 fuseAssignment') + + # Build block matrix in which block columns are stalks over activeCells and block rows are stalks over cells supported in assignment + # Blocks are ultimate restriction maps (not just one hops!) + + # Use (constrained?) least squares to solve for global assignment given existing assignment - if self.isNumeric(): + # Deserialize assignment + 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, From 3e55fce67c69b4204e935b367a6eafbe34369baa Mon Sep 17 00:00:00 2001 From: Michael Robinson Date: Tue, 19 Jun 2018 14:07:35 -0400 Subject: [PATCH 2/4] draft of linear-algebraic Sheaf.minimalExtend --- pysheaf/pysheaf.py | 29 +++++++++++++++++++++++++++-- 1 file changed, 27 insertions(+), 2 deletions(-) diff --git a/pysheaf/pysheaf.py b/pysheaf/pysheaf.py index a9b8e27..cdd62c0 100644 --- a/pysheaf/pysheaf.py +++ b/pysheaf/pysheaf.py @@ -878,14 +878,39 @@ def minimalExtend(self,assignment, activeCells=None, testSupport=None, method='n raise NotImplementedError('KernelProj only works for sheaves of vector spaces') if ord != 2: - warn('Kernel projection requires order 2 in fuseAssignment') + warn('Kernel projection requires order 2 in minimalExtend') # Build block matrix in which block columns are stalks over activeCells and block rows are stalks over cells supported in assignment # Blocks are ultimate restriction maps (not just one hops!) + mat=np.zeros((sum([self.cells[i].stalkDim for i in support]), + sum([self.cells[i].stalkDim for i in ac]))) + + # Compile dictionary of rows + rowstarts=dict() + rowidx=0 + for i in support: + rowstarts[i]=rowidx + rowidx+=self.cells[i].stalkDim + + colidx=0 # Current column + for i in activeCells: + if self.cells[i].stalkDim > 0: + 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, + colidx:colidx+self.cells[i].stalkDim]=cf.restriction.matrix + except ValueError: + pass + + colidx+=self.cells[i].stalkDim - # Use (constrained?) least squares to solve for global assignment given existing assignment + # Use least squares to solve for global assignment given existing assignment + asg=self.serializeAssignment(assignment,activeCells=support) # Confusingly, activeSupport here refers *only* to the support of the assignment + result=np.linalg.lstsq(mat,asg) # Deserialize assignment + return self.deserializeAssignment(result[0],ac,assignment) 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), From 6df3e4295172b08080b83ced6cd1b8b7dbf803bf Mon Sep 17 00:00:00 2001 From: Michael Robinson Date: Tue, 19 Jun 2018 14:16:46 -0400 Subject: [PATCH 3/4] it at least runs, but perhaps not correctly! --- pysheaf/pysheaf.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/pysheaf/pysheaf.py b/pysheaf/pysheaf.py index cdd62c0..b62d0fe 100644 --- a/pysheaf/pysheaf.py +++ b/pysheaf/pysheaf.py @@ -893,7 +893,7 @@ def minimalExtend(self,assignment, activeCells=None, testSupport=None, method='n rowidx+=self.cells[i].stalkDim colidx=0 # Current column - for i in activeCells: + for i in ac: if self.cells[i].stalkDim > 0: for cf in self.cofaces(i): # Iterate over all cofaces of this activeCell try: @@ -906,7 +906,8 @@ def minimalExtend(self,assignment, activeCells=None, testSupport=None, method='n colidx+=self.cells[i].stalkDim # Use least squares to solve for global assignment given existing assignment - asg=self.serializeAssignment(assignment,activeCells=support) # Confusingly, activeSupport here refers *only* to the support of the assignment + asg,bnds=self.serializeAssignment(assignment,activeCells=support) # Confusingly, activeSupport here refers *only* to the support of the assignment + print mat result=np.linalg.lstsq(mat,asg) # Deserialize assignment From f4fc242e49817f9a855b82a2c3dbc38417d0fba1 Mon Sep 17 00:00:00 2001 From: Michael Robinson Date: Tue, 19 Jun 2018 15:07:34 -0400 Subject: [PATCH 4/4] now passes regression tests --- pysheaf/pysheaf.py | 32 +++++++++++++++----------------- 1 file changed, 15 insertions(+), 17 deletions(-) diff --git a/pysheaf/pysheaf.py b/pysheaf/pysheaf.py index b62d0fe..15c1019 100644 --- a/pysheaf/pysheaf.py +++ b/pysheaf/pysheaf.py @@ -880,11 +880,6 @@ def minimalExtend(self,assignment, activeCells=None, testSupport=None, method='n if ord != 2: warn('Kernel projection requires order 2 in minimalExtend') - # Build block matrix in which block columns are stalks over activeCells and block rows are stalks over cells supported in assignment - # Blocks are ultimate restriction maps (not just one hops!) - mat=np.zeros((sum([self.cells[i].stalkDim for i in support]), - sum([self.cells[i].stalkDim for i in ac]))) - # Compile dictionary of rows rowstarts=dict() rowidx=0 @@ -892,26 +887,29 @@ def minimalExtend(self,assignment, activeCells=None, testSupport=None, method='n rowstarts[i]=rowidx rowidx+=self.cells[i].stalkDim - colidx=0 # Current column + 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, - colidx:colidx+self.cells[i].stalkDim]=cf.restriction.matrix + mat[rowstarts[supportidx]:rowstarts[supportidx]+self.cells[supportidx].stalkDim,:]=cf.restriction.matrix except ValueError: pass - - colidx+=self.cells[i].stalkDim - # Use least squares to solve for global assignment given existing assignment - asg,bnds=self.serializeAssignment(assignment,activeCells=support) # Confusingly, activeSupport here refers *only* to the support of the assignment - print mat - result=np.linalg.lstsq(mat,asg) - - # Deserialize assignment - return self.deserializeAssignment(result[0],ac,assignment) + # 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),