ratio_algorithm.py

import random
import itertools
import numpy as np

# counter example:
# a = np.array([25,  1, 46, 25, 45,  6, 25, 43,  9,  3, 21, 34, 41, 32, 36, 12, 33,
#         7, 42, 50, 22, 37, 16, 31,  3, 44, 28,  8,  4, 41, 15, 19, 29, 10,
#                31, 14,  2, 17, 31, 16,  9, 15, 15, 40, 37, 14, 10,  2, 19,  9])
# b = np.array([ 4, 37, 43, 11, 26, 17, 35, 42, 47, 31, 49, 39, 47, 39, 37, 15, 35,
#        26, 25, 43, 47, 49, 19, 42, 28, 49, 18, 19, 48, 35,  1, 15, 40, 15,
#               26, 14,  1,  9, 23, 15, 24, 34, 47, 15, 30, 38, 36, 49, 36, 12])


def maxsum_of_ratio(a,b,k):
    c = a/b
    sorted_c = np.argsort(-c)
    max_indexset = sorted_c[range(k)]
    ratio = (sum(a[max_indexset])+0.0)/sum(b[max_indexset])
    return ratio, max_indexset

def exhaustive(a,b,k):
    index = range(len(a))
    max = 0
    max_indexset = []
    for i in itertools.combinations(index, k):
        sigma_a = sum([a[j] for j in list(i)])
        sigma_b = sum([b[j] for j in list(i)])
        ratio = (sigma_a+0.0)/sigma_b
        if ratio > max:
            max = ratio
            max_indexset = list(i)
    return max, max_indexset

def alg(a,b,k):
    index = range(len(a))   
    sample = random.sample(index, k)

    lambda_prev = 0
    lambda_ = (sum(a[sample])+0.0)/sum(b[sample])
    c = a - lambda_* b
    sorted_c = np.argsort(-c)
   
    iters = 0
    while lambda_ > lambda_prev:
        iters += 1
        sample = sorted_c[range(k)]
        lambda_prev = lambda_
        lambda_ = (sum(a[sample])+0.0)/sum(b[sample])
        c = a - lambda_* b
        sorted_c = np.argsort(-c)

    sample = sorted_c[range(k)]
    return lambda_, sample, iters


# Algorithm with epsilon
def until_eps(b, indices, eps):
    sum_b = 0
    full_indices = []
    portion_index = (0,0)
    for i in indices:
        sum_b += b[i]
        full_indices.append(i)
        if sum_b == eps:
            break
        elif sum_b > eps:
            half_index = full_indices.pop()
            portion = (eps - sum(b[full_indices]) + 0.0)/b[half_index]
            portion_index = (portion, half_index)
            break

    return full_indices, portion_index

def alg_eps(a, b, eps):
    index = range(len(a))
    full_indices, portion_index = until_eps(b, index, eps)

    lambda_prev = 0
    sum_a = sum(a[full_indices]) + portion_index[0]*a[portion_index[1]]
    lambda_ = (sum_a+0.0)/eps
    c = a - lambda_* b

    iters = 0
    while lambda_ > lambda_prev:
        iters += 1
        indices = np.argsort(-c)
        full_indices, portion_index = until_eps(b, indices, eps)
        lambda_prev = lambda_
        sum_a = sum(a[full_indices]) + portion_index[0]*a[portion_index[1]]
        lambda_ = (sum_a+0.0)/eps
        c = a - lambda_* b
        sorted_c = np.argsort(-c)

    indices = [] + full_indices
    indices.append(portion_index[1])
    portion = portion_index[0]
    return lambda_, indices, portion, iters


if __name__ == '__main__':
    R = 50
    N = 100
    a = np.array([random.randint(1,R) for i in range(N)])
    b = np.array([random.randint(1,R) for i in range(N)])

    print a, b
    max, max_index =  exhaustive(a,b,3)
    print max, a[max_index], b[max_index]

    max, max_index, iters = alg(a,b,3)
    print max, a[max_index], b[max_index]

    print iters