Skip to content

Update Neutron Star Section#425

Open
Zhengjingyi0823 wants to merge 8 commits into
LSST-strong-lensing:mainfrom
Zhengjingyi0823:jingyi
Open

Update Neutron Star Section#425
Zhengjingyi0823 wants to merge 8 commits into
LSST-strong-lensing:mainfrom
Zhengjingyi0823:jingyi

Conversation

@Zhengjingyi0823
Copy link
Copy Markdown
Contributor

Adding function for Neutron stars.

@codecov
Copy link
Copy Markdown

codecov Bot commented May 19, 2026

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 98.61%. Comparing base (7fa85b4) to head (a17debf).
⚠️ Report is 2 commits behind head on main.

Additional details and impacted files
@@            Coverage Diff             @@
##             main     #425      +/-   ##
==========================================
- Coverage   98.61%   98.61%   -0.01%     
==========================================
  Files         104      106       +2     
  Lines        8030     8085      +55     
==========================================
+ Hits         7919     7973      +54     
- Misses        111      112       +1     
Files with missing lines Coverage Δ
slsim/Sources/Events/BNSMerger/BNSMerger_pop.py 100.00% <100.00%> (ø)
slsim/Sources/Events/event_pop.py 100.00% <100.00%> (ø)
slsim/Sources/Supernovae/supernovae_pop.py 100.00% <100.00%> (ø)

... and 1 file with indirect coverage changes

🚀 New features to boost your workflow:
  • ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

Copy link
Copy Markdown
Contributor

@sibirrer sibirrer left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

thank you very much @Zhengjingyi0823! I made some comments and requests to further polish documentation and code

else:
raise ValueError("model should be chosen from 'BNS' or 'SNIa'")

def calculate_event_rate(self, z):
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

please precisely document the specific return with units

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I have added the description of the outputs of functions and their units.

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this should be in the docstrings of this definition, not the init statement

Comment thread slsim/Sources/Events/event_pop.py Outdated
raise ValueError("model should be chosen from 'BNS' or 'SNIa'")

def calculate_event_rate(self, z):
if self.model_name == "BNS":
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

if you have the definitions of the two separate classes named the same, you don't need an if/else statement here

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I added a wrapper function called calculate_event_rate in supernovae_pop.py so that we don’t need to change old code which calls the calculate_SNIa_rate and can simplify code in event_pop.py at the same time.

ft_d = 1 / (t_d * (np.log(self.t_d_max / self.t_d_min)))
return ft_d

def z_from_time(self, t):
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this is a very generic definition that deserves to be at a different place where it can be re-used. This routine is not specific of what the BNSMerger_pop class is about

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I moved delay_time_destribution out of the class BNSMergerRate. Since it is normalized, I renamed it as norm_delay_time_distribution.

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this does not addresses your comment. This comment was about z_from_time(t)

self.t_d_max = 13.8 # Gyr, Hubble time approximately

def calculate_star_formation_rate(self, z):
"""Calculates the binary formation rate. (Eq 3 - Kuwahara et al. 2025)
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

is this star formation or binary formation?

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I originally kept the name for consistency with the supernova population code. Since this function is mainly used within the BNS merger model and only externally called in the BNS merger notebook, I renamed it to calculate_binary_formation_rate for clarity.

"""Calculates the binary formation rate. (Eq 3 - Kuwahara et al. 2025)

:param z: redshift (z>=0)
:param z_m: peak-related redshift parameter
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this does not seem to be an input

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I moved the fixed model parameters nu, a, b, and z_m out of the :param: section and rewrote their explanation in the correct format.


:param t_d: time delay (t_d>=0) in [Gyr]
:param t_d_min: minimum of t_d in [Gyr], assumed as 20 Myr
:param t_d_max: maximum of t_d in [Gyr], assumed as Hubble time
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this is not an input here. I would make this either a default-valued input or better document it in the class initialization what it means

(Eq 4 - Kuwahara et al. 2025)

:param t_d: time delay (t>=0)
:param t: time at final integrated redshift (t>=0)
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

is there a condition that t>t_d? or something like that?

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

  1. The condition t_d < t is already enforced by the integration range in calculate_event_rate. To make this clearer, I added this condition as a part of comment above the code.


:param z: redshift (z>=0)

:return: BNS merger rate R_m(z), following the unit of R_f
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

what units?

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The specific unit is added.

def calculate_event_rate(self, z):
"""Calculates the rate of events, such as BNS merger. (Eq 4 - Kuwahara et al. 2025)

:param z: redshift (z>=0)
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this is an array of redshifts, please specify. Do they need to be ascending/decending etc?

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I completed the description of z in function calculate_event_rate. Array z don't need to be sorted but in practice use of BNS merger model, it should be descending since it is converted from t by function z_from_time with the time-delay order from t_d_min to t_d_max.


rate_array = event_pop.calculate_event_rate(self.z_array)

npt.assert_almost_equal(rate_array[0], 0.03081, decimal=3)
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

here I would test against the specific class instance's routine (which should be the same) instead of against a numerical answer (that is a test for the individual classes test function

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I have modified test_event_pop.py to compare the output of EventPopulation with the corresponding specific population model class.

@Zhengjingyi0823 Zhengjingyi0823 requested a review from sibirrer May 26, 2026 04:10
Copy link
Copy Markdown
Contributor

@sibirrer sibirrer left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thank you @Zhengjingyi0823 ! I still have some comments that should be addressed

ft_d = 1 / (t_d * (np.log(self.t_d_max / self.t_d_min)))
return ft_d

def z_from_time(self, t):
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this does not addresses your comment. This comment was about z_from_time(t)

def __init__(self, model, cosmo, z_max):
"""
BNS model calls the function to calculate the rate of BNS merger. (Eq 4 - Kuwahara et al. 2025)
:return: BNS merger rate R_m(z), in (M_sol)yr^(-1)Gpc^(-3)
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

initialization do not have any returns.
You should make the two definitions to provide the same units. It's very dangerous to have different returns for different settings

else:
raise ValueError("model should be chosen from 'BNS' or 'SNIa'")

def calculate_event_rate(self, z):
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this should be in the docstrings of this definition, not the init statement

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants