RFI: Overhauling how Materials and Isotopes are handled.

[MicahGale]

As has been discussed for awhile the interface for Material, MaterialComponent, etc. needs to be overhauled. So here's the overall thoughts I have so far.

1. Preserve backwards compatibility for Material to avoid a major release.
2. Make Isotope parse all isotope formats that I like (see Add support for MCNP 6.3.1 materials #179). I'm thinking only the set:
   1. '1001'
   2. 1001
   3. 'H-1'
   4. 'H1'
   5. 'hydrogen-1'
   6. 'hydrogen1'
   7. 'h-1m1'...
   8. ('h', 1)
   9. (Element('H'), 2)
   10. ('hydrogen', 3)
3. Maintain Material.material_components as is.
4. Make Material a duck dictionary. This would accept any valid Isotope instance or Isotope string (above), and then return the fraction as a float.
   1. I just came up with the idea of treating element, isotopes as a 2d array. I feel like this is a good idea. Also could do 3D array of [element, isotope, isomer].

So I think the following code should all be valid:

mat = ....

#adding isotopes
mat["h1"] = 0.01
mat["hydrogen-2"] = 0.01
mat[Element("H"), 2] = 0.01

# deleting
del mat["H", :]

Thoughts @tjlaboss?

[MicahGale]

Other thought: there needs to be a way to update and view the libraries for every component, and also see the particle specific libraries. Maybe something like:

mat.def_libraries["neutron"] = "00c"
mat.iso_libraries[1001] = "80c"

[tjlaboss]

For the next minor release, this makes sense. A compatibility-breaking release that adds a friendly interface to materials and their components will be welcome in the future : )

For the cross section library, continuing to support the extension would be intuitive to users. At the very least, 1001.00c and H-1.00c (as of MCNP 6.3.1) must be supported. Allowing H1.00c would also be trivial. A third option in the tuple, such as (Element('H'), 2, '00c'), might be desirable.

[MicahGale]

Good point. Does MCNP support:

m1    1001.80c 0.5
      1001.00c 0.5

?

Intuitively, isotopes should be the unique key, not isotope + library, but I'm not sure if MCNP agrees. If this is the case then I would allow using mat['1001.80c] but just treat it as mat[1001].

[tjlaboss]

Yes, you can mix and match. MCNP gives a warning.

It will even let you do :

m1    1001.80c 0.5
      1001.80c 0.5

[MicahGale]

Let's check the docs. The one thing I'm not a fan of is: the main use case of libraries is changing them. So 1001.80c to 1001.00c rather than having both in the same material. Making the library part of the key makes the library change operation clunky, unless we provide methods for that.

[MicahGale]

As for keys. I think internally this should be stored as a 3- or 4-tuple of: element, A-number, isomeric state, library?. The isomeric state should always be included as it is with MCNP.

[tjlaboss]

The MCNP documentation gives no restrictions on what can coexist in the material definitions. We should assume that nothing is mutually exclusive (although at times it may be appropriate to raise warnings).

[MicahGale]

Reading section 5.6.1 (6.3 manual) there is not explicit guidance on if repeated ZAIDs are allowed.

Data-card Form	Mm z1 f1 z2 f2 . . . zK fK keyword = values(s)
zk	ZAID, Either a full ZZZAAA.abx or partial ZZZAAA element or nuclide identifier
	for each constituent k, where
	• ZZZZZZ represents the atomic number;
	• if AAA > 0, then AAA represents the atomic mass number; and
	• if AAA = 000, then AAA indicates a naturally occurring element ( 2 ).
	• ab is the alphanumeric library identifier; and
	• x is the class of data.

I'm interpreting z1... zK to mean that zK must be unique. However zK is the ZAID + library, so I think it should be the 4-tuple then that makes it unique.

[MicahGale]

See #504.

[MicahGale]

For atomic data which should be valid?

1. H
2. hydrogen
3. H0
4. hydrogen-0

Requiring an explicit 0 just feels jank and MCNP-like.

[MicahGale]

See e0b7b95 for example tests.

[MicahGale]

And 98a4ddc

[tjlaboss]

I don't mind the 0. Everyone used it until OpenMC switched to "-Nat"

[MicahGale]

The goal is to not perpetuate Stockholm syndrome.

[tjlaboss]

I want this on a T-shirt

[MicahGale]

Other thought: how should slices be handled?

I think should probably be the method

>>> mat["U",:]
{Isotope: ..., }

[tjlaboss]

That would be a neat feature. I'd use it for total uranium mass and fissile equivalent stuff.

[MicahGale]

I keep making things more and more complicated for myself :D

[MicahGale]

See these two gists for some thoughts I had:

• General material indexing

• Material slicing

Some initial feedback I have gotten from @pshriwise is that this isn't pythonic. It specifically violates the principle:

There should be one-- and preferably only one --obvious way to do it

Given this I think this should be greatly simplified, and only one (maybe) two ways should be allowed. ZAID is not an option, IMO.

I think symbol-A[meta][.library] notation should be that way. e.g., U-235m1.80c. There are questions of:

1. Should the hyphen be optional?
2. Should capitalization be enforced?
3. Should all MCNP nuclide notations be supported. (MCNP 6.3.1 will be expanding these options soon.)

As for slicing this creates a great justification for tuple notation. Once again simplification is needed. I think Element slicing should only support slicing by Z number. mat[33:92] is a lot clearer on Boron being excluded than mat["As":"U"].

Finally @gonuke was concerned with library slicing. It isn't the most intuitive.

As proposed mat["00c":"80p"] isn't allowed, but feels like it should be. The most intuitive interpretation is that it would be all strings between those when alphabetized. This probably is very similar to what I am proposing, but it isn't the clearest.

[MicahGale]

Next discussion: assuming #504 is the case. How should repeated nuclides be handled?

Current work

For the current design in mat_redesign I implemented a pointer data structure that then points to the right key in Material.material_components, which then returns just the fraction when using __getitem__ i.e., mat[idx].

This pointer is a nested dictionary that is structured as:

[Element][(A, meta_state)][Library]

So the current interface would be:

>>>mat["U-235.80c"]
0.05

However, if duplicate entries are supported this would need to return an iterable (at least sometimes).

One option is to always have it return a list:

>>>mat["U-235"]
[0.05]
>>>mat["U-238"]
[0.05, 0.10]

There's complicated in between options when it would sometimes return a scalar or an iterable. But that feels janky.

Using lists can lead to ambiguity. Imagine the two materials:

m1 1001.80c 0.01 1001.80c 0.01 1001.80c 0.01 1001.80c 0.01
m2 1001.00c 0.01 1001.00c 0.01 1001.80c 0.01 1001.80c 0.01

Both materials would return

>>>mat["H-1"]
[0.01, 0.01, 0.01, 0.01]

But those entries would have different ambiguous meanings.
This is in part because this includes an implicit library slice.

One option to avoid this would be to return a dictionary (or Material container thingy) whenever slicing occurs.

So in the above example:

>>> materials[1]["H-1"]
{Isotope("1001.80c"): [0.01] * 4}
>>> materials[2]["H-1"]
{Isotope("1001.00c"): [0.01, 0.01], Isotope("1001.80c"): [0.01, 0.01]}

This final interface feels the most correct because:

Explicit is better than Implicit

Though whenever setting values having to work with a list inside of a dictionary violates:

Simple is better than complex

[tjlaboss]

Don't use __getitem__; use a find() function

Rather than trying to do everything through a series[of][getitem][indices], implement Material.find() where the arguments control the behavior, similar to openmc.Tally.get_values().

Example:

>>> materials[2].find(element="H", A=1, xslib="00c", sum=True)
{Isotope("H-1.00c"): 0.03}
>>> materials[2].find(Z=1, H=1)
{
  Isotope("H-1.00c"): 0.01,   
  Isotope("H-1.00c"): 0.02, 
  Isotope("H-1.80c"): 0.03,
  Isotope("H-1.80c"): 0.04
}
>>> materials[2].find(zaid="1001", xslib="04p")
{}
>>> materials[2].find(Z=slice(92,118))
{Isotope("U-238.00c": 0.1, Isotope("Cf-252.01c"): 1e-5}

or the equivalent zip / tuple behavior.

[MicahGale]

This is because the dictionary is indexed by Isotope instances, with different memory addresses. The user can't just create a new Isotope to index the dictionary by (which, incidentally, led to quite a bit of confusion when using a deepcopy of a Materials collection). Which leads me to my next comment.

This can be solved by some combination of immutability, __hash__ and __eq__.

[MicahGale]

Don't use __getitem__; use a find() function

Rather than trying to do everything through a series[of][getitem][indices], implement Material.find() where the arguments control the behavior, similar to openmc.Tally.get_values().

Example:

>>> materials[2].find(element="H", A=1, xslib="00c", sum=True)
{Isotope("H-1.00c"): 0.03}
>>> materials[2].find(Z=1, H=1)
{
  Isotope("H-1.00c"): 0.01,   
  Isotope("H-1.00c"): 0.02, 
  Isotope("H-1.80c"): 0.03,
  Isotope("H-1.80c"): 0.04
}
>>> materials[2].find(zaid="1001", xslib="04p")
{}
>>> materials[2].find(Z=slice(92,118))
{Isotope("U-238.00c": 0.1, Isotope("Cf-252.01c"): 1e-5}

or the equivalent zip / tuple behavior.

After thinking about it for a few hours. I think you're right.

So:

1. store it in a List of Tuples. I think the optimal search that a dict would give isn't worth while. Most materials have < 20 isotopes, and probably 95% have < 100.
2. Implement __iter__/ __next__/__getitem__ in string way
3. Implement append
4. Implement find

I'm thinking:

def find(self, fancy_name=None, element=None, isotope=None, a=None, meta_state=None, library=None):
    """
    Filter stuff!


    :param fancy_name: filter by `Isotope.get_from_fancy_name`
    :type fancy_name: str, int
    :param element: Filter by element. Filter by Element, atomic symbol, or z number. Slicing is done by Z-number. Or can pass a function that accepts an Element and returns a bool.
    :type element: Element, str, int, slice, func
    ...
    :yields: an iterable of tuples of Isotopes and fraction 
    :rtype: generator

Then we could also implement find_values that only returns the fractions.

[MicahGale]

So far in my investigation the only way I have found to make a slice is either with foo[1:2] or slice. Are you aware of any syntax sugar that would make the function call easy?

[MicahGale]

Thought I just had: if we are using find: how do you update values? Obviously returning just a float as the fraction will not make that possible. We could return a wrapper object ala MaterialComponent, but better. Or we could just say find is read-only.

I was about suggest update(isotope, fraction). However this doesn't deal with the duplication problem.

Maybe the only way would be through __setitem__? Then users need to know the indices. So find might want to return:

[(idx, isotope, fraction), ....]

[MicahGale]

Implementing Dunders

Which dunders (__func__) functions should be implemented?

So far these will be implemented:

1. __setitem__, __getitem__, __delitem__ to provide a list like interface mat[1] = (Isotope(...), 0.1)
2. __hash__ , __eq__ (nvm on __hash__. Mutables should not hash).
3. __iter__
4. __contains__
5. __len__

New proposed ones. Specifically ones for doing easy material mixing. This works well with the submaterial concept that @dodu94 has:

1. __add__: allow adding two materials to create a mixture: mixture = water + boric_acid
2. __mul__: allow redefining, and diluting materials: mixture = water + boric_acid * 1e-6. This brings in questions of when to do normalization.
3. __rmul__: support: 1e-6 * boric_acid
4. __iadd__: allow in-place mixing: water += boric_acid * 1e-6.
5. __imul__: support in-place dilution. Not sure if I like this one: boric_acid *= 1e-6

Weirder ones that might be interesting:

1. __and__. This one feels weird. The only use case I can think of is trying to find isotopes in both of two materials. This would make it more set like, and I think duck-programming a set is a bad idea. Maybe we should have a method like: get_isotope_set? Then users could find: common_isos = mat1.get_isotope_set() & mat2.get_isotope_set()?
2. __div__: I don't see how this would make much sense: unless we do: mat* 1/div. Might be useful
3. __mod__: Not sure what mat % 2 means
4. __matmul__: while technically these are vectors they basically all exist in their own vector space so this feels like non-sense. Also what is the dot product of two materials?
5. __lshift__ and other bit tricks. I think these are non-sense
6. __reversed__. Could be good. Not sure if it's needed if users could just do: reversed(list(mat))
7. __call__. I don't have a use-case, but this just seems like a fun chaotic easter-egg.

Misc:

1. __bool__, probably should make truthiness based on the components list

[MicahGale]

but if you call it from an instance wouldn't you have to do mat1.mix(mat1...) still?

[dodu94]

Just to add my 2 cents on how I perform material mixing:

here you can find the method I am using

Basically the methods accepts a list of materials to mix, a list of the same size which indicates the fraction of each material in the new mix, a flag specifying if this percentage is to be intended as mass or atom fraction, and a library manager. The library manager is an object that I defined in F4Enix which has access to the user xsdir mostly to perform material translations from one nuclear data library to another. It also contains isotopic data that allow to switch from mass to atom fraction (which is also a material method worth having as a standalone).

You do not need densities in this way, which is MCNP-like

[tjlaboss]

That is true. You only need the original densities if you need to calculate the mixed density.

[tjlaboss]

but if you call it from an instance wouldn't you have to do mat1.mix(mat1...) still?

Yes, but it's still an error precursor; mat1.mix(mat2, mat3) would give a nasty surprise.

[MicahGale]

What if we made this a method of Materials? That way it's not a top level function but you avoid @classmethod errors.

You could then call: mat = problem.materials.mix(mat1, 0.1, mat2, 0.9).

This begs the questions:

1. Is this namespace appropriate?
2. Should it accept object instances or object numbers? I feel like numbers is a bad idea.
3. Should it accept *args with the every other pattern like above or: mix([mat1, mat2], [0.1, 0.9]). I think the latter would be clearer and is very pyplot.

This is the final feature "needed" for 1.0.0a1. Though we could always bump this to be a 1.1.0 feature.

I'm assuming material math is completely dead, right?

[MicahGale]

Other thought: should there be immutable Materials? I feel like some materials should be immutable. water is water ... probably. We're going to skip discussions on nuclear data. The use case of .add or .append building a material is important. I think though it would be nice to create a new immutable object from a mutable one that you can then trust to say: water is water.

[tjlaboss]

[MicahGale]

I got around to it eventually.

[tjlaboss]

The 1.0 Migration documentation should provide guidance on how to replace deprecated user-facing functionality; namely, material_component.

[MicahGale]

Do you think it's concrete enough? i was thinking of populating this during alpha testing.

[tjlaboss]

Not yet, but it should be there before, or no later than, the MontePy 1.0 release.

[MicahGale]

Agreed. I added it to my checklist for the alpha release. I'm thinking make the migration plan concrete with 1.0.0.a2

[MicahGale]

I started writing it up in #587.

[tjlaboss]

The migration guide has already helped me navigate deprecations during alpha testing.

[MicahGale]

The next question is: how should we provide information about the libraries that will actually be used?

As a reminder the precendence is:

1. nuclide library extensions: 1001.80c
2. material level default libraries nlib=80c
3. problem level defaults: m0 nlib=80c
4. xsdir ordering.

Reading Xsdir is out of scope right now.

How should users get this info?

is it nuclide.libraries[NEUTRON]?
Or material.get_library(i, NEUTRON)?
...