Wigner lacks precision

[GoogleCodeExporter]

What steps will reproduce the problem?
1. Calculate a Wigner function of a state in the high-dimensional Hilbert space:

xvec=linspace(-10,10,101)
yvec=xvec
X,Y = meshgrid(xvec, yvec)
psi=coherent(100,7)
contourf(X, Y, wigner(psi,xvec,yvec),100)

2. Observe the huge noise.

What is the expected output? What do you see instead?

There should be a nice coherent state on the righthand side of the image.
The iterative method produces noise on the left. The new laguerre method gives 
error.

What version of the product are you using? On what operating system?
QuTip 1.3, 2.0, 2.1. Mac OS 10.8, Linux Mint

Please provide any additional information below.

I guess the problem here consists of two parts. First part is likely related to 
the truncated Hilbert space which produces ghost features of the Wigner 
function for a certain number of dimensions. One can increase the number of 
dimensions or do something more clever to overcome it. Another part of the 
problem is the calculation precision. The noise here originates from it. I have 
changed the wigner procedure to use numpy.complex256 and float128, it helped a 
lot in this case. However, the outcome is still not perfect. It appears that on 
mac os and linux numpy.float128 is effectively only float96. The ultimate 
solution would be to use the arbitrary precision of mpmath.

The modified version of wigner.py is attached. It is from qutip 1.3.

Original issue reported on code.google.com by D.Vutshi on 15 Nov 2012 at 10:15

Attachments:

• wigner.py

[GoogleCodeExporter]

Original comment by nonhermitian on 19 Nov 2012 at 2:44

• Changed state: Accepted

[GoogleCodeExporter]

I did the necessary modification to the code. Now it uses gmpy2 library. 
Performance is much lower than the original version, however it is much faster 
than my first implementation with mpmath.

One can use it by importing the file (first install gmpy2 and all necessary 
libraries):
sys.path.append('/path to the file')
from wigner_gmpy2 import *
and now just call the function
wigner_precise_gmpy2(psi, xvec, yvec, g=sqrt(2), ddp=100)
ddp is the precision in bits, the standard double precision is 53 bits.

Original comment by D.Vutshi on 23 Nov 2012 at 2:53

Attachments:

• wigner_gmpy2.py

[GoogleCodeExporter]

Thank you for letting us know about this issue. I will add the code you sent to 
the new plugins page that we are making for the QuTiP website.  At some point I 
might add support for the mpmath module into the QuTiP Wigner function, however 
I am not aware of many people having an use for such a feature given that the 
goal is usually to get a small number of quanta.

Original comment by nonhermitian on 27 Nov 2012 at 1:51

• Changed state: Done

[GoogleCodeExporter]

mpmath is so slow that it is close to useless in this case. In my tests it 
takes about 5-10 minutes to calculate Wigner function for 20x20 points and 
Ndim=30 with mpmath. It takes the same time for gmpy2 implementation to produce 
Wigner function for 100x100 points Ndim=100. Mpmath would take hours in this 
case. Now I think that the best solution would be to use the fast arbitrary 
precision libraries (MPFR and MPC) directly from Cython or Fortran via f2py. 
Also a good compromise would be to use a fast quad precision library from C or 
Fortran.

Original comment by D.Vutshi on 27 Nov 2012 at 9:40

[GoogleCodeExporter]

mpmath will automatically detect gmpy, if it exists, and use it in its 
calculations.  Your current suggestions seems like a bit too much work for 
addressing the current problem.  Programming in mpmath, and using the parfor 
function, should be a good enough solution for anyone who runs into this 
particular problem provided they are not asking for too high a precision.  If 
someone needs a specialty tool for calculating these things a lot, then perhaps 
your method is the way to go.  I am currently writing the Laguerre based Wigner 
function in mpmath and have already added the option to run it in parallel.  
Hopefully, this will be good enough for now.

Original comment by nonhermitian on 5 Dec 2012 at 5:25

• Changed state: Started
• Added labels: Type-Enhancement
• Removed labels: Type-Defect

[GoogleCodeExporter]

This issue has been addressed by adding an FFT based Wigner method that does 
not use special functions.  This is available in the development version on 
github.

Original comment by nonhermitian on 11 Jul 2013 at 2:06

• Changed state: Fixed