How to build an adjacency matrix of stars connectivity from their physical kpc distances?

Question

I'm probing into the Illustris API, and gathering information from a specific cosmos simulation, for a given redshift value.

This is how I request the api:

import requests

baseUrl = 'http://www.tng-project.org/api/'
    
def get(path, params=None):
    # make HTTP GET request to path
    headers = {"api-key":"my_key"}
    r = requests.get(path, params=params, headers=headers)

    # raise exception if response code is not HTTP SUCCESS (200)
    r.raise_for_status()

    if r.headers['content-type'] == 'application/json':
        return r.json() # parse json responses automatically
    
    if 'content-disposition' in r.headers:
        filename = r.headers['content-disposition'].split("filename=")[1]
        with open(f'sky_dataset/simulations/{filename}', 'wb') as f:
            f.write(r.content)
        return filename # return the filename string
    return r

And this is how I get star coordinates for a given subhalo in this particular simulation. Note that -if I'm doing it right- distances have already been converted from ckpc/h to physical kpc:

import h5py
import numpy as np

simulation_id = 100
redshift = 0.57
subhalo_id = 99

scale_factor = 1.0 / (1+redshift)
little_h = 0.704

params = {'stars':'Coordinates,GFM_Metallicity'}

url = "http://www.tng-project.org/api/Illustris-1/snapshots/z=" + str(redshift) + "/subhalos/" + str(subhalo_id)
sub = get(url) # get json response of subhalo properties
saved_filename = get(url + "/cutout.hdf5",params) # get and save HDF5 cutout file

with h5py.File(f'sky_dataset/simulations/{saved_filename}') as f:
    # NOTE! If the subhalo is near the edge of the box, you must take the periodic boundary into account! (we ignore it here)
    dx = f['PartType4']['Coordinates'][:,0] - sub['pos_x']
    dy = f['PartType4']['Coordinates'][:,1] - sub['pos_y']
    dz = f['PartType4']['Coordinates'][:,2] - sub['pos_z']
    
    rr = np.sqrt(dx**2 + dy**2 + dz**2)
    rr *= scale_factor/little_h # ckpc/h -> physical kpc

    fig = plt.figure(figsize=(12,12))
    with mpl.rc_context(rc={'axes3d.grid': True}):
        ax = fig.add_subplot(projection='3d')

        # Plot the values
        ax.scatter(dx, dy, dz)
        ax.set_xlabel('X-axis')
        ax.set_ylabel('Y-axis')
        ax.set_zlabel('Z-axis')
    plt.show()

The above plots:

---

My aim is to build a connectivity network for this system, starting with an square (simetrical) adjacency matrix, whereby any two stars (or vertices) are connected if they lie within the linking length l of 1.2 Mpc, that is:

Aij = 1 if rij ≤ l, otherwise 0

where rij is the distance between the two vertices, i and j.

---

Any ideas on how to get this adjacency matrix, based on my linking length? Any help would be greatly appreciated.

Answer 1

Computationally, a brute force algorithm where you calculate $r_{ij}$ for each pair $(i,j)$ and set $A_{ij}$ if it is smaller than the cut-off distance runs in $O(N^2)$ time. While this was much too slow when I first did this for $\sim 700$ nearby stars back in the 1980s on a home computer, today even if you have a few thousand halos it might be acceptably fast since you likely only calculate $A_{ij}$ once.

However, there is a faster way of doing it. $$r_{ij}^2=(x_i-x_j)^2+(y_i-y_j)^2+(z_i-z_j)^2.$$ Note that if $(x_i-x_j)^2>l^2$, that is $|x_i-x_j|>l$, then the distance must be larger than $l$, so $A_{ij}=0$. You can hence start by calculating $\Delta x =x_i-x_j$, and if $|\Delta x|<l$ continue to calculate $\Delta y$, otherwise go to the next pair. If $|\Delta y|<l$ too, now you calculate $\Delta x^2+\Delta y^2 +\Delta z^2$ and check if it is $<l^2$: only then is $A_{ij}=1$ (and $A_{ji}=1$ too, no need to calculate that separately). This algorithm is technically also $O(N^2)$, but as long as $l\ll$ the diameter of the data set only the first comparison is needed for most points, making it really fast.

A trick is to sort the points by $x$ (takes $O(N \log N)$ time), and loop over $j$ from your starting point $i$ upwards and downwards: as soon as the $\Delta x$ test fails, you know you will not get any other neighbours for bigger/smaller $j$ and can increment $i$. This version is even faster.

Still, in a modern math environment parallelising the operations may make the brute force approach fast enough. When I write something like this in Matlab I typically run a loop over $i$, and then in parallel generate a vector $\Delta \mathbf{r}$ of the distances from point $i$ to all the others. Then I use the parallelised comparison operation to set the relevant elements of $\mathbf{A}$ to 1 (for big $N$ this has the advantage that I only keep a length $N$ vector in memory rather than a $N^2$ matrix; I also use a sparse matrix representation for $\mathbf{A}$). Generally in Matlab at least, parallelised operations are much, much faster than nested loops.

(If you want to go down deeper rabbit holes one can consider programming this for a GPU to really make it fly, but this only makes sense if you can hold your entire dataset in the GPU - the time it takes to send things there and back to memory can easily eat up the parallelism in the GPU.)

But remember, if you only calculate $\mathbf{A}$ once and save it, a running time of ten minutes may be totally fine. Optimizing algorithms is fun, but not always a good use of time.