Дано целое число n и 2 действительные последовательности { a_1 , ..., a_n } и { b_1 ,. .., b_n }, с a_i , b_i > 0, для всех i . Для данного фиксированного m <<em> n пусть { P_1 , ..., P_m } будет разделом множества {1, ..., n } как в P_1 U ... U P_n = {1, ..., n }, с попарно непересекающимся P_i (пустое пересечение). Я sh найду раздел размером м , который максимизирует выражение
$$ \ space \ max_ {P = \ {P_1, ..., P_n \}} \ sum_ {j = 1} ^ {m} \ frac {\ sum_ {i \ in P_j} a_i} {\ sum_ {i \ in P_j} b_i} $$
Количество разделов в наборе n выберите m , слишком большое, чтобы сделать его методом грубой силы. Есть ли итеративное или приблизительное решение, которое лучше?
Для понимания этой проблемы блок кода в конце решается с помощью грубой силы. Для проблем с разрешением c ( n ~ 1e6, k ~ 20) он непригоден как есть, но легко распространяется.
Редактировать : Предварительная сортировка a , b по значениям a ^ 2 / b всегда дает увеличивающиеся индексы секций:
a = rng.uniform(low=0.0, high=10.0, size=NUM_POINTS)
b = rng.uniform(low=0.0, high=10.0, size=NUM_POINTS)
ind = np.argsort(a/b)
(a,b) = (seq[ind] for seq in (a,b))
примерный прогон с
NUM_POINTS = 16
PARTITION_SIZE = 3
дает оптимальное разбиение
[[0, 1, 2, 3, 4, 5, 6, 7], [8, 9], [10, 11]]
, что является монотонным c в индексах. Я думаю, что могу доказать это. Если это так, поиск грубой силы может быть улучшен до n , выберите k -1 раз, все еще долго, но значительную экономию.
import numpy as np
import multiprocessing
import concurrent.futures
from functools import partial
from itertools import islice
rng = np.random.RandomState(55)
def knuth_partition(ns, m):
def visit(n, a):
ps = [[] for i in range(m)]
for j in range(n):
ps[a[j + 1]].append(ns[j])
return ps
def f(mu, nu, sigma, n, a):
if mu == 2:
yield visit(n, a)
else:
for v in f(mu - 1, nu - 1, (mu + sigma) % 2, n, a):
yield v
if nu == mu + 1:
a[mu] = mu - 1
yield visit(n, a)
while a[nu] > 0:
a[nu] = a[nu] - 1
yield visit(n, a)
elif nu > mu + 1:
if (mu + sigma) % 2 == 1:
a[nu - 1] = mu - 1
else:
a[mu] = mu - 1
if (a[nu] + sigma) % 2 == 1:
for v in b(mu, nu - 1, 0, n, a):
yield v
else:
for v in f(mu, nu - 1, 0, n, a):
yield v
while a[nu] > 0:
a[nu] = a[nu] - 1
if (a[nu] + sigma) % 2 == 1:
for v in b(mu, nu - 1, 0, n, a):
yield v
else:
for v in f(mu, nu - 1, 0, n, a):
yield v
def b(mu, nu, sigma, n, a):
if nu == mu + 1:
while a[nu] < mu - 1:
yield visit(n, a)
a[nu] = a[nu] + 1
yield visit(n, a)
a[mu] = 0
elif nu > mu + 1:
if (a[nu] + sigma) % 2 == 1:
for v in f(mu, nu - 1, 0, n, a):
yield v
else:
for v in b(mu, nu - 1, 0, n, a):
yield v
while a[nu] < mu - 1:
a[nu] = a[nu] + 1
if (a[nu] + sigma) % 2 == 1:
for v in f(mu, nu - 1, 0, n, a):
yield v
else:
for v in b(mu, nu - 1, 0, n, a):
yield v
if (mu + sigma) % 2 == 1:
a[nu - 1] = 0
else:
a[mu] = 0
if mu == 2:
yield visit(n, a)
else:
for v in b(mu - 1, nu - 1, (mu + sigma) % 2, n, a):
yield v
n = len(ns)
a = [0] * (n + 1)
for j in range(1, m + 1):
a[n - m + j] = j - 1
return f(m, n, 0, n, a)
def Bell_n_k(n, k):
''' Number of partitions of {1,...,n} into
k subsets, a restricted Bell number
'''
if (n == 0 or k == 0 or k > n):
return 0
if (k == 1 or k == n):
return 1
return (k * Bell_n_k(n - 1, k) +
Bell_n_k(n - 1, k - 1))
NUM_POINTS = 13
PARTITION_SIZE = 4
NUM_WORKERS = multiprocessing.cpu_count()
INT_LIST= range(0, NUM_POINTS)
REPORT_EACH = 10000
partitions = knuth_partition(INT_LIST, PARTITION_SIZE)
# Theoretical number of partitions, for accurate
# division of labor
num_partitions = Bell_n_k(NUM_POINTS, PARTITION_SIZE)
bin_ends = list(range(0,num_partitions,int(num_partitions/NUM_WORKERS)))
bin_ends = bin_ends + [num_partitions] if num_partitions/NUM_WORKERS else bin_ends
islice_on = list(zip(bin_ends[:-1], bin_ends[1:]))
# Have to consume it; can't split work on generator
partitions = list(partitions)
rng.shuffle(partitions)
slices = [list(islice(partitions, *ind)) for ind in islice_on]
return_values = [None] * len(slices)
futures = [None] * len(slices)
a = rng.uniform(low=0.0, high=10.0, size=NUM_POINTS)
b = rng.uniform(low=0.0, high=10.0, size=NUM_POINTS)
ind = np.argsort(a/b)
(a,b) = (seq[ind] for seq in (a,b))
def start_task():
print('Starting ', multiprocessing.current_process().name)
def _task(a, b, partitions, report_each=REPORT_EACH):
max_sum = float('-inf')
arg_max = -1
for ind,part in enumerate(partitions):
val = 0
for p in part:
val += sum(a[p])**2/sum(b[p])
if val > max_sum:
max_sum = val
arg_max = part
if not ind%report_each:
print('Percent complete: {:.{prec}f}'.
format(100*len(slices)*ind/num_partitions, prec=2))
return (max_sum, arg_max)
def reduce(return_values):
return max(return_values, key=lambda x: x[0])
task = partial(_task, a, b)
with concurrent.futures.ThreadPoolExecutor() as executor:
for ind,slice in enumerate(slices):
futures[ind] = executor.submit(task, slice)
return_values[ind] = futures[ind].result()
reduce(return_values)