Broadcasting Kets and Operators

[apkille]

@david-pl @Krastanov @amilsted Improves performance and reduces allocations in #16 by defining RecursiveArrayTools methods and Base.Broadcast._broadcast_get_index for QO types.

For the example discussed in #16 (copied below)

using QuantumOptics, DifferentialEquations, BenchmarkTools

ℋ = SpinBasis(100//1)
↓ = spindown(ℋ)
t₀, t₁ = (0.0, pi)
const σx = sigmax(ℋ)
const iσx = im*σx
schrod!(dψ, ψ, p, t) = QuantumOptics.mul!(dψ, iσx, ψ)
prob! = ODEProblem(schrod!, ↓, (t₀, t₁))

const ix = iσx.data
schrod_data!(dψ,ψ,p,t) = QuantumOptics.mul!(dψ, ix, ψ)
const u0 = (↓).data
prob_data! = ODEProblem(schrod_data!, u0, (t₀, t₁))

the performance is now as follows.

QO Ket:

julia> @benchmark sol = solve($prob!, DP5(); save_everystep=false)
BenchmarkTools.Trial: 2224 samples with 1 evaluation.
 Range (min … max):  2.230 ms …  2.409 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     2.243 ms              ┊ GC (median):    0.00%
 Time  (mean ± σ):   2.245 ms ± 13.415 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

      ▂▃▅▇█▄▇▇█▅▂                                             
  ▂▃▄▇████████████▇▅▅▄▄▄▄▄▃▃▃▃▃▃▂▃▂▃▂▂▂▂▂▂▂▂▁▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂ ▄
  2.23 ms        Histogram: frequency by time         2.3 ms <

 Memory estimate: 49.86 KiB, allocs estimate: 49.

Julia array:

julia> @benchmark sol = solve($prob_data!, DP5(); save_everystep=false)
BenchmarkTools.Trial: 2347 samples with 1 evaluation.
 Range (min … max):  2.107 ms …  2.274 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     2.127 ms              ┊ GC (median):    0.00%
 Time  (mean ± σ):   2.129 ms ± 12.497 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

          ▂▃▃▄▅█▅▇▆▅▄▅▃▄▁▁                                    
  ▂▁▃▃▄▅▇█████████████████▇█▆▆▄▅▅▄▃▃▃▃▃▃▃▂▃▂▃▂▂▂▂▂▂▂▁▂▂▁▂▂▂▂ ▄
  2.11 ms        Histogram: frequency by time        2.18 ms <

 Memory estimate: 56.34 KiB, allocs estimate: 47.

[Krastanov]

Just to make sure it is clarified: this is strictly an improvement over current capabilities, it only makes changes to the broadcast capabilities of our types, and it does not hook these new broadcasting capabilities in the time evolution solvers.

[apkille]

Looks like the manual piracy test is failing because of the new broadcasting interface. I'm not sure what this orgs policy on type piracy is.

[apkille]

I should also mention: if this PR is approved, a lot of the defined methods can probably be moved to QuantumInterface. I'm happy to do this, but it's up to the maintainers if we want to do this now or later (or at all 😄).

[david-pl]

@apkille @Krastanov the changes here look fine to me. Not sure what causes the piracy failure, but I'll leave that up to you.

I'd suggest we merge this PR (and by "we" I mean @Krastanov once you're happy with everything) and leave the move to QuantumInterface for another time (should be done though IMO). Do you guys agree?

[apkille]

Type piracies should be fine now.

I'd suggest we merge this PR (and by "we" I mean @Krastanov once you're happy with everything) and leave the move to QuantumInterface for another time (should be done though IMO). Do you guys agree?

Sure, that works for me.

[codecov]

Codecov Report

Attention: Patch coverage is 65.21739% with 24 lines in your changes missing coverage. Please review.

Project coverage is 92.13%. Comparing base (bb71f52) to head (145ddab).
Report is 12 commits behind head on master.

Files	Patch %	Lines
src/states.jl	59.52%	17 Missing ⚠️
src/operators_dense.jl	74.07%	7 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #172      +/-   ##
==========================================
- Coverage   93.00%   92.13%   -0.88%     
==========================================
  Files          25       26       +1     
  Lines        3104     3231     +127     
==========================================
+ Hits         2887     2977      +90     
- Misses        217      254      +37     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[apkille]

@Krastanov

[apkille]

@david-pl Is bumping the Julia compat to 1.10 OK with you? We are getting compat conflicts with RecursiveArrayTools otherwise.

[Krastanov]

FYI, 1.10 will be the new LTS in about a month or so, when 1.11 is released.

[Krastanov]

This looks good to me. I am requesting one minor change and it is good to merge.

TODOs:

• @david-pl , are you okay with bumping the minimal compat to the new Jula LTS?
• @apkille, you mentioned adding some broadcast benchmarks. I assume they will be in a separate PR and are not blocking for this
• @apkille, I suggested a minor import rewording to avoid polluting the namespace.

[Krastanov]

Oh, one more thing. I see a big part of the additions are not covered by the CI. Could you put your examples from the first comment in the test runner, maybe as a new file test_sciml_broadcast_interfaces.jl. I imagine some of this will be potentially added to a benchmark runner, but we would want this tested in the meanwhile.

[apkille]

@Krastanov should be good to go. About the benchmarks, I am adding them in a separate PR, which will be submitted later today or early tomorrow.

[Krastanov]

@apkille , could you check whether the new JET warnings come from something you have done? If not, could you bump the permitted number of warnings, to have this test pass?

@david-pl , @amilsted , @ChristophHotter , let me know if there are any worries about the bump to the new julia LTS as a lower bound. If not, I will merge this in about a week (do not hesitate to request a longer waiting period).

[apkille]

@Krastanov two of the JET warnings were because of me, and they were simple fixes. The rest were not from changes in this branch, so I bumped the permitted warnings from 24 to 28.

[Krastanov]

Merged! Thank you for contributing this!