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)

    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:
        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)

The above plots:

enter image description here

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.

  • $\begingroup$ Have you tried a brute force approach? That is, iterate pairwise through the entire list and check each distance? Note that $D=\sqrt{(x_i-x_j)^2+(y_i-y_j)^2+(z_i-z_j)^2}$ implies $D^2 = (x_i-x_j)^2+(y_i-y_j)^2+(z_i-z_j)^2$, so you can use $l^2 \geq (x_i-x_j)^2+(y_i-y_j)^2+(z_i-z_j)^2$ as your condition, without having to compute a square root each time. $\endgroup$
    – Connor Garcia
    Commented Apr 11, 2022 at 16:27
  • $\begingroup$ I imagine this question would get a better answer in a different forum since it's not specific to astronomy. You very likely want to look into spatial indexing algorithms like a quad-tree or R-Tree. $\endgroup$ Commented Apr 11, 2022 at 20:40

1 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.

  • 2
    $\begingroup$ I am discarding my very similar, partially written answer since I like yours better. If the axes dimensions of the OP's 3d plot are in Mpc, I would expect the adjacency matrix to be sparse (mostly zeroes). If you 3D bin the data into $l$-sized cubes, you can immediately assign zeroes in the adjacency matrix for any pair of stars in non-adjacent cubes (any cube index difference greater than 1). If your adjacency matrix was initialized to all ones, then, after this process, you only need to test the remaining ones. Of course, brute force is the best initial approach. $\endgroup$
    – Connor Garcia
    Commented Apr 11, 2022 at 20:48
  • 2
    $\begingroup$ @ConnorGarcia The cube assignment can be done in O(N) time, so this definitely beats my approach in time complexity (memory is another matter). Cool. $\endgroup$ Commented Apr 11, 2022 at 22:47
  • $\begingroup$ Thanks, I would say this binning is just a complement to your answer. The memory is also O(n) since each star has a single bin assignment. I am not sure why the OP selected but did not vote up this post. $\endgroup$
    – Connor Garcia
    Commented Apr 12, 2022 at 15:33

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .