Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Why is broadcasting restricted? #299

Open
Krastanov opened this issue Apr 3, 2021 · 5 comments
Open

Why is broadcasting restricted? #299

Krastanov opened this issue Apr 3, 2021 · 5 comments

Comments

@Krastanov
Copy link
Collaborator

I see that operations like f.(ket) are explicitly set to raise an error: https://github.com/qojulia/QuantumOpticsBase.jl/blob/201b6305440b456b4b1936356e82735cfe6b2f90/src/states.jl#L236

This seems much too heavy handed to me. From my reading on the discourse forum it seems the general advise is to never use broadcasting for anything but "objects with indices and axes". Thus, what I was expecting is that broadcasting would work directly on the raw .data field, without any overhead or other considerations. It feel like the current broadcasting implementation in QuantumOptics is in an uncomfortable in-between state: it is not a dumb "just give me the indexed object" implementation and it is not fully a "permit only meaningful quantum algebra". Given that broadcasting should explicitly not be the latter, I am asking whether the developers could consider removing much of these restrictions and just enable standard boring julia broadcasting?

One immediate benefit from this course of action would be drastically simplifying the DifferentialEquations interactions. I was surprised to see how there were layers of "recasting" and callbacks for saving. All of this code will become unnecessary (and the result would probably be a bit less overhead) if boring broadcasting is enabled.

A few other definitions that might be necessary for that goal:

Base.getindex(k::Ket, i) = getindex(k.data, i)
Base.setindex!(k::Ket, i, d) = setindex!(k.data, i, d)

Base.eltype(::Type{Ket{B,A}}) where {B,N,A<:AbstractVector{N}} = N
Base.vec(k::Ket) = vec(k.data)
Base.zero(k::Ket) = typeof(k)(k.basis, zero(k.data))

TLDR: Julia broadcasting is supposed to be dumb operation on indexed objects unaware of any algebraic properties. That would enable the DifferentialEquations library to do more of the heavy lifting.

@Krastanov
Copy link
Collaborator Author

These comments all stem from my exploration of the code after I started looking into implementing #298

@Krastanov
Copy link
Collaborator Author

My excitement over unrestricted broadcasting seems to be a bit too optimistic. The forum post gives a few more details and constraints https://discourse.julialang.org/t/what-is-the-interface-necessary-for-differentialequations-jl-if-your-data-type-does-not-subclass-abstractarray/58480/2

@david-pl
Copy link
Member

david-pl commented Apr 6, 2021

I agree that the broadcasting is in an awkward state, and in principle I wouldn't have a problem with substantially changing it/removing restrictions. The motivation behind implementing this was that otherwise in-place updates of operators would error due to Julia's scoping since 1.0. Other than that broadcasting was never really of concern because you could just always do it on the .data fields.

Julia broadcasting is supposed to be dumb operation on indexed objects unaware of any algebraic properties.

An important aspect for types in QuantumOptics is keeping track of the Hilbert spaces (bases). This is especially useful for PositionBasis and MomentumBasis which describe the same system but in two different bases. If you base broadcasting solely on the indexing, this will break down and you will be able to add operators defined on different Hilbert spaces.

The restriction is there to keep the broadcasting limited to functions that are sure to make sense when used together with operators, although that is far from pretty and should be changed.

@david-pl
Copy link
Member

david-pl commented Apr 6, 2021

One immediate benefit from this course of action would be drastically simplifying the DifferentialEquations interactions. I was surprised to see how there were layers of "recasting" and callbacks for saving. All of this code will become unnecessary (and the result would probably be a bit less overhead) if boring broadcasting is enabled.

This will not be so simple for a semiclassical.State consisting of a density operator (matrix) and a few classical variables (vector), where you actually need to rewrite the state at each time step. The whole "recasting" thing was introduced with the semiclassical module IIRC.

@Krastanov
Copy link
Collaborator Author

I have not looked into semiclassical, but I think the DifferentialEquations.jl team has already solved many of these issues with their special array types. Maybe that would be a fruitful approach.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants