added divide support, added normalization for Vector3d

[din14970]

Background

I'm just testing the waters with orix (big picture: I'm still trying to work out an improved template matching procedure), so hereby a small PR. I want to improve on the code I added in diffsims some time ago so that these functions use orix objects like Vector3d, Orientation, etc. I am also working on a stereographic/pole figure plotting function, so I don't have to take things over to MTEX all the time. For this, I felt I was lacking some basic functionality in the Vector3d class.

This PR

Divide did not seem to be supported for vectors, nor was there an easy "normalize" function. Now you can get normalized vectors in two ways with:

>>> normed = vectors.normalize()
>>> normed = vectors / vectors.norm

[hakonanes]

Hi @din14970,

great that you want to contribute!

Do you plan on using the normalize() method for something else than getting unit vectors? If not you can just use Vector3d.unit (https://orix.readthedocs.io/en/stable/reference.html#orix.vector.vector3d.Vector3d, see the third examples block), right?

The division is a really nice addition.

I am also working on a stereographic/pole figure plotting function

That would be great. I have a loose idea of where that would fit in orix, namely that we should have a orix.projections module with a base SphericalProjection class with StereographicProjection etc. subclasses. Similar to what MTEX has. This has been discussed previously here, and @dnjohnstone has already started on creating a direct Matplotlib transformation here.

Would your function(s) fit into such a class "framework"? If you don't want to write that framework, I might be able to do that either before or after you add your functionality.

[din14970]

Vector3d.unit

Aha, exactly what I was looking for basically. I figured it should be implemented somewhere but couldn't immediately find it so just quickly threw together normalize. Maybe normalize() can stay as an alias to unit?

and @dnjohnstone has already started on creating a direct Matplotlib transformation

Ah so it has been in the works for a while. Actually whatever was added in that PR looks much more advanced than what I have. Obviously I'm not aware of the big picture, for now I was just going with something quick and dirty to get some type of visual on vectors and orientations. I was just adding a module in the plot subpackage and tried to model a bit off of what is there. I converting the Vector3d to stereo projection simply with a function in this module (I split it between lower and upper hemisphere):

def _stereographic_coordinates(vector3d):
    """
    Converts Vector3d objects into stereographic coordinates

    The upper and lower hemisphere are separated so that everything
    can be plot within a unit circle.

    Parameters
    ----------
    vector3d : orix.vector.vector3d.Vector3d
        The list of vectors to convert to stereographic coordinates
    """
    vector3d = vector3d.normalize()
    vector3d_up = vector3d[vector3d.z >= 0]
    vector3d_down = vector3d[vector3d.z < 0]
    xsu = vector3d_up.x.data/(1+vector3d_up.z.data)
    ysu = vector3d_up.y.data/(1+vector3d_up.z.data)
    xsb = -vector3d_down.x.data/(-1+vector3d_down.z.data)
    ysb = -vector3d_down.y.data/(-1+vector3d_down.z.data)
    return xsu, ysu, xsb, ysb

Then one can just scatter plot xsu vs ysu. For now I was just going with full stereo plots:

Pole and inverse pole figures that take into consideration symmetry is something I didn't want to deal with.

[hakonanes]

Maybe normalize() can stay as an alias to unit?

I'm used to thinking of vectors, and unit vectors (not normalized vectors), so the Vector3d.unit is for me the most natural syntax. I'm not a fan of having two ways to the same answer, but would be happy to hear @pc494's thought on this.

Obviously I'm not aware of the big picture, for now I was just going with something quick and dirty to get some type of visual on vectors and orientations.

Something quick and dirty is a very good start, which we can build upon later.

I was just adding a module in the plot subpackage and tried to model a bit off of what is there.

What I hope we can build is a general projections module where

• the projection to stereographic coordinates is not directly coupled to plotting, so that it can be used in other situations
• the inverse projection is implemented at the same time, along the lines of how we've started on a kikuchipy.projections module in the kikuchipy package: https://github.com/pyxem/kikuchipy/blob/master/kikuchipy/projections/spherical_projection.py. We plan to move some of these projections to orix (kikuchipy depends on orix) in the future.

I converting the Vector3d to stereo projection simply with a function in this module (I split it between lower and upper hemisphere)

What you have here looks great. I would be happy to help get this into a orix.projections module with documentation etc.

[din14970]

Ok, removed normalize. Will look into the projection option but will make it a separate PR. Not sure why the notImplemented line is not getting covered by the tests in __rtruedivide__.

[hakonanes]

Will look into the projection option but will make it a separate PR.

Sounds good, we can discuss that then. (I think getting your stuff in first and then combined that with what's in #62 might be the simplest way forward.)

Not sure why the notImplemented line is not getting covered by the tests in __rtruedivide__.

I'll have a look later today or tomorrow.

[hakonanes]

Not sure why the notImplemented line is not getting covered by the tests in __rtruedivide__.

Don't know either, but you can run this test to cover it

@pytest.mark.xfail
def test_rdiv():
    v = Vector3d([1, 2, 3])
    other = "dracula"
    _ = other / v

By the way, have you run black code formatting? https://black.readthedocs.io/en/stable/installation_and_usage.html

[pc494]

I'm broadly happy with what's currently in this PR although I would like to ask about how unambigous we think:

Vector3d([1, 2,3]) / (1,2)

and

2 / Vector3d([1, 2,3])

are. I'm especially suspicious of the latter case.

[hakonanes]

>>> Vector3d([1, 2, 3]) / (1, 2)
Vector3d (2,)
[[1.  2.  3. ]
 [0.5 1.  1.5]]
>>> 2 / Vector3d([1, 2, 3])
Vector3d (1,)
[[2.     1.     0.6667]]

Actually, I find this broadcasting natural.

[din14970]

@pc494
I pretty much copied the multiplication rules.
The second one is the same behavior of a numpy array. I think defaulting to some element-wise operation is pretty natural.
The first one is a bit more weird, but I actually like this - often you want to repeat some operation on a 1D array for different values, and with numpy arrays it can be a mess.

[hakonanes]

Okay, I'm reconsidering the first example since it isn't supported by NumPy... The second case IS supported by NumPy. So, I'm ambivalent about the first case, but support the second one.

[din14970]

@hakonanes I've heard of black but haven't used it so far. I gave it a try and it looks like it reformats quite a lot of files in the repo...

[hakonanes]

@hakonanes I've heard of black but haven't used it so far. I gave it a try and it looks like it reformats quite a lot of files in the repo...

Ah yes, I don't think they are all formatted according to black, @pc494? I'm sorry if you now have to git restore <lots-of-files>. It suffices for you to run black orix/vector/vector3d.py orix/tests/test_vector3d.py.

[hakonanes]

Great, thanks a lot @din14970.

[hakonanes]

I haven't touched the Vector3d "ufuncs" before, so I'll let you, @pc494 and @din14970, decide which broadcasting we should support.

After this is decided, we should add these operations to the Vector3d docstring examples:

orix/orix/vector/vector3d.py

Line 67 in 08119d3

>>> 3 * v

Edit: And add a changelog entry. After that, I'm satisfied with this addition.

[pc494]

I still think Vector3d should was designed to and should behave as a vector rather than array. In matlab elementwise operations have a different syntax (./), but in python we don't have that. If you saw, written as a formulae

2 / vector

you would be rightly upset. It doesn't mean anything and I think for syntactic consistency we shouldn't offer this for the class that claims to be "vectors". Additionally, there is a workaround if you did want this (scalar divided by an array) ,which would be to use a Scalar of the correct shape.

Do either of you find this case convincing?

PS: Although it's perhaps not ideal it seems fairly obvious that a user writing v / 2 means the same thing as 0.5 * v so I would leave that in.

[hakonanes]

Do either of you find this case convincing?

Yes.

[din14970]

I don't have such a strong opinion about it, I just found it a bit awkward that division was not supported at all. I see your point that mathematically, right side division by a vector indeed doesn't make much sense. On the other hand, most users will be familiar with numpy and expect similar behavior; they may then find it annoying that you can't just do

1 / vector

to get element wise inverse, which might be a useful operation.

In addition, if one wants to stick to mathematical rigor, one could also ask questions about operations that are already supported like:

[1, 2, 3, 4] * vector

It makes more sense, but it's still stretching the * operator in my opinion.

I leave it up to you, but if it should be removed my opinion is that there should be an obvious way to act on the data in an element-wise fashion. To take the inverse of all the elements I think most users don't want to do:

Vector3d(1/vector.data)

have a different syntax (./), but in python we don't have that

The solution could be just that: implement element-wise behavior on the // operator?

[pc494]

This turned out to be a bit more tricky than I had anticipated. However the operation you described only works if the array and the shape of the vector are the broadcastable, which avoids some more unreasonable set ups. For example:

In [27]: v                                                                      
Out[27]: 
Vector3d (2, 2)
[[[0.736  0.2595 0.7069]
  [0.1947 0.9908 0.148 ]]

 [[0.18   0.4939 0.5601]
  [0.9818 0.6614 0.4525]]]

In [28]: [1,2,3,4] * v      
ValueError: operands could not be broadcast together with shapes (4,1) (2,2,3) 

Leaning towards making an issue and having a more comprehensive set of thoughts about how these ideas should work.

[din14970]

Maybe we can settle on the following:

• I remove the controversial right side divide by a vector __rtruedivide__ and just put in a not-implemented - division by a vector is undefined.
• I leave division of vectors by other things __truedivide__ so that the same rules as __mul__ are applied
• I make an issue about how to best deal with element-wise division, taking element-wise inverse etc.

[pc494]

Looks good to me, @hakonanes if you're happy could you merge this please.

[hakonanes]

Looks good.

I added your name to the package credits, @din14970. Feel free to remove it if you don't want it there.

Edit: I'll merge after you've said your OK with it or you've removed it.

[din14970]

@hakonanes Thanks a lot, although this is a super tiny addition. Hope to make some more meaningful contributions in the future.

[hakonanes]

@hakonanes Thanks a lot, although this is a super tiny addition. Hope to make some more meaningful contributions in the future.

I hope so as well!