docs: def ils(self, node, sorting_times=None, force_topology_change=True):
          """
          A constrained and approximation of ILS using nearest-neighbour interchange
  
          Process
          -------
          A node with at least three descendents is selected from an ultrametric tree
          (node '2', below)
  
                     ---0--...              ---0--...           ---0--...
                    |        |             |        |        --1--      |
                    |        R           --1--      R       |     |     R
         age        |                   |     |            -2-    |
          ^         |                   |     |           |   |   |
          |       --1--                -2-    |           |   |   |
          |      |     |       or     |   |   |     or    |   |   |
          |      |     |              |   |   |           |   |   |
          |     -2-    |              |   |   |           |   |   |
          |    |   |   |              |   |   |           |   |   |
          |    A   B   C              C   B   A           A   C   B
  
          Nodes 'A', 'B' and 'C' are rearranged into one of the three configurations
          [(A, B), C], [A, (B, C)], [(A, C), B]
  
          Nodes 1 and 2 are slid further up the tree, but no further than node 0
          (this is why it's a constrained version), by an amount drawn from a
          truncated exponential distribution.
  
          This is approximately corresponds to the case where A and B failed to
          coalesce in the branch 1->2, so they coalesce with C in the branch
          0 -> 1 instead
          """
          # node = '2', par = '1', gpar = '0' -- in above diagram
          n_2 = node
          n_1 = n_2.parent_node
          if n_1 == self.tree._tree.seed_node:
              logger.warn('Node 1 is the root - calling again on child')
              self.ils(n_2.child_nodes())
          n_0 = n_1.parent_node
          a, b = node.child_nodes()
          c, = node.sister_nodes()
  
          ages = [a.age, b.age, c.age, n_2.age, n_1.age, n_0.age]
  
          # Do topology changes
          if force_topology_change:
              swap_mode = random.choice([1, 2])
          else:
              swap_mode = random.choice([0, 1, 2])
  
          if swap_mode == 1:
              # Exchange 'a' and 'c'
              n_2.remove_child(a)
              n_1.remove_child(c)
              n_2.add_child(c)
              n_1.add_child(a)
  
          elif swap_mode == 2:
              # Exchange 'b' and 'c'
              n_2.remove_child(b)
              n_1.remove_child(c)
              n_2.add_child(c)
              n_1.add_child(b)
  
          # Do branch length adjustments
          # Bounds - between node 0 (upper) and node 1 (lower)
          min_unsorted_age = n_1.age
          max_unsorted_age = n_0.age
          if sorting_times is None:
              sorting_times = truncated_exponential(max_unsorted_age-min_unsorted_age,
                                                scale=0.1*(max_unsorted_age-min_unsorted_age),
                                                sample_size=2) # E(t) = n(n-1)/2, n = 3
              sorting_times += min_unsorted_age
              sorting_times = np.array([min_unsorted_age, ages[3]])
  
          # Adjust node 1 edge length
          new_n1_age = max(sorting_times)
          prev_age = ages[4]
          slide = (new_n1_age - prev_age)
          if slide < 1e-6:
              slide = 0
              new_n1_age = prev_age
          n_1.edge.length -= slide
          n_2.edge.length += slide
  
          # Adjust node 2 edge length
          new_n2_age = min(sorting_times)
          prev_age = ages[3]
          slide = (new_n2_age - prev_age)
          if slide < 1e-6:
              slide = 0
              new_n2_age = prev_age
          n_2.edge.length -= slide
  
          # Adjust a, b and c edge lengths
          if swap_mode == 0:
              a.edge.length = (new_n2_age - ages[0])
              b.edge.length = (new_n2_age - ages[1])
              c.edge.length = (new_n1_age - ages[2])
  
          elif swap_mode == 1:
              a.edge.length = (new_n1_age - ages[0])
              b.edge.length = (new_n2_age - ages[1])
              c.edge.length = (new_n2_age - ages[2])
  
          else:
              a.edge.length = (new_n2_age - ages[0])
              b.edge.length = (new_n1_age - ages[1])
              c.edge.length = (new_n2_age - ages[2])
  
          # used to be .reindex_taxa() before dendropy 4.
          # migrate_taxon_namespace is recommended migrated function,
          # but not sure if its even needed anymore.
          self.tree._tree.migrate_taxon_namespace(self.tree._tree.taxon_namespace)
  
          self.tree._tree.encode_bipartitions()
          self._validate()
          logger.debug(self.tree)