Comparing Flavors

[jpkneller]

I've been thinking about the problem Andrey identified with the new TwoFlavor, ThreeFlavor and FourFalvor classes we want to introduce in version 2.0. The issues is with comparing `flavors' across the cases with different numbers i.e. how do you compare NU_E_BAR from TwoFlavor with NU_E_BAR from ThreeFlavor or compare NU_X from TwoFlavor with NU_MU from FourFlavor, and return True in both cases even though the numerical values of the enumerations can be different. The solution I came up with is to overload the comparison operators for the *Flavor classes. To reduce the size of the code I re-assigned the numerical values and added a couple of extra short methods. Overloading the comparison operators allows a user to compare flavors from different schemes but it doesn't solve the problem that the numerical values are different. What do people think? Is this is a problem we need to fix? Anyway, the solution is below. Now this code returns True for all cases

f1 = TwoFlavor.NU_X
f2 = ThreeFlavor.NU_TAU

print(f2==f1)
print(f1==f2)

f1 = FourFlavor.NU_TAU
f2 = TwoFlavor.NU_X

print(f2==f1)
print(f1==f2)

Here's the relevant parts of neutrino.py

class TwoFlavor(IntEnum):
    """Enumeration of CCSN Neutrino flavors."""
    NU_E = 0
    NU_E_BAR = 1    
    NU_X = 2
    NU_X_BAR = 3
    
    def to_tex(self):
        """LaTeX-compatible string representations of flavor."""
        if '_BAR' in self.name:
            return r'$\overline{{\nu}}_{0}$'.format(self.name[3].lower())
        return r'$\{0}$'.format(self.name.lower())

    @property
    def is_electron(self):
        """Return ``True`` for ``TwoFlavor.NU_E`` and ``TwoFlavor.NU_E_BAR``."""
        return self.value in (TwoFlavor.NU_E.value, TwoFlavor.NU_E_BAR.value)

    @property
    def is_neutrino(self):
        """Return ``True`` for ``TwoFlavor.NU_E`` and ``TwoFlavor.NU_X``."""
        return self.value in (TwoFlavor.NU_E.value, TwoFlavor.NU_X.value)

    @property
    def is_antineutrino(self):
        return self.value in (TwoFlavor.NU_E_BAR.value, TwoFlavor.NU_X_BAR.value)
        
    def __eq__(self, other):
        if isinstance(other,ThreeFlavor) == True or isinstance(other,FourFlavor) == True:
            return other == self

            
class ThreeFlavor(IntEnum):
    """Enumeration of CCSN Neutrino flavors."""
    NU_E = 0
    NU_E_BAR = 1
    NU_MU = 2
    NU_MU_BAR = 3
    NU_TAU = 4
    NU_TAU_BAR = 5
    
    def to_tex(self):
        """LaTeX-compatible string representations of flavor."""
        if '_BAR' in self.name:
            return r'$\overline{{\nu}}_{0}$'.format(self.name[3].lower())
        return r'$\{0}$'.format(self.name.lower())

    @property
    def is_electron(self):
        """Return ``True`` for ``ThreeFlavor.NU_E`` and ``ThreeFlavor.NU_E_BAR``."""
        return self.value in (ThreeFlavor.NU_E.value, ThreeFlavor.NU_E_BAR.value)
        
    @property        
    def is_muon(self):
        """Return ``True`` for ``ThreeFlavor.NU_MU``, ``ThreeFlavor.NU_MU_BAR``."""
        return self.value in (ThreeFlavor.NU_MU.value, ThreeFlavor.NU_MU_BAR.value)

    @property
    def is_tauon(self):
        """Return ``True`` for ``ThreeFlavor.NU_TAU``, ``ThreeFlavor.NU_TAU_BAR``."""
        return self.value in (ThreeFlavor.NU_TAU.value, ThreeFlavor.NU_TAU_BAR.value)        

    @property
    def is_neutrino(self):
        """Return ``True`` for ``ThreeFlavor.NU_E`` and ``ThreeFlavor.NU_MU`` and ``ThreeFlavor.NU_TAU``."""
        return self.value in (ThreeFlavor.NU_E.value, ThreeFlavor.NU_MU.value, ThreeFlavor.NU_TAU.value)

    @property
    def is_antineutrino(self):
        return self.value in (ThreeFlavor.NU_E_BAR.value, ThreeFlavor.NU_MU_BAR.value, ThreeFlavor.NU_TAU_BAR.value)
        
    def __eq__(self, other):
        if isinstance(other,FourFlavor) == True:
            return other == self    
    
        f = self.value 

        if isinstance(other,TwoFlavor) == True and self.is_tauon == True:
             f -= 2
          
        if f == other.value:
            return True
        else:
            return False


class FourFlavor(IntEnum):
    """Enumeration of CCSN Neutrino flavors."""
    NU_E = 0
    NU_E_BAR = 1
    NU_MU = 2
    NU_MU_BAR = 3    
    NU_TAU = 4
    NU_TAU_BAR = 5    
    NU_S = 6
    NU_S_BAR = 7
    
    def to_tex(self):
        """LaTeX-compatible string representations of flavor."""
        if '_BAR' in self.name:
            return r'$\overline{{\nu}}_{0}$'.format(self.name[3].lower())
        return r'$\{0}$'.format(self.name.lower())

    @property
    def is_electron(self):
        """Return ``True`` for ``FourFlavor.NU_E`` and ``FourFlavor.NU_E_BAR``."""
        return self.value in (FourFlavor.NU_E.value, FourFlavor.NU_E_BAR.value)
        
       
    @property        
    def is_muon(self):
        """Return ``True`` for ``FourFlavor.NU_MU``, ``FourFlavor.NU_MU_BAR``."""
        return self.value in (FourFlavor.NU_MU.value, FourFlavor.NU_MU_BAR.value)

    @property
    def is_tauon(self):
        """Return ``True`` for ``FourFlavor.NU_TAU``, ``FourFlavor.NU_TAU_BAR``."""
        return self.value in (FourFlavor.NU_TAU.value, FourFlavor.NU_TAU_BAR.value)                

    @property
    def is_sterile(self):
        """Return ``True`` for ``FourFlavor.NU_S``, ``FourFlavor.NU_S_BAR``."""
        return self.value in (FourFlavor.NU_S.value, FourFlavor.NU_S_BAR.value)                
        
    @property
    def is_neutrino(self):
        """Return ``True`` for ``FourFlavor.NU_E`` and ``FourFlavor.NU_MU`` and ``FourFlavor.NU_TAU`` and ``FourFlavor.NU_S``."""
        return self.value in (FourFlavor.NU_E.value, FourFlavor.NU_MU.value, FourFlavor.NU_TAU.value, FourFlavor.NU_S.value)

    @property
    def is_antineutrino(self):
        return self.value in (FourFlavor.NU_E_BAR.value, FourFlavor.NU_MU_BAR.value, FourFlavor.NU_TAU_BAR.value, FourFlavor.NU_S_BAR.value)
        
        
    def __eq__(self, other):
        f = self.value 
        if isinstance(other,TwoFlavor) == True and self.is_tauon == True:
             f -= 2
            
        if f == other.value:
            return True
        else:
            return False

[Sheshuk]

(edited first post, just adding the syntax highlight)

[Sheshuk]

Thanks for bringing it up! I wanted to separate the Flavor changes from #308 and to have a discussion on separate PR, but discussing it in an issue is more productive.

My thoughts on the suggested solution:
I don't think that the problem with flavors equality is what we need to solve. Most probably the user wouldn't compare the flavors directly, but provide flavors as argument to access the desired flavor, like flux[Flavor.NU_E] - so maybe all checks can be done there, keeping the Flavor schemes cleaner.

Conversion of NU_X to be NU_MU+NU_TAU etc can be done via conversion matrices, allowing to easily convert object of one flavor scheme to another.

[Sheshuk]

As the flavor classes are used practically everywhere throughout the code, I suggest we summarize their desired behavior, before we discuss actual implementation.

My view of the problem:

Now we have several Flavor schemes. Whether they're Enums, or subranges or a single Enum, or general classes.

Classes containing the Flavor information

1. Transformation matrices of course, must explicitly know their flavor scheme.
2. Flux containers (they already have a flavor attribute, which lists the available flavors)
3. Conversion matrices (TwoFlavor->ThreeFlavor, for example, etc.) should keep information of two schemes.

I'm not sure if a every FlavorTransformation needs to know it (trivial example - NoTransformation class).

Also I don't think the SupernovaModel should keep it's output flavor scheme - it can just write it in the calculated flux.

Use cases

These Flavor schemes are to be accessed by user in two main cases:

1. Loop over all flavors in a given object:

for flavor in flux:
    plot_flavor(flux.energy, flux[flavor], label=flavor.to_tex())

2. Access the given flavor or flavors range:

prob_em = prob_matrix[Flavor.NU_E, Flavor.NU_MU]
flux_nux = flux[Flavor.NU_X]

and

3. Implicitly used when applying the transformation matrices, or making a flavor scheme conversion:

flux_2f = conversion_2f3f @ flux_3f #using an external conversion matrix
flux_oscillated = transformation_matrix @ flux_unoscillated

Requirements

Case 1:

We don't have any problems with using the wrong flavor scheme as long as it is stored in the object containers (and not just a global for flavor in Flavor).

❗ Requirement:
container class must provide an __iter__ method which returns its flavor list.

Case 2:

Here we want to avoid accessing the wrong flavor. So if user provides a wrong flavor we can either

• Raise an exception
• Make a conversion to the requested scheme (and write a warning that we did that)

❗ Requirement:
container class makes a check/conversion inside __getitem__ method

Case 3:

We need to check if the flavor schemes of the objects are compatible - similar as above.

❗ Requirement:
container class makes a check/conversion inside used method (like TransformationMatrix.__matmul__ method in this example)

[jpkneller]

I'm leaning towards making the conversion in limited circumstances when the flavor is being used as a key and the numerical value is immaterial. The circumstance where the conversion should be made is A) when there are equivalent keys in the different schemes (e.g. NU_E_BAR since it is common to all of them) and B) when the conversion is unambiguous e.g. the user passed in the flavor NU_MU from the ThreeFlavor scheme but the container uses TwoFlavor so the conversion is NU_MU -> NU_X. I haven't thought about whether this solution means a lot of new code but if it does, raise the exception.

[JostMigenda]

Before I look at the proposed implementation in #324, I wanted to write up my thinking at a very high level.

The key feature of snewpy is providing access to lots of different supernova models (and flavor transformations) with a single consistent interface:

• Users do not have to worry about what format some model’s input files are in. We provide the code to load data from different file formats and transform/interpolate as necessary.
• Users do not have to worry about whether some particular simulation distinguished between NU_X and NU_X_BAR or not. Right now, for models that don’t we automatically add a NU_X_BAR component while initializing the model, so any usage is independent of that implementation detail.
• Having this consistent interface that abstracts away implementation details also makes it easy to integrate snewpy as a dependency in other software packages like sntools, ASTERIA and others, which multiplies its impact.

Now, whether any particular simulation distinguished between NU_MU and NU_TAU or lumped them both together as NU_X is also an implementation detail. Users should not have to worry about that; and anytime some user code works for a TwoFlavor scheme but not for a ThreeFlavor scheme (or vice versa), I would consider that to be a serious bug.

In contrast, whether or not a simulation (or flavour transformation) includes sterile neutrinos is an actual physics difference; it is not an implementation detail. Ideally, we should still try to find a solution that ensures that code written for Two/ThreeFlavor schemes also works for FourFlavor schemes (and vice versa); but if that is not possible in some edge cases (or only with significant tradeoffs), I would be open to having some API differences here.

[Sheshuk]

I mostly agree about the final user perspective.

Except for this statement:

anytime some user code works for a TwoFlavor scheme but not for a ThreeFlavor scheme (or vice versa), I would consider that to be a serious bug

I think that's too extreme.
In some cases the user expects to have a TwoFlavor case, and writes code specifically for that. I agree that this should not be dictated by the implementation of the model though. So I think the best solution is to allow user to be declarative - if you want to work in ThreeFlavor basis, you explicitly convert the flux to that basis, and don't care about any implementation details.
In my solution in #324 it is as simple as:

flux = model.get_flux(...)>>ThreeFlavor
#now we know flux has 3 flavors

Note that it is somewhat similar to using astropy.units: when you get, say, an energy, no need to check its unit, just request to convert it:

energy = get_energy(...)<<u.GeV
#now we know energy is in GeV

Also there is another side: the implementation of flavor schemes should make our lives as developers easy - by making the code flexible and easier to extend.

So in summary I think:

• We should not impose the one and only flavor scheme
• We should allow flavor schemes conversion - in the user interfaces
• I don't care if we make our model.get_flux output the flux in ThreeFlavor by default - but the scheme will depend on the applied FlavorTransformation anyway, leading to steriles
• I'm strongly opposed to change the code in the existing models to use three flavors scheme, as it is a lot of changes to the code that already works and is tested. Let's use conversion matrices instead.

[JostMigenda]

[Note: I’m deliberately focussing on TwoFlavor vs. ThreeFlavor here. If models/transformations involve sterile neutrinos, that is a clear physical difference and I think it’s reasonable to require the user to explicitly keep track of that. The conversion matrices look like a good approach in that case. My point here is about the most common case, without steriles, and the best user interface for that.]

ThreeFlavor is preferred over TwoFlavor and should be the default

The TwoFlavor case an approximation that’s useful to save computing resources when running SN simulations; but it’s not physically meaningful. And of course different papers/codes have different conventions on whether NU_X means (a) NU_MU or NU_TAU, (b) the sum of both, or (c) the sum of those two and their respective antineutrinos. That alone clearly indicates to me that we all think in three flavours by default (and then map our thinking onto two flavours only if forced to).

SNOwGLoBES already enforces a ThreeFlavor scheme, with separate interaction channels (and associated cross section, etc. files) for NU_MU and NU_TAU.

In some cases the user expects to have a TwoFlavor case, and writes code specifically for that.

Users may “expect” a TwoFlavor case because that was the convention in some particular model paper, sure. But in that same sense, they might “expect” a “OnePointFiveFlavor” scheme (NU_E, NU_E_BAR and NU_X only), too. Or they might expect NU_X to mean either (b) or (c), while we use (a). So we will almost invariably break some initial user expectation.

The best thing we can do to help users recover from breaking these initial expectations is to provide a consistent interface. And as explained above, the ThreeFlavor scheme

If the TwoFlavor case is all we have, that’s one thing; but once we have a ThreeFlavor case, I see no reason why users would explicitly want to go back to the TwoFlavor case; just like there was no one who ever requested to go back to the OnePointFiveFlavor case.

Having optional support for TwoFlavor is harmful

I think there are several good reasons not to support the TwoFlavor case:

• The TwoFlavor case is a very common source of confusion or bugs:

  • Like you, I often have to double-check which of the conventions (a), (b) or (c) is the correct one in any given context.
  • There’s even an issue with the example code here: The quantity plotted as the TwoFlavor total is physically meaningless; since we use convention (a), the NU_X and NU_X_BAR flux should be weighted by a factor of 2 to get the same result as the ThreeFlavor total. (I’m not blaming you—this stuff is confusing and easy to forget. But if we developers of snewpy struggle with this, how do users stand any chance?)
    • Related: Almost every time I work with project students and show them FlavorTransformation.ipynb, I get questions about why transformations look like they don’t conserve neutrino number. It’s exactly because of this weighting factor.

• If we can convert all existing models to ThreeFlavor during initialisation (see “Note on implementation” below for how to do that), then we would not have to handle the 2to3 (or 3to2) case in conversion matrices. The code to generate conversion matrices to/from sterile schemes could become much simpler

  • And even across the API (beyond just the conversion matrices), it would eliminate a whole class of edge cases in all code that has to deal with flavors and make it a bit simpler and easier to test.

• “There should be one-- and preferably only one --obvious way to do it.” — If we give users the choice to convert back to a TwoFlavor scheme, they will have some extra cognitive load of having to decide whether to use that. Without a clear answer of when to use what option, choice can be harmful.

  • Related: Matplotlib has two different coding styles and I find it terribly frustrating that I’m never quite sure which one to use. (And then the docs use a different style than some legacy scripts I have, and I end up with a weird mix and ugh…)

Note on implementation

I'm strongly opposed to change the code in the existing models to use three flavors scheme, as it is a lot of changes to the code that already works and is tested. Let's use conversion matrices instead.

I agree with your point about wanting to minimize changes to the existing code; but am drawing a different conclusion from that: Generally, a best practice when using external data is to clean it up immediately after ingesting it, so all code that needs to understand details of the data format is in one place, and all downstream code is isolated from those details.

My preference would therefore be to convert to ThreeFlavours in the Model.__init__ function, right as we read data from the file. We would only need to change the mapping of which flavour corresponds to which column in the input file; all the code to read from the file and calculate the desired quantities would be unaffected. (This is also where we already did the OnePointFiveFlavor to TwoFlavor conversion.)
By applying this as soon as we read from the file, we only need to change one location and are guaranteed to have consistent data throughout our code. In contrast, applying the conversion matrix would have to happen at a later stage (get_flux, or maybe get_initial_spectra); but then e.g. self.luminosity (which is set in __init__) would remain hardcoded to be TwoFlavor.

[Sheshuk]

In general I agree that the TwoFlavor scheme is misleading, and should be discouraged, at least in the usage examples. However imposing a restriction means less freedom also for the developer side.

1) Our models are defined in TwoFlavor scheme

My preference would therefore be to convert to ThreeFlavours in the Model.__init__ function

The inner data representation is different for each model (loader), so that means writing different conversion code in each of those. While this can be done, I think this is a waste of time and effort.

We have code that is working and tested, and we can operate on its result, to get the desired behavior: like make conversion in get_flux method before returning the result.
It's the same for the final user, but we use standard conversion matrices on the same objects (flux.Containers) instead of implementing them for each model on whatever inner representation they have.

2) Our flavor transformations are defined in TwoFlavor scheme

While most of them can and should be switched to ThreeFlavor, as in #308, in some cases it can be more natural to use TwoFlavor (I'm no expert in this though).
I think having the freedom to define the flavor transformation in its easiest form and then just multiply by conversion matrix is much better than restrict the definition to ThreeFlavor case, and write large matrices by hand.

Right now we have the freedom to define any kind of neutrino basis, and conversion matrices, which might become useful to write implementation pretty close to the formulas in the papers, which means much easier implementation and cross-check.

[jpkneller]

Sorry for joining the discussion late, I've been out of the office the past couple of days.

"Generally, a best practice when using external data is to clean it up immediately after ingesting it, so all code that needs to understand details of the data format is in one place, and all downstream code is isolated from those details. My preference would therefore be to convert to ThreeFlavours in the Model.init function"

This is my preference also and I did that in #308 for the supernova models which inherit from PinchedModel. The init function for each model (e.g. Nakazato_2013) loads the data then the parent class (PinchedModel) checks whether turns the input data into three flavors. I split it this way to avoid duplicating the code to convert to three flavors. However I was lazy and didn't do the TwoFlavor>>ThreeFlavor conversion for either the Fornax models or the pre-supernova models so these require TwoFlavor >> ThreeFlavor conversion code in their get_initial_spectra functions. This can be easily fixed by copying the relevant block from PinchedModel. I would prefer not to add a method to SupernovaModel so that it remains an ABC.

I dislike TwoFlavor schemes for exactly the reasons indicated, you always have to check exactly what is meant by it. I don't have a strong preference for >> over << but I am reluctant to include the conversion from TwoFlavor to ThreeFlavor because, it is not unambiguous what this means: "In the face of ambiguity, refuse the temptation to guess." My preference is that if the user wants to do this conversion, they should specify it themselves. Even ThreeFlavor>>FourFlavor makes me shiver. The only conversion matrix we require for SNEWPY2 is FourFlavor>>ThreeFlavor which I buried inside neutrino.py and called Remove_Steriles.

[JostMigenda]

We have code that is working and tested, and we can operate on its result, to get the desired behavior: like make conversion in get_flux method before returning the result.

Doing the conversion only in get_flux is not sufficient; then users would still get the old format from get_initial_spectra and get_transformed_spectra and self.luminosity. But those three don’t contain instances of the NeutrinoFlux class, so we wouldn’t be able to use the conversion matrices as is.

I’m currently working on a demo showing how the changes to the __init__s would look like. I’ll be on annual leave Fri and Mon; I hope I can get this ready before then (if I don’t get pulled aside by that medical imaging project I work on as well).

Even ThreeFlavor>>FourFlavor makes me shiver.

I don’t have an issue with this: Every ThreeFlavor flux implicitly has zero sterile neutrinos; so explicitly adding this zeroed-out flux component feels very straightforward and intuitive to me. In this case, there are no conflicting conventions causing confusion.

[jpkneller]

Yes, we need to add the TwoFlavor>>ThreeFlavor conversion in the init functions for Fornax and preSN models. As I said, I was lazy (the init functions for them are big and messy).

[Sheshuk]

I think this is rather extreme.

This can be easily fixed by copying the relevant block from PinchedModel

IMO copying code is bad in general - it makes the software really hard to maintain and change.
And classes are good - they allow to put the logic in a single place, and reuse it everywhere. This is exactly what the conversion matrices do.

Copying code N times might be an easy fix, but it makes N times more probable for mistakes (and N times more time to check everything each time you need to touch it).

Even ThreeFlavor>>FourFlavor makes me shiver

If the model didn't have sterile neutrinos from the start, it will provide a ThreeFlavor output, and then you'll have to convert it to ThreeFlavor>>FourFlavor if you add the steriles in the FlavorTransformation. This requires conversion (as trivial as it is) - because it moves our arrays to different flavor basis. Do you suggest to manually add missing rows and columns instead of using matrix with defined and tested behaviour?

Also if we add more sterile flavors - will we, again, do more conversions by hand?