Generalize categorizes (keep state, etc.)

[zerothi]

Describe the feature

Currently the categories are quite strict in the sense that they provide the functionality to check individual atoms.

However, there are cases where it could be valuable to retain state information in the form of

• geometry from which it was categorized
• index in the above geometry
• neighbor information (when searching for neighbors)

A primary concern when trying to solve #202 is to reduce the overhead of doing the neighbor estimation twice, which is necessary when doing the categories.

Instead we might propose something like this:

1. a category is a base class that is only used to check stuff

2. when categorizing a Category, it should return a CategoryResult class specifying whether it succeeded (if res) and some information. I imagine something like:

   @dataclass
   class CategoryResult:
       category: Category
       info: PropertyDict # geometry, atom, [neighbors, ...]
       def __bool__(self): return not isinstance(self.category, NullCategory)

   This approach is not too different from scipy optimizers. My main worry is the optional variables associated with the categories.

3. Are there other ways to do this?

[zerothi]

Here it would also be nice if one could count the number of categories a single atom is existing in.
E.g. when extracting the rings (from #688), it might be useful to extract non-overlapping rings, in which case the initial categorization would return overlapping rings. Then one would need to figure out a single atom which is only in 1 ring. And then remove duplicates.

[zerothi]

Here are some more thoughts on how to generalize the categories.
Currently the categories are looped on individual atoms, and then checked.

However,

1. When we want to deal with categories involving multiple atoms, we have to know the full geometry before making a categorization.

2. When doing composite categories with categories involving multiple atoms (say rings), then it might be valuable to also restrict the ring to be a combination of another set of categories, i.e. 5 C and 1 N atom. This could possibly be done by having a category which fixes a preset number of atoms with another category, since there are no logical symbols involving *, we could overload the Category.__mul__ to accept a category and an integer? This could be errorprone as order of operations would matter

   # this would result in the (Ring & Z) * 5, which would never be true.
   cat = CategoryRing(...) & CategoryZ(Z=6) * 5

   Exactly how this would work I am not sure, it should be simple, and intuitive, ideas are welcome!

3. Categories should act on the complete geometry. So the interface would stay the same, yet the way it works would be changed.
   My idea would be that a categorization would always return a CategoryResult object, which stores the category it originates from, and a list of positive results. E.g. for a AtomZ it would return something like [[1], [2]]. The point of this would be that categories involving mulitple atoms could have a result index like [[1, 2, 3, 4], [7, 8, 9]]. The result object could then also hold additional information that is defined by the originating Category, this is partly inspired by the scipy, OptimizeResult which can be very versatile, but a bit hard to document. On the other hand, it might be simple enough for end-users since the additional information isn't always required knowledge (you just want the result arrays).
   This CategoryResult object would then have an __iter__ implementation that iterates the result index (as shown above).
   This should also do the equivalent for CompositeCategories, the result should encompass a result from the combined result of both categories. This might seem a bit dubious, however, I think it could be useful for the composites to contain the information on the nested categories. Consider for example the case where one wishes to extract the atoms only having 2 neigbors, and they must only be carbons, here it would be ideal if one could access a result from the neighbor result (it might contain the information on which atoms are its neighbors etc.) So we need a mechanism that lets one filter the result easily, since we can't push down the stored information from nested categories into a composite result (consider if there are two branches of categories, both with neighbor conditions, that would mean that there would be two (possibly different) additional information.
   Hence it might be ideal to do:

   neighbor = Neighbors(...)
   cat = neighbor & AtomZ(6)
   results = cat.categorize(geometry)
   for result in results.select(neighbor):
       ...

   this would only loop the positives of the composite category cat, but filter out the result array from the neighbor category (if the extra information is needed). If one is simply interested in the indices, one would simply do for result in results.

4. When dealing with composites we should, by default, force the categories to be fully expanded. I.e. each category in the Compostite should be parsed as though it didn't know of the other one. Hence one could do (same as above):

   for result in results.B: # if we should name it `A`/`B` (as is done now), or `lhs`/`rhs`.

   I am not sure if results.B would be too ambiguous with results.select(AtomZ(6))? Should it even be possible? I assume it could be useful to check for near misses with this functionality open, and it shouldn't contradict the methodology as they are separated.

5. Oh, and all of the above would deprecate the NullCategory which is actually kind of weird to deal with... So that would make things simpler.

Hmm, I might come up with some more, but here some discussion points, @pfebrer @ialcon comments?

[ialcon]

Hi @zerothi, I like very much the overall approach, and to me it looks quite logical and organized. I also like both the for result in results.select(neighbor): and for result in results.B: construction. I understand that "A" in your example would be the first category indicated in the composite creation, and "B" the second, right?

results.A --> neighbor
results.B --> AtomZ(6)

In this example what it is not fully clear to me is to what AtomZ(6) applies? To the targeted ia whose bonds are scrutinized (does ia have 2 bonds?) or to its neighbors? So, in other words, cat in cat = neighbor & AtomZ(6) would provide all carbon atoms having (say) 2 neighbors or it would provide all atoms having 2 carbon neighbors? I guess it is the first... right? There would be an easy way to do the second using a CompositeCategory?

Finally, returning to the iteration of results, which type of info would provide something like for result in results.B? In the case of neighbors, it would return the number of neighbors for each result, the indices of the neighbors for each result? Could this Category-specific results have attributes? Say:

for result in results.B :
    print(result.n_neighbors) # number of neighbors
    print(result.neighbors)     # indices of neighbors
    ...

And, super finally, wouldn't make more sense than this Category-specific information was obtained per result instead from the list of results ... :

for result in results :
    print(result.B.n_neighbors) # number of neighbors
    print(result.B.neighbors)     # indices of neighbors
    print(result.B.distances)     # distance to each neighbor
    ...

Of course I guess this also depends on how simple/complex you want each Category to be..

[pfebrer]

1. Agree

2/3. I don't understand the results you are proposing. Why do we need a 2D array of results? To me the result of categorizing should be an atom-wise property. I.e. an array of length na containing the category of each atom. I think anything else is overcomplicating the implementation and the usage. Perhaps we are trying to fit into the concept something that simply can't fit (ring identification).

4. Again I think this might be overcomplicating something as simple as categorization. Here I see the point that you are making, but do you have an example where it could be useful?

5. Agree, it is better that it disappears :)

I think a clean approach would be to first define atom-wise properties, maybe under an AtomProperty class: some properties could be:

• Z
• NNeighbors
• x, y, z
• ... basically all properties that are checked in the current categories and maybe some more.

These properties, when combined with a geometry, should return an array containing the property value for each atom.

Then, from them, using comparison operators, one could create a BooleanCategory (a subclass of AtomProperty maybe?), that could be used as a mask to select atoms for example:

• NNeighbors == 2: Boolean mask specifying if the atom has two neighbors or a different number.
• NNeighbors > 2: Boolean mask specifying if the atom has more than two neighbors.

These categories would be combined into composite categories as the current categories are now.

Then one could also use these properties to create IntCategory or StrCategory. Although here I don't see the need for special classes, they could be simple AtomProperty whose property is an integer or a string. But the point I want to make is that these would be used for classifying, rather than creating masks:

• NNeighbors is already a "classifier" since it returns integers.
• Clustered([Z, NNeighbors], some_clustering_algorithm): Here Clustered would take some other property (or multiple of them) and acording to some algorithm group atoms into clusters. Each atom would have an integer assigned which is the cluster index.

The idea is that these atom properties I think could be much more versatile and clean, as they would be defined as atom-wise arrays of data. So their shape is more predictable. And they can be reused for other things. E.g. to color or scale atoms in plots 😄

I would advocate for their returns being simple arrays, but maybe CategoryResult could be a wrapper that makes it look and work like an array while containing extra metadata.

I think the ring identification could be better off in a separate functionality unless we can frame it as an atom-wise problem. Otherwise it makes the concept of categories much more complicated.

[zerothi]

Hi @zerothi, I like very much the overall approach, and to me it looks quite logical and organized. I also like both the for result in results.select(neighbor): and for result in results.B: construction. I understand that "A" in your example would be the first category indicated in the composite creation, and "B" the second, right?

results.A --> neighbor
results.B --> AtomZ(6)

Indeed, that was the idea.

In this example what it is not fully clear to me is to what AtomZ(6) applies? To the targeted ia whose bonds are scrutinized (does ia have 2 bonds?) or to its neighbors? So, in other words, cat in cat = neighbor & AtomZ(6) would provide all carbon atoms having (say) 2 neighbors or it would provide all atoms having 2 carbon neighbors? I guess it is the first... right? There would be an easy way to do the second using a CompositeCategory?
The first, yes, AtomZ selects atoms with a given Z. Then the neighbors are determined on the smaller section.

Finally, returning to the iteration of results, which type of info would provide something like for result in results.B? In the case of neighbors, it would return the number of neighbors for each result, the indices of the neighbors for each result? Could this Category-specific results have attributes? Say:

for result in results.B :
    print(result.n_neighbors) # number of neighbors
    print(result.neighbors)     # indices of neighbors
    ...

Yeah, here I am a bit more undecided. Since you are already in the AtomZ category, you can't access the neighbors. That is a bit difficult on the composite categories.

And, super finally, wouldn't make more sense than this Category-specific information was obtained per result instead from the list of results ... :

for result in results :
    print(result.B.n_neighbors) # number of neighbors
    print(result.B.neighbors)     # indices of neighbors
    print(result.B.distances)     # distance to each neighbor
    ...

Of course I guess this also depends on how simple/complex you want each Category to be..
Agreed! This is a good suggestion, lets continue along that line.

1. Agree

2/3. I don't understand the results you are proposing. Why do we need a 2D array of results? To me the result of categorizing should be an atom-wise property. I.e. an array of length na containing the category of each atom. I think anything else is overcomplicating the implementation and the usage. Perhaps we are trying to fit into the concept something that simply can't fit (ring identification).

The main problem of the way things currently work is exactly that it is a category per atom.
If the category is per atom, then it requires additional checks to describe which atoms belongs to a ring or not? No?
Perhaps the 2D lists would only be interesting for those categories which deals with multiple atoms (such as the ring). But for AtomZ it could just as well be a 1D array.

4. Again I think this might be overcomplicating something as simple as categorization. Here I see the point that you are making, but do you have an example where it could be useful?

Another motivation of mine, is the re, it works quite similarly. You query some quantities, and loop the results.
A main use would be the rings, where looping the results is equivalent to looping the rings.

EDIT, sorry didn't do an actual example.

My motivation for this would be that one wanted to search for neighbors, but only with AtomZ, one could then check whether it was the lhs or the rhs which didn't match it by using this. It might not be necessary, so we could do without, and force the user to do a 2nd request. I guess this point could be left out!

5. Agree, it is better that it disappears :)

I think a clean approach would be to first define atom-wise properties, maybe under an AtomProperty class: some properties could be:

* `Z`

* `NNeighbors`

* `x`, `y`, `z`

* ... basically all properties that are checked in the current categories and maybe some more.

The problem with properties, is that one might be querying properties with manual tweaking values, and then it gets fairly complicated to have something defined generically, it is also rather hard to make more complex searches.

These properties, when combined with a geometry, should return an array containing the property value for each atom.

Then, from them, using comparison operators, one could create a BooleanCategory (a subclass of AtomProperty maybe?), that could be used as a mask to select atoms for example:

* `NNeighbors == 2`: Boolean mask specifying if the atom has two neighbors or a different number.

* `NNeighbors > 2`: Boolean mask specifying if the atom has more than two neighbors.

These categories would be combined into composite categories as the current categories are now.

Then one could also use these properties to create IntCategory or StrCategory. Although here I don't see the need for special classes, they could be simple AtomProperty whose property is an integer or a string. But the point I want to make is that these would be used for classifying, rather than creating masks:

* `NNeighbors` is already a "classifier" since it returns integers.

* `Clustered([Z, NNeighbors], some_clustering_algorithm)`: Here `Clustered` would take some other property (or multiple of them) and acording to some algorithm group atoms into clusters. Each atom would have an integer assigned which is the cluster index.

I think that this would just be difficult to extend, whereas my approach seems a bit more generic, and would allow one to search also based on other criteria. For instance, a ring with 5 Carbon and 1 Nitrogen.

The idea is that these atom properties I think could be much more versatile and clean, as they would be defined as atom-wise arrays of data. So their shape is more predictable. And they can be reused for other things. E.g. to color or scale atoms in plots 😄

I would advocate for their returns being simple arrays, but maybe CategoryResult could be a wrapper that makes it look and work like an array while containing extra metadata.
If wanting 1D arrays, it should be trivially handled, no? indices = np.concatenate([result.index for result in results])

I think the ring identification could be better off in a separate functionality unless we can frame it as an atom-wise problem. Otherwise it makes the concept of categories much more complicated.

That was why I suggested to make it based on the full geometry, and not necessarily an Atom property.

My main worry is that limiting its functionality basically makes it useless. Using a AtomZ category is just geometry.atoms.Z == 6, but nesting the logical operators becomes a head-ache, and here the categories would be useful.

@pfebrer see my edit!

[pfebrer]

The main problem of the way things currently work is exactly that it is a category per atom.

Isn't the problem that the category is computed atom by atom? I mean if you give the whole geometry to the property calculator and then it calculates an atom-wise property, that is not a limiting factor, no?

My motivation for this would be that one wanted to search for neighbors, but only with AtomZ, one could then check whether it was the lhs or the rhs which didn't match it by using this. It might not be necessary, so we could do without, and force the user to do a 2nd request. I guess this point could be left out!

I understand this, but maybe a cleaner solution would be to cache the results of categorization / property calculation, so that a second query is basically for free?

The problem with properties, is that one might be querying properties with manual tweaking values, and then it gets fairly complicated to have something defined generically, it is also rather hard to make more complex searches.

No, why? You can get at least to the same level of complexity. E.g. imagine that the NNeighbors property accepts a neighbor argument to filter. It would be exactly like the current neighbors category:

NNeighbor(R=3) # Gives you the number of neighbors up to a certain R.
NNeighbor(R=3, neighbor=Z > 1) # Gives you the number of neighbors without counting hydrogens.

The only difference is that you specify the properties, not the categories, and then you create "categories" by using the comparison operators. I'd say that if something you have more freedom to get as complex as you want easily. You can for example create "categories" by comparing two properties:

x > y # This would be a mask that selects atoms whose x coordinate is bigger than the y coordinate.

And actually I have thought of something that could work quite good with the ring finder. The returned property could be of shape (na, n_rings). Something like:

[
 [0, 0, 0], # Atom 0 belongs to [ring_0, ring_1, ring_2] ? 
 [0, 1, 1], # Atom 1 belongs to [ring_0, ring_1, ring_2] ? 
 [1, 0, 0], # Atom 2 belongs to [ring_0, ring_1, ring_2] ? 
 ...
]

Which could be a sparse matrix, of course. But the point is that this would still be an atom-wise property. And doing things like:

rings_property = RingFinder()
rings_count = np.sum(rings_property, axis=-1)

Would give you a property that computes the number of rings to which each atom belongs. Again you would then be free to create a category however complex you want using this property.

You could do many other cool things. E.g.:

# Property that returns a (na, nZ) array specifying how many neighbors of each type an atom has.
neighs = NNeighbors(dim1=Z) 

# Category that selects atoms that have 1 hydrogen neighbor and 2 carbon neighbors.
my_complex_category = neighs[:, 1] == 1 & neighs[:, 6] == 2

Wouldn't it be nice? 🥲

[pfebrer]

By the way all this lazy arithmetic (and the caching) could be supported by the nodes framework.

[ialcon]

I kind of also understand the point of @pfebrer and, somehow, agree with it. At the same time, I also understand that, as you say @zerothi, limiting to atom-wise categorization limits a lot what one could do (specially in terms of complexity and versatility).

Would it be possible to have Categories for atom-wise properties (that is, as @pfebrer advocates) and something like Classifiers (or a similar name, maybe Fragments ...) that classify atoms based on a particular property? In such a way Categories could be maintained clean and simple, whereas Classifiers could go crazier, without affecting the fundamental structure of Categories. Of course, Classifiers could internally use Categories to achieve their classification, but Categories in turn would not care at all about Classifiers and not being affected by them at all.

So, following the initial example by @zerothi, instead of doing:

cat = CategoryRing(...) & CategoryZ(Z=6) * 5

You could have:

get_rings = Ring([(AtomZ(6), 5), (AtomZ(7), 1)])
rings = Ring.classify(geometry)
for ring in rings :
    print(ring.na)
    print(ring.id_ring)
    print(ring.idx)
    print(ring.dihedral)
    ....

I am not sure if this makes sense to you, but to me is sort of trying to have some sort of hierarchy, so that we can keep separate the complexity of more versatile objects from the simplier atom-wise categories.

[zerothi]

The main problem of the way things currently work is exactly that it is a category per atom.

Isn't the problem that the category is computed atom by atom? I mean if you give the whole geometry to the property calculator and then it calculates an atom-wise property, that is not a limiting factor, no?

Yes, the problem is that it is an atom-wise property. That is exactly what I am trying to propose to go about. A Category should be able to encompass any number of atoms etc. It might be overkill for AtomZ to have a list of positives, an 1D array might work here. However, I think groups of positives can be made quite versatile (as always I could be wrong) ;)

My motivation for this would be that one wanted to search for neighbors, but only with AtomZ, one could then check whether it was the lhs or the rhs which didn't match it by using this. It might not be necessary, so we could do without, and force the user to do a 2nd request. I guess this point could be left out!

I understand this, but maybe a cleaner solution would be to cache the results of categorization / property calculation, so that a second query is basically for free?

Yeah, this might be overdoing it, perhaps we shouldn't worry here. ;)

The problem with properties, is that one might be querying properties with manual tweaking values, and then it gets fairly complicated to have something defined generically, it is also rather hard to make more complex searches.

No, why? You can get at least to the same level of complexity. E.g. imagine that the NNeighbors property accepts a neighbor argument to filter. It would be exactly like the current neighbors category:

NNeighbor(R=3) # Gives you the number of neighbors up to a certain R.
NNeighbor(R=3, neighbor=Z > 1) # Gives you the number of neighbors without counting hydrogens.

But what is the difference between this and a Category, aren't they the same?
Do you mean that categories are something that extracts properties, or?
I don't see a big difference between the two.

The only difference is that you specify the properties, not the categories, and then you create "categories" by using the comparison operators. I'd say that if something you have more freedom to get as complex as you want easily. You can for example create "categories" by comparing two properties:

x > y # This would be a mask that selects atoms whose x coordinate is bigger than the y coordinate.

I guess categories could do the same, no ?

And actually I have thought of something that could work quite good with the ring finder. The returned property could be of shape (na, n_rings). Something like:

[
 [0, 0, 0], # Atom 0 belongs to [ring_0, ring_1, ring_2] ? 
 [0, 1, 1], # Atom 1 belongs to [ring_0, ring_1, ring_2] ? 
 [1, 0, 0], # Atom 2 belongs to [ring_0, ring_1, ring_2] ? 
 ...
]

This just looks horrendous? Wouldn't these groupings be simpler to loop in the proposed result method? Kind of like pandas.groupby?

Which could be a sparse matrix, of course. But the point is that this would still be an atom-wise property. And doing things like:

rings_property = RingFinder()
rings_count = np.sum(rings_property, axis=-1)

Would give you a property that computes the number of rings to which each atom belongs. Again you would then be free to create a category however complex you want using this property.

You could do many other cool things. E.g.:

# Property that returns a (na, nZ) array specifying how many neighbors of each type an atom has.
neighs = NNeighbors(dim1=Z) 

# Category that selects atoms that have 1 hydrogen neighbor and 2 carbon neighbors.
my_complex_category = neighs[:, 1] == 1 & neighs[:, 6] == 2

Wouldn't it be nice? 🥲
I see, I agree that it could be nice if the methods are generic in some way. Hmm, this needs some thinking. I kind of don't like the way one should post-process things, it seems to be cumbersome when wanting to do more complicated things? By specifying constraints to what you request, you have less flexibility at query time, but simpler parsing.

I kind of also understand the point of @pfebrer and, somehow, agree with it. At the same time, I also understand that, as you say @zerothi, limiting to atom-wise categorization limits a lot what one could do (specially in terms of complexity and versatility).

Would it be possible to have Categories for atom-wise properties (that is, as @pfebrer advocates) and something like Classifiers (or a similar name, maybe Fragments ...) that classify atoms based on a particular property? In such a way Categories could be maintained clean and simple, whereas Classifiers could go crazier, without affecting the fundamental structure of Categories. Of course, Classifiers could internally use Categories to achieve their classification, but Categories in turn would not care at all about Classifiers and not being affected by them at all.
If we can get all functionality in 1 class, then I would really try and maintain it like that, I don't know if that will work, if not we might have to have something like your Fragments idea.

So, following the initial example by @zerothi, instead of doing:

cat = CategoryRing(...) & CategoryZ(Z=6) * 5

You could have:

get_rings = Ring([(AtomZ(6), 5), (AtomZ(7), 1)])
rings = Ring.classify(geometry)
for ring in rings :
    print(ring.na)
    print(ring.id_ring)
    print(ring.idx)
    print(ring.dihedral)
    ....

I am not sure if this makes sense to you, but to me is sort of trying to have some sort of hierarchy, so that we can keep separate the complexity of more versatile objects from the simplier atom-wise categories.

I somewhat agree, i just always find API's that have some kind of flexibility very nice to work with, it takes a bit of time to work it out. But...

That is why this discussion is important, what can be done, and how would it work ;)

[pfebrer]

Yes, the problem is that it is an atom-wise property. That is exactly what I am trying to propose to go about. A Category should be able to encompass any number of atoms etc. It might be overkill for AtomZ to have a list of positives, an 1D array might work here. However, I think groups of positives can be made quite versatile (as always I could be wrong) ;)

The way I see it, in the rings example you are creating multiple categories on the fly, one for each ring that you find: belongs_to_ring_0: (true, false), belongs_to_ring_1: (true, false), ... And then you are defining the value of each atom on that category. So effectively you are building a 2D boolean array like the one I propose. Then you choose to represent this sparse matrix as a list of lists storing the indices of the atoms that are True, because you think it is the most convenient way. But I believe if you think of it as an atom-wise property with an extra dimension (i_ring) it is very general because then you can reduce that result into whatever you want with common matrix operations.

# This would give you the number of rings each atom belongs to
np.sum(rings_property, axis=-1)

# This would give you some measure of "connectivity" (through rings) between atoms
rings_property @ rings_property.T

# This would give you the atoms that belong to ring 2
np.where(rings_property[:, 2])

# Etc...

This just looks horrendous? Wouldn't these groupings be simpler to loop in the proposed result method? Kind of like pandas.groupby?

The convenience of the shape of the result depends on what you want to achieve afterwards. It is not obvious to me that with the one that you propose it is simpler to achieve @ialcon's original goal (separating the two ribbons of an NPG).

I guess categories could do the same, no ?

Yes, but the point is that if you allow the user to build their own conditions you don't have to implement all the possible conditions as categories.

But what is the difference between this and a Category, aren't they the same?
Do you mean that categories are something that extracts properties, or?
I don't see a big difference between the two.

The main point that I want to make is that a category (as it is now) is the combination of a property + a condition. Internally inside categories we are implementing complex properties like number of neighbors, but we are not exposing that for users to use however they want. Also, we are choosing which conditions to implement, which limits the functionality.

Another thing is that what is now called Category is only a boolean category. But categorical variables are also typically integers or strings. An atom property could be of any type. And of course depending on what type it is, that property can be used for some things or for others.

Regarding the points of @ialcon, atom-wise properties could also have metadata, so I would say this is not an argument in favor of categories vs properties.