Burrows–Wheeler transform: Difference between revisions

Line 1,942:

An efficient implementation, based mainly on http://spencer-carroll.com/an-easy-to-understand-explanation-of-the-burrows-wheeler-transform/

Takes around about ten seconds to transform and invert a 100K string. Note: requires 0.8.0+

<lang Phix>--~~demo\rosetta\burrows_wheeler.exw~~

-- demo\rosetta\burrows_wheeler.exw

--/*

The traditional method:

7 banana$ $banana 6

6 $banana ~~===>~~ a$~~banan~~ 5

7 banana$ $banana 6

5 a$~~banan~~ ~~ana~~$~~ban~~ 3

6 $banana ===> a$banan 5

4 na$~~bana~~ ~~sort~~ ~~anana~~$b 1

5 a$banan ana$ban 3

3 ~~ana~~$~~ban~~ ~~banana~~$ 7

4 na$bana sort anana$b 1

2 ~~nana~~$ba ~~===>~~ na$~~bana~~ 4

3 ana$ban banana$ 7

1 ~~anana~~$b ~~nana~~$ba 2

2 nana$ba ===> na$bana 4

~~^ desired answer == "annb$aa"~~

1 anana$b nana$ba 2

^ desired answer == "annb$aa"

First ignore the numbers: the desired answer is found by creating a table of all

rotations of "banana$", sorting it, and then extracting the right-hand column.

However, there is no need to actually create such a table, which could be very

expensive for long strings, instead just number them logically (admittedly that

was somewhat arbitrarily chosen to get the indexes to work out nicely, I picked

the original index of the last character), and perform a custom sort on those.

The latter effectively just recreates the rotations one character at a time until

there is a mismatch (which there always will be since there is only one $).

The left hand column is my arbitrary numbering scheme and the right hand column

is those sorted into order, which is also the indexes to the original string of

the characters that we want.

The code below uses $ as the terminator, but eg 1 (== '\#01') should be fine,

except of course for the display of that on a console.

--*/

--*/

constant terminator = '$'

with javascript_semantics

constant terminator = '$'

function rot_sort(integer i,j, sequence s)

-- i,j are indexes of the last character, so bump before first compare.

function rot_sort(integer i,j, sequence s)

-- eg/ie rot_sort(i,j,s) should yield compare(rotate(s,i),rotate(s,j)),

-- i,j are indexes of the last character, so bump before first compare.

-- as in rot_sort(7,6,"banana$") == compare("banana$","$banana")

-- eg/ie rot_sort(i,j,s) should yield compare(rotate(s,i),rotate(s,j)),

-- - but one character at a time rather than constructing both.

-- as in rot_sort(7,6,"banana$") == compare("banana$","$banana")

integer l = length(s)

-- - but one character at a time rather than constructing both.

while true do

integer l = length(s)

i = mod(i,l)+1

while true do

j = mod(j,l)+1

i = mod(i,l)+1

integer c = compare(s[i],s[j])

j = mod(j,l)+1

if c!=0 then return c end if

integer c = compare(s[i],s[j])

end while

if c!=0 then return c end if

end function

end while

end function

function burrows_wheeler_transform(string s)

if find(terminator,s) then ?9/0 end if

function burrows_wheeler_transform(string s)

s &= terminator

if find(terminator,s) then return "error" end if

integer l = length(s)

s &= terminator

sequence t = custom_sort(routine_id("rot_sort"),tagset(l),{s})

integer l = length(s)

string res = repeat(' ',l)

sequence t = custom_sort(routine_id("rot_sort"),tagset(l),{s})

for i=1 to l do

string res = repeat(' ',l)

res[i] = s[t[i]]

for i=1 to l do

end for

res[i] = s[t[i]]

return res

end for

end function

return res

end function

--/*

Inversion. The traditional method is add column and sort, seven times,

--/*

to reconstruct the table above, then pick the entry that ends with the

Inversion. The traditional method is add column and sort, seven times,

marker. Showing that technique in full detail here is not helpful, and

to reconstruct the table above, then pick the entry that ends with the

like above that would be hideously inefficient for large strings.

marker. Showing that technique in full detail here is not helpful, and

like above that would be hideously inefficient for large strings.

$banana 1 $ (1 ) a 2

a$banan 2 a ( 1) n 6

~~ana~~$~~ban~~ 3 a ( 2) n 7

$banana 1 $ (1 ) a 2

~~anana~~$b 4 a ( 3) b 5

a$banan 2 a ( 1) n 6

~~banana~~$ 5 b ~~$ 1~~

ana$ban 3 a ( 2) n 7

na$~~bana~~ 6 n (2 ) a 3

anana$b 4 a ( 3) b 5

~~nana~~$ba 7 n (3 ) a 4

banana$ 5 b $ 1

^ ^ ~~^ ^ ^~~

na$bana 6 n (2 ) a 3

f l ~~f l t~~

nana$ba 7 n (3 ) a 4

^ ^ ^ ^ ^

f l f l t

However, we already have the last column, and the first is just that

sorted alphabetically, and with just those two, we have all possible

However, we already have the last column, and the first is just that

character pairings of the original message. The trick is in figuring

sorted alphabetically, and with just those two, we have all possible

out how to stitch them together in the right order. If you carefully

character pairings of the original message. The trick is in figuring

study the three that end in a, and the three that start in a, notice

~~the~~ ~~$banan,na$ban,nana$b~~ ~~parts~~ ~~are~~ ~~sorted~~ in the ~~same~~ order, ~~whether~~

out how to stitch them together in the right order. If you carefully

study the three that end in a, and the three that start in a, notice

they are prefixed with a or not. That is, the middle (parenthesised)

the $banan,na$ban,nana$b parts are sorted in the same order, whether

matching numbers are both 123, not 123 and say 231. It is quite hard

they are prefixed with a or not. That is, the middle (parenthesised)

to see that being useful, but eventually the penny should drop. The

matching numbers are both 123, not 123 and say 231. It is quite hard

right-hand 1 with an a rotated right gives the left-hand 1, and the

to see that being useful, but eventually the penny should drop. The

same goes for 2 and 3: they are in fact links to the prior pairing.

right-hand 1 with an a rotated right gives the left-hand 1, and the

In ~~other~~ ~~words~~ ~~the~~ ~~first~~ a in l ~~always corresponds~~ to the ~~first in~~ f,

same goes for 2 and 3: they are in fact links to the prior pairing.

the second to the second, and so on, and that (amazingly) forms the

In other words the first a in l always corresponds to the first in f,

order in which the pairings need to be daisy-chained together.

the second to the second, and so on, and that (amazingly) forms the

order in which the pairings need to be daisy-chained together.

Try following (1->)2a->6n->3a->7n->4a->5b->$, == reverse("banana"),

in the above f and t tables.

Try following (1->)2a->6n->3a->7n->4a->5b->$, == reverse("banana"),

in the above f and t tables.

The code below builds a queue of 'a' ({1,6,7}, built backwards) from

l (aka s), into which we pop into t the {2,3,4} of the 'a' in f, and

The code below builds a queue of 'a' ({1,6,7}, built backwards) then

likewise for all other letters, forming the links for each pairing.

we pop {2,3,4} into those slots in t as we find 'a' in f, likewise

See the trivial step 3 scan below, then go back and stare at f and

for all other letters, forming the links for each pairing as shown.

t as shown above, and once again, eventually the penny should drop.

See the trivial step 3 scan below, then go back and stare at f and

I will admit I had to read ten or so explanations before I got it.

t as shown above, and once again, eventually the penny should drop.

I will admit I had to read ten or so explanations before I got it.

--*/

--*/

function inverse_burrows_wheeler(string s)

function inverse_burrows_wheeler(string s)

if find('\0',s) then ?9/0 end if -- (doable, but needs some +1s)

if find('\0',s) then ?9/0 end if -- (doable, but needs some +1s)

integer l = length(s), c

integer l = length(s), c

string f = sort(s)

string f = sort(s)

sequence q = repeat(0,256), -- queue heads (per char)

sequence q = repeat(0,256), -- queue heads (per char)

x = repeat(0,l), -- queue links

x = repeat(0,l), -- queue links

t = repeat(0,l) -- reformed/pairing links

t = repeat(0,l) -- reformed/pairing links

-- Step 1. discover/build queues (backwards)

-- Step 1. discover/build queues (backwards)

for i=l to 1 by -1 do

for i=l to 1 by -1 do

c = s[i]

c = s[i]

x[i] = q[c]

x[i] = q[c]

q[c] = i

q[c] = i

end for

end for

-- Step 2. reform/pop char queues into pairing links

-- Step 2. reform/pop char queues into pairing links

for i=1 to l do

for i=1 to l do

c = f[i]

c = f[i]

t[q[c]] = i

t[q[c]] = i

q[c] = x[q[c]]

q[c] = x[q[c]]

end for

end for

-- Step 3. rebuild (backwards)

-- Step 3. rebuild (backwards)

c = find(terminator,f)

c = find(terminator,f)

string res = repeat(' ',l-1)

if c=0 then return "error" end if

for i=l-1 to 1 by -1 do

string res = repeat(' ',l-1)

c = t[c] -- (first time in, skip the end marker)

for i=l-1 to 1 by -1 do

res[i] = f[c]

c = t[c] -- (first time in, skip the end marker)

end for

res[i] = f[c]

return res

end for

end function

return res

end function

procedure test(string src)

string enc = burrows_wheeler_transform(src),

procedure test(string src)

dec = inverse_burrows_wheeler(enc)

string enc = burrows_wheeler_transform(src),

printf(1,"original: %s --> %s\n inverse: %s\n",{src,enc,dec})

dec = inverse_burrows_wheeler(enc)

end procedure

src = shorten(src,"characters")

enc = shorten(enc,"characters")

test("banana")

dec = shorten(dec,"characters")

test("dogwood")

printf(1,"original: %s --> %s\n inverse: %s\n",{src,enc,dec})

test("TO BE OR NOT TO BE OR WANT TO BE OR NOT?")

end procedure

test("SIX.MIXED.PIXIES.SIFT.SIXTY.PIXIE.DUST.BOXES")</lang>

test("banana")

test("dogwood")

test("TO BE OR NOT TO BE OR WANT TO BE OR NOT?")

test("SIX.MIXED.PIXIES.SIFT.SIXTY.PIXIE.DUST.BOXES")

<pre>