<p>Computationally, a brute force algorithm where you calculate <span class="math-container">$r_{ij}$</span> for each pair <span class="math-container">$(i,j)$</span> and set <span class="math-container">$A_{ij}$</span> if it is smaller than the cut-off distance runs in <span class="math-container">$O(N^2)$</span> time. While this was much too slow when I first did this for <span class="math-container">$\sim 700$</span> 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 <span class="math-container">$A_{ij}$</span> once.</p>
<p>However, there is a faster way of doing it. <span class="math-container">$$r_{ij}^2=(x_i-x_j)^2+(y_i-y_j)^2+(z_i-z_j)^2.$$</span>
Note that if <span class="math-container">$(x_i-x_j)^2&gt;l^2$</span>, that is <span class="math-container">$|x_i-x_j|&gt;l$</span>, then the distance must be larger than <span class="math-container">$l$</span>, so <span class="math-container">$A_{ij}=0$</span>. You can hence start by calculating <span class="math-container">$\Delta x =x_i-x_j$</span>, and if <span class="math-container">$|\Delta x|&lt;l$</span> continue to calculate <span class="math-container">$\Delta y$</span>, otherwise go to the next pair. If <span class="math-container">$|\Delta y|&lt;l$</span> too, now you calculate <span class="math-container">$\Delta x^2+\Delta y^2 +\Delta z^2$</span> and check if it is <span class="math-container">$&lt;l^2$</span>: only then is <span class="math-container">$A_{ij}=1$</span> (<em>and</em> <span class="math-container">$A_{ji}=1$</span> too, no need to calculate that separately). This algorithm is technically also <span class="math-container">$O(N^2)$</span>, but as long as <span class="math-container">$l\ll$</span> the diameter of the data set only the first comparison is needed for most points, making it really fast.</p>
<p>A trick is to sort the points by <span class="math-container">$x$</span> (takes <span class="math-container">$O(N \log N)$</span> time), and loop over <span class="math-container">$j$</span> from your starting point <span class="math-container">$i$</span> upwards and downwards: as soon as the <span class="math-container">$\Delta x$</span> test fails, you know you will not get any other neighbours for bigger/smaller <span class="math-container">$j$</span> and can increment <span class="math-container">$i$</span>. This version is even faster.</p>
<p>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 <span class="math-container">$i$</span>, and then in parallel generate a vector <span class="math-container">$\Delta \mathbf{r}$</span> of the distances from point <span class="math-container">$i$</span> to all the others. Then I use the parallelised comparison operation to set the relevant elements of <span class="math-container">$\mathbf{A}$</span> to 1 (for big <span class="math-container">$N$</span> this has the advantage that I only keep a length <span class="math-container">$N$</span> vector in memory rather than a <span class="math-container">$N^2$</span> matrix; I also use a sparse matrix representation for <span class="math-container">$\mathbf{A}$</span>). Generally in Matlab at least, parallelised operations are much, much faster than nested loops.</p>
<p>(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.)</p>
<p>But remember, if you only calculate <span class="math-container">$\mathbf{A}$</span> 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.</p>

CC BY-SA 4.0
2022-04-11T20:48:50.790Z
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.
Connor Garcia
CC BY-SA 4.0
2022-04-11T22:47:47.217Z
@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.
Anders Sandberg
CC BY-SA 4.0
2022-04-12T15:33:18.810Z
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.
Connor Garcia
CC BY-SA 4.0
2022-04-11T19:11:27.733Z
2022-04-11T19:11:27.733Z
Anders Sandberg
<p>I'm probing into the <a href="https://www.tng-project.org/" rel="nofollow noreferrer">Illustris</a> API, and gathering information from a specific cosmos simulation, for a given redshift value.</p>
<p>This is how I request the api:</p>
<pre><code>import requests

baseUrl = 'http://www.tng-project.org/api/'
    
def get(path, params=None):
    # make HTTP GET request to path
    headers = {&quot;api-key&quot;:&quot;my_key&quot;}
    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(&quot;filename=&quot;)[1]
        with open(f'sky_dataset/simulations/{filename}', 'wb') as f:
            f.write(r.content)
        return filename # return the filename string
    return r
</code></pre>
<p>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:</p>
<pre><code>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 = &quot;http://www.tng-project.org/api/Illustris-1/snapshots/z=&quot; + str(redshift) + &quot;/subhalos/&quot; + str(subhalo_id)
sub = get(url) # get json response of subhalo properties
saved_filename = get(url + &quot;/cutout.hdf5&quot;,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 -&gt; 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()
</code></pre>
<p>The above plots:</p>
<p><a href="https://i.sstatic.net/6MLzW.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/6MLzW.png" alt="enter image description here" /></a></p>
<hr />
<p>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 <em>l</em> of <em>1.2 Mpc</em>, that is:</p>
<pre><code>Aij = 1 if rij ≤ l, otherwise 0
</code></pre>
<p>where <em>rij</em> is the distance between the two vertices, <em>i</em> and <em>j</em>.</p>
<hr />
<p>Any ideas on how to get this adjacency matrix, based on my linking length? Any help would be greatly appreciated.</p>

CC BY-SA 4.0
2022-04-11T16:27:05.170Z
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.
Connor Garcia
CC BY-SA 4.0
2022-04-11T20:40:43.363Z
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.
Greg Miller
<galaxy><distances><python><galactic-halo>
How to build an adjacency matrix of stars connectivity from their physical kpc distances?
8-Bit Borges
# How to build an adjacency matrix of stars connectivity from their physical kpc distances?

(Asked by: username_0 on 2022-04-10)

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.

**Tags:** galaxy, distances, python, galactic-halo
> **Comments:**
> **username_1**: 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.  
> **username_2**: 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.  

---

## 1 answers

### Answer by username_3 on 2022-04-11. Score: 3 | Accepted answer

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| continue to calculate $\Delta y$, otherwise go to the next pair. If $|\Delta y| too, now you calculate $\Delta x^2+\Delta y^2 +\Delta z^2$ and check if it is $: 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.

> **Comments:**
> **username_1**: 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.  
> **username_3**: @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.  
> **username_1**: 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.  

---

