Why is broadcasting restricted?

[Krastanov]

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]

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

[Krastanov]

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]

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]

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]

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.

[Krastanov]

with PR 172 from QuantumOpticsBase this can be now closed

[apkille]

Note: broadcasting for semiclassical.State was dealt with in #404.