taskGNU Astronomy Utilities - Tasks: task #15803, Match program builds k-d tree and...


You are not allowed to post comments on this tracker with your current authentication level.

task #15803: Match program builds k-d tree and later read from it

Submitted by:  Mohammad Akhlaghi <makhlaghi>
Submitted on:  Thu 29 Oct 2020 10:51:11 PM UTC  
Should Start On:  Thu 29 Oct 2020 12:00:00 AM UTC Should be Finished on:  Thu 29 Oct 2020 12:00:00 AM UTC
Category:  Match Priority:  5 - Normal
Item Group:  Enhancement Status:  In Progress
Privacy:  Public Percent Complete:  90%
Assigned to:  Sachin Kumar Singh <sks_15> Open/Closed:  Open
Effort:  0.00

( Jump to the original submission)

Sun 18 Apr 2021 02:04:14 AM UTC, comment #9: 

In the process suggested in the previous comment, the coordinates of the two catalogs will be identical. Although that won't make any difference with regard to this particular test, a more robust test would be to add noise to one of the catalogs (or both) before the match.

Of course, the noise should be much smaller than the matching aperture size so you can still match (almost) everything.

I had been planning to add this feature (adding noise to table columns) in Gnuastro for a long time, but the motivation+time+free-mind never converged together until now! So I am happy to say that it has been implemented in Commit 5dc054e9d.

So if you install the latest Gnuastro (P.S.), you can run this command on the 'all-sub.fits' table to randomize the positions with a very small scatter of 0.00003 degrees! Feel free to change it as you like ;-)

asttable all-sub.fits \
         -cEDR3Name \
         -c'arith RAJ2000 0.00003 mknoise-sigma' \
         -c'arith DEJ2000 0.00003 mknoise-sigma' \
         -cGmag \
         --colmetadata=ARITH_1,RAJ2000,deg,"Right Ascension" \
         --colmetadata=ARITH_2,DEJ2000,deg,"Declination" \

Of course, Table also has the '--envseed' option to control the random-ness (and make it reproducible).

P.S. For those that can't build from the Git repository, the latest Gnuastro tarball and its corresponding PDF manual are available here:

Mohammad Akhlaghi <makhlaghi>
Project Administrator
Sat 17 Apr 2021 09:20:36 PM UTC, comment #8: 

Thanks a lot Sachin, I just had a fast look. Here are some points:

  • We won't be rebasing the commits into one commit for the final merge any more. We'll just rebase them all over the 'master' branch. So it is important that the history be well preserved and to have clear commit message. In this spirit I made some corrections in your commit message (Commit d9ee53845357b): 1) the title should start with the relevant component of Gnuastro (in this case "Library (match.h)"). 2) The first sentence of the body should start with "Until now, ..." and fully describe the problem that motivated this commit (before jumping into the solution). The solution should start in a separate paragraph starting with "With this commit ...". 3) Highlight codes with single quotes ('), not reverse quates (`). In Markdown, reverse quotes are indeed used to highlight code, but the commit messages aren't interpreted as markdown and a single quote is easier to read with no markdown. You can see the modified commit message in the P.S.
  • I noticed that the lines in 'during-dev-test-data/match-query.txt' ended with space characters. Git highlights them with red in a simple 'git diff' (or 'git log -p' after the commit). So before committing/pushing, be sure to run either of these two and make sure your commit never has such lines ;-).

After correcting the commit, I built Gnuastro over it and noticed that 'make check' failed because the test script wasn't executable, so I made the tiny Commit 44bedb4a over it. Besides making it executable, I also removed the '-match' part in the name: because the test is already prefixed with 'match/' and having the 'match' string repeated two times is redundant ;-).

About the tests, instead of random trails (which can also be hard to reproduce later), it may be easier to just use real data ;-). With this command you can download all Gaia stars around the celestial coordinates (0,0) within a radius of 1 degrees (and save them in 'all.fits'):

astquery vizier --dataset=gaiaedr3 --center=0,0 --radius=1 -cEDR3Name,RAJ2000,DEJ2000,Gmag -oall.fits

It results in 10056 stars. Now, you can select a random sub-set of the rows (for matching) using the '--rowrandom' option. For example with this command, you can extract 1000 rows of it randomly:

asttable all.fits --rowrandom=1000 -oall-sub.fits

You can now match this sub-set with the full table :-). Since its a sub-set you know a-priori the total number of successful matches (1000 in the case above), and things are easier to test.

To make the subset selection reproducible, see the Generating Random Numbers section of the book. Here is what you should do:

export GSL_RNG_TYPE=ranlxs1
export GSL_RNG_SEED=1619042801
asttable all.fits --envseed --rowrandom=1000

If you later want more stars (for more serious tests), you can simply increase the value to '--radius' in the 'astquery' command above :-).

The good thing with this process is that we won't need to keep data in the commits and at the very start of the benchmarking script, we can simply download the necessary dataset.

Finally, we will probably be using this benchmarking a lot in the future to improve/debug the '--kdtree=automatic' algorithm. So its good to write it in Bash and move it into the Match program directory. On the command-line, you can use the 'time' command before any command to get the resource usage of the command.

Note that for now, the benchmarking is only to see when Gnuastro's own sort-based algorithm should be switched to the k-d tree based algorithm ;-). For comparing with other k-d tree implementations, we'll need better control over the environment (using Maneage (by the way: Maneage's paper has been accepted into IEEE's CiSE)

P.S. Modified commit message (Commit d9ee53845357b)

Library (match.h): match_coordinate_ replaced by match_sort_based_

Until now, there was only a single match algorithm in Gnuastro, so the name
of the respective functions in the library had a 'match_coordinate_'
prefix. However, in this branch we are adding a new k-d tree based
matching, so that name for the initial algorithm could cause confusion
(because k-d tree also uses coordinates!).

With this commit, those same function names are now prefixed with
'match_sort_based_'. This change is done in all the relevant places:
bin/match/match.c, lib/gnuastro/match.h, lib/match.c and in

Furthermore, a 'make check' test has been added for the k-dtree based
matching which does simple matching based on the predefined inputs used for
the old matching algorithm.

Finally, some scripts have been added in the temporary directory for
benchmarking tests to see the efficiency of k-d tree over sort-based
matching. Two scripts are written in '/during-dev-test-data/scripts'
(namely 'kdtree-gen.py', which generates the pseudo random output tests,
and 'benchmark.py') which will calculate the run time for the scripts and
generate a graph for number of times the program is executed vs. time
taken. The script for benchmarking is simple and might take long time for
>10 executions. A possible efficient method is multithreading, which will
be implemented soon. After this we can implement the --kdtree=automatic
option based on the result.

Mohammad Akhlaghi <makhlaghi>
Project Administrator
Sat 17 Apr 2021 07:17:15 PM UTC, comment #7: 

I made the commit in which I've changed the name of `match_coordinate_*` functions which did sort-based matching to `match_sort_based_*`. Also, I added the tests to be performed during `make test` in tests/match/.

I also wrote the scripts to generate test files during benchmarking and test the efficiency of kdtree over sort-based matching. I've written briefly about that in the commit message.

Commit message:

Changes the name of `match_coordinate_*` functions which did sort based
matchin to `match_sort_based_*` for clarity. This change is done in
bin/match/match.c, lib/gnuastro/match.h, lib/match.c and in

Tests for the kdtree based matching is introduced which does simple
matching based on the predefined inputs given in tests/match/. Also,
added the test in Makefile.am to be included during `make test`.

Added scripts to generate test files during benchmarking and test the
efficiency of kdtree over sort-based matching. Two scripts are written
in /during-dev-test-data/scripts namely `kdtree-gen.py`, which generates
the pseudo-random output tests, and `benchmark.py` which will calculate
the run time for the scripts and generate a graph for a number of times
the program is executed vs. time taken. The script for benchmarking is
simple and might take a long time for >10 executions. A possible efficient
method is a multithreading, which will be implemented soon. After this, we
can implement the --kdtree=automatic option based on the result.

Sachin Kumar Singh <sks_15>
Project MemberIn charge of this item.
Mon 29 Mar 2021 01:49:26 AM UTC, comment #6: 

I done a full re-write in Commit a6c61757668 (commit message in P.S.) and look clear for the final steps now :-). In short, I noticed that all the troubles we were dealing with (one point matching multiple in the other) was already treated before in the sort-based matching.

As Ken Thompson says in The Art of Unix programming: "One of my most productive days was throwing away 1000 lines of code."! This fit so well here :-)!

Sachin, the job is almost done, I have put some next steps in the commit message, so if you can implement them in followup commits, we can finalize things and finish it off with the documentation. I am really eager to see the result of the test to see from what point the k-d tree becomes faster than the sort-based method and how it improves (maybe as a plot showing two curves).

P.S. Commit message:

Library (match.h): Using the old infra-structure for double-matches

Until now, we hadn't considered the fact that the main body of work for
cleaning the final match (removing multiple nearby matches and formatting
of the output) was already implemented in the sort-based match!

With this commit, instead of trying to debug the steps in the last few
commits, I just moved everything into a temporary file (deleted them from
'lib/match.c') and added everything modeled based on sort-based matching
implementation: using the same 'bina' array. Once that was implemented, it
was really easy to just use the final functions of the sort-based matching
and everything seems to work properly now, with no segmentation

In those functions of the sort-based matching that are now also used in the
k-d tree based matching, the '_coordinate_' part is removed and they have
been moved to a separate (generic) part of 'lib/match.c' to highlight that
they are low-level and shared between both matching methods.

Generally, in programming, its always good to avoid re-writing a function
and modularize/reuse previously written functions. So for example if a bug
is found (or a feature is added) in the sort-based matching, it will
automatically be implemented in the k-d tree based matching too and no
extra work is necessary.

Some other minor points:

- The 'inplace' argument to the sort-based matching was for the sorting
   step. It is not used in the k-d tree based matching, so I removed it
   from the arguments.

Things remaining:

- Find a good name for the sort-based matching.

- Complete all the expected features for k-d tree based matching and write
   tests for them during 'make check'.

- Test to see which matching methodology is faster for how many rows. It
   is expected that for very few rows (like 10 or 100) probably the
   sort-based matching is faster, but from a cetain number of rows and
   certain number of threads, the k-d tree based matching will be
   faster. Once we find this, we will be able to implement the 'automatic'

Mohammad Akhlaghi <makhlaghi>
Project Administrator
Sun 13 Dec 2020 02:53:07 AM UTC, comment #5: 

I made some further improvements in Commit 6b5e8ef1c7, please look into them carefully.

The logic of the program is clear (all the UI, like reading input columns, or setting the kdtree mode) is done in 'ui.c'. We shouldn't do UI-related things in the processing functions.

The main work is now done in that commit: the raw table that is filled in the multi-threaded function (p->c2match) contains the index in 'coords1' ('i1') that matches the respective index in 'coords2' ('i2'). Something like this: 'c2match[i2]=i1'.

You can use this to easily create an output (after the threading has finished) in the same format as 'gal_match_coordinates'.

Mohammad Akhlaghi <makhlaghi>
Project Administrator
Thu 10 Dec 2020 03:42:51 AM UTC, comment #4: 

A first implementation of the option reading and preparations has been done in Commit 61e7108b5b12 of my development branch on my development repository.

With the example code in the book for the two k-d tree related functions, it should be pretty streightforward to complete this task ;-).

Mohammad Akhlaghi <makhlaghi>
Project Administrator
Wed 04 Nov 2020 10:28:27 PM UTC, comment #3: 

After the first implementation is complete, we can add an internal/automatic decision to use k-d trees or not.

This is the scenario I have in mind: we test both matching scenarios by increasing the size of one (or both) catalog(s) gradually. At first (for example maybe when there are only 100 rows), the simple default matching should be faster.

But as the number increases, from some point, the k-d tree should become faster (accounting for the number of threads). We can then automatically decide to use a k-d tree instead of the basic matching.

In terms of implementation, we can add a new k-d tree mode: '--kdtree=automatic'.

In this scenario, the default value to '--kdtree' (in the main configuration file of Match) can be set to 'automatic'. However, in such a scenario, the user would have to have the option to disable k-d trees, so we should also add a "--kdtree=none" mode (which will definitely use the basic matching).

But this should be tested and done after we finish the initial implementation. So it isn't urgent now. I am just writing it here so we don't forget ;-).

Mohammad Akhlaghi <makhlaghi>
Project Administrator
Tue 03 Nov 2020 10:52:37 PM UTC, comment #2: 

Thanks for the nice review, its great to add this feature soon :-).

I think the main thing you missed is this: the k-d tree file should just include the k-d tree (the raw output of 'gal_kdtree_create', which is only two columns).

So when you want to search with that k-d tree, you need the actual table that the k-d tree was created with also.

The reason for this is that k-d trees usually become useful for large tables and large tables will take space to store (for example assume that 'large.fits' is a 100Mb table). There is no need to force the user to have an extra copy of the large table inside the 'kdtree.fits'.

Does this sound reasonable?

About the acceptable strings given to the '--kdtree' option, it just occurred to me that we may want a third alternative also:

'--kdtree=build': will build the k-d tree, and write it inside the file given to the '--output' option.

'--kdtree=internal': we will construct the k-d tree internally, use it in parallel, and free the k-d tree after the match.

'--kdtree=FILENAME': it will read the k-d tree from the given file, use it in parallel, and free it. Ofcourse, in this scenario, we need to check that the k-d tree in the given file has the same number of rows as 'large.fits'.

By the way, when matching in k-d tree mode, we should always assume that the k-d tree will be built/used over the first input catalog.

How does this sound?

Mohammad Akhlaghi <makhlaghi>
Project Administrator
Tue 03 Nov 2020 08:49:59 PM UTC, comment #1: 

I'm about to begin with the implementation of this task and just wanted to discuss the details of the implementation.

So, we simply have to add a --kdtree option in the match program which takes 2 parameters. If the parameter is `build`, we will make a KD tree using the X, Y coordinates in `large.fits` using the `gal_kdtree_create` function and we can output the kdtree in a file using the --output option.

Next for matching the coordinates in the `small.fits` it will use the `gal_kdtree_nearest_neighbour` to match each point in the `small.fits` to the nearest neighbor in the kd-tree.

If the above case is okay, I think the command should look as below, because now the call to `large.fits` is redundant as it is used to create kdtree which has already been done.

astmatch small.fits --ccol1=X,Y --kdtree=kdtree.fits

But I also understand that using the command below is more similar to how `astmatch` it is normally used and hence much more familiar to the users.

astmatch large.fits small.fits --ccol1=X,Y --ccol2=X,Y \

What do you think?

Sachin Kumar Singh <sks_15>
Project MemberIn charge of this item.
Thu 29 Oct 2020 10:51:11 PM UTC, original submission:  

Gnuastro's library now has a k-d tree library. With k-d trees, matching large catalogs can be greatly sped up.

So I propose the following strategy to match two catalogs:

Currently to match two catalogs, we call match like below. For the context here, let's call them 'large.fits' and 'small.fits', and let's assume we want to match them with their the 'X' and 'Y' columns, within an aperture of 0.5 pixels.

astmatch large.fits small.fits --ccol1=X,Y --ccol2=X,Y \

But to employ k-d trees as is available in the Gnuastro library, we can take the following approach. The user first runs 'astmatch' only on 'large.fits' with the '--kdtree=build' option like this (the k-d tree will be built in the value to '--output', or '-o'):

astmatch large.fits --ccol1=X,Y --kdtree=build -okdtree.fits

Later, when they want to match it with the small catalog they would run it like this (very similar to the first call):

astmatch large.fits small.fits --ccol1=X,Y --ccol2=X,Y \

Match would read the k-d tree from the given file, with 'large.fits' and 'small.fits', then it will use a multi-threaded function to call Gnuastro's k-d tree matching function in parallel. Besides the improved k-d tree search method, this would also benefit from parallel processing, thus doing the job very fast.

Mohammad Akhlaghi <makhlaghi>
Project Administrator


No files currently attached


Depends on the following items: None found

Items that depend on this one: None found


Carbon-Copy List
  • -email is unavailable- added by sks_15 (Posted a comment)
  • -email is unavailable- added by pedram
  • -email is unavailable- added by makhlaghi (Submitted the item)

    Do you think this task is very important?
    If so, you can add your encouragement to it.
    This task has 0 encouragements so far.

    Only logged-in users can vote.




    Follow 6 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2021-03-29 makhlaghi Percent Complete70% => 90%
    2020-12-13 makhlaghi Percent Complete20% => 70%
        SummaryMatch program should build k-d tree and later read from it => Match program builds k-d tree and later read from it
    2020-12-10 makhlaghi StatusNone => In Progress
        Percent Complete0% => 20%
    2020-10-30 pedram Carbon-Copy- => Added pedram

    Back to the top

    Powered by Savane 3.7