qojulia / QuantumOptics.jl

Library for the numerical simulation of closed as well as open quantum systems.
http://qojulia.org
Other
533 stars 104 forks source link

Return concise error message if Ket coefficients are not compatible with timeevolution methods #318

Closed dehond closed 3 years ago

dehond commented 3 years ago

Added a method check_psi0 to timeevolution_base that uses Julia's type matching to check whether the coefficients of the initial wave function in time evolution methods are complex. If they are not, a short error message is generated that tells the user the coefficients should be converted to complex.

I've added this check to the schroedinger, master, and mcwf methods. During some unit tests, schroedinger methods are called with non-Ket datatypes as initial wave functions, so there's a generic check_psi0 function that's called if the input is not of the Ket type. To make sure it passes the unit tests in test_timeevolution_abstractdata.jl I had to specify the input datatype to be

check_psi0(psi0::Ket{B, T} where B<:Basis where T<:Array)

instead of just

check_psi0(psi0::Ket)

In a sense the function is now checking the type twice, once to verify it's a Ket{B, T}, and then once more to make sure T is actually an array with complex coefficients. This seems suboptimal to me, but I couldn't immediately figure out a way to do it in one go. (For instance by having the value of the T that is matched in the function header available the function body.)

fixes #317

bastikr commented 3 years ago

I'm pretty sure this change is against what @david-pl tries to achieve - so please let's not merge this until he reviews this PR. Originally we took an approach more in the line of this PR, i.e. forcing many things to be of a specific type like states being stored as Vector{ComplexF64}. This has a few advantages like better error messages and guaranteeing that everything is reasonable fast since we can be sure the user does not use abstract types by accident (https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-abstract-container). However, the restriction to a specific type is something which is avoided in the Julia world since it harms flexibility.

ChristophHotter commented 3 years ago

@bastikr Okay, my first thought was also that it harms the flexibility but on the other hand it seemed to be very convenient to get a better error message in this special case. Anyway, I will discuss it with @david-pl.

david-pl commented 3 years ago

@dehond Sorry for the delay. As @bastikr said, I think implementing this change may harm flexibility, so I would argue against merging it. I will try to explain why:

Interoperability between Julia packages is based on using generic types. You already found one example of this, which is the test_timeevolution_abstractdata.jl, where Kets are used that have data which acts like an array but isn't a subtype to Array. While there may not be an obvious use-case for this, it might come in handy when working with special data structures. An example at the level of Operator would be that one can use a LinearMap (from LinearMaps.jl) as data which is useful sometimes.

While you took care of this properly by implementing the generic fallback method check_psi(x)=nothing, there could still be a case where one uses a special number type in the data array that is <:Number, but not <:Complex, while it still works with complex numbers using type promotion. In such a case, your check would falsely error. While I admit that this is quite a small edge case, which we could argue is not important, all we gain from this change is a different error message. That's why I would say it's probably not worth it.

Furthermore, the issue #317 all comes down to the following error

using LinearAlgebra
psi, dpsi = rand(5), rand(5)
H = rand(5,5)
mul!(dpsi,H,psi,-im,0)

which will reproduce your InexactError. The ugly stacktrace you report in #317 is because the error is thrown inside the time evolution. Long and incomprehensible stacktraces are an issue within Julia itself. If we start working around these ugly error messages, then where do we stop? Ugly stacktraces are thrown all the time due to comparably simple errors. I'd rather wait for a proper fix from the Julia developers to this issue.

I hope these explanations are sufficiently clear. In any case, thanks for the effort you put in this PR @dehond

dehond commented 3 years ago

@david-pl thanks for your detailed explanation. I was worried this could be something that throws off a new user of what I feel is otherwise a very user-friendly and intuitive package. But if, as you point out, the problem is rather Julia's incomprehensible stack traces in general, it doesn't make much sense indeed to code around this for every corner case that does something like this. I can see how the resulting code will be hard to understand and maintain. Thanks for having a look at my PR, though!