PostProcessing: Adds 'CHE_BBH' and 'NON_CHE_BBH' to selectable DCO types for ClassCOMPAS

[jeffriley]

PostProcessing:

Sometimes users may want to separate CHE/non-CHE BBHs - this change adds 'CHE_BBH' and 'NON_CHE_BBH' to selectable DCO types for ClassCOMPAS DCOmask (see setCOMPASDCOmask() in ClassCOMPAS.py).

[TomWagg]

Nice idea to add this! 😄

Left 2 small inline comments (one just a suggestion, the other a possible cleaner way of writing it)

More generally (and apologies if I am not thinking through the combination of conditions correctly) but I think there is an easier cleaner to make the masks than repeating the stuff (but the opposite conditions) in L86-92.

You could instead just find the CHE mask (currently L78-84). Then set the additional type masks after line 105 with

type_masks["CHE_BBH"] = np.logical_and(np.in1d(dco_seeds, che_seeds), type_masks["BBH"])
type_masks["NON_CHE_BBH"] = np.logical_and(np.logical_not(np.in1d(dco_seeds, che_seeds)), type_masks["BBH"])

I'm also not as familiar with CHE as you but could it also apply to other DCOs? If so it might be nice to just save the pure CHE mask that the user could combine with the others. Something like

self.CHE_mask = np.in1d(dco_seeds, che_seeds)

[TomWagg]

Possible bug: Did you mean to do np.logical_and(stellar_type_1_zams, 16) instead of just stellar_type_1_zams == 16? Apologies if I am misunderstanding.

Either way I think you could write these two lines as just which reads a bit more clearly to me.

che_mask = np.logical_and.reduce((stellar_type_1_zams == 16,
                                  stellar_type_2_zams == 16,
                                  che_ms_1 == True,
                                  che_ms_1 == True))

[jeffriley]

You're not misunderstanding. I have no idea why I did that :-) I think I must have been cleaning up the code and copy pasted without paying attention. That check (stellar_type_x_zams == 16) isn't strictly required - in the current model the che_ms_x flag only gets set for stars that were type 16 at zams (which is why I didn't notice - in the test I did that check made no (discernable) difference). Yes, your suggestion looks cleaner - I'll change it.

[TomWagg]

No strong opinion here but you don't need this condition here.

I guess you're doing this to speed up the function by not calculating the CHE stuff when you don't want to use it for the DCO mask? The alternate argument is that if you calculate it up front now you'll have all of the masks available in the class in case you want to use it later and then won't have to re-run this mask function (you could instead do self.DCOmask = self.CHE_BBHmask for example).

But as I say, no strong opinion, either works for me.

[jeffriley]

I guess you're doing this to speed up the function by not calculating the CHE stuff when you don't want to use it for the DCO mask?

Actually, no - the check is needed (or some other coding is required). The che_ms_x flags are not in the output files by default (maybe they should be now that CHE is enabled by default - I may open an issue for that), so get_COMPAS_variables() will fail if it is called and those columns aren't present in the file (we could put exceptions in get_COMPAS_variables(), but we'd still need to deal with those columns not being there). The check is so that we only try to get those columns if CHE is requested. The columns may still not be in the file, but there is more chance of them being there if the user is requesting CHE BBHs. Maybe we should be more defensive and test for all requested columns in get_COMPAS_variables() (via try/except) and return an error on failure...

[TomWagg]

Ahhh gotcha, I hadn't realised there weren't a default, makes sense in that case! We could definitely add some try/excepts, though I find the hdf5 errors are usually fairly clear about what it missing?

[jeffriley]

Oh yeah, the hdf5 errors are clear, but it's an inelegant way of failing. If we all had lots of time to devote to coding things correctly we'd control the code flow (i.e. trap the error, print an error message, and decide what to do - we might just fail, proceed without CHE, or maybe ask the user what they wanted to do (in this case, probably just fail)). It's a bit like 'cp' or 'ls' on Linux failing and printing a stack trace because a file doesn't exists - instead of just printing a simple error message. But we don't have lots of time, so as long as people understand what the problem is when it fails it's ok (just inelegant :-)

[jeffriley]

Thanks Tom! I'm sure you're right on all counts :-) I'll take a look later - busy this morning.

[jeffriley]

You could instead just find the CHE mask (currently L78-84). Then set the additional type masks after line 105 with

type_masks["CHE_BBH"] = np.logical_and(np.in1d(dco_seeds, che_seeds), type_masks["BBH"])
type_masks["NON_CHE_BBH"] = np.logical_and(np.logical_not(np.in1d(dco_seeds, che_seeds)), type_masks["BBH"])

The reason I separated them was that che_seeds is not available if "CHE_BBH" was not specified by the user (same reason for the addition of the "if types == "CHE_BBH" else np.repeat(False, len(dco_seeds))" construct to then end of some of the statements).

However, your suggestions are good and they have prompted me to rewrite the code (I think) more cleanly. I confess to not paying as much attention to clean and efficient code for the python post-processing code as I do to the C++ code - thanks for prompting me.

I'm also not as familiar with CHE as you but could it also apply to other DCOs? If so it might be nice to just save the pure CHE
mask that the user could combine with the others. Something like

self.CHE_mask = np.in1d(dco_seeds, che_seeds)

Well, I've never looked for anything other than BBH, but I don't know that there's any reason that they couldn't exist, so I've included your suggestion.

I think I have everything working as expected, but you should test it :-)

Thanks again for your suggestions - the code is better for them.

[TomWagg]

This looks great @jeffriley and worked perfectly for me!

Glad the suggestions were useful 😄