Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Absorption with cubepy #133

Open
wants to merge 26 commits into
base: master
Choose a base branch
from

Conversation

pawel21
Copy link
Collaborator

@pawel21 pawel21 commented Jan 10, 2023

Hi,

I have created PR for my modification of code for absorption using cubepy.
I added example notebook in agnpy/docs/tutorials

@review-notebook-app
Copy link

Check out this pull request on  ReviewNB

See visual diffs & provide feedback on Jupyter Notebooks.


Powered by ReviewNB

@pawel21
Copy link
Collaborator Author

pawel21 commented Mar 12, 2023

@cosimoNigro @jsitarek could you look on my PR? :)

@cosimoNigro
Copy link
Owner

I will @pawel21, sorry for the delay. I had little time in general and an issue to fix before this PR.
Will try in the next two weeks.

agnpy/absorption/absorption.py Outdated Show resolved Hide resolved
agnpy/absorption/absorption.py Outdated Show resolved Hide resolved
agnpy/absorption/absorption.py Show resolved Hide resolved
agnpy/absorption/absorption.py Outdated Show resolved Hide resolved
docs/tutorials/absorption_in_BLR_with_cubepy.ipynb Outdated Show resolved Hide resolved
"source": [
"SUM_RATIO_LINES_BLR = 30.652\n",
"\n",
"def get_absor_blrs(E, L_disk, R_line, r_blob_bh, z):\n",
Copy link
Collaborator

Choose a reason for hiding this comment

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

this function would be useful to be included in the main code of agnpy
XI_LINE and eps_abs should be then among the function parameters

@jsitarek
Copy link
Collaborator

sorry for a very late feedback on this
@pawel21 - I left a few comments on the PR
@cosimoNigro - how do you prefer to proceed - should this function substitute the current integration? or should it be added as a separate function (like Paweł did now)? In the latter case a test comparing the two methods should be added.
The former case might be more reasonable, because this method is more precise and do not have those crazy memore requirements of the standard integration scheme.

@jsitarek
Copy link
Collaborator

jsitarek commented Jul 6, 2023

digging through my old e-mails I realized that this PR is still pending.
@pawel21 I think that still 2 of my comments are not implemented/responded
@cosimoNigro any feedback about the general implementation and relation to the older code?

@cosimoNigro
Copy link
Owner

Sorry for letting it slip, will do the review later today.

Copy link
Owner

@cosimoNigro cosimoNigro left a comment

Choose a reason for hiding this comment

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

Thanks a lot for your work @pawel21, and apologies again if I got disconnected and took a while to review it.

I think once you fix mine and @jsitarek's comments we can merge it in.
Discussing how to handle the general approach towards integration (e.g. allowing to choose between trapezoidal and adaptive integration) is something we should discuss together at a certain point in the future, along with the re-organisation of the package, not in this PR.

# conversions
epsilon_1 = nu_to_epsilon_prime(nu, z)

def f(x):
Copy link
Owner

Choose a reason for hiding this comment

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

I think it's not a good idea to define a function within another function.
I think the nested function will be recompiled each time the outer function is called. Could you move it outside?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Done

epsilon_line,
R_line,
r,
mu=mu_to_integrate,
Copy link
Owner

Choose a reason for hiding this comment

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

the mu and phi arguments are not used in the function, you use directly the default mu_to_integrate and phi_to_integrate

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Done

"""

tau_blr_list = []
for freq in nu:
Copy link
Owner

Choose a reason for hiding this comment

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

I think the problem of using integrators as cubepy, that rely on functions, instead of trapz is that you do not get directly the output as an array but you are forced to loop over each frequency value. Anyhow, I think this is good for the moment.

Copy link
Collaborator

Choose a reason for hiding this comment

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

true, but since you do not have arrays you have much smaller memory usage, and the integration module automatically selects the points to evaluate the function in an optimal way, so a simple loop over frequencies is a small price to pay :-)

agnpy/absorption/absorption.py Show resolved Hide resolved
@@ -0,0 +1,275 @@
{
Copy link
Owner

@cosimoNigro cosimoNigro Jul 16, 2023

Choose a reason for hiding this comment

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

Line #1.    SUM_RATIO_LINES_BLR = 30.652

If this notebook has to be included in the docs I think it needs to be considerably improved, at the moment it has no commentary and it does not explain what its purpose is.


Reply via ReviewNB

Copy link

Choose a reason for hiding this comment

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

Let's remove the notebook for now so we can go further. I would also like to develop other features, but it will be easier if the code for calculation absorption with cubepy is in the main repo. In the next step, I will add notebooks with explanations for the multi-shell BLR model.

@cosimoNigro
Copy link
Owner

Another thing that is missing is an absorption function returning the actual attenuation factor in the cubepy case, now only tau is computed.

@pgliwny
Copy link

pgliwny commented Oct 2, 2024

Another thing that is missing is an absorption function returning the actual attenuation factor in the cubepy case, now only tau is computed.

What do you mean? Which function is missing? I tried to add a function similar to a function for another target.

@pgliwny
Copy link

pgliwny commented Oct 2, 2024

@jsitarek @cosimoNigro Could you please review my code?
I have addressed all of your suggestions and modified the code accordingly. I will be using this code to model absorption for new source in near future and to add a fit that allows for the consideration of multiple shells included in the BLR. For this purpose, it would be easier to have this code in the main repository for further improvements.

@jsitarek
Copy link
Collaborator

jsitarek commented Oct 3, 2024

Hi @pgliwny
The automatic CI checks fail:

FAILED agnpy/fit/tests/test_fit.py::TestFit::test_mrk_421_fit - IndexError: boolean index did not match indexed array along dimension 2; dimension is 1 but corresponding boolean dimension is 9
FAILED agnpy/fit/tests/test_fit.py::TestFit::test_pks_1510_fit - ValueError: non-broadcastable output operand with shape (6,1,1) doesn't match the broadcast shape (6,1,6)

maybe your PR changed something in the API

EDIT: I had a deeper look and the new code does not modify any of the existing functions so it should not be the reason for the fail of the tests.

Returns
-------
float
The value of the integrand for the given parameters.
Copy link
Collaborator

@jsitarek jsitarek Oct 3, 2024

Choose a reason for hiding this comment

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

it would be useful to specify in which units the integrand is specified (I think it should be cm)

**note** these are observed frequencies (observer frame)
z : float
redshift of the source
mu_s : float
Copy link
Collaborator

Choose a reason for hiding this comment

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

I think this function is following evaluate_tau_blr that is only valid for mu_s=1 (otherwise integration over path of the photon and over the distance from the black hole is not equivalent), so it should not have it as a parameter (@cosimoNigro can you also have a look at this?)

eps_abs: float, optional
Absolute precision for the calculation, default is 1e-6
"""
SUM_RATIO_LINES_BLR = 30.652
Copy link
Collaborator

Choose a reason for hiding this comment

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

hardcoding this is not safe, you should sum up the lines from the dictionary that you will use

XI_TARGET_LINE = self.target.xi_line
tau_z_all = np.zeros(nu.shape)

blr_shells = dict()
Copy link
Collaborator

@jsitarek jsitarek Oct 3, 2024

Choose a reason for hiding this comment

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

why do you define those dictionaries? I guess it is a remnant of a script where you were also plotting the effect of individual lines, but this is either way not possible this function.

here you are using the dictionaries only to sum up over all the lines, so you could instead created SphericalShellBLR and Absorption as temporary variables directly in the for loop

and applies the computation to all the available BLR lines (shells) using
their respective ratios.

Parameters
Copy link
Collaborator

Choose a reason for hiding this comment

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

how about adding a dictionary 'lines' to the parameters (and setting it by default to lines_dictionary)
this would allow to run the function with limited number of lines (selecting only the most important ones for speed) or with modified dictionary (if we have more accurate measurements for a particular source). It would also facilitate tests because you could then implement a quick test on a dictionary of 2 lines

Copy link

Choose a reason for hiding this comment

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

Great idea, but what do you think about doing it in the "second" iteration (after merging this code with all lines into the main branch)?

Copy link
Collaborator

Choose a reason for hiding this comment

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

if you prefer ok, but note that this is really minor change to the code, just adding a new parameter to the function, with a default value and then in the line 637 using this parameter instead of lines_dictionary.
And it would simplify implementation of my comment about tests

@@ -228,6 +228,31 @@ def test_abs_dt_vs_point_source(self):
# requires a 10% deviation from the two SED points
assert check_deviation(nu, tau_dt, tau_ps_dt, 0.1)

def test_abs_blr_cubepy(self):
Copy link
Collaborator

Choose a reason for hiding this comment

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

the test as it is written now does not test much: only if you do not have a crash because of not matching API, or some internal inconsistency in the code.
Instead you could run a test of comparing tau_blr_cubepy output with tau_blr

Copy link
Collaborator

@jsitarek jsitarek left a comment

Choose a reason for hiding this comment

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

I left some comments. In particular the test should be rewritten.

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.

4 participants