A sample by energy
Seismologists want to examine by hand k earthquakes from a catalogue, chosen at random but so that strong ones are picked more often than weak ones, in proportion to the energy released. The energy of an earthquake of magnitude M is proportional to 101.5·M: about 31.6 times more per unit of magnitude. The catalogue arrives as a stream, and storing it in full is undesirable.
They decided to use the Efraimidis–Spirakis weighted reservoir sampling. An event with weight w = 101.5·M and a random number u from the interval (0, 1) is given the key u1/w, and the sample consists of the k events with the largest keys. So that the result can be checked, the random numbers are already recorded in the catalogue next to the magnitudes.
Print the numbers of the chosen events.
Input
The first line holds n and k (1 ≤ k ≤ n ≤ 6·104). Then n lines: the magnitude M (from 5.00 to 9.50, two digits after the point) and the number u (0 < u < 1, nine digits after the point). Events are numbered from one in input order. It is guaranteed that the k-th and (k + 1)-th keys in decreasing order differ noticeably, and the answer does not depend on rounding errors if computed carefully.
Output
The numbers of the k chosen events in increasing order, separated by spaces.
Examples
Example 1
4 2 5.00 0.500000000 6.00 0.500000000 5.50 0.900000000 5.00 0.990000000
3 4
Example 2
3 1 5.00 0.999999999 7.00 0.000000001 6.00 0.500000000
1
Limits
- Time · Python
- 4000 ms
- Time · Go
- 1000 ms
- Memory
- 256 MB
- Judging
- 10 tests
Walkthrough
Open it step by step: each step gives away a little more. You can jump straight to the solution, but a problem solved with one hint stays with you longer.
1Hint 1
Compute u1/w for an earthquake of magnitude 9.5 and u = 0.5. How many decimal places are needed to tell it from one, and from the key of the same event with u = 0.500001?
2Hint 2
You can compare not the keys themselves but any increasing or decreasing function of them. Which function turns u1/w into a sum of two ordinary numbers?
3The idea
The key u1/w with M = 8 and u = 0.5 equals 1 − 6.9·10−13, and with M = 9.5, 1 − 3.9·10−15: for the strongest events the keys differ from one and from each other in the last bits of a floating-point number, and their order would be decided by rounding. You can compare not the keys but their logarithms: ln u1/w = ln u / w. These are small negative numbers, but their relative precision is the same as for any others. It is more convenient to take the logarithm once more, of the negation: ln(−ln u) − ln w = ln(−ln u) − 1.5·M·ln 10, and then the weight 101.5·M need not be computed at all. This quantity decreases as the key grows, so we need the k events with its smallest values. Keep them in a heap of size k, as in any “k best from a stream” problem.
4Walkthrough
The heap keeps the k best events, with the worst of them on top. A new event displaces it if it is better. After the pass, sort the numbers.
The first example shows that the randomness here is real. An M6.0 event with u = 0.5 is 31.6 times stronger than an M5.0 event, but loses to it when the latter has u = 0.99. The keys are compared via ln u/w: −0.693/31.6 = −0.0219 against −0.0101 (both in units of the weak event’s 1/w). In the second example an M5.0 event with u almost equal to one beats an M7.0 event that drew u = 10−9. On the real catalogue the sample consists almost entirely of strong events, because energies differ thousands of times, but not only of the very strongest.
Why the key u1/w does what is needed. For first place, the event with the largest key, it can be checked directly: the probability that event i takes it is wi/(w1 + … + wn), exactly in proportion to weight. Then the same repeats among the rest: that is how weighted sampling without replacement works, and Efraimidis and Spirakis proved this in 2006. Unlike the key “u times w”, which seems natural but does not give proportionality.
4+Solution
5Complexity
O(n log k) time and O(k) memory: one pass and a heap of size k.
6Common mistakes
- The key u·w. It too favours heavy events, but not in proportion to weight. With weights differing thousands of times, it simply picks the strongest earthquakes: the first example.
- Computing u^(1/w) directly. The keys of strong events differ from each other in the last bits, and the choice between them becomes a matter of rounding.
- Sorting the whole stream. Correct, but O(n log n) and the whole sample in memory; a heap of k elements manages in one pass.