FABLE-3DXRD / xfab

Other
6 stars 7 forks source link

ubi_to_u fails to produce a rotation matrix for 32 bit inputs #33

Closed jonwright closed 2 years ago

jonwright commented 2 years ago
import xfab.tools
print(repr(ubi))
xfab.tools.ubi_to_u( ubi )

So I think it fails because the argument was a 32 bit array ?

array([[ 3.5964367 , -1.7651347 , -0.7474788 ],
       [ 1.7880325 ,  2.5086281 ,  2.6700249 ],
       [-0.69801915, -2.6838543 ,  2.9860113 ]], dtype=float32)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
~\AppData\Local\Temp\ipykernel_18864\1475499514.py in <module>
      1 import xfab.tools
      2 print(repr(ubi))
----> 3 xfab.tools.ubi_to_u( ubi )

~\Anaconda3\envs\py39\lib\site-packages\xfab\tools.py in ubi_to_u(ubi)
    421     B = form_b_mat(unit_cell)
    422     U = n.transpose(n.dot(B, ubi))/(2*n.pi)
--> 423     if CHECKS.activated: checks._check_rotation_matrix(U)
    424 
    425     return U

~\Anaconda3\envs\py39\lib\site-packages\xfab\checks.py in _check_rotation_matrix(U)
     19     """
     20     if not np.allclose( np.dot(U.T, U), np.eye(3,3) ):
---> 21         raise ValueError("orientation matrix U is not unitary, np.dot(U.T, U)!=np.eye(3,3)")
     22 
     23     if not np.allclose( np.linalg.det(U), 1.0 ):

ValueError: orientation matrix U is not unitary, np.dot(U.T, U)!=np.eye(3,3)
AxelHenningsson commented 2 years ago

I get the same behaviour. Seems the absolute error is something like 1e-07. There's a lot of flops going on in these tools functions, so no surprise precision gets lost. I suggest to lower atol in the np.allclose( np.dot(U.T, U), np.eye(3,3), atol=1e-6) orthonormality check. Going to 1e-06 seems to support 32bit float safely but still fails on 16bit. But then again, one should probably not compute these things with anything below 32bit precision considering the amount of flops going on...

PR on way. 👍

jonwright commented 2 years ago

Pragmatically, I was thinking we could either:

Both of those seem plausible. It is more a question of whether the library should be double precision only (like javascript) or whichever precision people need (like C++ templates).

jonwright commented 2 years ago

For my own use, I would lean towards the templates approach. For serving a wider community it might be easier to just use doubles everywhere.

AxelHenningsson commented 2 years ago

I am not a big fan of casting at the start of every function, seems nice to be able to use whatever precision one has the need for. However, the question of how to do the checks then remains. One could cast the U matrix to the same type as UBI as you say, and then have an additional logic in the check and adapt the atolbased on the dtype, but that gets a bit dirty I think. Fundamentally, if we have lower precision we expect to also have to impose a looser bound on the checks though.

On another note: perhaps it would be better to check the ubi matrix directly in xfab.tools.ubi_to_u().

AxelHenningsson commented 2 years ago

Seems this would be a better check for a ubi matrix in general?

if np.dot(ubi[2,:], np.cross(ubi[0,:],ubi[1,:]))<0:
    raise ValueError("ubi matrix must hold lattice vectors forming a feasible unit cell")
AxelHenningsson commented 2 years ago

Seem the question that needs answering before any new PR can be suggested is if we want to do checks for dtypes in every single xfab.tools function and then react to those types in the checks module and adapt tolerances accordingly.

Another option is to have a tolerance that can be set, i.e xfab.CHECKS.atol, which defaults to 1e-8 but can be set depending on the precision needed. This would sort of suggest to people to use 64bit but still allow other types if needed. What do you think of that?

jonwright commented 2 years ago

I think the code in checks was fine? numpy.allclose is usually doing the right thing by default for floating point.

For the left handed ubi matrices the test seems to be np.linalg.det(ubi)>0. Not sure whether zero should be allowed (e.g. truncation rods).

There is indeed a decision to make about the API. Is the returned type a "python float" precision numpy array, or does the returned type match the input type?

I guess the intention is for this to work for typical python numbers (tuple/list of floats), so probably it should go with doubles everywhere and cast arguments at the start of functions? This does rule out using float128 however.

AxelHenningsson commented 2 years ago

I think the code in the checks will not adapt depending on dtypes? For 32bit floats it seems atol must be lowered. The problem is not that the input gets casted to 64bit, the problem is that the computation made in xfab.tools actually loses the precision and produces a poor rotation matrix.

jonwright commented 2 years ago

You are right - numpy.allclose is not doing what I was expecting. It seems like the atol for a rotation matrix could be looking at something like np.finfo( dtype ).resolution ?

For the numerics, I did not get any improvement using a QR, so I assume the algorithm in xfab is pretty good ?

def ubi_to_u_qr( ubi ):    # ubi   == inv(UB) ; ubi.T ==  inv(u).T inv(b).T
    q, r = np.linalg.qr( ubi.T )
    s = np.diag( np.sign( np.diag(r) ) )
    return np.dot(q, s)

It seems like casting the argument avoids almost all of these problems xfab.tools.ubi_to_u(ubi.astype(float)) works fine. This would be a good reason for putting something like arg=np.asarray(arg, dtype=float) near the top of some functions?

AxelHenningsson commented 2 years ago

Yeah, I tried QR as well, seems to be a very minor improvement (although a lot nicer implementation code wise I think).

I know numpy.linalg runs the policy of "only 32bit or better". If you try 16bit with numpy.linalg.qr for instance you will get a type error like: TypeError: array type float16 is unsupported in linalg

This might be nicer for xfab than secretly casting peoples stuff into other stuff? :slightly_smiling_face:

So we can add a dtype check (perhaps as a decorator?) to xfab.tools

Cheers Axel

jonwright commented 2 years ago

The formal way to do typing seems to be annotations: https://numpy.org/devdocs/reference/typing.html ... but I did not learn to use those yet.

With decorators then things like introspection can get confusing.

The casting makes a private copy inside the function:

def ubi_to_u(ubi_input):
    ubi = n.asarray( ubi_input, float )
    unit_cell = ubi_to_cell(ubi)
    B = form_b_mat(unit_cell)
    U = n.transpose(n.dot(B, ubi))/(2*n.pi)
    if CHECKS.activated: checks._check_rotation_matrix(U)
    return U

... it gives back a valid 64bit rotation with no side effects. Currently, without the casting it gives back a 32bit rotation which is stored in a 64bit array.

Changing the API to make it send back a 32bit array seems like it will be more problematic than just computing a correct answer in 64 bit?

AxelHenningsson commented 2 years ago

Yeah, lets go with the casting copy stuff suggestion then I would say. Annotations probably will break python 2.7 anyways.

Will draft a new PR asap I get some time on my hands.

Cheers Axel

AxelHenningsson commented 2 years ago

Unfortunately casting back and forth between 64bit and 32bit will make an orientation matrix not pass a np.allclose()

import numpy as np
u, _ = np.linalg.qr( np.random.standard_normal((3, 3)) )
u = np.asarray(u, np.float32)
u = np.asarray(u, float)
assert np.allcose(u.T.dot(u), np.eye(3,3))

This means that functions such as xfab.tools.u_to_ubi() and xfab.tools.u_to_rod() simply cannot handle anything but 64bit floats, or else they will fail in the xfab.checks. This reflects in some sense the fact that 32bit floats simply gives poor precision for rotation matrix operations.

ps: @jonwright Awaiting review for updated PR.