2 tests fail in python 3.4 on Debian 8 about at HOT 14 CLOSED

Hi,

The warnings at the end describe why the tests failed:

RuntimeWarning: invalid value encountered in sqrt
    emit3 = numpy.sqrt(numpy.array([det(r66[s, s]) for s in _submat]))

RuntimeWarning: invalid value encountered in sqrt
    [det(r44[s, s], check_finite=False) for s in _submat[:2]]))

As for the cause of the warnings, I am not entirely sure. I previously encountered the warnings but put it at the bottom of my todo list as I didn't think it affected anything; I left the following note:

• "Diamond's lattice configuration causes emit3[1] to fail an evaluation as it is a numpy.nan because one of the values trying to be square-rooted is very small and so accidentally negative. - Do we use the abs value, do we check for nans and evaluate accordingly, is it supposed to be -ve/very small, should I test this on other lattices?"

I first encountered this on RHEL 7 with python 2.7.3 and numpy 1.11.1, I will see what other environments I can replicate it in; is there any chance that you could try either python 3.4 on Debian 9 or python 3.5 on Debian 8 so that we might establish whether operating system/python version plays a role. Could you also please let me know what numpy version you are using.
Perhaps @lfarv might be able to shed light on this? In the meantime, I will take another look at it.

from at.

Ok, I think I have found the cause. in the calculation of emit3 numpy.array([det(r66[s, s]) for s in _submat]) gave me [ 9.35924096e-18 1.47527700e-75 9.72176963e-12], now 1.47527700e-75 is below 1.18e-38 and so is denormal. It seems that numpy's handling of denormal numbers leaves much to be desired. For me, the following seemed to fix the problem:

hold = numpy.array([det(r66[s, s]) for s in _submat])
hold[hold<1.2e-38] = 0.0
emit3 = numpy.sqrt(hold)
#emit3 = numpy.sqrt(numpy.array([det(r66[s, s]) for s in _submat]))

I had great trouble consistently replicating the issue, I have two environments with the same python and numpy versions on the same machine and yet only one raised the warnings. That said, the above fix did stop the warnings, in that environment, and didn't seem to break anything else, in any environment. Please let me know if this works for you, especially if you are easily able to consistently get warnings.
Perhaps someone that understands the calculations in ohmi_envelope could come up with a more elegant solution, or explain, to me, what is supposed to be going on in the emit3 calculation.

from at.

It seems like my original comment about it being a negative value may have been correct. My fix above also, accidentally, corrects for negative values, but using

hold = numpy.array([det(r66[s, s]) for s in _submat])
bugged_indices = numpy.intersect1d(numpy.where(hold>0.0)[0], numpy.where(hold<1.2e-38)[0])
if bugged_indices:
    hold[bugged_indices] = 0.0
emit3 = numpy.sqrt(hold)

the warnings return. Sure enough from a print statement I get [9.35975251e-18, -4.27088931e-75, 9.72176961e-12], maybe the positive value that I got previously was a fluke? I have however found a very elegant possible fix, as interestingly if you specify the precision of the array as numpy.float32 it works.

emit3 = numpy.sqrt(numpy.array([det(r66[s, s]) for s in _submat], dtype=numpy.float32))

This is because it converts any values below the float32 precision threshold to 0.0 and since our troublesome values always seem to be below that threshold they are converted and so our problem is resolved.
Someone who understands the calculation, is this a viable fix or would the loss of precision matter, also why are we getting such small and or negative values?

Edit: it appears that in my working environments I am consistently getting 1.47e-75, but in the environment that produces warnings the value is -4.27e-75, why could this be?

from at.

Hello @T-Nicholls, I may have helped write the code for the eigenemittance algorithm in the Matlab AT. I was never totally satisfied with it. I thought I understood well how to go from a sigma matrix to the three eigenemittances, but there are some ambiguities that I never totally clarified.
Unfortunately, I haven't had much time to follow all the great develpments with pyat here.
Can you provide me a link to where this calculation is done in the repository, and I can try to have a look?

update:
I found it. I'll have a look.

from at.

Hello all. @bnash, it's nice to hear from you and see you are still following what's happening!
We are dealing here with a lattice without coupling. So the true result is a perfect zero. Whatever comes out, be it positive or negative, results from numerical errors. And a range of x.e-75 looks perfectly reasonable. So for me 0.0 or Nan does not make a big difference. I do not like the idea of cheating with the results:

• rounding to 32 bits unnecessarily spoils any result, I would avoid,
• taking the absolute value is even worse: it's wrong. Only truncating negative values to zero is acceptable, since we know that the square of the emittance must be positive or null.

@bnash, you are welcome to look back at the algorithm, but I am not sure that we can get rid of so small numerical errors. Otherwise, I propose either doing nothing (no vertical emittance for a lattice without coupling is not so stupid), or truncating negative values to zero.

More generally, I wonder why we do not get exactly the same results for any test, at least for those running on 64-bit Intel-like processors (probably all of us). It would be interesting to see at which point it diverges. Tracking through a drift involves only addition and multiplication (done by the processor, no library involved, I guess). Then going through focusing elements uses trigonometric functions (any difference?). @carmignani, apart from a different OS version (which should not have any role here), do you know if it's a different numpy version?

Update: @bnash, if I remember well, the ambiguity was in attributing the three Eigen emittance to the closest plane (x, z or s), and to avoid attributing the same one to two planes. The search could be improved by restricting the choice once a plane has been attributed, but this should not help for numerical errors…

from at.

Thanks @lfarv! I've wanted to stay more involved, but didn't find the time. Thanks for the clarification about this issue here.
I agree that taking the absolute value of the emittance doesn't make sense.

I'm curious how the tests have been set up here in pyat. I've been working on Zgoubi recently and we set up a set of tests to make sure we don't change the physics as we work on the code. The tests are based on ndiff, where you can give a level of tolerance and still consider two results the same. It seems like something like this is necessary since running on different platforms will always result in slightly different numerical values for an algorithm.

from at.

Hello,
I was using the same numpy version in the two cases: 1.16.2.

from at.

Hi @bnash, we mainly use the numpy.testing module for comparing results, it does support tolerances and in this case the numerical errors in the calculations would be well within the limits. However, because numpy's square root outputs nan for negative values the evaluation fails. From the test output:

E           AssertionError: 
E           Arrays are not almost equal to 7 decimals
E           
E           x and y nan location mismatch:
E            x: array([1.3252791e-10,           nan])
E            y: array([1.3252791e-10, 0.0000000e+00])

Thanks all for the explanations/clarifications; if I understand Laurent correctly we are not particularly concerned with the value of that second cell, so one option for fixing the tests would be to only check the first cell. Another option would be to use a try except to check if it is nan if it fails the original comparison. Or the third option is to fix the cause of the warnings and replace negative values with 0 before they go into the square root. I lean more towards the last, as I don't particularly like warnings for what is expected behaviour, but leave it up to the more knowledgable to decide. Perhaps @lfarv could add the fix for this onto the endo of #99?

from at.

@T-Nicholls: the second item is the vertical equilibrium emittance, which is theoretically zero for the 2 examples run in the test. However this test is meaningful for the 'err.mat' lattice where the vertical emittance is non-zero (this is why I am pushing to have this lattice in the test).
So if everyone agrees, I'll replace negative values with zeros in #99

from at.

Sounds good @lfarv, the timing discrepancy for real lattices with the Matlab engine is still on my list; hopefully, I'll get a chance to take another look at it soon.

from at.

Could you put the fix on a new branch so that @carmignani can test it separate to the plotting changes?

from at.

Done in #105

from at.

This should now be fixed on the master branch, could you please try again @carmignani and @abdomit?

from at.

Hello,
yes, I confirm. Now all the tests pass in both debian8+python3.4 and debian9+python3.5.

from at.