Speeding up `CellComplex.construct()`

[chenzhaiyu]

The current construct() is slow. Polyhedron.intersection(other) takes most of the time, as profiled with _lineprofiler:

Timer unit: 1e-06 s

Total time: 27.422 s
File: /Users/zhaiyu/opt/anaconda3/envs/abspy/lib/python3.8/site-packages/abspy/complex.py
Function: construct at line 390

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
   390                                               @profile
   391                                               def construct(self, exhaustive=False):
   392                                                   """
   393                                                   Construct cell complex.
   394                                           
   395                                                   Two-stage primitive-in-cell predicate.
   396                                                   (1) bounding boxes of primitive and existing cells are evaluated
   397                                                   for possible intersection. (2), a strict intersection test is performed.
   398                                           
   399                                                   Generated cells are stored in self.cells.
   400                                                   * query the bounding box intersection.
   401                                                   * optional: intersection test for polygon and edge in each potential cell.
   402                                                   * partition the potential cell into two. rewind if partition fails.
   403                                           
   404                                                   Parameters
   405                                                   ----------
   406                                                   exhaustive: bool
   407                                                       Do exhaustive partitioning if set True
   408                                                   """
   409         1        266.0    266.0      0.0          logger.info('constructing cell complex')
   410         1          4.0      4.0      0.0          tik = time.time()
   411                                           
   412         1      13862.0  13862.0      0.1          pbar = range(len(self.bounds)) if self.quiet else trange(len(self.bounds))
   413        42      16878.0    401.9      0.1          for i in pbar:  # kinetic for each primitive
   414                                                       # bounding box intersection test
   415                                                       # indices of existing cells with potential intersections
   416        41     292416.0   7132.1      1.1              indices_cells = self._intersect(self.bounds[i], self.planes[i], exhaustive)
   417        41        134.0      3.3      0.0              assert len(indices_cells), 'intersection failed! check the initial bound'
   418                                           
   419                                                       # half-spaces defined by inequalities
   420                                                       # no change_ring() here (instead, QQ() in _inequalities) speeds up 10x
   421                                                       # init before the loop could possibly speed up a bit
   422        82      56790.0    692.6      0.2              hspace_positive, hspace_negative = [Polyhedron(ieqs=[inequality]) for inequality in
   423        41      35014.0    854.0      0.1                                                  self._inequalities(self.planes[i])]
   424                                           
   425                                                       # partition the intersected cells and their bounds while doing mesh slice plane
   426        41        153.0      3.7      0.0              indices_parents = []
   427                                           
   428      1013       3274.0      3.2      0.0              for index_cell in indices_cells:
   429       972    1727068.0   1776.8      6.3                  cell_positive = hspace_positive.intersection(self.cells[index_cell])
   430       972    1697056.0   1745.9      6.2                  cell_negative = hspace_negative.intersection(self.cells[index_cell])
   431                                           
   432       972      53076.0     54.6      0.2                  if cell_positive.dim() != 3 or cell_negative.dim() != 3:
   433                                                               # if cell_positive.is_empty() or cell_negative.is_empty():
   434                                                               """
   435                                                               cannot use is_empty() predicate for degenerate cases:
   436                                                                   sage: Polyhedron(vertices=[[0, 1, 2]])
   437                                                                   A 0-dimensional polyhedron in ZZ^3 defined as the convex hull of 1 vertex
   438                                                                   sage: Polyhedron(vertices=[[0, 1, 2]]).is_empty()
   439                                                                   False
   440                                                               """
   441       105        113.0      1.1      0.0                      continue
   442                                           
   443                                                           # incrementally build the adjacency graph
   444       795       1518.0      1.9      0.0                  if self.graph is not None:
   445                                                               # append the two nodes (UID) being partitioned
   446       795       7322.0      9.2      0.0                      self.graph.add_node(self.index_node + 1)
   447       795       3261.0      4.1      0.0                      self.graph.add_node(self.index_node + 2)
   448                                           
   449                                                               # append the edge in between
   450       795       5137.0      6.5      0.0                      self.graph.add_edge(self.index_node + 1, self.index_node + 2)
   451                                           
   452                                                               # get neighbours of the current cell from the graph
   453       795      30508.0     38.4      0.1                      neighbours = self.graph[list(self.graph.nodes)[index_cell]]  # index in the node list
   454                                           
   455       795       1630.0      2.1      0.0                      if neighbours:
   456                                                                   # get the neighbouring cells to the parent
   457       794     119186.0    150.1      0.4                          cells_neighbours = [self.cells[self._index_node_to_cell(n)] for n in neighbours]
   458                                           
   459                                                                   # adjacency test between both created cells and their neighbours
   460                                                                   # todo:
   461                                                                   #   Avoid 3d-3d intersection if possible. Unsliced neighbours connect with only one child
   462                                                                   #   - reduce computation by half - can be further reduced using vertices/faces instead of
   463                                                                   #   polyhedron intersection. Sliced neighbors connect with both children
   464                                           
   465      8090      13708.0      1.7      0.0                          for n, cell in enumerate(cells_neighbours):
   466                                           
   467      7296   13415562.0   1838.8     48.9                              interface_positive = cell_positive.intersection(cell)
   468                                           
   469      7296     218415.0     29.9      0.8                              if interface_positive.dim() == 2:
   470                                                                           # this neighbour can connect with either or both children
   471      5100      69081.0     13.5      0.3                                  self.graph.add_edge(self.index_node + 1, list(neighbours)[n])
   472      5100    9149388.0   1794.0     33.4                                  interface_negative = cell_negative.intersection(cell)
   473      5100     144452.0     28.3      0.5                                  if interface_negative.dim() == 2:
   474      2894      39336.0     13.6      0.1                                      self.graph.add_edge(self.index_node + 2, list(neighbours)[n])
   475                                                                       else:
   476                                                                           # this neighbour must otherwise connect with the other child
   477      2196      30033.0     13.7      0.1                                  self.graph.add_edge(self.index_node + 2, list(neighbours)[n])
   478                                           
   479                                                               # update cell id
   480       795       1257.0      1.6      0.0                      self.index_node += 2
   481                                           
   482       795       1197.0      1.5      0.0                  self.cells.append(cell_positive)
   483       795        849.0      1.1      0.0                  self.cells.append(cell_negative)
   484                                           
   485                                                           # incrementally cache the bounds for created cells
   486       795     137263.0    172.7      0.5                  self.cells_bounds.append(cell_positive.bounding_box())
   487       795     114524.0    144.1      0.4                  self.cells_bounds.append(cell_negative.bounding_box())
   488                                           
   489       795       1246.0      1.6      0.0                  indices_parents.append(index_cell)
   490                                           
   491                                                       # delete the parent cells and their bounds. this does not affect the appended ones
   492       836       1132.0      1.4      0.0              for index_parent in sorted(indices_parents, reverse=True):
   493       795       1029.0      1.3      0.0                  del self.cells[index_parent]
   494       795        959.0      1.2      0.0                  del self.cells_bounds[index_parent]
   495                                           
   496                                                           # remove the parent node (and subsequently its incident edges) in the graph
   497       795        807.0      1.0      0.0                  if self.graph is not None:
   498       795      16774.0     21.1      0.1                      self.graph.remove_node(list(self.graph.nodes)[index_parent])
   499                                           
   500         1          3.0      3.0      0.0          self.constructed = True
   501         1        305.0    305.0      0.0          logger.info('cell complex constructed: {:.2f} s'.format(time.time() - tik))

[chenzhaiyu]

One way is to avoid 3D-3D intersections wherever possible. Another is to introduce multi-processing/threading into the workflow (e.g., when looking at which neighbors to link to).

[chenzhaiyu]

Attempted to avoid 3D-3D intersection by introducing an unsliced flag for neighboring cell, in which case the neighboring cell can only connect to one child of the partition. This actually made it even slower, possibly due to the overhead of 2D-3D intersection for predicting the flag.

[chenzhaiyu]

Partly addressed in #20 by introducing OBB to reduce 3D-3D intersections.