hemantaph / ler

Gravitational waves lensing rate calculator
MIT License
2 stars 2 forks source link

Time delay sorting #3

Closed ohannuks closed 2 months ago

ohannuks commented 9 months ago

@hemantaph

Code:

import ler
from ler.rates import LeR
ler = LeR(verbose=False, npool=32)
unlensed_param = ler.unlensed_cbc_statistics(resume=True)
lensed_param = ler.lensed_cbc_statistics(size=100000, resume=True)
rate_ratio, unlensed_param_detectable, lensed_param_detectable =ler.rate_comparision_with_rate_calculation()
lensed_param_detectable.keys()
print(lensed_param_detectable['time_delays'])

Output:

array([[0.00000000e+00, 2.21840881e+06, 3.21975840e+06, 8.12671915e+06],
       [0.00000000e+00, 6.69836256e+06, 6.94011945e+06, 7.60169732e+06],
       [1.89326954e+07, 1.58780987e+06, 3.97358452e+05, 0.00000000e+00],
       ...,
       [0.00000000e+00, 2.03573957e+04, 2.56343643e+04, 1.08990879e+06],
       [0.00000000e+00, 9.27605219e+07, 0.00000000e+00, 0.00000000e+00],
       [0.00000000e+00, 3.36367791e+06, 3.36431948e+06, 3.61677533e+06]])

Most of the time delays are sorted, but the entry [2] seems to be reverse sorted. Is this expected?

Additional test:

import numba
@numba.jit
def is_sorted(a, n_elements):
    array = np.zeros(len(a))
    for i in range(len(a)):
         for j in range(n_elements[i]-1):
            if a[i][j] < a[i][j+1] :
                array[i] = 1
            else:
                array[i] = 0
    return array
# Test if time delays are sorted
is_sorted(lensed_param_detectable['time_delays'], lensed_param_detectable['n_images'])

Output:

array([1., 1., 0., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 0., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 0., 1., 1., 1.,
       1., 1., 1., 1., 0., 1., 1., 1., 0., 1., 1., 1., 0., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1., 1., 0., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 0., 1., 0., 1., 1., 1., 0.,
       0., 1., 1., 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1., 0., 1., 1.,
       1., 1., 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       0., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 0., 1., 1., 0., 1., 1., 1., 1., 1.,
       1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 0.,
       1., 1., 1., 1., 1., 1., 0., 1., 1., 1., 1., 0., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 0., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 0., 1., 1., 1., 1., 1., 1.,
       1.])

Most time delays are sorted, but a subset of them seems not to be.

Here are some of the time delays that are not sorted:

idx=np.array(is_sorted(lensed_param_detectable['time_delays'], lensed_param_detectable['n_images']), dtype=bool)
lensed_param_detectable['time_delays'][~idx]

Output:

array([[1.89326954e+07, 1.58780987e+06, 3.97358452e+05, 0.00000000e+00],
       [1.35931493e+07, 1.35098053e+07, 1.30916784e+07, 0.00000000e+00],
       [4.24965433e+07, 7.45334306e+06, 7.31561501e+06, 0.00000000e+00],
       [1.08477077e+08, 6.56249809e+07, 6.39096024e+07, 0.00000000e+00],
       [8.21647148e+06, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
       [9.46242094e+05, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
       [6.79474908e+05, 5.47475601e+04, 5.46652766e+04, 0.00000000e+00],
       [1.54439684e+07, 6.55248053e+06, 1.35104910e+05, 0.00000000e+00],
       [1.26731315e+07, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
       [2.12395774e+06, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
       [1.94328929e+07, 1.49791481e+06, 1.20741069e+05, 0.00000000e+00],
       [8.50901798e+06, 2.12890688e+06, 1.93639510e+06, 0.00000000e+00],
       [1.96050868e+08, 8.15432375e+06, 1.60148969e+06, 0.00000000e+00],
       [1.34020097e+07, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
       [1.11721941e+07, 7.42542614e+06, 7.42480866e+06, 0.00000000e+00],
       [5.50277841e+06, 5.37987667e+06, 5.29881519e+06, 0.00000000e+00],
       [3.06205471e+06, 6.03839696e+05, 4.76103751e+05, 0.00000000e+00],
       [1.56547619e+06, 1.26649689e+06, 1.03858835e+06, 0.00000000e+00],
       [4.90648396e+07, 7.75617389e+06, 2.45764948e+06, 0.00000000e+00],
       [5.32693453e+05, 2.62352514e+05, 1.39790283e+04, 0.00000000e+00],
       [1.04435157e+06, 9.74113143e+05, 9.63079875e+05, 0.00000000e+00],
       [7.33391520e+07, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
       [4.91482521e+05, 4.76539279e+05, 1.99022612e+05, 0.00000000e+00],
       [7.49013164e+06, 4.24963782e+06, 3.64685617e+06, 0.00000000e+00],
       [2.02259991e+05, 1.40748296e+04, 0.00000000e+00, 0.00000000e+00],
       [5.63403539e+08, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
       [2.37005483e+08, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
       [3.55995143e+06, 2.86064043e+06, 5.10503702e+05, 0.00000000e+00],
       [5.01127015e+06, 1.40572513e+06, 1.26516599e+06, 0.00000000e+00]])
hemantaph commented 9 months ago

Thanks for raising the issue, I will check and get back to you. But it should be all sorted out.

On Fri, 2 Feb 2024 at 10:03 AM, Otto @.***> wrote:

Code:

import ler from ler.rates import LeR ler = LeR(verbose=False, npool=32) unlensed_param = ler.unlensed_cbc_statistics(resume=True) lensed_param = ler.lensed_cbc_statistics(size=100000, resume=True) rate_ratio, unlensed_param_detectable, lensed_param_detectable =ler.rate_comparision_with_rate_calculation() lensed_param_detectable.keys() print(lensed_param_detectable['time_delays'])

Output:

array([[0.00000000e+00, 2.21840881e+06, 3.21975840e+06, 8.12671915e+06], [0.00000000e+00, 6.69836256e+06, 6.94011945e+06, 7.60169732e+06], [1.89326954e+07, 1.58780987e+06, 3.97358452e+05, 0.00000000e+00], ..., [0.00000000e+00, 2.03573957e+04, 2.56343643e+04, 1.08990879e+06], [0.00000000e+00, 9.27605219e+07, 0.00000000e+00, 0.00000000e+00], [0.00000000e+00, 3.36367791e+06, 3.36431948e+06, 3.61677533e+06]])

Most of the time delays are sorted, but the entry [2] seems to be reverse sorted. Is this expected?

Additional test:

import numba @numba.jit def is_sorted(a, n_elements): array = np.zeros(len(a)) for i in range(len(a)): for j in range(n_elements[i]-1): if a[i][j] < a[i][j+1] : array[i] = 1 else: array[i] = 0 return array

Test if time delays are sorted

is_sorted(lensed_param_detectable['time_delays'], lensed_param_detectable['n_images'])

Output:

array([1., 1., 0., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1., 0., 1., 1., 1., 0., 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 0., 1., 0., 1., 1., 1., 0., 0., 1., 1., 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 0., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 0., 1., 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1.])

— Reply to this email directly, view it on GitHub https://github.com/hemantaph/ler/issues/3, or unsubscribe https://github.com/notifications/unsubscribe-auth/AMW7ZWWOG6E245D4NY5IZFLYRRCNJAVCNFSM6AAAAABCV5KDIWVHI2DSMVQWIX3LMV43ASLTON2WKOZSGEYTGOBXHAYTENI . You are receiving this because you are subscribed to this thread.Message ID: @.***>

ohannuks commented 9 months ago

Great! I think it might just be that we're not using the sort_time_delays option in lestronomy when we solve for the image positions (that may be my bad!)

hemantaph commented 9 months ago

I have changed the lens equation solver as given below, but the problem persists. This happens in 1% of the of the events. Or should I reorder this 1% sample?

Note: This has nothing to do with whether it is analytic or not.

(
    x0_image_position,
    x1_image_position,
) = lens_eq_solver.image_position_from_source(
    sourcePos_x=x_source,
    sourcePos_y=y_source,
    kwargs_lens=kwargs_lens,
    solver="analytical",
    magnification_limit=1.0 / 1000.0,
    arrival_time_sort=True,
)