Update to the HarmonicOscillatorWaveFunction with a wave function recurrence that uses Cython. #2553
Conversation
…urrence that uses Cython.
|
I will resolve the issues! Matheus. |
There was a problem hiding this comment.
Thank you for the contribution.
Could you add a towncrier entry?
Some tests would also be nice. I can't easily tell if you are computing the same thing looking at the code (not that I can confirm there is no bug in the original version.)
Also it would be better to have the cython code in qutip/_distribution.pyx. A new folder for one file feels overkill.
qutip/distributions.py
Outdated
|
|
||
| self.data += k * psi.data_as()[n, 0] | ||
| self.data += ( | ||
| wc.psi_n_single_fock_multiple_position_complex( |
There was a problem hiding this comment.
We have a scaling factor omega. It should be kept (or removed from every Distribution.)
There was a problem hiding this comment.
I created the _distribution.pyx method, I will provide tests and towncrier entries.
qutip/distributions.py
Outdated
| """ | ||
| This class uses a function from the wavefunction_cython | ||
| module, which is a compiled version for the Numpy and | ||
| Cython version of QuTip, from the wavefunction_cython | ||
| module of the Fast_Wave <https://github.com/fobos123deimos/fast-wave> | ||
| package, which is released under | ||
| the BSD license, with the following copyright notice: | ||
| CORDEIRO, MATHEUS; BEZERRA, I. P. ; VASCONCELOS, H. H. M. . Efficient | ||
| Computation of the Wave Function ψn(x) using Hermite Coefficient Matrix in | ||
| Python. In: 7º Workshop Escola de Computação e Informação Quântica | ||
| (7ª WECIQ),2024, Rio de Janeiro.ANAIS DO 7º WORKSHOP-ESCOLA DE COMPUTAÇÃO | ||
| E INFORMAÇÃO QUÂNTICA. Rio de Janeiro: CEFET/RJ, 2024. p. 56-60. | ||
| """ |
There was a problem hiding this comment.
This docstring is not added to the Distribution one, but replace it.
So it should have a description of what it is, it's parameters etc.
It's fine to add reference to your paper, but that not as the only description of the class.
| # BSD 3-Clause License | ||
| # | ||
| # Copyright (c) 2024, Matheus Gomes Cordeiro | ||
| # All rights reserved. | ||
| # | ||
| # Redistribution and use in source and binary forms, with or without | ||
| # modification, are permitted provided that the following conditions are met: | ||
| # | ||
| # 1. Redistributions of source code must retain the above copyright notice, this | ||
| # list of conditions and the following disclaimer. | ||
| # | ||
| # 2. Redistributions in binary form must reproduce the above copyright notice, | ||
| # this list of conditions and the following disclaimer in the documentation | ||
| # and/or other materials provided with the distribution. | ||
| # | ||
| # 3. Neither the name of the copyright holder nor the names of its | ||
| # contributors may be used to endorse or promote products derived from | ||
| # this software without specific prior written permission. | ||
| # | ||
| # THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" | ||
| # AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE | ||
| # IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE | ||
| # DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE | ||
| # FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL | ||
| # DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR | ||
| # SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER | ||
| # CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, | ||
| # OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE | ||
| # OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. | ||
|
|
There was a problem hiding this comment.
Would you be comfortable to be simply included as "QuTiP developers and contributors" of the main BSD 3-Clause License of Qutip?
Having you as the author of the file give the feeling that no one else should touch it and we don't want that.
If merged, we will eventually need to modify this file, maybe merge it or split it etc.
| @cython.cfunc | ||
| @cython.locals(x_size=np.npy_intp, j=int, i=int, k=int, temp1=complex, temp2=complex) | ||
| @cython.boundscheck(False) | ||
| cpdef np.ndarray[np.complex128_t, ndim=1] psi_n_single_fock_multiple_position_complex(int n, np.ndarray[np.complex128_t, ndim=1] x): |
There was a problem hiding this comment.
This is the only function used...
Can you remove the others?
| cimport numpy as np | ||
| import numpy as np | ||
| from libc.math cimport exp, sqrt, pi | ||
| from cmath import exp as cexp, sqrt as csqrt, pi as cpi |
There was a problem hiding this comment.
To import the c++ versions:
| from cmath import exp as cexp, sqrt as csqrt, pi as cpi | |
| cdef extern from "<complex>" namespace "std" nogil: | |
| double complex cexp "exp"(double complex x) | |
| double complex csqrt "sqrt"(double complex x) |
There was a problem hiding this comment.
I had trouble solving this, so I left it as is.
| @cython.nogil | ||
| @cython.cfunc | ||
| @cython.locals(x_size=np.npy_intp, j=int, i=int, k=int, temp1=complex, temp2=complex) | ||
| @cython.boundscheck(False) |
There was a problem hiding this comment.
You could add @cython.cdivision(True) to run in pure cython.
| for j in range(x_size): | ||
| result[0, j] = (cpi ** (-0.25)) * cexp(-(x[j] ** 2) / 2) |
There was a problem hiding this comment.
If you want to optimize as much as possible, you could pre-compute the power temp1 = cpi ** -0.25 / 2.
But I wouldn't be surprised the compiler do it with the optimization we use (-O3).
| @@ -20,6 +20,7 @@ | |||
|
|
|||
| from . import isket, ket2dm, state_number_index | |||
| from .wigner import wigner, qfunc | |||
| from _distributions import psi_n_single_fock_multiple_position_complex | |||
There was a problem hiding this comment.
| from _distributions import psi_n_single_fock_multiple_position_complex | |
| from ._distributions import psi_n_single_fock_multiple_position_complex |
There was a problem hiding this comment.
Thank you! Could you provide your test code so I can use it as a reference to find the problem with this issue?
There was a problem hiding this comment.
I managed to identify the error.
There was a problem hiding this comment.
The problem was the decorator @cython.cdivision(True). I removed it.
There was a problem hiding this comment.
With cdivision(True), we need to be more careful about types.
Still having test would be good.
There was a problem hiding this comment.
I will do the tests, but is it ok to remove cdivision(True)?
Checklist
Thank you for contributing to QuTiP! Please make sure you have finished the following tasks before opening the PR.
You can use pycodestyle to check your code automatically
docfolder, and the notebook. Feel free to ask if you are not sure.doc/changes/<PR number>.<type>'type' can be one of the following: feature, bugfix, doc, removal, misc, or deprecation (see here for more information).Delete this checklist after you have completed all the tasks. If you have not finished them all, you can also open a Draft Pull Request to let the others know this on-going work and keep this checklist in the PR description.
Description
I am adding an efficient module for calculating the wave function that uses a recurrence for it in Cython, and I am applying one of its functions in the update method of the HarmonicOscillatorWaveFunction class in the distribution.py module. This module is an adaptation of one of the modules from the package I developed for efficient wave function calculation: Fast Wave.
Related issues or PRs
Please mention the related issues or PRs here. If the PR fixes an issue, use the keyword fix/fixes/fixed followed by the issue id, e.g. fix #1184