Arf layout (#5910)

* added arf_layout

* reference to docstring and comparison to spring layout

* rebase to origin main

* black re-format

* Left aligned docstring text

* Cleaned up computation and update variables to new docstring

* Updated naming tests. Added input check on arf_layout parameter `a`

* Fixed Linter issues for py38 target

* Fixed Linter issues for target p38

* linter issue fixed

2 files changed, 139 insertions, 1 deletions

diff --git a/networkx/drawing/layout.py b/networkx/drawing/layout.py
index b6d2afe5..a4029942 100644
--- a/networkx/drawing/layout.py
+++ b/networkx/drawing/layout.py
@@ -32,6 +32,7 @@ __all__ = [
     "fruchterman_reingold_layout",
     "spiral_layout",
     "multipartite_layout",
+    "arf_layout",
 ]
 
 
@@ -1110,6 +1111,118 @@ def multipartite_layout(G, subset_key="subset", align="vertical", scale=1, cente
     return pos
 
 
+def arf_layout(
+    G,
+    pos=None,
+    scaling=1,
+    a=1.1,
+    etol=1e-6,
+    dt=1e-3,
+    max_iter=1000,
+):
+    """Arf layout for networkx
+
+    The attractive and repulsive forces (arf) layout [1]
+    improves the spring layout in three ways. First, it
+    prevents congestion of highly connected nodes due to
+    strong forcing between nodes. Second, it utilizes the
+    layout space more effectively by preventing large gaps
+    that spring layout tends to create. Lastly, the arf
+    layout represents symmmetries in the layout better than
+    the default spring layout.
+
+    Parameters
+    ----------
+    G : nx.Graph or nx.DiGraph
+        Networkx graph.
+    pos : dict
+        Initial  position of  the nodes.  If set  to None  a
+        random layout will be used.
+    scaling : float
+        Scales the radius of the circular layout space.
+    a : float
+        Strength of springs between connected nodes. Should be larger than 1. The greater a, the clearer the separation ofunconnected sub clusters.
+    etol : float
+        Graduent sum of spring forces must be larger than `etol` before succesful termination.
+    dt : float
+        Time step for force differential equation simulations.
+    max_iter : int
+        Max iterations before termination of the algorithm.
+
+    References
+    .. [1] "Self-Organization Applied to Dynamic Network Layout", M. Geipel,
+            International Jounral of Modern Physics C, 2007, Vol 18, No 10, pp. 1537-1549.
+            https://doi.org/10.1142/S0129183107011558 https://arxiv.org/abs/0704.1748
+
+    Returns
+    -------
+    pos : dict
+        A dictionary of positions keyed by node.
+
+    Examples
+    --------
+    >>> G = nx.grid_graph((5, 5))
+    >>> pos = nx.arf_layout(G)
+
+    """
+    import warnings
+
+    import numpy as np
+
+    if a <= 1:
+        msg = "The parameter a should be larger than 1"
+        raise ValueError(msg)
+
+    pos_tmp = nx.random_layout(G)
+    if pos is None:
+        pos = pos_tmp
+    else:
+        for node in G.nodes():
+            if node not in pos:
+                pos[node] = pos_tmp[node].copy()
+
+    # Initialize spring constant matrix
+    N = len(G)
+    # No nodes no computation
+    if N == 0:
+        return pos
+
+    # init force of springs
+    K = np.ones((N, N)) - np.eye(N)
+    node_order = {node: i for i, node in enumerate(G)}
+    for x, y in G.edges():
+        if x != y:
+            idx, jdx = (node_order[i] for i in (x, y))
+            K[idx, jdx] = a
+
+    # vectorize values
+    p = np.asarray(list(pos.values()))
+
+    # equation 10 in [1]
+    rho = scaling * np.sqrt(N)
+
+    # looping variables
+    error = etol + 1
+    n_iter = 0
+    while error > etol:
+        diff = p[:, np.newaxis] - p[np.newaxis]
+        A = np.linalg.norm(diff, axis=-1)[..., np.newaxis]
+        # attraction_force - repulsions force
+        # suppress nans due to division; caused by diagonal set to zero.
+        # Does not affect the computation due to nansum
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore")
+            change = K[..., np.newaxis] * diff - rho / A * diff
+        change = np.nansum(change, axis=0)
+        p += change * dt
+
+        error = np.linalg.norm(change, axis=-1).sum()
+        if n_iter > max_iter:
+            break
+        n_iter += 1
+    return {node: pi for node, pi in zip(G.nodes(), p)}
+
+
 def rescale_layout(pos, scale=1):
     """Returns scaled position array to (-scale, scale) in all axes.
 
diff --git a/networkx/drawing/tests/test_layout.py b/networkx/drawing/tests/test_layout.py
index f24d0038..faf734cf 100644
--- a/networkx/drawing/tests/test_layout.py
+++ b/networkx/drawing/tests/test_layout.py
@@ -66,6 +66,7 @@ class TestLayout:
         nx.kamada_kawai_layout(G)
         nx.kamada_kawai_layout(G, dim=1)
         nx.kamada_kawai_layout(G, dim=3)
+        nx.arf_layout(G)
 
     def test_smoke_string(self):
         G = self.Gs
@@ -80,6 +81,7 @@ class TestLayout:
         nx.kamada_kawai_layout(G)
         nx.kamada_kawai_layout(G, dim=1)
         nx.kamada_kawai_layout(G, dim=3)
+        nx.arf_layout(G)
 
     def check_scale_and_center(self, pos, scale, center):
         center = np.array(center)
@@ -175,6 +177,10 @@ class TestLayout:
         pos = nx.circular_layout(self.Gi)
         npos = nx.fruchterman_reingold_layout(self.Gi, pos=pos)
 
+    def test_smoke_initial_pos_arf(self):
+        pos = nx.circular_layout(self.Gi)
+        npos = nx.arf_layout(self.Gi, pos=pos)
+
     def test_fixed_node_fruchterman_reingold(self):
         # Dense version (numpy based)
         pos = nx.circular_layout(self.Gi)
@@ -242,6 +248,8 @@ class TestLayout:
         assert vpos == {}
         vpos = nx.kamada_kawai_layout(G, center=(1, 1))
         assert vpos == {}
+        vpos = nx.arf_layout(G)
+        assert vpos == {}
 
     def test_bipartite_layout(self):
         G = nx.complete_bipartite_graph(3, 5)
@@ -402,7 +410,6 @@ class TestLayout:
         for k, v in expectation.items():
             assert (s_vpos[k] == v).all()
         s_vpos = nx.rescale_layout_dict(vpos, scale=2)
-
         expectation = {
             0: np.array((-2, -2)),
             1: np.array((2, 2)),
@@ -411,6 +418,24 @@ class TestLayout:
         for k, v in expectation.items():
             assert (s_vpos[k] == v).all()
 
+    def test_arf_layout_partial_input_test(self):
+        """
+        Checks whether partial pos input still returns a proper position.
+        """
+        G = self.Gs
+        node = nx.utils.arbitrary_element(G)
+        pos = nx.circular_layout(G)
+        del pos[node]
+        pos = nx.arf_layout(G, pos=pos)
+        assert len(pos) == len(G)
+
+    def test_arf_layout_negative_a_check(self):
+        """
+        Checks input parameters correctly raises errors. For example,  `a` should be larger than 1
+        """
+        G = self.Gs
+        pytest.raises(ValueError, nx.arf_layout, G=G, a=-1)
+
 
 def test_multipartite_layout_nonnumeric_partition_labels():
     """See gh-5123."""