generalize hadamard(n) with a lookup table of known Hadamard matrices (fix #2)

stevengj · stevengj · commit 17a076ce6f69 · 2014-07-07T14:46:46.000-04:00
diff --git a/README.md b/README.md
@@ -48,10 +48,22 @@ and sums them to recover the signal, with no `n` factor.
   functions take the same arguments as `fwht` and `ifwht` and have the
   same normalizations, respectively.
 
-* Also provided is a function `hadamard(n)` which returns the
-  (natural-order) Hadamard matrix of order `n`, similar to the Matlab
-  function of the same name.  Currently, `n` must be a power of two,
-  in which case this function is equivalent to `ifwht_natural(eye(n), 1)`.
+## Hadamard matrices
+
+We also provide a a function `hadamard(n)` which returns a Hadamard
+matrix of order `n`, similar to the Matlab function of the same name.
+The known Hadamard matrices up to size 256 are currently supported
+(via a lookup table), along with products of these sizes and powers of
+two.
+
+The return value of `hadamard(n)` is a matrix of `Int8` values.  If
+you are planning to do matrix computations with this matrix, you may
+want to convert to `Int` or `Float64` first via `int(hadamard(n))` or
+`float(hadamard(n))`, respectively.
+
+For many sizes, the Hadamard matrix is not unique; the `hadamard`
+function returns an arbitrary choice.  For power-of-two sizes, the
+choice is equivalent to `ifwht_natural(eye(n), 1)`.
 
 ## Author
 
diff --git a/src/Hadamard.jl b/src/Hadamard.jl
@@ -201,23 +201,86 @@ for f in (:ifwht_natural, :ifwht_dyadic, :ifwht)
 end
 
 ############################################################################
-# create (floating-point) Hadamard matrix, similar to Matlab function,
-# in natural order.
+# Utilities to work with a precomputed cache of known Hadamard matrices
+# of various sizes, produced by util/fetchhadamard.jl from Sloane's web page
+# and stored as BitMatrices
+
+function readcache(cachefile::String)
+    B = BitMatrix[]
+    open(cachefile, "r") do io
+        while !eof(io)
+            k = int(ntoh(read(io, Int64)))
+            b = BitArray(k, k)
+            # Hack: use internal binary data from BitArray for efficiency
+            bits = read(io, eltype(b.chunks), length(b.chunks))
+            for i = 1:length(bits)
+                b.chunks[i] = ntoh(bits[i])
+            end
+            push!(B, b)
+        end
+    end
+    return B
+end
+readcache() = readcache(joinpath(Pkg.dir("Hadamard"), "src", "cache.dat"))
+
+function printsigns{T<:Real}(io::IO, A::AbstractMatrix{T})
+    m, n = size(A)
+    println(io, m, "x", n, " sign matrix from ", typeof(A))
+    for i = 1:m
+        for j = 1:n
+            print(io, A[i,j] > 0 ? "+" : "-")
+        end
+        println(io)
+    end
+end
+printsigns(A) = printsigns(STDOUT, A)
+
+function frombits(B::BitMatrix)
+    A = convert(Matrix{Int8}, B)
+    for i = 1:length(A)
+        A[i] = A[i] == 0 ? -1 : 1
+    end
+    return A
+end
+
+const hadamards = BitMatrix[] # read lazily below
+
+############################################################################
+# Create Int8 order-n Hadamard matrix (in natural order), by factorizing
+# n into a product of known Hadamard sizes (if possible)
+
+const H2 = Int8[1 1; 1 -1]
 
 function hadamard(n::Integer)
-    if power_of_two(n)
-        return ifwht_natural(eye(n), 1)
-    else
-        # Punt.  Matlab supports power-of-two multiples of 12 or 20,
-        # but lots of other sizes are known as well (see e.g. Sloane's
-        # web page http://neilsloane.com/hadamard/) and there is no
-        # particularly efficient way to construct them other than
-        # a look-up-table of known factors.  [Given Hadamard matrices Hm
-        # and Hn of sizes m and n, a Hadamard matrix of size mn is kron(Hm,Hn).]
-        # It is not really clear what factors are important to support,
-        # nor which of the (generally nonunique) Hadamard matrices to pick.
-        throw(ArgumentError("non power-of-two Hadamard matrices aren't supported"))
+    n < 1 && throw(ArgumentError("size n=$n should be positive"))
+    n0 = n
+    H = reshape(Int8[1], 1,1)
+    if !ispow2(n)
+        isempty(hadamards) && append!(hadamards, readcache())
+        for i = length(hadamards):-1:1
+            k = size(hadamards[i], 1)
+            if rem(n, k) == 0
+                Hk = frombits(hadamards[i])
+                while true
+                    H = kron(H, Hk)
+                    n = div(n, k)
+                    rem(n, k) != 0 && break
+                end
+            end
+        end
+    end
+    if !ispow2(n)
+        n >>= trailing_zeros(n)
+        throw(ArgumentError("unknown Hadamard factor $n of $n0"))
+    end
+    # Note: it would be faster to do a "power-by-squaring" like algorithm
+    # here, where we repeatedly double the order via kron(H, H), but
+    # it's not clear to me that we care much about performance here.
+    while n > 1
+        H = kron(H, H2)
+        n >>= 1
     end
+    return H
 end
 
 ############################################################################
diff --git a/src/cache.dat b/src/cache.dat
diff --git a/test/runtests.jl b/test/runtests.jl
@@ -1,4 +1,5 @@
 using Hadamard
+using Base.Test
 
 H8 = [   1     1     1     1     1     1     1     1
          1     1     1     1    -1    -1    -1    -1
@@ -29,19 +30,19 @@ H8d = [  1     1     1     1     1     1     1     1
 
 ieye(n) = int8(eye(n))
 
-@assert ifwht(eye(8),1) == H8
-@assert ifwht(eye(8),2)' == H8
-@assert ifwht_natural(eye(8),1) == H8n
-@assert ifwht_natural(eye(8),2)' == H8n
-@assert ifwht_dyadic(eye(8),1) == H8d
-@assert ifwht_dyadic(eye(8),2)' == H8d
+@test ifwht(eye(8),1) == H8
+@test ifwht(eye(8),2)' == H8
+@test ifwht_natural(eye(8),1) == H8n
+@test ifwht_natural(eye(8),2)' == H8n
+@test ifwht_dyadic(eye(8),1) == H8d
+@test ifwht_dyadic(eye(8),2)' == H8d
 
-@assert ifwht(ieye(8),1) == H8
-@assert ifwht(ieye(8),2)' == H8
-@assert ifwht_natural(ieye(8),1) == H8n
-@assert ifwht_natural(ieye(8),2)' == H8n
-@assert ifwht_dyadic(ieye(8),1) == H8d
-@assert ifwht_dyadic(ieye(8),2)' == H8d
+@test ifwht(ieye(8),1) == H8
+@test ifwht(ieye(8),2)' == H8
+@test ifwht_natural(ieye(8),1) == H8n
+@test ifwht_natural(ieye(8),2)' == H8n
+@test ifwht_dyadic(ieye(8),1) == H8d
+@test ifwht_dyadic(ieye(8),2)' == H8d
 
 H32 = [  1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
  1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1
@@ -76,23 +77,38 @@ H32 = [  1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
  1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1
  1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 ]
 
-@assert ifwht(eye(32),1) == H32
-@assert ifwht(eye(32),2)' == H32
+@test ifwht(eye(32),1) == H32
+@test ifwht(eye(32),2)' == H32
 
-@assert ifwht(ieye(32),1) == H32
-@assert ifwht(ieye(32),2)' == H32
+@test ifwht(ieye(32),1) == H32
+@test ifwht(ieye(32),2)' == H32
 
 X = reshape(sin([1:1024*32]), 1024,32);
 norminf(A) = maximum(abs(A))
 
 for f in (:fwht, :fwht_natural, :fwht_dyadic)
     fi = symbol(string("i", f))
     @eval begin
-        @assert norminf(X - $fi($f(X))) < 1e-14
-        @assert norminf(X - $fi($f(X,1),1)) < 1e-14
-        @assert norminf(X - $fi($f(X,2),2)) < 1e-14
-        @assert norminf($f($f(X,1),2) - $f(X)) < 1e-14
-        @assert norminf($f(X',1)' - $f(X,2)) < 1e-14
+        @test norminf(X - $fi($f(X))) < 1e-14
+        @test norminf(X - $fi($f(X,1),1)) < 1e-14
+        @test norminf(X - $fi($f(X,2),2)) < 1e-14
+        @test norminf($f($f(X,1),2) - $f(X)) < 1e-14
+        @test norminf($f(X',1)' - $f(X,2)) < 1e-14
     end
 end
 
+@test hadamard(8) == H8n
+@test ifwht_natural(eye(32), 1) == hadamard(32)
+@test ifwht_natural(eye(2), 1) == hadamard(2)
+
+for i = 4:4:1000
+    H = int(try
+        hadamard(i)
+    catch
+        Int[]
+    end)
+    if !isempty(H)
+        println("checking unitarity of hadamard($i)")
+        @test vecnorm(H'*H - size(H,1)*I) == 0
+    end
+end
diff --git a/util/fetchhadamard.jl b/util/fetchhadamard.jl
@@ -0,0 +1,137 @@
+# Fetch a bunch of precomputed Hadamard matrices from neilsloane.com/hadamard/
+# and convert them into BitArrays.  Then print the sizes and bitarray 64-bit
+# chunks to give us a compact representation of the data for inlining.
+
+const sloane = "http://neilsloane.com/hadamard/had." # URL prefix
+
+# There are multiple inequivalent Hadamard matrices for many sizes;
+# we just pick one, semi-arbitrarily.
+
+const files = [
+# powers of 2 (and 1) are handled specially
+#               1 => "1",
+#               2 => "2",
+               12 => "12",
+               20 => "20.1",
+               28 => "28.pal2",
+               36 => "36.pal2",
+               44 => "44.pal",
+               52 => "52.will",
+               60 => "60.pal",
+               68 => "68.pal",
+               72 => "72.pal",
+               76 => "76.pal2",
+               84 => "84.pal",
+               88 => "88.tpal",
+               92 => "92.will",
+               96 => "96.tpal",
+               100 => "100.will",
+               104 => "104.pal",
+               108 => "108.pal",
+               116 => "116.will",
+               120 => "120.tpal",
+               124 => "124.pal2",
+               132 => "132.pal",
+               136 => "136.tpal",
+               140 => "140.pal",
+               148 => "148.pal2",
+               152 => "152.pal",
+               156 => "156.will",
+               164 => "164.pal",
+               168 => "168.pal",
+               172 => "172.will",
+               180 => "180.pal",
+               184 => "184.twill",
+               188 => "188.tur",
+               196 => "196.pal2",
+               204 => "204.pal2",
+               212 => "212.pal",
+               216 => "216.tpal",
+               220 => "220.pal2",
+               228 => "228.pal",
+               232 => "232.twill",
+               236 => "236.od",
+               244 => "244.will",
+               248 => "248.twill",
+               252 => "252.pal",
+               ]
+
+# the inverse of the printsigns function: parse +-+-... data (Sloane's format)
+# into a bitarray.
+function parsesigns(io::IO)
+    lines = BitVector[]
+    linelen = -1
+    for s in eachline(io)
+        s = chomp(s)
+        line = BitVector(length(s))
+        for (i,c) in enumerate(s)
+            line[i] = c=='+' ? true : c=='-' ? false : error("invalid char $c")
+        end
+        if linelen == -1
+            linelen = length(line)
+        elseif linelen != length(line)
+            error("line ", length(lines)+1, " has length ", length(line),
+                  " != ", linelen)
+        end
+        push!(lines, line)
+        length(lines) == linelen && break
+    end
+    B = BitArray(length(lines), linelen)
+    for i = 1:length(lines)
+        for j = 1:linelen
+            B[i,j] = lines[i][j]
+        end
+    end
+    B
+end
+parsesigns(s::String) = parsesigns(IOBuffer(s))
+
+open("cache.dat", "w") do io
+    f = tempname()
+    for k in sort!(collect(keys(files)))
+        println("FETCHING order $k...\n")
+        B = try
+            download(sloane*files[k]*".txt", f)
+            open(parsesigns, f, "r")
+        catch
+            println("# ERROR for size $k: ", sloane*files[k]*".txt", " -> ", f)
+            rethrow()
+        end
+        size(B) != (k,k) && error("size $k mismatch for ", files[k])
+        write(io, hton(int64(k)), map(hton, B.chunks))
+    end
+
+    # size 428 is in a different format, grr
+    k = 428
+    println("FETCHING order $k...\n")
+    B = try
+        download(sloane*"428.txt", f)
+        s = readall(f)
+        s = replace(s, r" +1", "+")
+        s = replace(s, r" *-1", "-")
+        parsesigns(s)
+    catch
+        println("# ERROR for size $k: ", sloane*"428.txt", " -> ", f)
+        rethrow()
+    end
+    size(B) != (k,k) && error("size $k mismatch for 428.txt")
+    write(io, hton(int64(k)), map(hton, B.chunks))
+
+    # others not on Sloane's page can be found at http://www.iasri.res.in/webhadamard/
+    #= ... the site seems to be broken (truncated output) at the moment ...
+    for k in [260,268,276,284,288,292,300,308,316,324,332,340,348,356,364,372,380,388,396,404,412,420,436,444,452,460,468,476,484,492,500,508,516,520,524,532,536,540,548,552,556,564,568,572,576,580,584,588,596,600,604,612,616,620,628,632,636,644,648,652,660,664,676,680,684,692,696,700,708,712,724,728,732,740,744,748,756,760,764,772,776,780,788,792,796,804,808,812,820,824,828,836,840,844,852,860,868,872,884,888,900,904,908,916,920,924,932,936,940,948,952,956,964,968,972,980,984,988,996,1000]
+        B = try
+            download("http://www.iasri.res.in/webhadamard/Validate.dll?MfcISAPICommand=Hadamard&Order=$k&Normalization=3&output=2", f)
+            s = readall(f)
+            s = replace(s, r"</p><p> *", "\n")
+            s = replace(s, r"<html>.*<p> *"s, "")
+            parsesigns(s)
+        catch
+            println("# ERROR for size $k from www.iasri.res.in/webhadamard")
+            rethrow()
+        end
+        size(B) != (k,k) && error("size $k mismatch for webhadamard")
+    end
+    write(io, hton(int64(k)), map(hton, B.chunks))
+    =#
+end
