Get the critical compositions along the tieline between two
        compositions. I.e. where the decomposition products change.
        The endpoints are also returned.
        Args:
            comp1, comp2 (Composition): compositions that define the tieline
        Returns:
            [(Composition)]: list of critical compositions. All are of
                the form x * comp1 + (1-x) * comp2

def get_critical_compositions(self, comp1, comp2):
        """
        Get the critical compositions along the tieline between two
        compositions. I.e. where the decomposition products change.
        The endpoints are also returned.
        Args:
            comp1, comp2 (Composition): compositions that define the tieline
        Returns:
            [(Composition)]: list of critical compositions. All are of
                the form x * comp1 + (1-x) * comp2
        """

        n1 = comp1.num_atoms
        n2 = comp2.num_atoms
        pd_els = self.elements

        # the reduced dimensionality Simplexes don't use the
        # first element in the PD
        c1 = self.pd_coords(comp1)
        c2 = self.pd_coords(comp2)

        # none of the projections work if c1 == c2, so just return *copies*
        # of the inputs
        if np.all(c1 == c2):
            return [comp1.copy(), comp2.copy()]

        intersections = [c1, c2]
        for sc in self.simplexes:
            intersections.extend(sc.line_intersection(c1, c2))
        intersections = np.array(intersections)

        # find position along line
        l = (c2 - c1)
        l /= np.sum(l ** 2) ** 0.5
        proj = np.dot(intersections - c1, l)

        # only take compositions between endpoints
        proj = proj[np.logical_and(proj > -self.numerical_tol,
                                   proj < proj[1] + self.numerical_tol)]
        proj.sort()

        # only unique compositions
        valid = np.ones(len(proj), dtype=np.bool)
        valid[1:] = proj[1:] > proj[:-1] + self.numerical_tol
        proj = proj[valid]

        ints = c1 + l * proj[:, None]
        # reconstruct full-dimensional composition array
        cs = np.concatenate([np.array([1 - np.sum(ints, axis=-1)]).T,
                             ints], axis=-1)
        # mixing fraction when compositions are normalized
        x = proj / np.dot(c2 - c1, l)
        # mixing fraction when compositions are not normalized
        x_unnormalized = x * n1 / (n2 + x * (n1 - n2))
        num_atoms = n1 + (n2 - n1) * x_unnormalized
        cs *= num_atoms[:, None]
        return [Composition((c, v) for c, v in zip(pd_els, m)) for m in cs]