Now affine transforms seem to work...

notebook/quantics.py

@@ -0,0 +1,98 @@
+import numpy as np
+import scipy.interpolate
+
+
+def decompose_int(A, n, tol=None):
+    A = np.asarray(A)
+    Ushape = A.shape[:n]
+    VTshape = A.shape[n:]
+
+    Aflat = A.reshape(np.prod(Ushape), np.prod(VTshape))
+    Uflat, VTflat = decompose_int_matrix(Aflat, tol)
+    U = Uflat.reshape(*Ushape, -1)
+    VT = VTflat.reshape(-1, *VTshape)
+    return U, VT
+
+
+def decompose_int_matrix(A, tol=None):
+    A = np.asarray(A)
+    if np.issubdtype(A.dtype, np.integer):
+        A = A.astype(np.float32)
+    if tol is None:
+        tol = max(2, np.sqrt(A.shape[0])) * np.finfo(A.dtype).eps
+
+    # Truncated SVD
+    U, s, VT = np.linalg.svd(A, full_matrices=False)
+    s_mask = (s > tol * s[0])
+    U = U[:, s_mask]
+    s = s[s_mask]
+    VT = VT[s_mask, :]
+
+    # Truncate small elements
+    U_mask = np.abs(U) > tol
+    VT_mask = np.abs(VT) > tol
+    U[~U_mask] = 0
+    VT[~VT_mask] = 0
+
+    # Undo scaling
+    U_count = U_mask.sum(0)
+    VT_count = VT_mask.sum(1)
+    U *= np.sqrt(U_count)
+    VT *= np.sqrt(VT_count)[:, None]
+    s /= np.sqrt(U_count * VT_count)
+    np.testing.assert_allclose(s, 1, atol=tol, rtol=0)
+
+    # Round stuff
+    U_int = np.round(U).astype(int)
+    VT_int = np.round(VT).astype(int)
+    np.testing.assert_allclose(U_int, U, atol=tol, rtol=0)
+    np.testing.assert_allclose(VT_int, VT, atol=tol, rtol=0)
+
+    # Regauge values
+    U_sgn = np.sign(U_int.sum(0))
+    np.testing.assert_equal(np.abs(U_sgn), 1)
+    U_int *= U_sgn
+    VT_int *= U_sgn[:, None]
+
+    # We want a bit matrix
+    np.testing.assert_equal(U_int, U_int.astype(bool))
+    np.testing.assert_equal(VT_int, VT_int.astype(bool))
+
+    # Return
+    return U_int, VT_int
+
+
+def interleave_bits(A, K):
+    A = np.asarray(A)
+    R = A.ndim // K
+    if A.shape != (2,) * (K*R):
+        raise ValueError("not of proper quantics shape")
+
+    order = np.hstack([np.arange(i, K*R, R) for i in range(R)])
+    return A.transpose(*order)
+
+
+def quantize(A, interleave=False):
+    A = np.asarray(A)
+    bits = np.log2(A.shape)
+    if not (bits == bits.astype(int)).all():
+        raise ValueError("invalid shape: " + str(A.shape))
+    bits = bits.astype(int)
+    A = A.reshape((2,) * bits.sum())
+
+    if interleave:
+        if not (bits == bits[0]).all():
+            raise ValueError("unable to interleave")
+        A = interleave_bits(A, len(bits))
+    return A
+
+
+def print_nonzeros(A, headers=None):
+    A = np.asarray(A)
+    N = A.ndim
+    if headers:
+        headers = tuple(headers)
+        print(("%2s " * N) % headers)
+    fmt = " ".join(("%2i",) * N)
+    for x in zip(*A.nonzero()):
+        print(fmt % x)

src/affine.jl

@@ -13,7 +13,7 @@ function affine_transform_mpo(
         throw(ArgumentError("vector is not correctly dimensioned"))
 
     # get the tensors so that we know how large the links must be
-    tensors = affine_transform_matrices(
+    tensors = affine_transform_tensors(
                     R, SMatrix{M, N, Int}(A), SVector{M, Int}(b))
 
     # Create the links
@@ -35,115 +35,75 @@ function affine_transform_mpo(
     return mpo
 end
 
-function affine_transform_matrices(
+function affine_transform_tensors(
             R::Int, A::SMatrix{M, N, Int}, b::SVector{M, Int}
             ) where {M, N}
     # Checks
-    0 <= R <= 62 ||
-        throw(ArgumentError("invalid value of the length R"))
-
-    # Matrix which maps out which bits to add together. We use 2's complement,
-    # -x == ~x + 1, which means we use two masks: A1 operating on set bits, and
-    # A0 operating on unset bits. We group the +1 together with the shift.
-    A1 = @. clamp(A, 0, typemax(Int))
-    A0 = @. clamp(-A, 0, typemax(Int))
-    b = b .+ vec(sum(Bool.(A0), dims=2))
-    #println(A1, A0, b)
-
-    # Amax is the maximum value that can be reached by multiplying it with set
-    # or unset bits.
-    Amax = vec(sum(A0 + A1, dims=2))
+    0 <= R ||
+        throw(DomainError(R, "invalid value of the length R"))
 
     # The output tensors are a collection of matrices, but their first two
     # dimensions (links) vary
-    tensors = Array{Bool, 4}[]
-    sizehint!(tensors, R)
+    tensors = Vector{Array{Bool, 4}}(undef, R)
 
-    # We have to carry all but the least siginificant bit to the adjacent
-    # tensor. α is the index of the "outgoing" carrys, which at the beginning
-    # is simply all zeros.
-    maxcarry_α = zeros(typeof(b))
-    ncarry_α = prod(maxcarry_α .+ 1)
-    base_α = _indexbase_from_size(maxcarry_α, 1)
+    # The initial carry is zero
+    carry = [zero(SVector{M, Int})]
     for r in 1:R
-        # Copy the outgoing carry α of the previous tensor to the "incoming"
-        # carry β of the current tensor.
-        maxcarry_β = maxcarry_α
-        ncarry_β = ncarry_α
-
         # Figure out the current bit to add from the shift term and shift
         # it out from the array
-        bcurr = @. b & 1
-        b = @. b >> 1
+        bcurr = @. Bool(b & 1)
 
-        # Determine the maximum outgoing carry. It is determined by the maximum
-        # previous carry plus the maximum value of the mapped indices plus the
-        # current shift, all divided by two. In the case of the last tensor, we
-        # discard all carrys.
+        # Get tensor. For the last tensor, we ignore the carry and just
+        # or together all possible combinations.
+        new_carry, new_tensor = affine_transform_tensor(A, bcurr, carry)
         if r == R
-            maxcarry_α = zeros(typeof(maxcarry_α))
-            ncarry_α = 1
-            base_α = _indexbase_from_size(maxcarry_α, 0)
+            all_values = reduce(.|, eachslice(new_tensor, dims=1))
+            tensors[r] = reshape(all_values, 1, size(all_values)...)
         else
-            maxcarry_α =  @. (maxcarry_β + Amax + bcurr) >> 1
-            ncarry_α = prod(maxcarry_α .+ 1)
-            base_α = _indexbase_from_size(maxcarry_α, 1)
+            tensors[r] = new_tensor
         end
 
-        values = zeros(Bool, ncarry_α, ncarry_β, 1 << M, 1 << N)
-        #println("\n r=$r $(bcurr)")
-        # Fill values array. The idea is the following: we iterate over all
-        # possible values of the carry from the previous tensor (c_β). Each of
-        # those corresponds to a bond index β.
-        allcarry_β = map(i -> 0:i, maxcarry_β)
-        for (β, _c_β) in enumerate(Iterators.product(allcarry_β...))
-            c_β = SVector{M}(_c_β)
-
-            # Now we iterate over all possible indices of the input legs, which
-            # are all combination of values of N input bits
-            allspin_j = ntuple(_ -> 0:1, N)
-            for (j, _σ_j) in enumerate(Iterators.product(allspin_j...))
-                σ_j = SVector{N, Bool}(_σ_j)
-
-                # We apply the transformation y = Ax + b to the current bits
-                # σ_j and adding the incoming carry c_β. The least significant
-                # bits are then the output σ_i while the higher-order bits
-                # form the outgoing carry c_α
-                ifull = A1 * σ_j + A0 * (.~σ_j) + bcurr + c_β
-                c_α = @. ifull >> 1
-                σ_i = SVector{M, Bool}(@. ifull & 1)
-
-                # Map the outgoing carry to an index and store
-                α = mapreduce(*, +, c_α, base_α, init=1)
-                i = _spins_to_number(σ_i) + 1
-                #println("$(c_α) * 2 + $(σ_i) <- $(c_β) * 2 + $(σ_j)")
-                values[α, β, i, j] = true
-            end
-        end
-        push!(tensors, values)
+        # Set carry to the next value
+        carry = new_carry
+        b = @. b >> 1
     end
     return tensors
 end
 
-"""
-    get_int_inverse(A::AbstractMatrix{<:Integer})
+function affine_transform_tensor(
+            A::SMatrix{M, N, Int}, b::SVector{M, Bool},
+            carry::AbstractVector{SVector{M, Int}}) where {M, N}
+    # The basic idea here is the following: we compute r = A*x + b + c for all
+    # "incoming" carrys d and all possible bit vectors, x ∈ {0,1}^N.  Then we
+    # decompose r = 2*c + y, where y is again a bit vector, y ∈ {0,1}^M, and
+    # c is the "outgoing" carry, which may be negative.  We then store this
+    # as something like out[d][c, x, y] = true.
+    out = Dict{SVector{M, Int}, Array{Bool, 3}}()
+    for (c_index, c) in enumerate(carry)
+        for (x_index, x) in enumerate(Iterators.product(ntuple(_ -> 0:1, N)...))
+            r = A * SVector{N, Bool}(x) + b + SVector{N, Int}(c)
+            y = Bool.(r .& 1)
+            d = r .>> 1
+            y_index = _spins_to_number(y) + 1
+
+            d_mat = get!(out, d) do
+                return zeros(Bool, length(carry), 1 << M, 1 << N)
+            end
+            d_mat[c_index, x_index, y_index] = true
+        end
+    end
 
-Return inverse matrix to integer A, ensuring that it is indeed integer.
-"""
-function get_int_invserse(A::AbstractMatrix{<:Integer})
-    Ainv = inv(A)
-    Ainv_int = map(eltype(A) ∘ round, Ainv)
-    isapprox(Ainv, Ainv_int, atol=size(A,1)*eps()) ||
-        throw(DomainError(Ainv, "inverse is not integer"))
-    return typeof(A)(Ainv)
+    # We translate the dictionary into a vector of carrys (which we can then
+    # pass into the next iteration) and a 4-way tensor of output values.
+    carry_out = Vector{SVector{M, Int}}(undef, length(out))
+    value_out = Array{Bool, 4}(undef, length(out), length(carry), 1<<M, 1<<N)
+    for (p_index, p) in enumerate(pairs(out))
+        carry_out[p_index] = p.first
+        value_out[p_index, :, :, :] = p.second
+    end
+    return carry_out, value_out
 end
 
-_indexbase_from_size(v::SVector{M,T}, init::T) where {M,T} =
-    SVector{M,T}(_indexbase_from_size(Tuple(v), init))
-_indexbase_from_size(v::Tuple{}, init::T) where {T} = ()
-_indexbase_from_size(v::Tuple{T, Vararg{T}}, init::T) where {T} =
-    init, _indexbase_from_size(v[2:end], init * v[1])...
-
 _spins_to_number(v::SVector{<:Any, Bool}) = _spins_to_number(Tuple(v))
 _spins_to_number(v::Tuple{Bool, Vararg{Bool}}) =
     _spins_to_number(v[2:end]) << 1 | v[1]