michitaro / healpix

An implementation of HEALPix in JavaScript / TypeScript
MIT License
13 stars 4 forks source link

norder limitation #25

Open triple7 opened 2 years ago

triple7 commented 2 years ago

Hi,

I'm using healpix for a different project, and have been looking around for implementations and came across yours which seems interesting. However, the limitation is at norder <= 15. Current norder doesn't allow high resolution images of nebulas for example when overlaying different frequency image survey maps (infrared, nir, xray etc). Can you tell me what the exact reason for this is? I might be able to help somehow.

I downloaded healpix 3.8 and the doc provides a healpix2 cpp that allows for much higher resolutions, which indicates that the norder can be increased?

Any pointers to the reasons of the limitations greatly appreciated as I can further investigate.

michitaro commented 2 years ago

The direct root of the limitation is here. https://github.com/michitaro/healpix/blob/6282f3e0da97264e71fed468044daaf44e8296bb/src/index.ts#L546 This function is for converting scheme from RING to NEST.

If the function could handle values of order > 16, the limitation can be eliminated. (But, we have to care about the range of values being not bigger than Number.MAX_SAFE_INTEGER) https://developer.mozilla.org/en-US/docs/Web/JavaScript/Reference/Global_Objects/Number/MAX_SAFE_INTEGER

michitaro commented 2 years ago

BTY, this library's license is MIT and HEALPix's license is GPL, so we cannot port HEALPix's implementation into this library.

triple7 commented 2 years ago

Thanks for these, I will study them in parallel to the healpix primer I have been taking notes from. Here are some thoughts:

I'll get back when I find something interesting

Thanks again

On 27 Dec 2021, at 11:39 am, Koike Michitaro @.***> wrote:

The direct root of the limitation is here. https://github.com/michitaro/healpix/blob/6282f3e0da97264e71fed468044daaf44e8296bb/src/index.ts#L546 https://github.com/michitaro/healpix/blob/6282f3e0da97264e71fed468044daaf44e8296bb/src/index.ts#L546 This function is for converting scheme from RING to NEST.

If the function could handle values of order > 16, the limitation can be eliminated. (But, we have to care about the range of values being not bigger than Number.MAX_SAFE_INTEGER) https://developer.mozilla.org/en-US/docs/Web/JavaScript/Reference/Global_Objects/Number/MAX_SAFE_INTEGER https://developer.mozilla.org/en-US/docs/Web/JavaScript/Reference/Global_Objects/Number/MAX_SAFE_INTEGER — Reply to this email directly, view it on GitHub https://github.com/michitaro/healpix/issues/25#issuecomment-1001283353, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAJHQFOUDG673PWHDHBGOWTUS67TRANCNFSM5KYW7M4Q. Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub. You are receiving this because you authored the thread.

michitaro commented 2 years ago

Sorry, I forgot there was a similar discussion. See this issue. https://github.com/michitaro/healpix/issues/8

JavaScript's bit-operation supports only 32bit integers, so bit_{,de}combine cannot handle large integers. If the function could handle values without bit-operations, the limitation would be reduced.

triple7 commented 2 years ago

Yes, that's what I thought as we are talking powers of 2 to subdivide each unit pixel of the sphere. Have you considered breaking the bit to smaller parts? like taking the overflow bits and shifting them across to another 32-bit integer? I'm going to try these to see how far it can be taken. I'm looking at getting norder to up to 23.

On 27 Dec 2021, at 4:44 pm, Koike Michitaro @.***> wrote:

Sorry, I forgot there was a similar discussion. See this issue.

8 https://github.com/michitaro/healpix/issues/8

JavaScript's bit-operation supports only 32bit integers, so bit_{,de}combine cannot handle large integers. If the function could handle values without bit-operations, the limitation would be reduced.

— Reply to this email directly, view it on GitHub https://github.com/michitaro/healpix/issues/25#issuecomment-1001369587, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAJHQFLCFTGAJCLQCH2LMO3UTADMXANCNFSM5KYW7M4Q. Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub. You are receiving this because you authored the thread.

triple7 commented 2 years ago

I found a thin wrapper for the native BIGINT type which provides all bit operations, arithmetic operations and others. I think these are sufficient to provide larger npix integers for healpix. I will keep reading your code and understand it better and try a version using some larger powers of 2.

On 27 Dec 2021, at 4:44 pm, Koike Michitaro @.***> wrote:

Sorry, I forgot there was a similar discussion. See this issue.

8 https://github.com/michitaro/healpix/issues/8

JavaScript's bit-operation supports only 32bit integers, so bit_{,de}combine cannot handle large integers. If the function could handle values without bit-operations, the limitation would be reduced.

— Reply to this email directly, view it on GitHub https://github.com/michitaro/healpix/issues/25#issuecomment-1001369587, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAJHQFLCFTGAJCLQCH2LMO3UTADMXANCNFSM5KYW7M4Q. Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub. You are receiving this because you authored the thread.

michitaro commented 2 years ago

Have you considered breaking the bit to smaller parts? like taking the overflow bits and shifting them across to another 32-bit integer?

I think this is a good idea and more feasible than introducing BigInt.

I wonder introducing BigInt needs a lot of changes.

triple7 commented 2 years ago

I started converting operators to bigInt versions and am adding sqrt using newton’s method. Will see if it will work on all test cases properly

Sent from sputnik7

On 28/12/2021, at 5:04 PM, Koike Michitaro @.***> wrote:

 Have you considered breaking the bit to smaller parts? like taking the overflow bits and shifting them across to another 32-bit integer? I think this is a good idea and more feasible than introducing BigInt.

I wonder introducing BigInt needs a lot of changes.

— Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications on the go with GitHub Mobile for iOS or Android. You are receiving this because you authored the thread.

michitaro commented 2 years ago

The current test case generator generates only cases for norder <= 16, so tweak the file before your tests. https://github.com/michitaro/healpix/blob/master/spec/generate-testcase.py The case generator generates cases using random numbers, so the cases might not be comprehensive (edge cases might be missing). Keep in mind that.

I wonder BigInt and Newton Method bring some disadvantages on performance, so a flag to turn off BigInt mode will be helpful.

triple7 commented 2 years ago

I’m using test cases on the Hubble space telescope hips surveys. I checked the potential performance issues, however I’m not sure if 12 million bigint operations a second would be an issue in applications that require healpix. I’m writing a separate set of functions, wrapped into a HealpixIndex object. I’m just filtering the function chain to see where the bigints come in to play. Nside can go up to 2 to 26 as bigint’s max stable is 2 to 52. So far the biggest possible value would be 12nsidenside which is 2 to 53. It’s the possible trig functions or others i have not clearly identified which may be an issue but I think it’s worth it.

Sent from sputnik7

On 28/12/2021, at 5:34 PM, Koike Michitaro @.***> wrote:

 The current test case generator generates only cases for norder <= 16, so tweak the file before your tests. https://github.com/michitaro/healpix/blob/master/spec/generate-testcase.py The case generator generates cases using random numbers, so the cases might not be comprehensive (edge cases might be missing). Keep in mind that.

I wonder BigInt and Newton Method bring some disadvantages on performance, so a flag to turn off BigInt mode will be helpful.

— Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications on the go with GitHub Mobile for iOS or Android. You are receiving this because you authored the thread.

michitaro commented 2 years ago

I’m not sure if 12 million bigint operations a second would be an issue in applications that require healpix.

Wow, bigint operations are so fast? If so, they won't be bottleneck.

I’m writing a separate set of functions, wrapped into a HealpixIndex object.

😲

triple7 commented 2 years ago

I've read through all the functions you wrote, and compared them with my notes, and your header and nomenclature for the variables, projection methods and representations. it seems as though p and q come from the bit_combine and bit_decombine with some propagation into a few more functions so it shouldn't be too bad. Reconverting to number in crucial areas with pi divisors and others should isolate most of the bigInt requirements.

one thing of note, since a lot of transforms have a definite pattern on the 3 ring region cases, there might be a way to just generate some lookup tables. nside, jj, pp and some others. This can probably boost performance slightly.

On 28 Dec 2021, at 11:39 pm, Koike Michitaro @.***> wrote:

I’m not sure if 12 million bigint operations a second would be an issue in applications that require healpix.

Wow, bigint operations are so fast? If so, they won't be bottleneck.

I’m writing a separate set of functions, wrapped into a HealpixIndex object.

😲

— Reply to this email directly, view it on GitHub https://github.com/michitaro/healpix/issues/25#issuecomment-1002110650, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAJHQFNF34OSS4OQVVZI2R3UTG4YXANCNFSM5KYW7M4Q. Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub. You are receiving this because you authored the thread.

triple7 commented 2 years ago

Just wanted to double check with you in terminology::

Are both p and q meant to be ints or floats?

Thanks

On 28 Dec 2021, at 11:39 pm, Koike Michitaro @.***> wrote:

I’m not sure if 12 million bigint operations a second would be an issue in applications that require healpix.

Wow, bigint operations are so fast? If so, they won't be bottleneck.

I’m writing a separate set of functions, wrapped into a HealpixIndex object.

😲

— Reply to this email directly, view it on GitHub https://github.com/michitaro/healpix/issues/25#issuecomment-1002110650, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAJHQFNF34OSS4OQVVZI2R3UTG4YXANCNFSM5KYW7M4Q. Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub. You are receiving this because you authored the thread.

triple7 commented 2 years ago

Hi,

So I took the tests suite and customised them as I wasn't able to run npm run my-test. I extracted the json test suite and wrote my own routine to compare the various expected/resulting outputs. I got the below:

vec2pix_nest mismatch 0 vec2pix_ring mismatch 0 ang2pix_nest mismatch 0 ang2pix_ring mismatch 0 nest2Ring mismatch 0 ring2nest mismatch 0 pix2vec_nest mismatch 0 pix2vec_ring mismatch 0 nside2pixarea mismatch 0 nside2resol mismatch 0 max_pixrad expected 0.48146047147606696 got 0.4814604714760669 mismatch 13 corners_nest mismatch 0 corners_ring mismatch 0 orderpix2uniq mismatch 0 uniq2orderpix mismatch 0 nside2order mismatch 0

For some reason maxrad has a distance of abs(a-b) of 5.xxxxe-17, which I suspect is due to conversion from BigInt to Number before passing it through some divisor.

Coincidentally, 10e-17 is probably the limit of base 10 digits on an integer, so the conversion issue may be there. I even tried abs(a-b) > 10e-5 then false but no luck so far. Will dig deeper until all have 0 mismatch.

On 28 Dec 2021, at 11:39 pm, Koike Michitaro @.***> wrote:

I’m not sure if 12 million bigint operations a second would be an issue in applications that require healpix.

Wow, bigint operations are so fast? If so, they won't be bottleneck.

I’m writing a separate set of functions, wrapped into a HealpixIndex object.

😲

— Reply to this email directly, view it on GitHub https://github.com/michitaro/healpix/issues/25#issuecomment-1002110650, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAJHQFNF34OSS4OQVVZI2R3UTG4YXANCNFSM5KYW7M4Q. Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub. You are receiving this because you authored the thread.

triple7 commented 2 years ago

Hi again, and HNY :)

So here is the result of testing for norder up to 24. note: the comparisons are slowing the process down as some have to go through arrays of n elements. note2: The example mismatch really isn't one if you think about it. There's a single 10e-17th decimal change between the python output and the js result.

vec2pix_nest mismatch 0 vec2pix_ring mismatch 0 ang2pix_nest mismatch 0 ang2pix_ring mismatch 0 nest2ring mismatch 0 ring2nest mismatch 0 pix2vec_nest mismatch 0 pix2vec_ring mismatch 0 nside2pixarea mismatch 0 nside2resol mismatch 0 max_pixrad expected example 0.48146047147606696 got 0.4814604714760669 mismatch 21 corners_nest mismatch 0 corners_ring mismatch 0 orderpix2uniq mismatch 0 uniq2orderpix mismatch 0 nside2order mismatch 0 made 278096 comparisons in 3.011 seconds

On 28 Dec 2021, at 11:39 pm, Koike Michitaro @.***> wrote:

I’m not sure if 12 million bigint operations a second would be an issue in applications that require healpix.

Wow, bigint operations are so fast? If so, they won't be bottleneck.

I’m writing a separate set of functions, wrapped into a HealpixIndex object.

😲

— Reply to this email directly, view it on GitHub https://github.com/michitaro/healpix/issues/25#issuecomment-1002110650, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAJHQFNF34OSS4OQVVZI2R3UTG4YXANCNFSM5KYW7M4Q. Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub. You are receiving this because you authored the thread.

michitaro commented 2 years ago

Oh, sorry for the delayed response and Happy new year! Wow! Super progress!

michitaro commented 2 years ago

If you can make a PR, I will tweak the test codes to suit your codes.

triple7 commented 2 years ago

Hi,

I'm just having some issues with testing the function query_disc_inclusive_nest. I'm not sure what it's meant to return, normally a list of ipix?

Once I got this one fixed, I can send you the repo.

On 5 Jan 2022, at 11:50 am, Koike Michitaro @.***> wrote:

If you can make a PR, I will tweak the test codes to suit your codes.

— Reply to this email directly, view it on GitHub https://github.com/michitaro/healpix/issues/25#issuecomment-1005312879, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAJHQFIESCITYIJ76XDSNYLUUOPW5ANCNFSM5KYW7M4Q. Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub. You are receiving this because you authored the thread.

michitaro commented 2 years ago

query_disc_inclusive_nest does not return any list but accepts a callback function whose only parameter is ipix. like this: https://github.com/michitaro/healpix/blob/7c0e9d8f21cf1126b255ac37b90865fece321413/docs/example/src/query_disc/main.ts#L26

This is to avoid allocating a large amount of memory when the resultant ipxs are too many. If you want to get the list of ipxs, you can do something like this:

const a = []

for (const ipix of query_disc_inclusive_nest(...)) {
  a.push(ipix)
}
michitaro commented 2 years ago

How the example of healpix.query_disc_inclusive_nest work can be seen here. https://michitaro.github.io/healpix/query_disc.html

triple7 commented 2 years ago

thanks for the explanation below. However, I'm not sure what function to callback.

You wrote:

const a = []

for (const ipix of query_disc_inclusive_nest(...)) { a.push(ipix) }

the last argument of query_inclusive_nest is meant to be a callback. Do I just call function (ipix) as the last argument?

Sorry, not very familiar with the ES6 js semantics.

On 11 Jan 2022, at 12:09 pm, Koike Michitaro @.***> wrote:

query_disc_inclusive_nest does not return any list but accepts a callback function whose only parameter is ipix. like this: https://github.com/michitaro/healpix/blob/7c0e9d8f21cf1126b255ac37b90865fece321413/docs/example/src/query_disc/main.ts#L26 https://github.com/michitaro/healpix/blob/7c0e9d8f21cf1126b255ac37b90865fece321413/docs/example/src/query_disc/main.ts#L26 This is to avoid allocating a large amount of memory when the resultant ipxs are too many. If you want to get the list of ipxs, you can do something like this:

const a = []

for (const ipix of query_disc_inclusive_nest(...)) { a.push(ipix) } — Reply to this email directly, view it on GitHub https://github.com/michitaro/healpix/issues/25#issuecomment-1009534697, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAJHQFP3LT3CUS6BUJPEZFDUVOGNZANCNFSM5KYW7M4Q. Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub. You are receiving this because you authored the thread.

triple7 commented 2 years ago

No all good, I did some tests for length of arrays when pushing ipix through the callbacks and they match.

I'll run some more tests then show you the results on a high order HST survey image

On 11 Jan 2022, at 12:11 pm, Koike Michitaro @.***> wrote:

How the example of healpix.query_disc_inclusive_nest work can be seen here. https://michitaro.github.io/healpix/query_disc.html https://michitaro.github.io/healpix/query_disc.html — Reply to this email directly, view it on GitHub https://github.com/michitaro/healpix/issues/25#issuecomment-1009535519, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAJHQFPDSB4ZW7JMQ5TVBUTUVOGU3ANCNFSM5KYW7M4Q. Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub. You are receiving this because you authored the thread.

triple7 commented 2 years ago

Hi,

So I've gotten the pixlist working from query_disc_inclusive_nest, however the results don't match. I tried adding values to the healpy module function ipix = hp.query_disc(nside=8, vec=[1, 0, 0], radius=0.008328107069661798)

it returns an empty array. The results I get from my modification matches another healpix module, which are:

nside=8, v=[1, 0, 0] radius=0.008328107069661798 ​ 304n ​ 293n ​ 282n ​ 271n

Would it be possible to create some test cases of returned ipix lists from the generate-testcase.py?

On 11 Jan 2022, at 12:11 pm, Koike Michitaro @.***> wrote:

How the example of healpix.query_disc_inclusive_nest work can be seen here. https://michitaro.github.io/healpix/query_disc.html https://michitaro.github.io/healpix/query_disc.html — Reply to this email directly, view it on GitHub https://github.com/michitaro/healpix/issues/25#issuecomment-1009535519, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAJHQFPDSB4ZW7JMQ5TVBUTUVOGU3ANCNFSM5KYW7M4Q. Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub. You are receiving this because you authored the thread.

michitaro commented 2 years ago

Hi, Sorry for the delayed response. Use hp.query_disc(nside=8, vec=[1, 0, 0], radius=0.008328107069661798, inclusive=True, nest=True) instead. It also returns the same results.

import healpy

ipix = healpy.query_disc(nside=8, vec=[1, 0, 0], radius=0.008328107069661798, inclusive=True, nest=True)
print(ipix)
# [271 282 293 304]