Merge pull request #194 from grahamgower/cli

updates for release

CHANGELOG.md

@@ -1,9 +1,11 @@
 # Changelog
 
-## unreleased
+## [0.4.0] - 2023-05-03
 
-* Added `scale_bar` option to `demesdraw.tubes()` to draw a scale bar
-  that indicates the population size.
+* Documented the layout algorithm used in tubes plots.
+* Added `scale_bar` option to `demesdraw.tubes()`, and `--scale-bar` option
+  to the `demesdraw tubes` CLI, to draw a scale bar that indicates the
+  population size.
 * Removed dependency on `cvxpy` by going back to `scipy`
   for optimisation of the deme positions. The trust-constr
   method with linear constraints seems to work well.

demesdraw/__main__.py

@@ -66,10 +66,23 @@ class TubesCommand(Command):
     def __init__(self, subparsers):
         super().__init__(subparsers, "tubes")
 
+        self.parser.add_argument(
+            "--scale-bar",
+            action="store_true",
+            default=False,
+            help="Draw a scale bar that indicates population size.",
+        )
+
     def __call__(self, args):
         graph = demes.load(args.input_file)
         fig, ax = demesdraw.utils.get_fig_axes(aspect=args.aspect, scale=args.scale)
-        demesdraw.tubes(graph, ax=ax, log_time=args.log_time, title=args.title)
+        demesdraw.tubes(
+            graph,
+            ax=ax,
+            log_time=args.log_time,
+            title=args.title,
+            scale_bar=args.scale_bar,
+        )
         if args.output_file is not None:
             fig.savefig(args.output_file)
         else:

demesdraw/tubes.py

@@ -536,19 +536,18 @@ def random_migration_time(migration, log_scale: bool) -> float:
             ax.transData, ax.transAxes
         )
         x_offset = min(min(t.size1) for t in tubes.values())
-        # x_offset = max(max(t.size2) for t in tubes.values()) - width
-        ax_sb = ax.inset_axes([x_offset, -0.15, width, 0.01], transform=transform)
+        ax_sb = ax.inset_axes([x_offset, -0.16, width, 1e-6], transform=transform)
 
         ax_sb.yaxis.set_major_locator(matplotlib.ticker.NullLocator())
         ax_sb.spines[["left", "right", "top"]].set_visible(False)
-        ax_sb.xaxis.set_ticks_position("bottom")
+        ax_sb.xaxis.set_ticks_position("top")
         locator = matplotlib.ticker.AutoLocator()
         locator.set_params(integer=True)
         ax_sb.xaxis.set_major_locator(locator)
         ax_sb.xaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator())
         ax_sb.set_xlim(0, width)
         ax_sb.set_xlabel("deme size (individuals)")
-        ax_sb.xaxis.set_label_position("top")
+        ax_sb.xaxis.set_label_position("bottom")
 
     # Status bar text when in interactive mode.
     def format_coord(x, y):

demesdraw/utils.py

@@ -283,7 +283,7 @@ def _intersect(
     return not (dm_j.end_time >= dm_k.start_time or dm_k.end_time >= dm_j.start_time)
 
 
-def coexistence_indices(graph: demes.Graph) -> List[Tuple[int, int]]:
+def _coexistence_indices(graph: demes.Graph) -> List[Tuple[int, int]]:
     """Pairs of indices of demes that exist simultaneously."""
     contemporaries = []
     for j, deme_j in enumerate(graph.demes):
@@ -293,7 +293,7 @@ def coexistence_indices(graph: demes.Graph) -> List[Tuple[int, int]]:
     return contemporaries
 
 
-def successors_indices(graph: demes.Graph) -> Dict[int, List[int]]:
+def _successors_indices(graph: demes.Graph) -> Dict[int, List[int]]:
     """Graph successors, but use indices rather than deme names"""
     idx = {deme.name: j for j, deme in enumerate(graph.demes)}
     successors = dict()
@@ -303,7 +303,7 @@ def successors_indices(graph: demes.Graph) -> Dict[int, List[int]]:
     return successors
 
 
-def interactions_indices(graph: demes.Graph, *, unique: bool) -> List[Tuple[int, int]]:
+def _interactions_indices(graph: demes.Graph, *, unique: bool) -> List[Tuple[int, int]]:
     """Pairs of indices of demes that exchange migrants (migrations or pulses)."""
     idx = {deme.name: j for j, deme in enumerate(graph.demes)}
     interactions = []
@@ -530,9 +530,9 @@ def optimise_positions(
     :return:
         A dictionary mapping deme names to positions.
     """
-    contemporaries = coexistence_indices(graph)
-    successors = successors_indices(graph)
-    interactions = interactions_indices(graph, unique=unique_interactions)
+    contemporaries = _coexistence_indices(graph)
+    successors = _successors_indices(graph)
+    interactions = _interactions_indices(graph, unique=unique_interactions)
     if len(contemporaries) == 0:
         # There are no constraints, so stack demes on top of each other.
         return {deme.name: 0 for deme in graph.demes}
@@ -556,33 +556,31 @@ def optimise_positions(
         C.append(c)
     constraints = [scipy.optimize.LinearConstraint(C, lb=sep, ub=np.inf)]
 
-    with warnings.catch_warnings():
-        res = scipy.optimize.minimize(
-            _optimise_positions_objective,
-            x0,
-            args=(successors, interactions),
-            # The objective function returns a 2-tuple of (f_x, g_x), where f_x
-            # is the objective evalutated at x, and g_x is the Jacobian of f
-            # evaluated at x.
-            jac=True,
-            # Don't use the quasi-Newton Hessian approximation (the default for
-            # method="trust-constr" if nothing is specified). This performs
-            # poorly, and produces warnings, for some of the test cases.
-            # Writing the Hessian function manually would be tedious,
-            # and I'd certainly get it wrong. But since I did manually
-            # implement the Jacobian, scipy can use a finite-difference
-            # approximation for the Hessian, by specifying either "2-point",
-            # "3-point", or "cs". The "2-point" option seems to work fine.
-            hess="2-point",
-            method="trust-constr",
-            constraints=constraints,
-            bounds=scipy.optimize.Bounds(lb=np.min(x0) - sep, ub=np.max(x0) + sep),
-        )
+    res = scipy.optimize.minimize(
+        _optimise_positions_objective,
+        x0,
+        args=(successors, interactions),
+        # The objective function returns a 2-tuple of (f_x, g_x), where f_x
+        # is the objective evalutated at x, and g_x is the Jacobian of f
+        # evaluated at x.
+        jac=True,
+        # Don't use the quasi-Newton Hessian approximation (the default for
+        # method="trust-constr" if nothing is specified). This performs
+        # poorly, and produces warnings, for some of the test cases.
+        # Writing the Hessian function manually would be tedious,
+        # and I'd certainly get it wrong. But since I did manually
+        # implement the Jacobian, scipy can use a finite-difference
+        # approximation for the Hessian, by specifying either "2-point",
+        # "3-point", or "cs". The "2-point" option seems to work fine.
+        hess="2-point",
+        method="trust-constr",
+        constraints=constraints,
+        bounds=scipy.optimize.Bounds(lb=np.min(x0) - sep, ub=np.max(x0) + sep),
+    )
 
     if not res.success:
         warnings.warn(
             f"Failed to optimise: {res}\n\n"
-            "Demesdraw can probably do better. "
             "Please report this in the issue tracker at "
             "https://github.com/grahamgower/demesdraw/issues"
         )

tests/test_cli.py

@@ -73,6 +73,9 @@ def test_output_file_is_created(self, subcommand, params):
             assert output_file.exists()
             assert output_file.stat().st_size > 0
 
+    def test_tubes_scale_bar(self):
+        self.test_output_file_is_created("tubes", "--scale-bar")
+
     @pytest.mark.parametrize("subcommand", ["tubes", "size_history"])
     def test_input_from_stdin(self, subcommand):
         input_file = tests.example_files()[0]

tests/test_utils.py

@@ -137,7 +137,7 @@ def test_no_interations_isolated_demes(self, n, unique):
         for j in range(n):
             b.add_deme(f"deme{j}")
         graph = b.resolve()
-        interactions = utils.interactions_indices(graph, unique=unique)
+        interactions = utils._interactions_indices(graph, unique=unique)
         assert len(interactions) == 0
 
     @pytest.mark.parametrize("unique", (True, False))
@@ -147,7 +147,7 @@ def test_no_interactions_common_ancestor(self, unique):
         b.add_deme("b", ancestors=["a"])
         b.add_deme("c", ancestors=["a"])
         graph = b.resolve()
-        interactions = utils.interactions_indices(graph, unique=unique)
+        interactions = utils._interactions_indices(graph, unique=unique)
         assert len(interactions) == 0
 
     @pytest.mark.parametrize("n_migrations", (1, 2, 3))
@@ -166,12 +166,12 @@ def test_asymmetric_migrations(self, n_migrations):
                 rate=1e-5,
             )
         graph = b.resolve()
-        interactions = utils.interactions_indices(graph, unique=False)
+        interactions = utils._interactions_indices(graph, unique=False)
         assert (
             interactions == [(1, 2)] * n_migrations
             or interactions == [(2, 1)] * n_migrations
         )
-        interactions = utils.interactions_indices(graph, unique=True)
+        interactions = utils._interactions_indices(graph, unique=True)
         assert interactions == [(1, 2)] or interactions == [(2, 1)]
 
     @pytest.mark.parametrize("n_migrations", (1, 2, 3))
@@ -186,13 +186,13 @@ def test_symmetric_migrations(self, n_migrations):
                 demes=["b", "c"], start_time=start_time, end_time=end_time, rate=1e-5
             )
         graph = b.resolve()
-        interactions = utils.interactions_indices(graph, unique=False)
+        interactions = utils._interactions_indices(graph, unique=False)
         assert len(interactions) == 2 * n_migrations
         counts = collections.Counter(interactions)
         assert len(counts) == 2
         assert counts[(1, 2)] == n_migrations
         assert counts[(2, 1)] == n_migrations
-        interactions = utils.interactions_indices(graph, unique=True)
+        interactions = utils._interactions_indices(graph, unique=True)
         assert interactions == [(1, 2)] or interactions == [(2, 1)]
 
     @pytest.mark.parametrize("n_pulses", (1, 2, 3))
@@ -206,11 +206,11 @@ def test_pulses(self, n_pulses):
                 sources=["b"], dest="c", time=100 * j / n_pulses + 1, proportions=[0.1]
             )
         graph = b.resolve()
-        interactions = utils.interactions_indices(graph, unique=False)
+        interactions = utils._interactions_indices(graph, unique=False)
         assert (
             interactions == [(1, 2)] * n_pulses or interactions == [(2, 1)] * n_pulses
         )
-        interactions = utils.interactions_indices(graph, unique=True)
+        interactions = utils._interactions_indices(graph, unique=True)
         assert interactions == [(1, 2)] or interactions == [(2, 1)]
 
 
@@ -743,7 +743,7 @@ def test_n_immortal_demes(self, n):
         for j in range(n):
             b.add_deme(f"deme{j}")
         graph = b.resolve()
-        idx = utils.coexistence_indices(graph)
+        idx = utils._coexistence_indices(graph)
         assert set(idx) == set(itertools.combinations(range(n), 2))
 
     @pytest.mark.parametrize("n", (1, 2, 3))
@@ -756,7 +756,7 @@ def test_succesive_descendents(self, n):
                 ancestors = [f"deme{j-1}"]
             b.add_deme(f"deme{j}", ancestors=ancestors, epochs=[dict(end_time=100 - j)])
         graph = b.resolve()
-        idx = utils.coexistence_indices(graph)
+        idx = utils._coexistence_indices(graph)
         assert len(idx) == 0
 
     @pytest.mark.parametrize("unique", (True, False))
@@ -771,7 +771,7 @@ def test_tree_like1(self, unique):
         b.add_deme("b", ancestors=["a"], start_time=100)
         b.add_deme("c", ancestors=["b"], start_time=50)
         graph = b.resolve()
-        idx = utils.coexistence_indices(graph)
+        idx = utils._coexistence_indices(graph)
         assert len(idx) == 3
         assert (0, 1) in idx
         assert (0, 2) in idx
@@ -791,7 +791,7 @@ def test_tree_like2(self, unique):
         b.add_deme("c", ancestors=["b"], start_time=50)
         b.add_deme("d", ancestors=["c"], start_time=10)
         graph = b.resolve()
-        idx = utils.coexistence_indices(graph)
+        idx = utils._coexistence_indices(graph)
         assert len(idx) == 3
         assert (1, 2) in idx
         assert (1, 3) in idx
@@ -809,7 +809,7 @@ def test_ordering_doesnt_change(self, graph, unique):
         positions2 = utils.optimise_positions(
             graph, positions=positions1, sep=sep, unique_interactions=unique
         )
-        for j, k in utils.coexistence_indices(graph):
+        for j, k in utils._coexistence_indices(graph):
             if positions1[graph.demes[j].name] < positions1[graph.demes[j].name]:
                 assert positions2[graph.demes[j].name] < positions2[graph.demes[j].name]
             elif positions1[graph.demes[j].name] > positions1[graph.demes[j].name]:
@@ -826,7 +826,7 @@ def test_separation_distance(self, graph, unique, sep):
             graph, positions=positions1, sep=sep, unique_interactions=unique
         )
         epsilon = 1e-3  # Small amount of wiggle room for numerical error.
-        for j, k in utils.coexistence_indices(graph):
+        for j, k in utils._coexistence_indices(graph):
             assert (
                 abs(positions2[graph.demes[j].name] - positions2[graph.demes[k].name])
                 >= sep - epsilon
@@ -851,7 +851,7 @@ def test_island_model(self, n, unique, migrations):
             graph, positions=positions1, sep=sep, unique_interactions=unique
         )
         epsilon = 1e-3  # Small amount of wiggle room for numerical error.
-        for j, k in utils.coexistence_indices(graph):
+        for j, k in utils._coexistence_indices(graph):
             assert (
                 abs(positions2[graph.demes[j].name] - positions2[graph.demes[k].name])
                 >= sep - epsilon
@@ -869,8 +869,8 @@ def test_optimise_gradients(self, graph, unique):
         positions = utils.minimal_crossing_positions(
             graph, sep=sep, unique_interactions=unique
         )
-        successors = utils.successors_indices(graph)
-        interactions = utils.interactions_indices(graph, unique=unique)
+        successors = utils._successors_indices(graph)
+        interactions = utils._interactions_indices(graph, unique=unique)
 
         x0 = np.array([positions[deme.name] for deme in graph.demes])
         # Place the first deme at position 0.