Merge pull request #44 from RexWzh/rex-dev

Gauss elimination for linear equation

LICENSE

@@ -1,4 +1,4 @@
-Copyright (c) 2019 Jérémie Gillet <[email protected]>
+Copyright (c) 2019 Jérémie Gillet <[email protected]> and contributors
 
 Permission is hereby granted, free of charge, to any person obtaining a copy
 of this software and associated documentation files (the "Software"), to deal

Project.toml

@@ -20,9 +20,11 @@ julia = "1.3"
 
 [extras]
 ImageTransformations = "02fcd773-0e25-5acc-982a-7f6622650795"
+QRDecoders = "d4999880-6331-4276-8b7d-7ee1f305cff8"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TestImages = "5e47fb64-e119-507b-a336-dd2b206d9990"
 
 [targets]
-test = ["ImageTransformations", "Test", "Random", "TestImages"]
+test = ["ImageTransformations", "QRDecoders", "Test", "Random", "StatsBase", "TestImages"]

src/QRCoders.jl

@@ -49,4 +49,7 @@ include("styles/style.jl")
 # Generate and export QR code
 include("export.jl")
 
+# Tools for equation solving
+include("styles/equation.jl")
+
 end # module

src/polynomial.jl

@@ -34,7 +34,7 @@ In our notations, the generator matrix is rotate by 180 degrees.
 
 module Polynomial
 
-export Poly, geterrcode, mult, divide, generator_matrix
+export Poly, geterrcode, mult, divide, generator_matrix, gfinv, encodepoly
 
 import Base: length, iterate, ==, <<, +, *, ÷, %, copy, zero, eltype
 
@@ -311,6 +311,16 @@ Create the Generator Polynomial of degree `n`.
 generator(n::Int) = generator(UInt8, n)
 generator(::Type{T}, n::Int) where T = prod([Poly{T}([powtable[i], one(T)]) for i in 1:n])
 
+"""
+    encodepoly(msgpoly::Poly{T}, n::Int) where T
+
+Encode the message polynomial to received polynomial.
+"""
+function encodepoly(msgpoly::Poly{T}, n::Int) where T
+    f = msgpoly << n
+    f + f % generator(T, n)
+end
+
 """
     geterrcode(f::Poly, n::Int)
 
@@ -368,7 +378,8 @@ e.g. a_0, ..., a_n.
 In this sense, we still have G⋅x = c, where x is the message polynomial 
 and c is the received polynomial.
 
-To get an ordinary generator matrix, just rotate it by 180 degrees.
+To get an ordinary generator matrix, just rotate it by 180 degrees,
+i.e. @view(G[end:-1:1, end:-1:1])
 """
 generator_matrix(msglen::Int, necwords::Int) = generator_matrix(UInt8, msglen, necwords)
 function generator_matrix(::Type{T}, msglen::Int, necwords::Int) where T

src/styles/equation.jl

@@ -0,0 +1,77 @@
+# Tools for linear equation
+
+"""
+    gauss_elimination(A::AbstractMatrix, b::AbstractVector)
+
+Solve the linear equations `Ax=b` by Gauss elimination.
+"""
+gauss_elimination(A::AbstractMatrix, b::AbstractVector) = gauss_elimination!(copy(A), copy(b))
+function gauss_elimination!(A::AbstractMatrix, b::AbstractVector)
+    x = gauss_elimination!(A, reshape(b, :, 1))
+    reshape(x, :) # convert to vector
+end
+
+"""
+    gauss_elimination(A::AbstractMatrix, B::AbstractMatrix)
+
+Solve the linear equations `Ax=B` by Gauss elimination.
+In particular, `Ax=I` returns the inverse of `A`.
+"""
+gauss_elimination(A::AbstractMatrix, B::AbstractMatrix) = gauss_elimination!(copy(A), copy(B))
+function gauss_elimination!(A::AbstractMatrix, B::AbstractMatrix)
+    # check the size of A and B
+    m, n = size(A)
+    m == n || throw(DimensionMismatch("A must be square"))
+    size(B, 1) == m || throw(DimensionMismatch("A and B have different size"))
+    # Gauss elimination
+    @inbounds for i in 1:n
+        # Find the pivot
+        pivot = findfirst(!iszero, @view(A[i:end, i]))
+        pivot === nothing && throw(ArgumentError("A is singular"))
+        pivot += i - 1
+        # Swap the pivot row with the current row
+        currow = @views divide.(A[pivot, :], A[pivot, i]) # set leading entry to 1
+        curval = @views divide.(B[pivot, :], A[pivot, i])
+        A[pivot, :] = @view(A[i, :]) # `=` is faster than `.=` in general
+        B[pivot, :] = @view(B[i, :])
+        A[i, :], B[i, :] = currow, curval
+        # Eliminate the current column
+        for j in 1:n
+            j == i && continue
+            # rowⱼ = rowⱼ - Aⱼᵢ * rowᵢ
+            if !iszero(A[j, i])    
+                B[j, :] .⊻= mult.(A[j, i], curval)
+                A[j, :] .⊻= mult.(A[j, i], currow)
+            end
+        end
+    end
+    return B
+end
+
+"""
+    gfinv(A::AbstractMatrix)
+
+Return the inverse of `A` in GF(256).
+"""
+Polynomial.gfinv(A::AbstractMatrix) = gauss_elimination(A, one(A))
+
+"""
+    fillblank( block::AbstractVector{<:Integer}
+             , validinds::AbstractVector{<:Integer}
+             , necwords::Int)
+
+Use Gauss elimination to fill the blank entries in `block`
+with valid indices `misinds`.
+"""
+function fillblank( block::AbstractVector{<:Integer}
+                  , validinds::AbstractVector{<:Integer}
+                  , necwords::Int)
+    reclen = length(block)
+    length(validinds) + necwords == reclen || throw(DimensionMismatch(
+        "The number of missing indices must be equal to `necwords`"))
+    msglen = reclen - necwords
+    G = generator_matrix(eltype(block), msglen, necwords)
+    G = @view(G[end:-1:1, end:-1:1]) # rotate by 180 degree
+    msg = @views gauss_elimination(G[validinds, :], block[validinds])
+    return mult(G, msg) # encode message
+end

src/styles/locate.jl

@@ -41,6 +41,12 @@ function getindexes(v::Int)
 end
 
 ## 2.2 split indexes into several segments(de-interleave)
+"""
+    getecinfo(v::Int, eclevel::ErrCorrLevel)
+
+Get the error correction information.
+"""
+getecinfo(v::Int, eclevel::ErrCorrLevel) = @views ecblockinfo[eclevel][v, :]
 
 """
     getsegments(v::Int, mode::Mode, eclevel::ErrCorrLevel)
@@ -53,7 +59,7 @@ The procedure is similar to `deinterleave` in `QRDecoders.jl`.
 function getsegments(v::Int, eclevel::ErrCorrLevel)
     # initialize
     ## get information about error correction
-    necwords, nb1, nc1, nb2, nc2 = ecblockinfo[eclevel][v, :]
+    necwords, nb1, nc1, nb2, nc2 = getecinfo(v, eclevel)
     ## initialize blocks
     expand(x) = (8 * x - 7):8 * x
     segments = vcat([Vector{Int}(undef, 8 * nc1) for _ in 1:nb1],

test/runtests.jl

@@ -5,6 +5,8 @@ using Random
 using ImageCore
 using ImageTransformations
 using TestImages
+using StatsBase
+using QRDecoders.Syndrome: fillerasures!
 
 using QRCoders:
     # build
@@ -24,7 +26,8 @@ using QRCoders:
     # style
     unicodeplot, getindexes, getsegments,
     imagebyerrcor, animatebyerrcor,
-    pickcodewords
+    pickcodewords, getecinfo,
+    gauss_elimination, fillblank
 
 using QRCoders.Polynomial:
     # operator for GF(256) integers
@@ -33,7 +36,8 @@ using QRCoders.Polynomial:
     # operator for polynomials
     iszeropoly, degree, zero, unit,
     rpadzeros, rstripzeros, generator, 
-    geterrcode, euclidean_divide
+    geterrcode, euclidean_divide,
+    encodepoly
 
 # random polynomial
 randpoly(::Type{T}, n::Int) where T = Poly{T}([rand(0:255, n-1)..., rand(1:255)])
@@ -62,5 +66,8 @@ include("tst_overall.jl")
 # style
 include("tst_style.jl")
 
+# equations
+include("tst_equation.jl")
+
 # final message
 unicodeplotbychar("https://github.com/JuliaImages/QRCoders.jl") |> println

test/tst_equation.jl

@@ -0,0 +1,93 @@
+# Test for Gauss elimination
+function fillblankbyFA( block::AbstractVector{<:Integer}
+                      , misinds::AbstractVector{<:Integer}
+                      , necwords::Int)
+    n = length(block) # number of bytes
+    recpoly = Poly(reverse(block))
+    orgpoly = fillerasures!(recpoly, n .- misinds, necwords)
+    return @view(orgpoly.coeff[end:-1:1])
+end
+
+@testset "Fill blank -- Gauss elimination vs Forney algorithm" begin
+    # test for fillblank
+    msglen = rand(1:100)
+    necwords = rand(1:255 - msglen)
+    reclen = msglen + necwords
+    block = rand(0:255, reclen)
+    validinds = sample(1:msglen, msglen, replace=false)
+    misinds = setdiff(1:reclen, validinds)
+    # Note that Gauss elimination takes O(N^3)
+    # while Forney Algorithm takes O(N^2)
+    validblock = fillblank(block, validinds, necwords)
+    blockbyFA = fillblankbyFA(block, misinds, necwords)
+    @test validblock == blockbyFA
+    @test validblock[validinds] == block[validinds]
+    # test for UInt8
+    msglen = rand(1:100)
+    necwords = rand(1:255 - msglen)
+    reclen = msglen + necwords
+    block = rand(UInt8, reclen)
+    validinds = sample(1:msglen, msglen, replace=false)
+    misinds = setdiff(1:reclen, validinds)
+    validblock = fillblank(block, validinds, necwords)
+    blockbyFA = fillblankbyFA(block, misinds, necwords)
+    @test validblock == blockbyFA
+    @test validblock[validinds] == block[validinds]
+end
+
+## use Gauss elimination to find the inverse of the matrix
+@testset "Generator matrix -- inverse of submatrix" begin
+    # Test for matrix inverse
+    ## Int
+    A = [1 2 3; 4 5 6; 7 8 9]
+    A1 = gfinv(A)
+    @test mult(A, A1) |> isone
+    A = [1 2 3; 4 5 6; 5 7 9] # note that this matrix is invertible in GF(256)
+    A1 = gfinv(A)
+    @test mult(A, A1) |> isone
+    ## UInt8
+    A = UInt8[1 2 3; 4 5 6; 7 8 9]
+    A1 = gfinv(A)
+    @test mult(A, A1) |> isone
+    ## test throw
+    @test_throws ArgumentError gfinv([1 2; 1 2])
+    A = [1 2 3; 4 5 6; 7 8 9]
+    A[3, :] = A[1, :] .⊻ A[2, :]
+    @test_throws ArgumentError gfinv(A)
+    ## any `msglen` rows of a generator matrix are linearly dependent
+    msglen = 3
+    necwords = 5
+    reclen = msglen + necwords
+    A = generator_matrix(msglen, necwords)
+    rows = sample(1:reclen, msglen, replace=false)
+    @test mult(gfinv(A[rows, :]), A[rows, :]) |> isone
+
+    necwords, _, msglen = getecinfo(1, Medium())
+    A = generator_matrix(msglen, necwords)
+    rows = sample(1:msglen + necwords, msglen, replace=false)
+    @test mult(gfinv(A[rows, :]), A[rows, :]) |> isone
+    ### random test
+    for eclevel in eclevels, v in 1:40
+        necwords, _, msglen, _, msglen2 = getecinfo(v, eclevel)
+        A = generator_matrix(msglen, necwords)
+        rows = sample(1:msglen + necwords, msglen, replace=false)
+        @test mult(gfinv(A[rows, :]), A[rows, :]) |> isone
+        A = generator_matrix(msglen2, necwords)
+        rows = sample(1:msglen2 + necwords, msglen2, replace=false)
+        @test mult(gfinv(A[rows, :]), A[rows, :]) |> isone
+    end
+
+end
+
+@testset "Linear equations" begin
+    A = [1 2 3; 4 5 6; 7 8 9]
+    b = [1, 2, 3]
+    x = gauss_elimination(A, b)
+    @test mult(A, x) == b
+    b = reshape([1, 2, 3], :, 1)
+    x = gauss_elimination(A, b)
+    @test mult(A, x) == b
+    B = rand(0:255, 3, rand(1:10))
+    x = gauss_elimination(A, B)
+    @test mult(A, x) == B
+end

test/tst_operation.jl

@@ -177,17 +177,18 @@ end
     @test mult(A, mult(B, C)) == mult(mult(A, B), C)
     @test mult(A, mult(B, b)) == mult(mult(A, B), b)
 
-    encode(f::Poly, n::Int) = (f << n) + (f << n) % generator(n)
     G = generator_matrix(3, 4)
     f = randpoly(3)
-    @test mult(G, f.coeff) == encode(f, 4).coeff
+    @test mult(G, f.coeff) == encodepoly(f, 4).coeff
 
     m = rand(1:254)
     n = rand(1:255-m)
     G = generator_matrix(m, n)
     @test size(G) == (m+n, m)
     f = randpoly(m)
-    @test mult(G, f.coeff) == encode(f, n).coeff
+    @test mult(G, f.coeff) == encodepoly(f, n).coeff
+    # @btime mult(G, f.coeff); # 8.311 μs (86 allocations: 10.66 KiB)
+    # @btime encodepoly(f, n).coeff; # 1.401 μs (25 allocations: 1.88 KiB)
 
     @test powx(Int, 2) == Poly([0, 0, 1])
     @test powx(0) == Poly{UInt8}([1])