Monday, January 31, 2011

Binary Indexed Tree.

Let us discuss an interesting problem today related with Cumulative frequency. Consider a hotel database which stores timings and number of phone-calls arrived during the interval.

Time => Number of calls

10 - 12 => 30

12 - 14 => 5

14 - 16 => 12

16 - 18 => 25

...

In this example, cumulative sum of calls made before 18 means Number of calls made in the interval 16-18 + (14-16) ... (10-12). It could be considered as a Range-sum. Understanding of this problem is very essential to realize the need for Binary indexed Tree.

Let us assume that this information had been stored in an array, which has interval mapped as index, and number of calls as values. Mapping from time interval to array index is not very important for this discussion.

Now the problem is to find a suitable data structure which could make updates, insertions and queries on this information in an efficient way. Let us think of naive solutions first.

Recursion Approach:

If we have to find out number of calls made before 18, then we will have to add up the values prior to that interval.

Mathematically, CumulativeSumTill(n) = valueAt(n) + CumulativeSumTill(n-1), a recursive relation. We know this could be implemented with recursion, with complexity O(n), as it needs traversal of all elements occur before the given element. Finding cumulative sum at an index is not efficient in this way.

If we have to add up a new interval information, we could easily add in source array in constant time as addition would generally happen at the end of the array. So with updating any interval's info also. So, updates happen at constant time, very efficient.

Iterative or Dynamic Programming Approach:

With Dynamic programming, iterative approach, if we could memoize the intermediate results in a separate table, we could make complexity of Find or Query, O(1), a constant. In the below given code, we had memoized the cumulative sum results in a separate table. So, getting a CS at an index is very efficient. But think of an update at an index '0', which would render all the stored results invalid. Every updates would be of complexity O(n), very inefficient.

void MemoizeCumulativeSum()
{
CumulativeSum[0] = SourceArray[0];

for( int i = 1; i < lengthOfSource; i++)
{
CumulativeSum[i] = CumulativeSum[i-1] + SourceArray[i];
}
}

Summary of Naive Solution:
Dynamic programming approach or storing intermediate results approach makes query fast, but updates slow. On the fly calculation, recursion based approach, makes updates fast, but query slow. This is a perfect instance of Space-Time trade-off. So, in order to have both query and update happen in an efficient way, we need to think of a different data-structure which combines above-said two approaches in a perfect way.

Binary Indexed Tree:
We are going to discuss a data structure which helps to query and update in logarithmic time.

Binary indexed tree or Fenwick tree, could be realized using an array which holds cumulative sum or sub-cumulative sum as its values. Every index of this array would hold partial cumulative sum or complete cumulative sum according to the BIT algorithm. This array makes query and update happen in logarithm time as described below.

Basics:

Any number could be expressed in 2's powers. So, number of digits of any number in 2's power form would be log n of base 2. Presto, this is what exactly the property being made use in BIT.

Initially let us assume an auxiliary table as in iterative approach, contains all zeroes. While updating some value in source array, we need to make sure all the indices in the auxiliary table which contains partial or cumulative sum with that value, updated. This may look like naive iterative approach; but the beauty of this data-structure is that it doesn't need to update all the values as in naive memoization. So, this could be considered as an educated or intelligent memoization. On the average, it may need log n updates in auxiliary table which makes update efficient. Same way, query needs sum of certain values in auxiliary table unlike naive recursion approach, which needs to process all the values in source array. BIT needs to process only log n values for a query. So, query and updates are efficient in this data-structure.

BIT Core Algorithm:

A complete implementation of BIT which supports update and query, has been given at the end of this blog.

While the start of the algorithm, prepare an auxiliary array of source array's size and initialize all the elements with zero. Now loop through all the elements in the source array and update the auxiliary table according to below code. Below algorithm calculates the indices to be updated, as specified by the logic, nextIndex += (nextIndex & -nextIndex). ANDing a value and its negated isolates a non-zero-least significant bit. If an index is "xxx1000" then, index & ~index would result in "0001000". Now the new index is "xxx1000" + "0001000" would unset the least significant non-zero bit in "xxx1000". Hence this algorithm, in the worst case would update log2 n indices in the auxiliary table.


Below algorithm would be used to reflect any modifications made in source array; if value in some index had been incremented in the source array, this must be called in order to keep the BIT updated. The same argument with regards to complexity of initialization holds good for updates also. If you couldn't get a hold of algorithm, no worries. We will discuss the algorithm intuitively later.


Same way, query algorithm is as follows:



Intuitive Analysis:

Actually, after all updates, BIT or auxiliary table will contain either partial or complete cumulative frequency depends on the index.

Indices which are powers of two contain complete cumulative frequency of that index. Others would contain partial.

'6' can be written as sum of powers of two as 4 + 2. Now take the value at 4, value at index 2 starting from 4 and add up to get Cumulative Frequency of 6.

Same way '7' can be written as 4 + 2 + 1. Take value at 4, value at index 4+2 and value at index 4+2+1 and add up those. Intuitively, this is what exactly happens in Query algorithm.

Same could be applied for Update algorithm also. Only complexity of this algorithm is to understand the procedure. :)

BIT Class implementation:

class BinaryIndexedTree
{
public:
BinaryIndexedTree(vector<int> inputArray)
{
internalArray = NULL;
InitializeBIT(inputArray);
}
void IncrementValue(int value, int index)
{
int indexToModify = index + 1;

while(indexToModify < arraySize)
{
internalArray[indexToModify] += value;
indexToModify += (indexToModify & -indexToModify);
}
}
int Query(int index)
{
int indexToRetrieve = index + 1;
int retValue = 0;
while(indexToRetrieve)
{
retValue += internalArray[indexToRetrieve];
indexToRetrieve -= (indexToRetrieve & -indexToRetrieve);
}
return retValue;
}
private:
void InitializeBIT(vector<int> inputArray)
{
// Initialize internal array;
// Zeroth index is Sentinel.
arraySize = inputArray.size() + 1;
internalArray = new int[arraySize];
for(int i = 0; i < arraySize; i++)
internalArray[i] = 0;
for(int i = 1; i < arraySize; i++)
{
int valueToBeAdded = inputArray[i - 1];
internalArray[i] += valueToBeAdded;
int k = i;
while( k < arraySize)
{
k += (k & -k);
internalArray[k] += valueToBeAdded;
}
}
}
int *internalArray;
int arraySize;
};

Some more links for this topic



Thanks for Reading.

Sunday, January 30, 2011

Random Samples from a huge set - Intuitive reservoir sampling

Consider a set of cardinality N. We are given a problem of collecting 'n' random samples from this set. This can be easily done with standard random number generators.

Please review my another blog, if you are not comfortable with Fisher-Yates algorithm for Random shuffling. http://karthikpresumes.blogspot.com/2011/06/random-shuffling-intuitive-fisher-yates.html

Now a bit more complex sampling. What if the number 'N' is not known a priori. So, we couldn't generate random numbers, as we don't know the range.
Interestingly, what if the parent set is keep growing; N is time-varying.

Two Pass algorithm:

Let us consider a situation where 'N' is not varying, but not known a priori. In this case, N could be calculated with a single pass of the set. After this pass, problem boils down to initial one, we discussed. Cool, but number of passes needed is two.

Single Pass Reservoir Sampling algorithm:

Is there any way to make it one pass? Yes. Reservoir sampling.


In this post, I will try to explain this algorithm in an intuitive way.

Base case:

Let us assume the samples set as reservoir of 'n' numbers. If n is equal to N, problem is base case, which is trivial. Reservoir in this case would have all the samples, perhaps in a different order.

Simple case:

Now, assume N = n+1, we need n samples out of n + 1.

Initially fill the reservoir with first n elements. After processing element at index n+1, reservoir's elements should be having the probability n / (n + 1) of being placed in reservoir. It means, n + 1th element would be placed with probability n / (n+1) . This could be easily achieved by generating a random number between 1 .. n + 1. If generated one is above n, then place the element in reservoir.

Placing an element in filled reservoir, needs some element else gets removed. Please remember, after all, probability of elements being in Reservoir should be n /N.

Probability of an element, x in reservoir to get removed is selecting n + 1 th element for being in Reservoir and selecting the element x for replacing; which is n / (n + 1) * 1 / n => 1 / (n + 1). So, probability of staying of x is n / (n + 1), 1 - 1 /(n + 1), which is the required probability. This could be easily generalized for any index above n.

Generalization:

For general index i, probability of being placed in Reservoir is n / i. Probability of just removing a random element from Reservoir is 1 / n. So, combined probability of a random element to get replaced is 1 / i. Probability of the same element to stay in Reservoir is (i - 1) / i, which was there in Reservoir with Probability n / (i - 1). So, combined probability is n / i.

Luckily, this will work for time-varying 'n' also.

Thanks for reading.

Thursday, January 27, 2011

Dealing with Really big numbers.

Yesterday evening, I was trying to solve a problem which deals with multiplying really huge integers. As we know, in most of the programming languages integer data type has limit on the size of integer. Generally, it depends on the processor word size. In 32 bit machines, it would be 32 bits.
So, biggest unsigned integer that could be held in, is 4294967296, just 10 digits. What if we have to deal with 100 digit numbers???

Problem: Find the sum of natural numbers till N, N will be really huge, say 80 digits.

Obviously, we couldn't work with basic data type Integer. We can implement our own Abstract Data Type for this. Still, choosing perfect algorithm for processing is difficult.

There are some languages, which do support big integers naturally, like Haskell, Pike. Some languages like Java, supports big integers using different data type rather than usual Integer type. Haskell's limit for integer depends only on Main memory. Pike also supports integers of huge size.

I evaluated Pike yesterday. Certain characteristics of Pike are really interesting. Syntax in Pike is almost equivalent to C /C++. It is garbage collected and Object oriented language. I am just showing a Pike program for the above problem here.

int main()
{
string countStr; int count;
countStr = Stdio.stdin->gets();
sscanf(countStr, "%d", count);
while( count--)
{
int realInt;
string strInt;
strInt = Stdio.stdin->gets();
sscanf(strInt, "%d", realInt);
int res = realInt * (realInt + 1) / 2;
string outStr = sprintf("%d\n", res);
write(outStr);
}

return 0;

}

You could find there is no much difference in the syntax, if you are a C/ C++ programmer. It is very fast in handling really big integers. I tested and played with some other library functions in Pike. Hope you will also try and play with this innovative language.

In next blog, I will explain my experiences with Haskell.

Thanks for reading.

Tuesday, January 25, 2011

Storing partial results, Dynamic Programming

Today, I tried to resolve a programming puzzle; which is a bit different than usual problems. Generally a problem like Finding fibonacci series Nth number, has two variants of solutions, Recursive and Iterative. Recursive solutions needs solution for sub-problem first and makes up the complete solution in bottom-up manner.

Fib(n) = Fib(n - 1) + Fib( n -2 )

Ignoring base cases, this would be complete recursive solution. Calculating in this way, would be easier for us to implement, but computation complexity would be huge. For instance, not only Fib(n) calculates Fib(n-2), Fib(n-1) also calculates the same. Same calculations would recur several times. Cause of this is not storing the computed result. For this kind of problems, obvious way to efficiently solve is to store intermediate results. This leads to Iterative solution.

Fib[n] = Fib[n - 1] + Fib[n - 2]
Note: Array is being accessed, not any more recursive function calls. All calculations happen only once.

Computation complexity comes down to O(n). This is one way of viewing Space Time complexity trade off. We successfully traded off space for Time.

There are some problems, don't need all partial results. In that case calculating all the intermediate values would be waste of time. Think about space complexity also.

One such example, taken from CodeChef

One strange bank will provide you three possible coins for a coin of denomination $N, as follows:

$N/ 2, $N/ 3, $N/ 4. N is an integer.

Changing coin would be profitable, for denominations like $12, $24. But for 3, this would cause a loss of $1. (3 /2 + 3 / 3 + 3/4) = $2.

What would be the maximum profit possible with $X.

We can't solve this problem iteratively as we do need result for X, not for 1, 2 ... X. Since X can be considerably big, we can't store all intermediate values also. At last, there is no need for all intermediate calculations for calculating Xth value.

This suggests us to go for Recursion based one. But we know storing intermediate results would help to avoid same calculations several times. So, storing only a part of intermediate results in a Map, or a Hash Table or a Simple array would be very useful.