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

New framework for model atmospheres #1007

Open
wants to merge 27 commits into
base: feature-blending
Choose a base branch
from

Conversation

aprsa
Copy link
Contributor

@aprsa aprsa commented Feb 15, 2025

This huge PR incorporates a new model atmosphere framework. It is largely but not completely finished, and needs further polish. Several tests are failing but the main functionality has been restored. The rest will go in from spinoff PRs.

Reviewers: please don't refrain to comment your hearts out -- best to work out any framework issues now.

aprsa and others added 21 commits January 8, 2025 22:00
This is a major rewrite of the logic that supports model atmosphere definitions. The parent class, ModelAtmosphere, has been thoroughy documented. Several derived classes are implemented that employ the new concept.

Things to do still:
* final code clean-up
* update node list on impute
Imputation of stellar atmospheres should be out of scope for phoebe itself. The capability to impute model atmospheres has been migrated to a standalone code that should be used to update model atmosphere files, not phoebe internals.
For missing values that lie within the convex hull, it is sometimes necessary to add nodes to model atmosphere axes (for example, teff=5000K for phoenix). An add_axis_node() method was added to the ModelAtmosphere class that does this.
* The logic that pertains to model atmospheres is now in phoebe/atmospheres/models.py.
* Model atmosphere instantiation from fits files is moved to a class method from_path().
* Class attributes of ModelAtmosphere and its subclasses now suffice for passband handling.
* A master list of model atmospheres is in _atmtable.
* pb.save() has been updated to use the new ModelAtmosphere classes.
I'm not too happy with how associated axes are treated and may need to revisit it further. The next step is to update all passband functions to use atm.basic_axes instead of predefined teff/logg/abun arrays. That will likely shed additional light on how to best treat associated axes.
Providing associated axes to class-level meta makes little sense because it is table-specific. Thus, all associated axes have been removed from model atmosphere definitions.
* Dave reported that limb extrapolation sometimes failed; this was due to Inorms being extrapolated to negative values. This is now fixed.
* passband load/save methods fixed to properly account for the removal of associated axes from the model atmosphere class.
* Tweaks to pb.save/load and adding TMAP atm_models

* Added wls_file option to define name of npy wavelength file
Ndpolator meta tables were switched from tuples to dictionaries; this commit updates the code to conform to the upstream change. Version 1.2.2 is explicitly requested.
This is a simplified version of @Raindogjones' PR in addition to the number of other bugfixes.
* model atmospheres for 'extern_planckint' and 'extern_atmx' implemented;
* save() method simplified to export a pre-2.5-compatible version on demand;
* code updated further to the new model atmosphere framework.

After these changes, 9/61 tests still fail and bugs still need to be tracked down and fixed. But we're getting there.
This causes issues with tests because overloaded query points are not handled yet in passbands.py so expect breakage.
* Fix to photon-weighted ptf_area

* Fix to energy-weighted ptf_area
… new framework

The bug in Inorm caused log-intensities instead of intensities to be returned; fixed. This also exposed a difference between the original and the new blackbody treatment that needs to be further investigated.

The same difference causes the discrepancy between interpolated and rigorously computed Sun-Earth intensities to creep up a little bit, so the tolerance had to be increased from 0.1% to 1%, at least temporarily.

Several code tweaks had been done to further adapt the code to the new model atmosphere framework.
@aprsa aprsa self-assigned this Feb 15, 2025
@kecnry
Copy link
Member

kecnry commented Feb 16, 2025

looks like all tests are failing at the same point - is that an easy fix?

@aprsa
Copy link
Contributor Author

aprsa commented Feb 16, 2025

looks like all tests are failing at the same point - is that an easy fix?

Not really, no. It will take further work to update everything to the new framework. This was a trade-off between a large and a larger PR.

__version__ = '2.4.17.dev+feature-blending'
__version__ = '2.5.0.alpha'
Copy link
Member

Choose a reason for hiding this comment

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

remind me why this is necessary? Can we not stick with the existing version tagging scheme for feature branches?

Copy link
Contributor Author

@aprsa aprsa Feb 16, 2025

Choose a reason for hiding this comment

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

because 2.4 and 2.5 store passband data differently, and version triggers different pathways in the code (i.e., how passband files are loaded). See here for an example:

https://github.com/aprsa/phoebe2/blob/ac7e009daf9fec24e0bd9f736a5863f112d55cfc/phoebe/atmospheres/passbands.py#L555

Copy link
Member

Choose a reason for hiding this comment

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

Can that check for the version tag until we know this will officially be in 2.5 and merge into the release branch?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Fixed by commit aprsa@a7ec893.

Copy link
Member

@kecnry kecnry left a comment

Choose a reason for hiding this comment

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

Started with a quick skim through the resulting changes to univers/bundle/compute and luckily those are quite minimal and make sense to me. I'll have to dig through the new infrastructure separately.

Copy link
Contributor

@Raindogjones Raindogjones left a comment

Choose a reason for hiding this comment

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

Difficult to properly review while tests are failing, but I don't see any issues with the new framework. I still think that the major hurdle to any new atmosphere is preparing the required input files rather than processing them, but this should ultimately make it easier to support atmospheres with more/fewer than 3 basic axes.

if f'{atm.name}:ext' in self.content:
if export_to_pre25 and atm.name == 'blackbody':
data.append(fits.ImageHDU(self.ndp[atm.name].table['ext@energy']['grid'], name=f'{atm.prefix.upper()}EGRID'))
data.append(fits.ImageHDU(self.ndp[atm.name].table['ext@photon']['grid'], name=f'{atm.prefix.upper()}PGRID'))
Copy link
Contributor

Choose a reason for hiding this comment

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

Not sure I understand this. Aren't blackbody extinction grids also called BBXEGRID and BBXPGRID?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

they are now, but not in the pre-2.4.17 passbands.

if ld_func == 'interp':
if atm == 'blackbody' and 'blackbody:Inorm' in self.content and ldatm in ['ck2004', 'phoenix', 'tmap_sdO', 'tmap_DA', 'tmap_DAO', 'tmap_DO', 'tremblay']:
if atm.name == 'blackbody' and 'blackbody:Inorm' in self.content and ldatm.name in ['ck2004', 'phoenix', 'tmap_sdO', 'tmap_DA', 'tmap_DAO', 'tmap_DO', 'tremblay']:
Copy link
Contributor

Choose a reason for hiding this comment

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

Shouldn't we try to avoid these lengthy lists of atmospheres? Perhaps just removing the unwanted atms from _atms_supported?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

yep, this is WIP. This simplifies to:

if ... and hasattr(ldatm, 'mus') ...

prefix = 'bb'

basic_axis_names = ['teffs']
teffs = 10**np.linspace(2.5, 5.7, 97) # this corresponds to the 316K-501187K range.
Copy link
Contributor

Choose a reason for hiding this comment

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

Increase the number of points to try and maintain consistency with previous implementation

Copy link
Member

Choose a reason for hiding this comment

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

Can we use np.geomspace or np.logspace?

Copy link
Contributor

@MichaelAbdul-Masih MichaelAbdul-Masih left a comment

Choose a reason for hiding this comment

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

all looks good from my side. I wasn't too familiar with the previous setup, but this current way of doing things seems straight forward enough to me.

Copy link
Member

@kecnry kecnry left a comment

Choose a reason for hiding this comment

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

I like the overall clean design of the atmosphere subclasses, which is a sign that the design is at least in the right direction 😉

I do think we shouldn't invent a new versioning system here, if at all possible, but I think the rest of my comments are just minor suggestions that I'll let you decide whether to adopt or not (or in some cases defer to follow-up tickets - just please create them if you agree we should do them eventually so we don't forget).

Comment on lines +9 to +10
A parent class for handling model atmosphere data. Please note that only
derived classes should be instantiated.
Copy link
Member

Choose a reason for hiding this comment

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

could we enforce this? Perhaps with a NotImplementedError in the init or something?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Not directly, because subclasses call the parent's init method.

Copy link
Collaborator

Choose a reason for hiding this comment

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

There are AbstractBaseClasses which throw errors when instantiating

Copy link
Member

Choose a reason for hiding this comment

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

Easiest way would be to check __class__ and raise an error if it's still the base class, any subclasses would then not raise an error

Copy link
Member

Choose a reason for hiding this comment

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

(but also not a big deal, we just don't expose the API for this in the docs and it won't be a problem)

Comment on lines +30 to +32
* `units` (float): intensity unit conversion factor. The intensities are
usually given in erg/s/cm^2/A, which is converted to W/m^3 by
multiplying with this factor.
Copy link
Member

Choose a reason for hiding this comment

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

would we ever want to just support inputting quantities and determining this automatically?

# default axes:
basic_axis_names = ['teffs', 'loggs', 'abuns']

def __init__(self, basic_axes=None, from_path=False):
Copy link
Member

Choose a reason for hiding this comment

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

might be able to get around the from_path thing by defining __new__ with most of the logic below and sharing that between __init__ and from_path. Since this is internal though, that probably could be done later down the road if you want to create a ticket and defer until later.

# function, basic_axes will be automatically determined from
# axis definitions via class attributes.

if hasattr(self, 'intensity') and callable(self.intensity):
Copy link
Member

Choose a reason for hiding this comment

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

this might be cheaper (I'm actually not sure - take it or leave it 🤷‍♂️ )

Suggested change
if hasattr(self, 'intensity') and callable(self.intensity):
if callable(getattr(self, 'intensity', None)):

Comment on lines +95 to +97
else:
if basic_axes is None:
raise ValueError('basic_axes must be defined.')
Copy link
Member

Choose a reason for hiding this comment

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

Suggested change
else:
if basic_axes is None:
raise ValueError('basic_axes must be defined.')
elif basic_axes is None:
raise ValueError('basic_axes must be defined.')

Comment on lines +168 to +169

return NotImplementedError
Copy link
Member

Choose a reason for hiding this comment

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

Suggested change
return NotImplementedError
return NotImplementedError("parse_rules must be implemented by the subclass")

Comment on lines +175 to +176
spun by the current axes. If you do not know why you would
use this method, you likely should not use it.
Copy link
Member

Choose a reason for hiding this comment

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

😂

prefix = 'bb'

basic_axis_names = ['teffs']
teffs = 10**np.linspace(2.5, 5.7, 97) # this corresponds to the 316K-501187K range.
Copy link
Member

Choose a reason for hiding this comment

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

Can we use np.geomspace or np.logspace?

Comment on lines +252 to +254
def __init__(self, *args, **kwargs):
super().__init__(*args, **kwargs)

Copy link
Member

Choose a reason for hiding this comment

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

Doesn't hurt, but could save some lines by removing unless you think we'll want to override in the near future? (Same goes for other subclasses)

Suggested change
def __init__(self, *args, **kwargs):
super().__init__(*args, **kwargs)

Comment on lines -158 to +151
Step #3: compute Castelli & Kurucz (2004) intensities. To do this, the
tables/ck2004 directory needs to be populated with non-filtered
intensities available for download from %static%/ck2004.tar.
Step #3: instantiate a model atmosphere object and compute intensities:

```py atmdir =
os.path.abspath(os.path.join(os.path.dirname(os.path.abspath(__file__)),
'tables/ck2004')) pb.compute_ck2004_response(atmdir) ```
```py atm = CK2004ModelAtmosphere('path/to/ck2004')
pb.compute_intensities(atm) ```

Step #4: repeat step #3 for other model atmospheres.

Step #4: -- optional -- import WD tables for comparison. This can only
Step #5: -- optional -- import WD tables for comparison. This can only
Copy link
Member

Choose a reason for hiding this comment

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

do we have a follow-up ticket to update the passbands tutorial?

It turns out that requiring the version to be < 2.4.17.dev+feature-blending does not suffice because 2.4.17 fails that test. Thus, "<" is now changed to "!=", and will be reverted to 2.5.0 once released.
Instead of passing query_pts, the framework now relies on query_table that packs query_pts alongside column names. Each function in passbands.py knows what columns it needs to do its thing, so we can pack anything and everything necessary into query_tables, and respective functions will pick out what they need. This allows us to combine different atmospheres and it greatly simplifies the logic of existing functionality.

A multitude of bugs has been fixed. The only tests currently failing are extinction tests. That, and the recomputation of blackbody passband tables, is next.
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.

5 participants