JMK no matter what

Linear Fractional Programming

After ACM ICPC Seoul Regional 2009, I hear problem D could easily modeled as a Linear Fractional Programming problem; lewha0 pointed me to the solution of this kind of problem before - here's a brief review of the solution.

Problem

Given a sequence of 0 and 1's, return a continuous subsequence with maximal average. The subsequence must not be shorter than L.

Solution

I am aware of two approaches.

  1. Try every length up to 2L: given a subsequence of length 2L, by pigeonhole principle, one of the two L-length halves has the same or better average.
  2. Generalize the problem and use LFP.

LFP can solve this generalized version of the problem: Given A[] and B[], return max interval (i, j) where (A[i]+..+A[j])/(B[i]+..+B[j]). In this problem, A[] is either 0 or 1, and B[] is always 1.

Let's start with a guess; will the final solution be smaller than x? Let us assert x is an upper bound of sum(A[i..j])/sum(B[i..j]), for all possible interval S,


x > \frac{\sum_{i\in S}{A_i}}{\sum_{i \in S}{B_i}}

Rearranging, we get


\sum{xB_i}-\sum{A_i} > 0

So x is an upper bound if and only if for every possible interval S, the above equation is positive. We can easily check this by calculating the vector x\vec{B}-\vec{A} and find the interval with minimum sum - and see if it's less than zero. If it is greater than zero, x is indeed an upper bound.

Now, we can easily do binary searching over possible ranges of x. source code follows

lang:py
from random import randint
def brute(a,b,l):
    ret = 0
    for i in range(len(a)):
        asum, bsum = 0, 0
        for j in range(i,len(a)):
            asum += a[j]
            bsum += b[j]
            if j-i+1 >= l:
                ret = max(ret, float(asum) / bsum)
    return ret

def lfp(a,b,l):
    lo, hi = 0.0, 100*1000
    for iter in range(100):
        mid = (lo + hi) / 2.0
        c = [b[i] * mid - a[i] for i in range(len(a))]
        psum = [0] * len(a)
        psum[0] = c[0]
        for i in range(1,len(a)):
            psum[i] = psum[i-1] + c[i]
        mx = 0
        subret = 1e200
        for i in range(l-1,len(a)):
            # min range ending at i
            subret = min(subret, psum[i] - mx)
            mx = max(mx, psum[i-l+1])

        if subret <= 0:
            lo = mid
        else:
            hi = mid
    return lo

for iter in xrange(1000):
    N = 1000
    a = [randint(1,100) for i in range(N)]
    b = [randint(50,100) for i in range(N)]
    l = randint(5,N)
    print brute(a,b,l) - lfp(a,b,l)
2009-11-07 02:23:39 | JM | /logs/library/ | 1 Comments
ltdtl
2009-11-10 05:43:55
아 파이선 코드만 보면 눈이 아픕니다.
왜 그럴까요?

Leave a comment