From 52592bb84ed19d745f00e95b018b04d85eec4104 Mon Sep 17 00:00:00 2001
From: Skylar Chan <54919210+schance995@users.noreply.github.com>
Date: Tue, 5 Aug 2025 21:11:25 -0400
Subject: [PATCH] fix: choose partition centers and sizes randomly.

This commit implements the following update from
bct (https://sites.google.com/site/bctnet/all-help-headers):

> Dec 2016: Updated code so that both partition centers and partition
> sizes are chosen at random. Also added in a constraint on partition
> placement that prevents boxes from being located outside the edges of
> the network. This helps prevent skewed results due to boundary effects
> arising from the finite size of the network. (Lia Papadopoulos)
---
 bct/algorithms/physical_connectivity.py | 14 ++++++++++++--
 1 file changed, 12 insertions(+), 2 deletions(-)

diff --git a/bct/algorithms/physical_connectivity.py b/bct/algorithms/physical_connectivity.py
index 55cf7cf..afcf319 100644
--- a/bct/algorithms/physical_connectivity.py
+++ b/bct/algorithms/physical_connectivity.py
@@ -141,7 +141,17 @@ def rentian_scaling(A, xyz, n, seed=None):
     # and the number of edges traversing the boundary of the partition (e)
     while count < n:
         # define cube endpoints
-        randx = np.sort((1 + nmax - nmin) * rng.random_sample((2,)))
+        randx = np.sort((nmax - nmin) * rng.random((2,)))
+        randy = np.sort((nmax - nmin) * rng.random((2,)))
+        randz = np.sort((nmax - nmin) * rng.random((2,)))
+
+        # ensure cube is within network boundaries
+        randx[0] = max(randx[0], nmin)
+        randx[1] = min(randx[1], nmax)
+        randy[0] = max(randy[0], nmin)
+        randy[1] = min(randy[1], nmax)
+        randz[0] = max(randz[0], nmin)
+        randz[1] = min(randz[1], nmax)
 
         # find nodes in cube
         l1 = xyzn[:, 0] > randx[0]
@@ -151,7 +161,7 @@ def rentian_scaling(A, xyz, n, seed=None):
         l5 = xyzn[:, 2] > randx[0]
         l6 = xyzn[:, 2] < randx[1]
 
-        L, = np.where((l1 & l2 & l3 & l4 & l5 & l6).flatten())
+        (L,) = np.where((l1 & l2 & l3 & l4 & l5 & l6).flatten())
         if np.size(L):
             # count edges crossing at the boundary of the cube
             E[count] = np.sum(A[np.ix_(L, np.setdiff1d(range(m), L))])