Bump REQUIRE to v0.7-beta2 (#230)

* Bump REQUIRE to v0.7-beta2

* Fix deprecations

* Fix deprecations

* Remove Compat imports

* Change WignerSymbols version

* First successful build

* Fock tests pass

* Some more tests pass

* All tests except printing pass

* Replace Complex128 with ComplexF64 in tests

* Rename Complex128 to ComplexF64

* Add FFTW requirement

* Update printing

* Add Arpack requirement

* Add rounding to printing

* Update tests to use 0.7

* Fix compilation deprecation warnings

* Fix printing

* Fix deprecations

* Rename full to dense

* Fix some more deprecations

* More deprecations

* Fix operator deprecations

* Fix particle deprecations

* Fix all deprecations occurring in tests

* Update REQUIRE

* Update appveyor

* Fix a bug in subspace; export norm

* Fix silly bug in subspace

.travis.yml

@@ -3,7 +3,7 @@ os:
     - osx
     - linux
 julia:
-    - 0.6
+    - 0.7
     - nightly
 matrix:
     allow_failures:

REQUIRE

@@ -1,7 +1,8 @@
-julia 0.6
-Compat 0.64.0
-OrdinaryDiffEq 3.19.1
-DiffEqCallbacks 1.1
-StochasticDiffEq 4.4.5
-RecursiveArrayTools
-WignerSymbols 0.1.1
+julia 0.7-beta2
+OrdinaryDiffEq 4.7.1
+DiffEqCallbacks 2.1
+StochasticDiffEq 5.6.0
+RecursiveArrayTools 0.17.2
+WignerSymbols 0.2.0
+FFTW
+Arpack

appveyor.yml

@@ -1,7 +1,7 @@
 environment:
   matrix:
-  - JULIA_URL: "https://julialang-s3.julialang.org/bin/winnt/x86/0.6/julia-0.6-latest-win32.exe"
-  - JULIA_URL: "https://julialang-s3.julialang.org/bin/winnt/x64/0.6/julia-0.6-latest-win64.exe"
+  - JULIA_URL: "https://julialang-s3.julialang.org/bin/winnt/x86/0.7/julia-0.7-latest-win32.exe"
+  - JULIA_URL: "https://julialang-s3.julialang.org/bin/winnt/x64/0.7/julia-0.7-latest-win64.exe"
   - JULIA_URL: "https://julialangnightlies-s3.julialang.org/bin/winnt/x86/julia-latest-win32.exe"
   - JULIA_URL: "https://julialangnightlies-s3.julialang.org/bin/winnt/x64/julia-latest-win64.exe"
 

src/QuantumOptics.jl

@@ -2,11 +2,14 @@ __precompile__()
 
 module QuantumOptics
 
+using SparseArrays, LinearAlgebra
+
 export bases, Basis, GenericBasis, CompositeBasis, basis,
         tensor, ⊗, permutesystems,
-        states, StateVector, Bra, Ket, basisstate,
+        states, StateVector, Bra, Ket, basisstate, norm,
                 dagger, normalize, normalize!,
-        operators, Operator, expect, variance, identityoperator, ptrace, embed,
+        operators, Operator, expect, variance, identityoperator, ptrace, embed, dense, tr,
+            sparse,
         operators_dense, DenseOperator, projector, dm,
         operators_sparse, SparseOperator, diagonaloperator,
         operators_lazysum, LazySum,

src/bases.jl

@@ -8,8 +8,6 @@ export Basis, GenericBasis, CompositeBasis, basis,
 
 import Base: ==, ^
 
-using Compat
-
 """
 Abstract base class for all specialized bases.
 
@@ -103,8 +101,8 @@ tensor(b1::Basis, b2::Basis) = CompositeBasis(Int[prod(b1.shape); prod(b2.shape)
 tensor(b1::CompositeBasis, b2::CompositeBasis) = CompositeBasis(Int[b1.shape; b2.shape], Basis[b1.bases; b2.bases])
 function tensor(b1::CompositeBasis, b2::Basis)
     N = length(b1.bases)
-    shape = Vector{Int}(N+1)
-    bases = Vector{Basis}(N+1)
+    shape = Vector{Int}(undef, N+1)
+    bases = Vector{Basis}(undef, N+1)
     for i=1:N
         shape[i] = b1.shape[i]
         bases[i] = b1.bases[i]
@@ -115,8 +113,8 @@ function tensor(b1::CompositeBasis, b2::Basis)
 end
 function tensor(b1::Basis, b2::CompositeBasis)
     N = length(b2.bases)
-    shape = Vector{Int}(N+1)
-    bases = Vector{Basis}(N+1)
+    shape = Vector{Int}(undef, N+1)
+    bases = Vector{Basis}(undef, N+1)
     for i=1:N
         shape[i+1] = b2.shape[i]
         bases[i+1] = b2.bases[i]

src/fock.jl

@@ -5,6 +5,7 @@ export FockBasis, number, destroy, create, displace, fockstate, coherentstate
 import Base: ==
 
 using ..bases, ..states, ..operators, ..operators_dense, ..operators_sparse
+using SparseArrays
 
 
 """
@@ -34,7 +35,7 @@ Number operator for the specified Fock space.
 """
 function number(b::FockBasis)
     diag = complex.(0.:b.N)
-    data = spdiagm(diag, 0, b.N+1, b.N+1)
+    data = spdiagm(0 => diag)
     SparseOperator(b, data)
 end
 
@@ -45,7 +46,7 @@ Annihilation operator for the specified Fock space.
 """
 function destroy(b::FockBasis)
     diag = complex.(sqrt.(1.:b.N))
-    data = spdiagm(diag, 1, b.N+1, b.N+1)
+    data = spdiagm(1 => diag)
     SparseOperator(b, data)
 end
 
@@ -56,7 +57,7 @@ Creation operator for the specified Fock space.
 """
 function create(b::FockBasis)
     diag = complex.(sqrt.(1.:b.N))
-    data = spdiagm(diag, -1, b.N+1, b.N+1)
+    data = spdiagm(-1 => diag)
     SparseOperator(b, data)
 end
 
@@ -65,7 +66,7 @@ end
 
 Displacement operator ``D(α)`` for the specified Fock space.
 """
-displace(b::FockBasis, alpha::Number) = expm(full(alpha*create(b) - conj(alpha)*destroy(b)))
+displace(b::FockBasis, alpha::Number) = exp(dense(alpha*create(b) - conj(alpha)*destroy(b)))
 
 """
     fockstate(b::FockBasis, n)
@@ -82,7 +83,7 @@ end
 
 Coherent state ``|α⟩`` for the specified Fock space.
 """
-function coherentstate(b::FockBasis, alpha::Number, result=Ket(b, Vector{Complex128}(length(b))))
+function coherentstate(b::FockBasis, alpha::Number, result=Ket(b, Vector{ComplexF64}(undef, length(b))))
     alpha = complex(alpha)
     data = result.data
     data[1] = exp(-abs2(alpha)/2)

src/manybody.jl

@@ -9,6 +9,7 @@ import ..fock: number, create, destroy
 import ..nlevel: transition
 
 using ..bases, ..states, ..operators, ..operators_dense, ..operators_sparse
+using SparseArrays
 
 """
     ManyBodyBasis(b, occupations)
@@ -26,7 +27,7 @@ mutable struct ManyBodyBasis <: Basis
     occupations_hash::UInt
 
     function ManyBodyBasis(onebodybasis::Basis, occupations::Vector{Vector{Int}})
-        new([length(occupations)], onebodybasis, occupations, hash(occupations))
+        new([length(occupations)], onebodybasis, occupations, hash(hash.(occupations)))
     end
 end
 
@@ -63,8 +64,8 @@ Return a ket state where the system is in the state specified by the given
 occupation numbers.
 """
 function basisstate(basis::ManyBodyBasis, occupation::Vector{Int})
-    index = findfirst(basis.occupations, occupation)
-    if index == 0
+    index = findfirst(isequal(occupation), basis.occupations)
+    if isa(index, Nothing)
         throw(ArgumentError("Occupation not included in many-body basis."))
     end
     basisstate(basis, index)
@@ -227,7 +228,7 @@ where ``X`` is the N-particle operator, ``x`` is the one-body operator and
 ``|u⟩`` are the one-body states associated to the
 different modes of the N-particle basis.
 """
-function manybodyoperator(basis::ManyBodyBasis, op::T)::T where T <: Operator
+function manybodyoperator(basis::ManyBodyBasis, op::T)::T where T<:Operator
     @assert op.basis_l == op.basis_r
     if op.basis_l == basis.onebodybasis
         result =  manybodyoperator_1(basis, op)
@@ -303,7 +304,7 @@ function manybodyoperator_2(basis::ManyBodyBasis, op::SparseOperator)
         value = values[j]
         for m=1:N, n=1:N
             # println("row:", row, " column:"column, ind_left)
-            index = ind2sub((S, S, S, S), (column-1)*S^2 + row)
+            index = Tuple(CartesianIndices((S, S, S, S))[(column-1)*S^2 + row])
             C = coefficient(occupations[m], occupations[n], index[1:2], index[3:4])
             if C!=0.
                 result.data[m,n] += C*value

src/master.jl

@@ -8,7 +8,7 @@ using ...bases, ...states, ...operators
 using ...operators_dense, ...operators_sparse
 
 
-const DecayRates = Union{Vector{Float64}, Matrix{Float64}, Void}
+const DecayRates = Union{Vector{Float64}, Matrix{Float64}, Nothing}
 
 """
     timeevolution.master_h(tspan, rho0, H, J; <keyword arguments>)
@@ -20,7 +20,7 @@ Further information can be found at [`master`](@ref).
 function master_h(tspan, rho0::DenseOperator, H::Operator, J::Vector;
                 rates::DecayRates=nothing,
                 Jdagger::Vector=dagger.(J),
-                fout::Union{Function,Void}=nothing,
+                fout::Union{Function,Nothing}=nothing,
                 kwargs...)
     check_master(rho0, H, J, Jdagger, rates)
     tmp = copy(rho0)
@@ -43,7 +43,7 @@ function master_nh(tspan, rho0::DenseOperator, Hnh::Operator, J::Vector;
                 rates::DecayRates=nothing,
                 Hnhdagger::Operator=dagger(Hnh),
                 Jdagger::Vector=dagger.(J),
-                fout::Union{Function,Void}=nothing,
+                fout::Union{Function,Nothing}=nothing,
                 kwargs...)
     check_master(rho0, Hnh, J, Jdagger, rates)
     tmp = copy(rho0)
@@ -86,7 +86,7 @@ non-hermitian Hamiltonian and then calls master_nh which is slightly faster.
 function master(tspan, rho0::DenseOperator, H::Operator, J::Vector;
                 rates::DecayRates=nothing,
                 Jdagger::Vector=dagger.(J),
-                fout::Union{Function,Void}=nothing,
+                fout::Union{Function,Nothing}=nothing,
                 kwargs...)
     isreducible = check_master(rho0, H, J, Jdagger, rates)
     if !isreducible
@@ -130,7 +130,7 @@ at [`master_dynamic`](@ref).
 """
 function master_nh_dynamic(tspan, rho0::DenseOperator, f::Function;
                 rates::DecayRates=nothing,
-                fout::Union{Function,Void}=nothing,
+                fout::Union{Function,Nothing}=nothing,
                 kwargs...)
     tmp = copy(rho0)
     dmaster_(t, rho::DenseOperator, drho::DenseOperator) = dmaster_nh_dynamic(t, rho, f, rates, drho, tmp)
@@ -164,7 +164,7 @@ operators:
 """
 function master_dynamic(tspan, rho0::DenseOperator, f::Function;
                 rates::DecayRates=nothing,
-                fout::Union{Function,Void}=nothing,
+                fout::Union{Function,Nothing}=nothing,
                 kwargs...)
     tmp = copy(rho0)
     dmaster_(t, rho::DenseOperator, drho::DenseOperator) = dmaster_h_dynamic(t, rho, f, rates, drho, tmp)
@@ -181,13 +181,13 @@ master_nh_dynamic(tspan, psi0::Ket, f::Function; kwargs...) = master_nh_dynamic(
 
 
 # Recasting needed for the ODE solver is just providing the underlying data
-function recast!(x::Array{Complex128, 2}, rho::DenseOperator)
+function recast!(x::Array{ComplexF64, 2}, rho::DenseOperator)
     rho.data = x
 end
-recast!(rho::DenseOperator, x::Array{Complex128, 2}) = nothing
+recast!(rho::DenseOperator, x::Array{ComplexF64, 2}) = nothing
 
 function integrate_master(tspan, df::Function, rho0::DenseOperator,
-                        fout::Union{Void, Function}; kwargs...)
+                        fout::Union{Nothing, Function}; kwargs...)
     tspan_ = convert(Vector{Float64}, tspan)
     x0 = rho0.data
     state = DenseOperator(rho0.basis_l, rho0.basis_r, rho0.data)
@@ -206,7 +206,7 @@ end
 # or a matrix.
 
 function dmaster_h(rho::DenseOperator, H::Operator,
-                    rates::Void, J::Vector, Jdagger::Vector,
+                    rates::Nothing, J::Vector, Jdagger::Vector,
                     drho::DenseOperator, tmp::DenseOperator)
     operators.gemm!(-1im, H, rho, 0, drho)
     operators.gemm!(1im, rho, H, 1, drho)
@@ -257,7 +257,7 @@ function dmaster_h(rho::DenseOperator, H::Operator,
 end
 
 function dmaster_nh(rho::DenseOperator, Hnh::Operator, Hnh_dagger::Operator,
-                    rates::Void, J::Vector, Jdagger::Vector,
+                    rates::Nothing, J::Vector, Jdagger::Vector,
                     drho::DenseOperator, tmp::DenseOperator)
     operators.gemm!(-1im, Hnh, rho, 0, drho)
     operators.gemm!(1im, rho, Hnh_dagger, 1, drho)

src/mcwf.jl

@@ -6,13 +6,15 @@ using ...bases, ...states, ...operators
 using ...operators_dense, ...operators_sparse
 using ..timeevolution
 using ...operators_lazysum, ...operators_lazytensor, ...operators_lazyproduct
+using Random, LinearAlgebra
 import OrdinaryDiffEq
+
 # TODO: Remove imports
 import DiffEqCallbacks, RecursiveArrayTools.copyat_or_push!
 import ..recast!
-Base.@pure pure_inference(fout,T) = Core.Inference.return_type(fout, T)
+Base.@pure pure_inference(fout,T) = Core.Compiler.return_type(fout, T)
 
-const DecayRates = Union{Vector{Float64}, Matrix{Float64}, Void}
+const DecayRates = Union{Vector{Float64}, Matrix{Float64}, Nothing}
 
 """
     mcwf_h(tspan, rho0, Hnh, J; <keyword arguments>)
@@ -112,7 +114,7 @@ function mcwf(tspan, psi0::Ket, H::Operator, J::Vector;
         kwargs...)
     else
         Hnh = copy(H)
-        if typeof(rates) == Void
+        if typeof(rates) == Nothing
             for i=1:length(J)
                 Hnh -= 0.5im*Jdagger[i]*J[i]
             end
@@ -255,11 +257,11 @@ function integrate_mcwf(dmcwf::Function, jumpfun::Function, tspan,
                         kwargs...)
 
     tmp = copy(psi0)
-    as_ket(x::Vector{Complex128}) = Ket(psi0.basis, x)
+    as_ket(x::Vector{ComplexF64}) = Ket(psi0.basis, x)
     as_vector(psi::Ket) = psi.data
     rng = MersenneTwister(convert(UInt, seed))
     jumpnorm = Ref(rand(rng))
-    djumpnorm(x::Vector{Complex128}, t, integrator) = norm(as_ket(x))^2 - (1-jumpnorm[])
+    djumpnorm(x::Vector{ComplexF64}, t, integrator) = norm(as_ket(x))^2 - (1-jumpnorm[])
 
     if !display_beforeevent && !display_afterevent
         function dojump(integrator)
@@ -280,7 +282,7 @@ function integrate_mcwf(dmcwf::Function, jumpfun::Function, tspan,
     else
         # Temporary workaround until proper tooling for saving
         # TODO: Replace by proper call to timeevolution.integrate
-        function fout_(x::Vector{Complex128}, t::Float64, integrator)
+        function fout_(x::Vector{ComplexF64}, t::Float64, integrator)
             recast!(x, state)
             fout(t, state)
         end
@@ -322,7 +324,7 @@ function integrate_mcwf(dmcwf::Function, jumpfun::Function, tspan,
                          save_positions = (false,false))
         full_cb = OrdinaryDiffEq.CallbackSet(callback,cb,scb)
 
-        function df_(dx::Vector{Complex128}, x::Vector{Complex128}, p, t)
+        function df_(dx::Vector{ComplexF64}, x::Vector{ComplexF64}, p, t)
             recast!(x, state)
             recast!(dx, dstate)
             dmcwf(t, state, dstate)
@@ -344,7 +346,7 @@ function integrate_mcwf(dmcwf::Function, jumpfun::Function, tspan,
 end
 
 function integrate_mcwf(dmcwf::Function, jumpfun::Function, tspan,
-                        psi0::Ket, seed, fout::Void;
+                        psi0::Ket, seed, fout::Nothing;
                         kwargs...)
     function fout_(t, x)
         psi = copy(x)
@@ -366,7 +368,7 @@ Default jump function.
 * `J`: List of jump operators.
 * `psi_new`: Result of jump.
 """
-function jump(rng, t::Float64, psi::Ket, J::Vector, psi_new::Ket, rates::Void)
+function jump(rng, t::Float64, psi::Ket, J::Vector, psi_new::Ket, rates::Nothing)
     if length(J)==1
         operators.gemv!(complex(1.), J[1], psi, complex(0.), psi_new)
         psi_new.data ./= norm(psi_new)
@@ -409,7 +411,7 @@ The non-hermitian Hamiltonian is given in two parts - the hermitian part H and
 the jump operators J.
 """
 function dmcwf_h(psi::Ket, H::Operator,
-                 J::Vector, Jdagger::Vector, dpsi::Ket, tmp::Ket, rates::Void)
+                 J::Vector, Jdagger::Vector, dpsi::Ket, tmp::Ket, rates::Nothing)
     operators.gemv!(complex(0,-1.), H, psi, complex(0.), dpsi)
     for i=1:length(J)
         operators.gemv!(complex(1.), J[i], psi, complex(0.), tmp)
@@ -488,13 +490,13 @@ corresponding set of jump operators is calculated.
 """
 function diagonaljumps(rates::Matrix{Float64}, J::Vector{T}) where T <: Operator
     @assert length(J) == size(rates)[1] == size(rates)[2]
-    d, v = eig(rates)
+    d, v = eigen(rates)
     d, [sum([v[j, i]*J[j] for j=1:length(d)]) for i=1:length(d)]
 end
 
 function diagonaljumps(rates::Matrix{Float64}, J::Vector{T}) where T <: Union{LazySum, LazyTensor, LazyProduct}
     @assert length(J) == size(rates)[1] == size(rates)[2]
-    d, v = eig(rates)
+    d, v = eigen(rates)
     d, [LazySum([v[j, i]*J[j] for j=1:length(d)]...) for i=1:length(d)]
 end