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

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)