Compressed frequencies: Representing frequencies with less bits

In a cache memory you usually store the most frequently used data that is currently in a larger but slower memory. For example, you keep your most frequently accessed files cached in RAM rather than on your hard drive. Since you can't fit all the contents of your hard disk in RAM, you keep only the most frequently used files that can fit. You will still need to access your hard disk once in a while in order to access your less frequently used files but the average file access time will now be greatly reduced.

The problem is how to keep count of the number of times each file is being used in order to know which is the most frequently used. The obvious solution is to associate each file with a number and increment that number each time it is used. But numbers take space as well, and sometimes this becomes a significant problem. You might not afford to waste 4 or 8 bytes of memory worth of frequency integers for every item. Is there a way to bring down the number of bytes used by the frequency integers without losing their usefulness?

Here is an example of a 4 byte int integer number in memory representing the number 45723:

00000000 00000000 10110010 10011011

The most obvious thing you can do is to tear away the most significant bits (the ones on the left which have a larger value) by using smaller range number types such as the 2 byte short, which gives us 10110010 10011011. If the frequency is sometimes, but rarely, larger than the ranges provided by these smaller types, then you can just cap it off by stopping incrementation once the maximum number is reached. For example, if we're using a two byte short, then the maximum number this can be is 65535. Once the frequency reaches this number, then it gets frozen and never incremented again. Many frequences follow a zipfian distribution, meaning that the vast majority of items will have a small frequency, followed by a handful of very frequent items. An example of this is words in a document where most words will only occur once and only a few words such as "the" and "of" will occur frequently. If this is the case then you will be fine with capping off your frequencies since only a few items will have a large frequency and it might not be important to order these high frequency items among themselves.

It might seem more useful instead to tear away the least significant bits (the ones on the right which have a smaller value) instead, since these are less useful. The way you do this is to divide the frequency by a constant and keep only the whole number part. For example, if we divide the above number by 256, we'd be shifting the bits by one byte to the right, which gives us 00000000 00000000 00000000 10110010. The least significant byte has been removed which means that we can use less bytes to store the frequency. But in order to do that you need to first have the actual frequency which defeats the purpose. So what we can do is to simulate the division by incrementing the frequency only once every 256 times. If we do that then the resulting number will always be a 256th of the actual frequency which is the frequency without the least significant byte. But how do you know when to increment the frequency next? If you keep a separate counter which counts to 256 in order to know when to increment next then you lose the space you would have saved. Instead we can do it stochastically using random numbers. Increment the frequency with a probability of 1 in 256 and the frequency will be approximately a 256th of the actual frequency.

By combining these two ideas together we can reduce an 8 byte frequency into a single byte and that byte will be one of the original 8 bytes of the actual frequency. Here is a Python function that increments an integer with a compressed frequency that is a certain number of bytes long and with a certain number of least significant bytes torn off.

def compressed_increment(frequency, bytes_length, bytes_torn):
    if frequency < 256**bytes_length: #cap the frequency to the maximum number that can be contained in bytes_length bytes
        if random.randint(1, 256**bytes_torn) == 1: #increment with a probability of 1 in 256^bytes_length (number to divide by to shift the frequency by that number of bytes)
            return frequency + 1

Of course this is a lossy compression. Information is lost. This means that the compressed frequency is not useful in certain situations, such as when you want to also decrement the frequency or when approximate frequencies are inadequate.

Translating an arbitrary integer into a circular array index

A circular array is an array which is connected at the edges, that is, it has no beginning or end and traversing the array will eventually get you back to where you started. Of course in practice a normal array is used and the given index is mapped into a valid array index. This is usually done using modulo operations (the remainder after dividing the index by the array length). But what if you need to also allow negative indexes?

Let's say you have an array of length 5 and you want to use it as a circular array.

0

1

2

3

4

If you start at index 2 and move 1 to the right then you end up in index 3. But if you move 3 to the right you end up in index 0. On the other hand if you move 3 to the left then you end up in index 4. Here are some other examples of this translation:

Starting index	Add to it	Resultant	Translated index
2	+1	3	3
2	-1	2	2
4	+1	5	0
0	-1	-1	4
2	+3	4	0
2	-3	-1	4
2	+7	9	4
2	-7	-5	0
2	+15	17	2
2	-15	-13	2

We need a general formula to map arbitrary resultant integers into corresponding indexes. The modulo operator will not map negative numbers correctly:

6 % 5 = 1
5 % 5 = 0
4 % 5 = 4
3 % 5 = 3
2 % 5 = 2
1 % 5 = 1
0 % 5 = 0
-1 % 5 = -1
-2 % 5 = -2
-3 % 5 = -3
-4 % 5 = -4
-5 % 5 = 0
-6 % 5 = -1

This is because if the whole number division of a negative number N divided by a positive number P is D, then the remainder would be the number X such that D*P + X = N holds. For example, 4/5 = 0 remainder 4, because 0*5 + 4 = 4. Another example, -4/5 = 0 remainder -4 because 0*5 + -4 = -4.

Now in order to get a mapping from arbitrary resultant integers to corresponding indexes in a circular array we need to use the following formula:

Given a resultant R and a length of array L, the corresponding index is

(R%L + L)%L

Fractions and decimals

Let's take a break from irrational numbers and focus a bit on the rational ones. Rational numbers can be either be of a fixed number of decimal digits such as 0.123, called terminating decimals, or have a part of its decimal digits which repeat forever such as 0.272727..., called recurring decimals. Both terminating and recurring decimals can be represented as fractions. Let's see how to convert one to the other.

Fraction to decimal

Terminating decimals

The way to convert a fraction to a decimal is of course through the familiar long division algorithm, which was already shown in a previous blog post of mine. The way I show it here is how I learned it at school, which is a fast method but which leaves nearly nothing explained in its working. Let's convert 43/5 to decimal form:

          
5 )4 3

We start by seeing how many times the denominator 5 goes into the first digit of the numerator, 4. It goes 0 times into it and leaves a remainder of 4. We write the integer quotient as the first digit at the top. The remainder we write in front of the next digit in the numerator.

   0         
5 )4 43

We now see how many times 5 goes into 43. It goes 8 times into it and leaves 3 as a remainder. We write the integer quotient as the second digit at the top. We now have no more digits in the numerators, so we add a decimal point and a zero after it to create more digits. We also add a decimal point to the quotient at the top.

   0  8 .     
5 )4 43 . 30

And repeat.

   0  8 .  6  
5 )4 43 . 30 00

Since the last remainder was 0, if we had to continue from here onwards we'd be adding nothing by zeros to the top quotient which would be pointless. So instead we declare the number as terminating decimal and stop there. The answer the top quotient, that is, 43/5 = 8.6

What's happening here is that we're first trying to find the tens digit of the quotient by seeing how many times 50 goes into 43 which is 0 (tens) remainder 43. Then we're trying to find the units digit of the quotient by seeing how many times 5 goes into 43 which is 8 (units) remainder 3. Then we're trying to find the tenths digit of the quotient by seeing how many times 0.5 goes into 3, or equivalently, how many times 5 goes into 30, which is 6 (tenths) remainder 0.

Recurring decimals

Let's use the previous method to convert 1/3 to decimal form:

          
3 )1

We start by seeing how many times the denominator 3 goes into the first digit of the numerator, 1. It goes 0 times into it and leaves a remainder of 1.

   0 .       
3 )1 . 10

We now see how many times 3 goes into 10. It goes 3 times into it and leaves 1 as a remainder.

   0 .  3    
3 )1 . 10 10

And repeat.

   0 .  3  3 
3 )1 . 10 10 10

As you can see, this process will continue indefinitely, meaning that the 3 at the top will keep on repeating itself. This makes the quotient a recurring decimal, the proof of which is that one of the remainders after the decimal point was reached appeared twice and hence the same number must come out again.

Decimal to fraction

Terminating decimals

To convert a terminating decimal to a fraction you simply multiply it by a power of 10 that is large enough to make it a whole number, then put that same power of 10 as the denominator under the whole number. For example, if you have 12.3456 to convert to fraction form:

12.3456 = 12.3456 x 10000 / 10000
        = 123456/10000
        = 7716/625

Recurring decimals

Of course the previous method cannot be used when the fractional part is infinitely long. But there is a trick we can use. The fraction 1/9 has the following decimal expansion:

   0 .  1  1 ...
9 )1 . 10 10 ...

In other words, it gives a decimal number with an infinite sequence of 1s. We can use this to our advantage so that we can obtain infinite sequences of any digit by simply multiplying the digit by 1/9.

1/9 = 0.111...
2/9 = 0.222...
3/9 = 0.333...
etc

The interesting thing about this is that when you try to get an infinite sequence of 9s, you get the whole number 1. This is one of the proofs that 0.999... = 1.

But this is only good for single digit recurrences. For two digit recurrences we can use 1/99:

    0 .  0  1  0  1 ...
99 )1 . 10 100 10 100 ...

This is useful. We can replace every 1 with any digit by multiplying 1/99 by that digit.

1/99 = 0.0101...
2/99 = 0.0202...
3/99 = 0.0303...
etc

We can also shift those digits one place to the left by multiplying the digit by 10:

10/99 = 0.1010...
20/99 = 0.2020...
30/99 = 0.3030...
etc

So we can control both digits by adding these two fractions together:

10/99 + 3/99 = 0.1010... + 0.0303... = 0.1313...

Which of course simplifies to

13/99 = 0.1313...

That is, you just multiply the 2 digit number you want to be repeated by 99.

Can this be extended to any number of digits? 1/999 gives the following expansion:

     0 .  0  0   1  0  0   1 ...
999 )1 . 10 100 1000 10 100 1000 ...

So with 1/999 we can repeat any 3 digit recurrences. With each new 9 added to the denominator we are postponing the number of times the remainder must be moved before it is big enough to be divided by the denominator. That postponement results in another 0 added to the recurring decimal.

So in short, to convert a recurring decimal number to a fraction, just take the repeating part of the decimal and put it over a denominator consisting of as many 9s as there are digits in the repeated number:

0.01230123... = 0123/9999

But this will only work if the decimal number consists of nothing but repeating numbers. What if you have the following decimal:

210.67801230123...

In this case, the repeating part is the "0123" but it starts with other non-repeating digits. Not to worry, we just separate the two parts of the number into a sum:

210.67801230123... = 210.678 + 0.00001230123...

The recurring term can be multiplied by a power of 10 that will make it lose its leading zeros. Of course we can't just multiply it by a power of 10 without also dividing it by the same number in order to avoid changing its value:

210.67801230123... = 210.678 + (1000 x 0.00001230123... / 1000)
210.67801230123... = 210.678 + (0.01230123... / 1000)

We can now convert each term separately:

210.67801230123... = 210.678 + (0.01230123... / 1000)
                   = 210678/1000 + ((0123/9999) / 1000)
                   = 140437963/666600

Finding the nth row in every group in SQL

Let's say you have the following log table which stores the dates of each access:

TimeStamp	IP	UserName
2005-10-30 10:45:03	172.16.254.10	jlor
2005-10-30 10:46:31	172.16.254.12	kpar
2005-10-31 09:14:13	172.16.254.14	jlor
2005-10-31 09:25:42	172.16.254.16	kpar
2005-10-31 12:41:14	172.16.254.19	jlor
2005-11-01 07:15:15	172.16.254.20	kpar

You are asked to make a report of the last time each user has accessed the system using SQL.

At first you try using GROUP BY but then realize that it's not so simple to include the IP field along with the TimeStamp and UserName. GROUP BY works when you're interested in aggregating every field that is not used to group the records. In other words, you can easily do this:

SELECT UserName, MAX(TimeStamp)
FROM log
GROUP BY UserName

USERNAME	MAX(TIMESTAMP)
jlor	October, 31 2005 12:41:14+0000
kpar	November, 01 2005 07:15:15+0000

But if you also want to show the corresponding IP address of the access with the latest time stamp, you'd have a problem using simple SQL. If you add the IP field in the SELECT statement, you'd end up with the first IP in the table that belongs to the corresponding user, rather than the IP of the latest time stamp.

SELECT UserName, IP, MAX(TimeStamp)
FROM log
GROUP BY UserName

USERNAME	IP	MAX(TIMESTAMP)
jlor	172.16.254.10	October, 31 2005 12:41:14+0000
kpar	172.16.254.12	November, 01 2005 07:15:15+0000

The way to do this is to simulate the GROUP BY statement using more expressiveness methods.

MS SQL Server

In MS SQL Server, this is achieved using the ROW_NUMBER function. This function gives a number for each row (1, 2, 3, ...) which can be used inside a SELECT statement. The cool thing about this function is that the numbering can be made to restart for every different value in a field. So if we used it on the UserName field we'd have the following:

SELECT UserName, IP, TimeStamp, ROW_NUMBER() OVER(PARTITION BY UserName ORDER BY TimeStamp DESC)
FROM log

USERNAME	IP	TIMESTAMP	COLUMN_3
jlor	172.16.254.10	October, 30 2005 10:45:03+0000	1
jlor	172.16.254.14	October, 31 2005 09:14:13+0000	2
jlor	172.16.254.19	October, 31 2005 12:41:14+0000	3
kpar	172.16.254.12	October, 30 2005 10:46:31+0000	1
kpar	172.16.254.16	October, 31 2005 09:25:42+0000	2
kpar	172.16.254.20	November, 01 2005 07:15:15+0000	3

It even orders the rows by user name and it lets you say how you want the rows of each user to be ordered so that you can say how you want the numbering. Using the SQL above, the row with the latest time stamp of each user has a 1 in the last column. This allows us to select it. It will have to be inside a nested query however in order to be used in a WHERE statement.

SELECT UserName, IP, TimeStamp
FROM (
  SELECT UserName, IP, TimeStamp, ROW_NUMBER() OVER(PARTITION BY UserName ORDER BY TimeStamp DESC) AS rank
  FROM log
) AS t
WHERE rank = 1

USERNAME	IP	TIMESTAMP	COLUMN_3
jlor	172.16.254.19	October, 31 2005 12:41:14+0000
kpar	172.16.254.20	November, 01 2005 07:15:15+0000

Notice that you can even find when the second to last time an access was made by changing the 1 in the WHERE statement to a 2.

You can experiment with this in this SQL Fiddle.

MySQL

Unfortunately MySQL doesn't have a function as nifty as ROW_NUMBER so instead we'll have to simulate that using variables. In MySQL you can create variables using the SET statement and then update them within a SELECT statement so that they change for each row, like this:

SET @row_number := 0;
SELECT UserName, IP, TimeStamp, @row_number := @row_number + 1
FROM log

USERNAME	IP	TIMESTAMP	@ROW_NUMBER := @ROW_NUMBER + 1
jlor	172.16.254.10	October, 30 2005 10:45:03+0000	1
kpar	172.16.254.12	October, 30 2005 10:46:31+0000	2
jlor	172.16.254.14	October, 31 2005 09:14:13+0000	3
kpar	172.16.254.16	October, 31 2005 09:25:42+0000	4
jlor	172.16.254.19	October, 31 2005 12:41:14+0000	5
kpar	172.16.254.20	November, 01 2005 07:15:15+0000	6

This is only half the story of course. We want the numbering to restart for every user and we also want this to happen after sorting the rows by user name. We also want the rows belonging to each user to be sorted by time stamp. A simple ORDER BY statement can handle the sorting part:

SET @row_number := 0;
SELECT UserName, IP, TimeStamp, @row_number := @row_number + 1
FROM log
ORDER BY UserName, TimeStamp DESC

USERNAME	IP	TIMESTAMP	@ROW_NUMBER := @ROW_NUMBER + 1
jlor	172.16.254.19	October, 31 2005 12:41:14+0000	1
jlor	172.16.254.14	October, 31 2005 09:14:13+0000	2
jlor	172.16.254.10	October, 30 2005 10:45:03+0000	3
kpar	172.16.254.20	November, 01 2005 07:15:15+0000	4
kpar	172.16.254.16	October, 31 2005 09:25:42+0000	5
kpar	172.16.254.12	October, 30 2005 10:46:31+0000	6

The restarting of numbering is a little less simple. We have to keep track of what the previous value was using another variable and we have to also choose between setting row_number to 1 or to increment it by 1. Here is the code:

SET @row_number := 0;
SET @prev_username := NULL;
SELECT UserName, IP, TimeStamp, @row_number := CASE WHEN UserName = @prev_username THEN @row_number + 1 ELSE 1 END, @prev_username := UserName
FROM log
ORDER BY UserName, TimeStamp DESC

USERNAME	IP	TIMESTAMP	@ROW_NUMBER := CASE WHEN USERNAME = @PREV_USERNAME THEN @ROW_NUMBER + 1 ELSE 1 END	@PREV_USERNAME := USERNAME
jlor	172.16.254.19	October, 31 2005 12:41:14+0000	1	jlor
jlor	172.16.254.14	October, 31 2005 09:14:13+0000	2	jlor
jlor	172.16.254.10	October, 30 2005 10:45:03+0000	3	jlor
kpar	172.16.254.20	November, 01 2005 07:15:15+0000	1	kpar
kpar	172.16.254.16	October, 31 2005 09:25:42+0000	2	kpar
kpar	172.16.254.12	October, 30 2005 10:46:31+0000	3	kpar

The CASE statement selects a value to set row_number. If the current row's user name is the same as the previous one's then the value will be one more than row_number it currently is. Otherwise it is set to 1. After that variable is set, the prev_username variable is set to the current row's user name.

Finally we can now use this to select the latest access for each user.

SET @row_number := 0;
SET @prev_username := NULL;
SELECT UserName, IP, TimeStamp
FROM (
  SELECT UserName, IP, TimeStamp, @row_number := CASE WHEN UserName = @prev_username THEN @row_number + 1 ELSE 1 END AS rank, @prev_username := UserName
  FROM log
  ORDER BY UserName, TimeStamp DESC
) AS t
WHERE rank = 1

USERNAME	IP	TIMESTAMP
jlor	172.16.254.19	October, 31 2005 12:41:14+0000
kpar	172.16.254.20	November, 01 2005 07:15:15+0000

Notice that you can even find when the second to last time an access was made by changing the 1 in the WHERE statement to a 2.

You can experiment with this in this SQL Fiddle.

The Lempel Ziv Welch (LZW) compression algorithm

The Lempel Ziv Welch algorithm (LZW) is a classic compression algorithm published in 1984. It's a simple but practical algorithm that should be under every geek's belt and is often used in combination with other techniques.

The basic idea

Let's start with a plain English description of how this algorithm works. Let's say that we want to compress the following input string:

xxxxyyyyxxxxxxxxxxxx

---

Compression

If we're lazy, all we have to do to produce a valid output is to represent each letter using 2 characters.

0x0x0x0x0y0y0y0y0x0x0x0x0x0x0x0x0x0x0x0x

According to the compressed language of LZW, when a character pair starts with a "0", that means that the second character is the original letter. This has doubled the size of the sequence, but hopefully this is not the actual output. When the first character is something other than "0", the character pair becomes a reference to some prior long sequence. These references are what will compress the sequence.

---

In order to compress, we need to use a special table called a "dictionary" which maps 2 character values to the string they represent. Right off the bat, the dictionary will start with the following values:

Actual string	2 character value
a	0a
b	0b
...	...
x	0x
y	0y
z	0z

We start scanning through the input string from the first letter and find the longest sequence of letters which are in the dictionary. At this point, that would obviously be the first letter, since the dictionary only have single letters.

xxxxyyyyxxxxxxxxxxxx

After finding the longest string from the start which is in the dictionary, "x", including also the next letter in the input will make the shortest string which is not in the dictionary. This unregistered string is "xx".

We replace the longest string found in the dictionary with the corresponding 2 character value:

0xxxxyyyyxxxxxxxxxxxx

The next shortest string not in the dictionary is then added to the dictionary.

Actual string	2 character value
a	0a
b	0b
...	...
z	0z
xx	1a

Now we continue scanning the input string from after the last replacement.

---

Again, we look for the longest string that is in the dictionary. That would be "xx", which we added in the previous step.

0xxxxyyyyxxxxxxxxxxxx

We replace this string with the corresponding 2 character value and add into the dictionary this string plus the next letter ("xxx").

0x1axyyyyxxxxxxxxxxxx

Actual string	2 character value
a	0a
...	...
z	0z
xx	1a
xxx	1b

Repeat the process from right after the last replacement.

---

0x1axyyyyxxxxxxxxxxxx

This time it was "x" on its own that was the longest string in the dictionary since "xy" was not in the dictionary.

0x1a0xyyyyxxxxxxxxxxxx

Actual string	2 character value
...	...
z	0z
xx	1a
xxx	1b
xy	1c

Repeat.

---

0x1a0xyyyyxxxxxxxxxxxx

Longest string in dictionary was "y", shortest string not in dictionary was "yy".

0x1a0x0yyyyxxxxxxxxxxxx

Actual string	2 character value
...	...
z	0z
xx	1a
xxx	1b
xy	1c
yy	1d

Repeat.

---

0x1a0x0yyyyxxxxxxxxxxxx

Longest string in dictionary was "yy", shortest string not in dictionary was "yyy".

0x1a0x0y1dyxxxxxxxxxxxx

Actual string	2 character value
...	...
z	0z
xx	1a
xxx	1b
xy	1c
yy	1d
yyy	1e

Repeat.

---

0x1a0x0y1dyxxxxxxxxxxxx

Longest string in dictionary was "y", shortest string not in dictionary was "yx".

0x1a0x0y1d0yxxxxxxxxxxxx

Actual string	2 character value
...	...
xx	1a
xxx	1b
xy	1c
yy	1d
yyy	1e
yx	1e

Repeat.

---

0x1a0x0y1d0yxxxxxxxxxxxx

Longest string in dictionary was "xxx", shortest string not in dictionary was "xxxx".

0x1a0x0y1d0y1bxxxxxxxxx

Actual string	2 character value
...	...
xx	1a
xxx	1b
xy	1c
yy	1d
yyy	1e
yx	1e
xxxx	1f

Repeat.

---

0x1a0x0y1d0y1bxxxxxxxxx

Longest string in dictionary was "xxxx", shortest string not in dictionary was "xxxxx".

0x1a0x0y1d0y1b1fxxxxx

Actual string	2 character value
...	...
xx	1a
xxx	1b
xy	1c
yy	1d
yyy	1e
yx	1e
xxxx	1f
xxxxx	1g

Repeat.

---

0x1a0x0y1d0y1b1fxxxxx

Longest string in dictionary was "xxxxx" which resulted in the whole input string being consumed, hence there is nothing left to add to the dictionary.

0x1a0x0y1d0y1b1f1g

---

0x1a0x0y1d0y1b1f1g

This has resulted in an output string which is 2 characters shorter. Obviously the string was engineered to be compressed. Had it been a longer string then the dictionary would have contained longer strings which would lead to more compression.

---

Decompression

Compression is quite straightforward and so should decompression. Except that it's a little less straightforward because of a special case that can sneak up on you if you don't know about it (some online sources don't mention it).

If we knew what the dictionary contains then we can simply replace each 2 character value with its corresponding string. But in order to know from the start what the dictionary is we would have to include it with the output string, which would be a considerable amount of extra bytes. Instead, decompression basically consists of guessing what the dictionary contained one 2 character value at a time.

---

We start with the obvious. The dictionary surely contained all the single letters.

Actual string	2 character value
0a	a
0b	b
...	...
0x	x
0y	y
0z	z

The 2 character value is always one of the above initial dictionary (since that is the first shortest string not in the dictionary), so we go ahead and take care of that.

0x1a0x0y1d0y1b1f1g

x1a0x0y1d0y1b1f1g

From here on, if a 2 character value is in the dictionary then we just replace it with its corresponding string. After each replacement we use the replaced string to update the dictionary (as will be shown in an examples further down). The problem is if the 2 character value is not in the dictionary, as is the case now. This is the special case.

---

The next 2 character value is not in the dictionary. But notice that it is the very next 2 character value that will enter the dictionary, "1a". When this is the case, the following scenario must have taken place:

The input string has been determined to start with an "x", but the following letters are unknown.

x _ _ _

We know that the following letters were replaced with the next available 2 character value (from the compressed string), which means that it must be the shortest string that was not in the dictionary after "x" was replaced with "0x".

This shortest string must have been the last string that was replaced with a 2 character value ("x"), plus the letter after it. So the dictionary must have looked something like this:

Actual string	2 character value
...	...
0x	x
0y	y
0z	z
1a	x_

(Notice that the blank is an unknown letter)

So then the first letter of this unknown string being referred to by the 2 character value "1a" is "x". In that case then we know what the second letter is in the input string: an "x", according to the dictionary we're constructing.

x x _ _

But wait, since the 2 character value "1a" is referring to the first "x" followed by the next letter, and since we have determined that the next letter was "x", then "1a" must be referring to "xx".

Actual string	2 character value
...	...
0x	x
0y	y
0z	z
1a	xx

In general, every time we encounter this situation, where 2 character value is not in the dictionary and is the next value to be added to the dictionary, the string being referred to is the previous replaced string followed by the same string's first letter. In this case, the previous replaced string was "x" whose first letter is "x", so the string referred to by "1a" is "xx".

x1a0x0y1d0y1b1f1g

xxx0x0y1d0y1b1f1g

After every replacement after the very first (this is the second), we need to update the dictionary with a new 2 character value. The update needs to reflect what was added during compression when the previous replacement (the first in this case) took place. This is because you need to know what letter follows the replacement in order to know which string was added to the dictionary.

Remember that during compression we were adding to the dictionary the shortest string that was not in the dictionary, which consisted of the longest string found in the dictionary (the replaced string) followed by the next letter. After the very first replacement we made, the "x", the next letter in the string is the first letter of the second replacement we made, the "xx". So we add to the dictionary the previous replacement followed by the first letter of the current replacement.

Actual string	2 character value
...	...
0x	x
0y	y
0z	z
1a	xx

Finally we've covered all the steps needed to get repeatin'. Let's continue.

---

xxx0x0y1d0y1b1f1g

This one is easy as it is already in the dictionary.

xxxx0y1d0y1b1f1g

The previous replacement was "xx", the current replacement was "x". So to the dictionary we add "xxx".

Actual string	2 character value
...	...
0x	x
0y	y
0z	z
1a	xx
1b	xxx

Repeat.

---

xxxx0y1d0y1b1f1g

In the dictionary.

xxxxy1d0y1b1f1g

The previous replacement was "x", the current replacement was "y". So to the dictionary we add "xy".

Actual string	2 character value
...	...
0x	x
0y	y
0z	z
1a	xx
1b	xxx
1c	xy

Repeat.

---

xxxxy1d0y1b1f1g

Next value is the next 2 character value to be added to the dictionary. Add previous replacement, "y", followed by its own first letter. So add "yy".

xxxxyyy0y1b1f1g

The previous replacement was "y", the current replacement was "yy". So to the dictionary we add "yy".

Actual string	2 character value
...	...
0x	x
0y	y
0z	z
1a	xx
1b	xxx
1c	xy
1d	yy

Repeat.

---

xxxxyyy0y1b1f1g

In the dictionary.

xxxxyyyy1b1f1g

The previous replacement was "yy", the current replacement was "y". So to the dictionary we add "yyy".

Actual string	2 character value
...	...
0x	x
0y	y
0z	z
1a	xx
1b	xxx
1c	xy
1d	yy
1e	yyy

Repeat.

---

xxxxyyyy1b1f1g

In the dictionary.

xxxxyyyyxxx1f1g

The previous replacement was "y", the current replacement was "xxx". So to the dictionary we add "yx".

Actual string	2 character value
...	...
0x	x
0y	y
0z	z
1a	xx
1b	xxx
1c	xy
1d	yyy
1e	yx

Repeat.

---

xxxxyyyyxxx1f1g

Next value is the next 2 character value to be added to the dictionary. Add previous replacement, "xxx", followed by its own first letter. So add "xxxx".

xxxxyyyyxxxxxxx1g

The previous replacement was "xxx", the current replacement was "xxxx". So to the dictionary we add "xxxx".

Actual string	2 character value
...	...
0x	x
0y	y
0z	z
1a	xx
1b	xxx
1c	xy
1d	yyy
1e	yx
1f	xxxx

Repeat.

---

xxxxyyyyxxxxxxx1g

Next value is the next 2 character value to be added to the dictionary. Add previous replacement, "xxxx", followed by its own first letter. So add "xxxxx".

xxxxyyyyxxxxxxxxxxxx

Complete.

---

xxxxyyyyxxxxxxxxxxxx

In practice

In practice, we do not use 2 character values as that would have a large amount of overhead which reduces the amount of compression possible. Instead we work at the bit level and use 12 bit values, a little over 1 byte. The more bits are used, the bigger the dictionary can be and the longer the strings that are added will be, but this needs to be compromised with the overhead.

You can find code for different programming languages here.

Testing for uniformly distributed random events

Uniformly distributed randomness is when the probability of each individual outcome of the randomness is equal. A fair coin toss is uniformly distributed because the probability of getting a heads is equal to the probability of getting a tails.

Say you developed a program which uses randomness. How can you test that the randomness is actually uniformly distributed? An example is testing that a mutation function for a genetic algorithm mutates each bit with equal probability. Another example is that a random shuffle of a list results in each item occupying any position in the list with equal probability.

Creating the test

The way to test for uniform distribution is by running the random process multiple times and recording the output each time. Count the number of times each output is given. After many trials, the frequencies should be similar. For example, in order to test if a coin flip is uniformly distributed, count the number of times of each heads and tails. The counts should be similar. In the case of shuffling a list, make a table with a row for each position in the list and a column for each item in the list and after each shuffle, count the number of times each item ended up in each position. The counts across the whole table should be similar.

Here is an example unit test which tests if a dice simulator is uniformly distributed. Assume that there is a function called "similarity" which tells you how similar to each other a bunch of frequencies are.

import unittest

class TestDice(unittest.TestCase):

    def test_uniform(self):
        trials = 1000000
        best_similarity = similarity({ 1 : trials//6, 2 : trials//6, 3 : trials//6, 4 : trials//6, 5 : trials//6, 6 : trials//6 }.values())
        threshold = 0.9*best_similarity

        frequencies = { 1 : 0, 2 : 0, 3 : 0, 4 : 0, 5 : 0, 6 : 0 }
        for i in range(trials):
            outcome = roll_dice()
            frequencies[outcome] += 1

        self.assertTrue(similarity(frequencies.values()) >= threshold)

This program starts by creating a threshold against which to compare the similarity of the frequency list. This threshold is 90% of the similarity measure when the frequencies are all the same. The percentage can be tweaked of course. The similarity between the actual frequencies must be at least as high as this threshold for the test case to pass.

It could be that case that instead of a similarity function there is a "distance" function which is 0 when all frequencies are the same and greater otherwise. In this case the code should be a little different.

import unittest

class TestDice(unittest.TestCase):

    def test_uniform(self):
        trials = 1000000
        worst_distance = distance({ 1 : trials, 2 : 0, 3 : 0, 4 : 0, 5 : 0, 6 : 0 }.values())
        threshold = 0.1*best_distance

        frequencies = { 1 : 0, 2 : 0, 3 : 0, 4 : 0, 5 : 0, 6 : 0 }
        for i in range(trials):
            outcome = roll_dice()
            frequencies[outcome] += 1

        self.assertTrue(distance(frequencies.values()) <= threshold)

In this case instead of comparing against the best similarity, we compare against the worst similarity which is when only one of the outcomes gets all of the frequency. The threshold would then be 10% of this worst measure. The distance between the actual frequencies must not exceed this threshold for the test case to pass.

Measuring the distance/similarity

The problem is how to measure the similarity between a list of frequencies, in such a way that the similarity is quantified into a single number by a program (a unit test) in order to check if it is greater than or smaller than some threshold. Here are three methods to quantify the similarity between a list of frequencies.

Standard Deviation

The average of a bunch of frequencies is a number which is nearest to every frequency. But an average is only half of the story, it is also important to know how close the frequencies are to this average. To measure how close all the frequencies are to their average you need to find their standard deviation, which is basically the average difference between every frequency and their average. It is found in the following way:

distance(xs) = sqrt( average( (x - average(xs))^2 for x in xs ) )

The closer they are to their average, the more similar they are. So the smaller the standard deviation, the more similar the frequencies are; hence, it is a distance measure rather than a similarity measure.

Chi-squared test

Let's say that you're conducting an experiment to measure the effects of smoking cigarettes on a person's health. You take two groups of people, one group consisting of smokers and another group consisting of non-smokers, and you measure their collective health. I have no idea how to measure collective health but let's say that you have two numbers measure how healthy each group is. How do you quantify how similar the two numbers are in order to determine how significant the effect of smoking is on a person's health?

A common way to do this in science is by using Pearson's Chi-squared test which is used to measure how similar a list of pairs of numbers are. In order to apply this to similarity between frequencies, we can measure how similar each frequency is to their average.

distance(xs) = sum( (x - average(xs))^2 / average(xs) for x in xs )

This number will be smaller the more similar the frequencies are, making this a distance measure rather than a similarity measure.

Min max ratio

Of course there's a less fancy way to measure how similar a bunch of frequencies are, and that is by measuring the difference between the largest and smallest frequency. If the largest and smallest numbers are the same then all the numbers are the same. If they are different then that is the largest difference between any two numbers.

The difference between maximum and minimum isn't a useful measure of similarity since the similarity between 10 and 5 is not the same as the similarity between 1000 and 995. Instead it is more useful to divide the minimum by the maximum. But this will result in a value of zero whenever one of the frequencies is zero (and hence the minimum is zero), which is not measuring anything useful either. A solution for this is to add one to both the minimum and maximum frequency in order to avoid having zeros but still have a maximum similarity of 1.

similarity(xs) = (min(xs) + 1)/(max(xs) + 1)

Apart from this function giving a larger number the more similar the frequencies are, the difference between this measure and the others is that if all the frequencies are the same except for one, this method will identify the frequencies as different, whereas the others will be affected less. This might be considered an advantage when checking if a random process is uniformly distributed.

Dijkstra's Algorithm, Uniform Cost Search, and A* Search (graph search series)

Let's say that you're developing a game which involves moving a hero around in order to avoid a monster that is chasing him or her. It's easy to program the hero to just move around wherever the arrow keys on the keyboard indicate, but how do you make the monster automatically follow a sensible path towards the hero?

If there are no obstacles in the way then the monster just has to move in the direction of the hero, but you don't need me to tell you that. In order to get from point A to point B with obstacles in the way you need to use some path finding algorithms, which ideally find the shortest path. Let's start from Dijkstra's Algorithm.

Dijkstra's Algorithm

Dijkstra's Algorithm finds the shortest path from a point to every other point. In other words it is not ideal for finding the shortest path between two points. This might not be what you need but it's a good basis to understand the more focused algorithms.

It is a dynamic programming algorithm which exploits the following observation about shortest paths. Consider the following diagram which shows the shortest path between two crosses with an obstacle in the middle:



If you take any point in the path, such as the circle, the shortest path from this point and one of the crosses is the same path as the shortest path between the two crosses. In other words, imagine that you are searching for the shortest path between the two crosses, and you meet a wizard who tells you that he knows what this shortest path is but that he will not tell you how to take it; however, he will tell you that the circle is on this path. With this information you can simplify your search by looking for the shortest path from one cross to the circle and from the circle to the other cross. These two paths can then be combined into a single shortest path between the two crosses.

Obviously if the circle is not on the shortest path between the two crosses then this strategy will not work, but it is useful for reducing the number of different paths to check. In general we probably will not have access to such a wizard, but we can still take advantage of this fact by looking for the shortest path between the first cross and every point along the way to the second cross. The shortest path to the second cross must begin with one of these paths. But how can we find the shortest path to these intermediate points? But using intermediate points that are along the way to these intermediate points of course. This strategy allows us to find the shortest path between two points bit by bit, by finding other shorter shortest paths along the way.

This is how Dijkstra's algorithm works. Here is the algorithm:

distance_from_start = { starting_point : 0 }
previous_point = { starting_point : null }
next_points = [ starting_point ]

while next_points is not empty do
    next_point = point in next_points where distance_from_start[point] is minimum
    remove next_point from next_points

    for each neighbouring_point around next_point do
        new_distance = distance_from_start[next_point] + distance from next_point to neighbouring_point
        if neighbouring_point is not in distance_from_start or new_distance < distance_from_start[neighbouring_point] then
            add neighbouring_point to next_points if not in distance_from_start

            distance_from_start[neighbouring_point] = new_distance
            previous_point[neighbouring_point] = next_point

The algorithm works by starting at the starting point and checking every neighbouring point around it. The path from the starting point to each one of these points is considered to be the shortest path to the points. Each of these points is checked for its neighbouring points also and if the point is a previously visited point, then the distance of the new path found is compared to the distance of the previous path that was used. If the new distance is shorter than the previous one then the new path is used instead of the old one. This is repeated until every point has been checked. At the end, every point is associated with some previous point which if followed to the starting point will form the shortest path to the starting point.

Uniform Cost Search

Uniform Cost Search is Dijkstra's Algorithm which is focused on finding a single shortest path to a single finishing point rather than a shortest path to every point. It does this by stopping as soon as the finishing point is found.

Here is the algorithm:

distance_from_start = { starting_point : 0 }
previous_point = { starting_point : null }
next_points = [ starting_point ]

while next_points is not empty do
    next_point = point in next_points where distance_from_start[point] is minimum
    if next_point = finishing_point then exit
    remove next_point from next_points

    for each neighbouring_point around next_point do
        new_distance = distance_from_start[next_point] + distance from next_point to neighbouring_point
        if neighbouring_point is not in distance_from_start or new_distance < distance_from_start[neighbouring_point] then
            add neighbouring_point to next_points if not in distance_from_start

            distance_from_start[neighbouring_point] = new_distance
            previous_point[neighbouring_point] = next_point

It is not enough to check if a neighbouring point is the finishing point to exit as there might be other points which are the finishing point's neighbours that are better. In order to be sure that no other neighbouring point is better, you have to wait until it is the closest point available in the list of next points to explore. In this way you can be sure that any other point that can lead to the finishing point is further away from the starting point (and thus have a longer path) than the finishing node is.

The problem with this approach is that although it is guaranteed to give you the shortest path from the starting point to the finishing point, it does so because it checks every possible path blindly, which would be too long for practical use in a game. A faster algorithm exists.

A* Search

A* Search is the informed version of Uniform Cost Search. It differs in that you have to give it a way to estimate how close any point is to the finishing point which it will use to make informed decisions on which point it should follow next. This takes the "blindly" part out of the Uniform Cost Search. The estimate is used in order to estimate what the total distance of the path from the starting point to the finishing point is, provided that a particular point is included in the path. The way the path distance is estimated is by adding together the distance of the point from the starting point (which is done as usual) and the estimated distance from the point to the finishing point. In this way, the selection of the next point to explore is not made just on the point's distance from the starting point, but on the whole path length.

The important thing is that this estimate is never more than the actual distance to the finishing point, although it can be shorter. If it is longer then it will not give the shortest path. The reason for this is because as soon as an actual path is found (when the finishing point is discovered), if the actual distance of this path is shorter than the underestimated distances of the other paths, then it can be safely assumed that the others paths actual distances will be more than the actual path found. On the other hand, if the estimated distances are overestimated, then they cannot be compared to an actual path distance since their actual distance might be less.

Another thing to notice is that the closer the estimate is to zero, the more the search will behave like Uniform Cost Search and the more likely the path given is to be optimal. A good estimate would be the Euclidian distance (straight line distance) to the finishing point. This estimate is used to select the most promising point to process from all the ones discovered.

Here is the algorithm:

distance_from_start = { starting_point : 0 }
previous_point = { starting_point : null }
next_points = [ starting_point ]

while next_points is not empty do
>>> next_point = point in next_points where distance_from_start[point] + estimate to finishing point is minimum
    if next_point = finishing_point then exit
    remove next_point from next_points

    for each neighbouring_point around next_point do
        new_distance = distance_from_start[next_point] + distance from next_point to neighbouring_point
        if neighbouring_point is not in distance_from_start or new_distance < distance_from_start[neighbouring_point] then
            add neighbouring_point to next_points if not in distance_from_start

            distance_from_start[neighbouring_point] = new_distance
            previous_point[neighbouring_point] = next_point