Merge pull request #170 from Jutho/jh/randisometry

restore randisometry, add tests and deprecations

Author: Jutho

src/TensorKit.jl

@@ -56,7 +56,7 @@ export spacetype, sectortype, storagetype, scalartype, tensormaptype
 export blocksectors, blockdim, block, blocks
 
 # random methods for constructor
-export randuniform, randnormal, randisometry, randhaar
+export randisometry, randisometry!, rand, rand!, randn, randn!
 
 # special purpose constructors
 export zero, one, one!, id, isomorphism, unitary, isometry
@@ -123,7 +123,7 @@ using SparseArrays: SparseMatrixCSC, sparse, nzrange, rowvals, nonzeros
 
 import Base.Meta
 
-using Random: Random
+using Random: Random, rand!, randn!
 
 using PackageExtensionCompat
 

src/auxiliary/deprecate.jl

@@ -36,6 +36,24 @@ for f in (:rand, :randn, :zeros, :ones)
     end
 end
 
+Base.@deprecate(randuniform(dims::Base.Dims), rand(dims))
+Base.@deprecate(randuniform(T::Type{<:Number}, dims::Base.Dims), rand(T, dims))
+Base.@deprecate(randnormal(dims::Base.Dims), randn(dims))
+Base.@deprecate(randnormal(T::Type{<:Number}, dims::Base.Dims), randn(T, dims))
+Base.@deprecate(randhaar(dims::Base.Dims), randisometry(dims))
+Base.@deprecate(randhaar(T::Type{<:Number}, dims::Base.Dims), randisometry(T, dims))
+
+for (f1, f2) in ((:randuniform, :rand), (:randnormal, :randn), (:randisometry, :randisometry), (:randhaar, :randisometry))
+    @eval begin
+        Base.@deprecate TensorMap(::typeof($f1), T::Type, P::HomSpace) $f2(T, P)
+        Base.@deprecate TensorMap(::typeof($f1), P::HomSpace) $f2(P)
+        Base.@deprecate TensorMap(::typeof($f1), T::Type, cod::TensorSpace, dom::TensorSpace) $f2(T, P, cod, dom)
+        Base.@deprecate TensorMap(::typeof($f1), cod::TensorSpace, dom::TensorSpace) $f2(cod, dom)
+        Base.@deprecate Tensor(::typeof($f1), T::Type, space::TensorSpace) $f2(T, space)
+        Base.@deprecate Tensor(::typeof($f1), space::TensorSpace) $f2(space)
+    end
+end
+
 Base.@deprecate EuclideanProduct() EuclideanInnerProduct()
 
 #! format: on

src/auxiliary/random.jl

@@ -1,24 +1,3 @@
-"""
-    randuniform([::Type{T}=Float64], dims::Dims{N}) -> Array{T,N}
-
-Create an array of size `dims` with random entries uniformly distributed in the allowed
-values of `T`.
-
-See also [`randnormal`](@ref), [`randisometry`](@ref) and[`randhaar`](@ref).
-"""
-randuniform(dims::Base.Dims) = randuniform(Float64, dims)
-randuniform(::Type{T}, dims::Base.Dims) where {T<:Number} = rand(T, dims)
-
-"""
-    randnormal([::Type{T}=Float64], dims::Dims{N}) -> Array{T,N}
-
-Create an array of size `dims` with random entries distributed according to the standard normal distribution.
-
-See also [`randuniform`](@ref), [`randisometry`](@ref) and[`randhaar`](@ref).
-"""
-randnormal(dims::Base.Dims) = randnormal(Float64, dims)
-randnormal(::Type{T}, dims::Base.Dims) where {T<:Number} = randn(T, dims)
-
 """
     randisometry([::Type{T}=Float64], dims::Dims{2}) -> Array{T,2}
     randhaar([::Type{T}=Float64], dims::Dims{2}) -> Array{T,2}
@@ -29,9 +8,18 @@ See also [`randuniform`](@ref) and [`randnormal`](@ref).
 """
 randisometry(dims::Base.Dims{2}) = randisometry(Float64, dims)
 function randisometry(::Type{T}, dims::Base.Dims{2}) where {T<:Number}
-    return dims[1] >= dims[2] ?
-           MatrixAlgebra.leftorth!(randnormal(T, dims), QRpos(), 0)[1] :
-           throw(DimensionMismatch("cannot create isometric matrix with dimensions $dims; isometry needs to be tall or square"))
+    return randisometry(Random.default_rng(), T, dims)
+end
+function randisometry(rng::Random.AbstractRNG, ::Type{T},
+                      dims::Base.Dims{2}) where {T<:Number}
+    return randisometry!(rng, Matrix{T}(undef, dims))
 end
 
-const randhaar = randisometry
+randisometry!(A::AbstractMatrix) = randisometry!(Random.default_rng(), A)
+function randisometry!(rng::Random.AbstractRNG, A::AbstractMatrix)
+    dims = size(A)
+    dims[1] >= dims[2] ||
+        throw(DimensionMismatch("cannot create isometric matrix with dimensions $dims; isometry needs to be tall or square"))
+    Q, = MatrixAlgebra.leftorth!(Random.randn!(rng, A), QRpos(), 0)
+    return copy!(A, Q)
+end

src/tensors/tensor.jl

@@ -179,73 +179,84 @@ for (fname, felt) in ((:zeros, :zero), (:ones, :one))
     end
 end
 
-for randfun in (:rand, :randn, :randexp)
-    randfun! = Symbol(randfun, :!)
+for randf in (:rand, :randn, :randexp, :randisometry)
     _docstr = """
-        $randfun([rng=default_rng()], [T=Float64], codomain::ProductSpace{S,N₁},
+        $randf([rng=default_rng()], [T=Float64], codomain::ProductSpace{S,N₁},
                      domain::ProductSpace{S,N₂}) where {S,N₁,N₂,T} -> t
-        $randfun([rng=default_rng()], [T=Float64], codomain ← domain) -> t
+        $randf([rng=default_rng()], [T=Float64], codomain ← domain) -> t
         
-        Generate a tensor `t` with entries generated by `$randfun`.
+    Generate a tensor `t` with entries generated by `$randf`.
+    
+    See also [`($randf)!`](@ref).
     """
     _docstr! = """
-        $randfun!([rng=default_rng()], t::AbstractTensorMap) -> t
+        $(randf)!([rng=default_rng()], t::AbstractTensorMap) -> t
         
-    Fill the tensor `t` with entries generated by `$randfun!`.
+    Fill the tensor `t` with entries generated by `$(randf)!`.
+    
+    See also [`($randf)`](@ref).
     """
 
+    if randf != :randisometry
+        randfun = GlobalRef(Random, randf)
+        randfun! = GlobalRef(Random, Symbol(randf, :!))
+    else
+        randfun = randf
+        randfun! = Symbol(randf, :!)
+    end
+
     @eval begin
-        @doc $_docstr Random.$randfun(::Type, ::HomSpace)
-        @doc $_docstr! Random.$randfun!(::Type, ::HomSpace)
+        @doc $_docstr $randfun(::Type, ::HomSpace)
+        @doc $_docstr! $randfun!(::Type, ::HomSpace)
 
         # converting `codomain` and `domain` into `HomSpace`
-        function Random.$randfun(codomain::TensorSpace{S},
-                                 domain::TensorSpace{S}) where {S<:IndexSpace}
-            return Random.$randfun(codomain ← domain)
+        function $randfun(codomain::TensorSpace{S},
+                          domain::TensorSpace{S}) where {S<:IndexSpace}
+            return $randfun(codomain ← domain)
         end
-        function Random.$randfun(::Type{T}, codomain::TensorSpace{S},
-                                 domain::TensorSpace{S}) where {T,S<:IndexSpace}
-            return Random.$randfun(T, codomain ← domain)
+        function $randfun(::Type{T}, codomain::TensorSpace{S},
+                          domain::TensorSpace{S}) where {T,S<:IndexSpace}
+            return $randfun(T, codomain ← domain)
         end
-        function Random.$randfun(rng::Random.AbstractRNG, ::Type{T},
-                                 codomain::TensorSpace{S},
-                                 domain::TensorSpace{S}) where {T,S<:IndexSpace}
-            return Random.$randfun(rng, T, codomain ← domain)
+        function $randfun(rng::Random.AbstractRNG, ::Type{T},
+                          codomain::TensorSpace{S},
+                          domain::TensorSpace{S}) where {T,S<:IndexSpace}
+            return $randfun(rng, T, codomain ← domain)
         end
 
         # accepting single `TensorSpace`
-        Random.$randfun(codomain::TensorSpace) = Random.$randfun(codomain ← one(codomain))
-        function Random.$randfun(::Type{T}, codomain::TensorSpace) where {T}
-            return Random.$randfun(T, codomain ← one(codomain))
+        $randfun(codomain::TensorSpace) = $randfun(codomain ← one(codomain))
+        function $randfun(::Type{T}, codomain::TensorSpace) where {T}
+            return $randfun(T, codomain ← one(codomain))
         end
-        function Random.$randfun(rng::Random.AbstractRNG, ::Type{T},
-                                 codomain::TensorSpace) where {T}
-            return Random.$randfun(rng, T, codomain ← one(domain))
+        function $randfun(rng::Random.AbstractRNG, ::Type{T},
+                          codomain::TensorSpace) where {T}
+            return $randfun(rng, T, codomain ← one(domain))
         end
 
         # filling in default eltype
-        Random.$randfun(V::TensorMapSpace) = Random.$randfun(Float64, V)
-        function Random.$randfun(rng::Random.AbstractRNG, V::TensorMapSpace)
-            return Random.$randfun(rng, Float64, V)
+        $randfun(V::TensorMapSpace) = $randfun(Float64, V)
+        function $randfun(rng::Random.AbstractRNG, V::TensorMapSpace)
+            return $randfun(rng, Float64, V)
         end
 
         # filling in default rng
-        function Random.$randfun(::Type{T}, V::TensorMapSpace) where {T}
-            return Random.$randfun(Random.default_rng(), T, V)
+        function $randfun(::Type{T}, V::TensorMapSpace) where {T}
+            return $randfun(Random.default_rng(), T, V)
         end
-        Random.$randfun!(t::AbstractTensorMap) = Random.$randfun!(Random.default_rng(), t)
+        $randfun!(t::AbstractTensorMap) = $randfun!(Random.default_rng(), t)
 
         # implementation
-        function Random.$randfun(rng::Random.AbstractRNG, ::Type{T},
-                                 V::TensorMapSpace) where {T}
+        function $randfun(rng::Random.AbstractRNG, ::Type{T},
+                          V::TensorMapSpace) where {T}
             t = TensorMap{T}(undef, V)
-            Random.$randfun!(rng, t)
+            $randfun!(rng, t)
             return t
         end
 
-        function Random.$randfun!(rng::Random.AbstractRNG, t::AbstractTensorMap)
+        function $randfun!(rng::Random.AbstractRNG, t::AbstractTensorMap)
             for (_, b) in blocks(t)
-                Random.$randfun!(rng, b)
+                $randfun!(rng, b)
             end
             return t
         end

test/tensors.jl

@@ -308,12 +308,25 @@ for V in spacelist
                 @test HrA12array ≈ convert(Array, HrA12)
             end
         end
+        @timedtestset "Multiplication of isometries: test properties" begin
+            W2 = V4 ⊗ V5
+            W1 = W2 ⊗ (oneunit(V1) ⊕ oneunit(V1))
+            for T in (Float64, ComplexF64)
+                t1 = randisometry(T, W1, W2)
+                t2 = randisometry(T, W2 ← W2)
+                @test t1' * t1 ≈ one(t2)
+                @test t2' * t2 ≈ one(t2)
+                @test t2 * t2' ≈ one(t2)
+                P = t1 * t1'
+                @test P * P ≈ P
+            end
+        end
         @timedtestset "Multiplication and inverse: test compatibility" begin
             W1 = V1 ⊗ V2 ⊗ V3
             W2 = V4 ⊗ V5
             for T in (Float64, ComplexF64)
                 t1 = rand(T, W1, W1)
-                t2 = rand(T, W2, W2)
+                t2 = rand(T, W2 ← W2)
                 t = rand(T, W1, W2)
                 @test t1 * (t1 \ t) ≈ t
                 @test (t / t2) * t2 ≈ t
@@ -331,9 +344,9 @@ for V in spacelist
                 W1 = V1 ⊗ V2 ⊗ V3
                 W2 = V4 ⊗ V5
                 for T in (Float32, Float64, ComplexF32, ComplexF64)
-                    t1 = rand(T, W1, W1)
+                    t1 = rand(T, W1 ← W1)
                     t2 = rand(T, W2, W2)
-                    t = rand(T, W1, W2)
+                    t = rand(T, W1 ← W2)
                     d1 = dim(W1)
                     d2 = dim(W2)
                     At1 = reshape(convert(Array, t1), d1, d1)