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

Updated piecewise orbital model #1545

Merged
merged 13 commits into from
Sep 28, 2023

Conversation

poneill129
Copy link
Contributor

Tried merging these changes into the existing PR but looks like its opened as a new one, apologies!
Tagged reviewers and commenters so you don't miss it. (@scottransom @aarchiba @dlakaplan)

Had a chance to address Anne's comments on the piecewise orbital model.

I will add the responses to the comments in the original PR #1378 as well as here in a follow-up comment.

@poneill129 poneill129 changed the title Updated to reflect comments Updated piecewise orbital model Mar 31, 2023
@codecov
Copy link

codecov bot commented Mar 31, 2023

Codecov Report

Patch coverage: 67.96% and no project coverage change.

Comparison is base (25f38ca) 68.27% compared to head (484508e) 68.28%.

Additional details and impacted files
@@           Coverage Diff            @@
##           master    #1545    +/-   ##
========================================
  Coverage   68.27%   68.28%            
========================================
  Files         100      101     +1     
  Lines       23111    23596   +485     
  Branches     4010     4160   +150     
========================================
+ Hits        15779    16112   +333     
- Misses       6342     6427    +85     
- Partials      990     1057    +67     
Files Changed Coverage Δ
src/pint/models/binary_bt.py 66.79% <64.75%> (-21.21%) ⬇️
...nt/models/stand_alone_psr_binaries/BT_piecewise.py 71.07% <71.07%> (ø)
src/pint/models/__init__.py 100.00% <100.00%> (ø)

... and 3 files with indirect coverage changes

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

@dlakaplan
Copy link
Contributor

Can you share in the docstring and/or here how such a model is defined?

@scottransom
Copy link
Member

This looks great, @poneill129. Once the docstring is added to as requested by David K, and the minor conflict in src/pint/models/__init__.py is fixed, do you think this is basically ready to merge? My student is wanting to test/use it and it would be better to do that (I think) if it is merged than trying to work in a side branch.

@poneill129
Copy link
Contributor Author

This looks great, @poneill129. Once the docstring is added to as requested by David K, and the minor conflict in src/pint/models/__init__.py is fixed, do you think this is basically ready to merge? My student is wanting to test/use it and it would be better to do that (I think) if it is merged than trying to work in a side branch.

Hi, the docstring turned into a more involved problem than anticipated. Trying to replicate the convention for other standalone binary models I found myself in a rabbit hole.

As part of the model declaration an arbitrary piece and associated range are added. Then attempting to add a piece directly to the model (not via the pint facing binary_piecewise module) is bugged in that the dict that is used in the update_binary_object() causes an error when adding a piece.

I will need to find a way of producing the bug to report it, I will emphasise that this does not happen when piecwise model is declared via the wrapper module. Effectively it’s a timing model that I can only make work via the pint interface.

@scottransom
Copy link
Member

@poneill129 Coming back to this again... Is the current status still accurate based on your last comment above with an error when you try to add a piece? If so, might you be able to email me with some example code as to how/when that is happening? And also, what is your "wrapper" module that allows things to work? Is there an example usage of that? I'm going to start playing with this this week as two of my students really do need it starting now-ish.

@poneill129
Copy link
Contributor Author

Hi @scottransom, yes I think I've cracked it! I'll email you a demo script asap and stage my changes.

Just to keep it clear, when I refer to a stand alone model I'm talking about BT_model.py/BT_piecewise.py, and when I say timing model I mean "obtained from get_model()".

So my issue was in documenting the piecewise stand alone model in the style of the BT model, I tried following the steps albeit with the piecewise stand alone model, not realising update_object(t,updates) is not how the arguments are passed by the timing model (anymore?). I needed to step away at that point because I could not see how the binary_instance was "noticing" the toas while the piecewise variables worked fine during a "get_model" or residual calculation. Can think of this as the stand alone model is in a state of waiting for toas to allocate to its predefined pieces.

I thought this was an indication that there was some foundational oversight that meant that I couldn't operate on it as a stand alone model. But looking over it with fresh eyes I've found that this was not the case, I did catch a minor bug in the toa-piece allocation step (never was relevant when calling using a timing model but popped up during this exploration). Just folding this into the documentation now, along with the recommended changes to "XR_0000" to denote piece ranges. I suppose if others agree that the steps included in the documentation for BT_model.py are for an old update_input method there's something else that can be changed.

self.remove_param("SINI")
self.T0.value = 1
self.A1.value = 1
self.add_group_range(50000, 51000, frozen=False, j=0)
Copy link
Contributor

Choose a reason for hiding this comment

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

Can this be initialized with None instead of random MJDs?

@dlakaplan
Copy link
Contributor

Hi, @poneill129 . Great that you made progress. I started to look at your code but I realize you may have some updates waiting to go. Should I keep reviewing things now or wait?

@poneill129
Copy link
Contributor Author

Hi @dlakaplan, The proposed change is only very small, affecting only the "toa_belongs_in_group" function in BT_piecewise.py. I'm still trying to pull things together. Thank you for the comments you've already given, if you want to keep going (avoiding toa_belongs_in_groups) I'll pull the changes in while I go!

@@ -21,6 +21,7 @@
# Import all standard model components here
from pint.models.astrometry import AstrometryEcliptic, AstrometryEquatorial
from pint.models.binary_bt import BinaryBT
from pint.models.binary_piecewise import BinaryBTPiecewise
Copy link
Contributor

Choose a reason for hiding this comment

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

This is fine, but you could also consider defining BinaryBTPiecewise in binary_bt.py (like how dispersion_model.py has DispersionDM and DispersionDMX

Copy link
Contributor

Choose a reason for hiding this comment

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

(the standalone model piece should be in a separate file)

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I've shifted the binary_piecewise class into BTmodel now, made changes in init to accommodate

)
)

def lock_groups(self):
Copy link
Contributor

Choose a reason for hiding this comment

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

Is this needed? Is it used?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

So I found myself needing this function as validates/setups cannot be run before all pieces and corresponding group boundaries are declared. Otherwise there are floating piecewise parameters that don't have a group to go with them (which validate does not like). The function is intended to be be used after all add_group_range and add_piecewise_param have been called

for t0n in T0X_mapping.values():
T0Xs[t0n] = getattr(self, t0n).quantity

for t0_name, t0_value in T0Xs.items():
Copy link
Contributor

Choose a reason for hiding this comment

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

would it make more sense to do a single loop over the PLB/PUB mappings? And then for each index check whether there was a A1X or T1X value and set them if so?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Implemented the sorting algorithm described


if j_plb != j_pub:
raise ValueError(
f"Group boundary mismatch error. Number of detected upper bounds: {j_plb}. Number of detected upper bounds:{j_pub}"
Copy link
Contributor

Choose a reason for hiding this comment

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

probably include a space before the values

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Added extra space

@dlakaplan
Copy link
Contributor

You might look at toa_select.py. That has tools to do some of what you are talking about. You can see how it's used by DMX (for instance).

str elements, look like ['0000','0001'] for two TOAs where one refences T0X_0000 or T0X_0001.
"""
# asks the model what group a single TOA/list of TOAs
barycentric_toa = self._parent.get_barycentric_toas(toa)
Copy link
Contributor

Choose a reason for hiding this comment

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

I think passing a TOAs object here would be fine. then you don't have to worry about re-barycentering (which can be slow). Again, look at how dmx_dm does things.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

reduced to barycentric_toa = toa (this is one of a few convenience functions which are going to change with testing)

)
return group_indexes

def get_group_indexes_in_four_digit_format(self):
Copy link
Contributor

Choose a reason for hiding this comment

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

Is this needed? If you have the method above to return integers can't you just use:

[f"{int(...):04d" for x in ...]

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Yeah, this was an one of the functions I would use in the testing script. Does nothing essential, scrapped

)
return group_indexes

def get_T0Xs_associated_with_toas(self, toas):
Copy link
Contributor

Choose a reason for hiding this comment

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

maybe instead of having specialized methods for TOX and A1X, you can have 1 method to see which group(s) TOAs are associated with? Then the TOX and A1X can just be gotten using the existing mappings.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

implemented new function paramX_per_toa('T0',toa) which does as was described returns same as get_T0Xs_associated_with_toas. This function changes the structure of the testing suite since now a lot of convenience functions can be consolidated into this

from .binary_generic import PSR_BINARY


def extract(lst):
Copy link
Contributor

Choose a reason for hiding this comment

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

Is this needed?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Nope, and its gone now

"""Query whether a TOA/list of TOAs belong(s) to a specific group
Parameters
----------
toas :
Copy link
Contributor

Choose a reason for hiding this comment

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

try to be explicit with the types/class references. So this would be of type pint.toa.TOAs.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I left some of these unchanged in the interest of time, since some of these functions are being removed

@dlakaplan
Copy link
Contributor

Hi, @poneill129 . Overall this looks good, although I think you could simplify some of the smaller functions. But I notice that you don't have any tests yet. It would be great to have some tests like:

  • make a model with upper/lower ranges for the whole data-set, and compare the binary delays with existing BT
  • make a models divided into two ranges. make 2 BT models and compare appropriately
  • test the various setup/adding/removing functions. Make sure that errors are raised appropriately for the conditions you give (like overlapping ranges)
    You might look at some of the tests for DMX, WaveX or similar to see various examples.

@poneill129
Copy link
Contributor Author

poneill129 commented Sep 7, 2023

Hi, @poneill129 . Overall this looks good, although I think you could simplify some of the smaller functions.

Cheers again for these, and I'm working my way through them now

But I notice that you don't have any tests yet. It would be great to have some tests like:

I don't know if it's been missed or I'm misunderstanding but I should have a test file in the form of "test_BT_piecewise.py". It contains tests that might be similar to:

  • make a models divided into two ranges. make 2 BT models and compare appropriately
  • make a model with upper/lower ranges for the whole data-set, and compare the binary delays with existing BT

I'm fairly sure there's a test in there for testing the error raised when overlapping groups are declared

  • test the various setup/adding/removing functions. Make sure that errors are raised appropriately for the conditions you give (like overlapping ranges)

Regardless, I'm going to look into the test file again and identify which tests closely resemble what you've outlined, and I'll adapt what is there if I can't find an equivalent test.

You might look at some of the tests for DMX, WaveX or similar to see various examples.

Will do, thanks for the pointer!
Cheers for these again!

@dlakaplan
Copy link
Contributor

Great - I missed the tests since the browser didn't show me the code initially.

@poneill129
Copy link
Contributor Author

I've folded in a lot of the changes discussed, though in doing so, some of the functions in what was binary_piecewise.py have been made redundant. Those functions were used in the test suite which I am now working on fixing.

@poneill129
Copy link
Contributor Author

This PR is a little untested. This should be functional, and I'm trialing some demo notebooks which I will circulate for you to play around with while I address the testing suite.

Currently dealing a bug due to my current setup that means binary_piecewise still exists and get_model() calls are trying to use it rather than going via binary_bt like it should now.

@poneill129
Copy link
Contributor Author

Hi all, more updates on the piecewise model (this time on the testing suite side of things). The tests appear to pass in my local version of this, it is only in the CI tests that tests for other functions begin failing. I haven't found a way to narrow what it is that's causing it, but I have a feeling there is something lingering after the piecewise tests are done.

I don't know if its possible for this to be the reason, but if its not I cannot see why the test is failing since it shouldn't be affected by any of the changes I've made.

And to ease nervous minds, I've painstakingly looked through each of the tests, and done some quick fitting using pulsar data, to verify the model and tests are still behaving as intended (i.e. everything is not broken!)

@dlakaplan
Copy link
Contributor

Hi, Patrick. The CI failure that I saw was not something you did. The test was for something else and was close to passing, but occasionally it fails because the tolerances may be too tight. I just tried to re-run it.

@abhisrkckl
Copy link
Contributor

Hi @poneill129, from reading the comments above, it seems that this PR is pretty much ready. Are you planning to do anything more in this PR?

@poneill129
Copy link
Contributor Author

Not at the moment, happy with it. But if anything else needs clearing up I'm able to make changes.

@abhisrkckl abhisrkckl merged commit 01d3153 into nanograv:master Sep 28, 2023
8 checks passed
@abhisrkckl abhisrkckl mentioned this pull request Oct 4, 2023
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