Найдите общее кратное простых степеней (2, 3, 5, 7) › N и минимизируйте показатели

Причина, по которой я это делаю, состоит в том, чтобы найти ближайшее число, превышающее N, которое является общим кратным степеней простых чисел, чтобы иметь возможность использовать FFTW.

Насколько я понимаю, это проблема оптимизации/линейного программирования. И мне удалось сформулировать математически:

minimize(a log(2) + b log(3) + c log(5) + ...)

Constraints::
  a log(2) + b log(3) + c log(5) + ..... >= log (N)  (known parameter: N)
  {a, b, c ...} >= 0      (i.e. are positive)
  {a, b, c ...} % 1 == 0  (i.e. are integers)

где a, b, c — показатели степени и целые числа.

Я хотел бы реализовать то же самое численно. Предпочтительно использовать scipy.optimize.

EDIT: я немного изменил уравнения на основе комментариев. Даже если реализации не существует, алгоритм был бы полезен.


python scipy optimization linear-programming integer-programming
person jadelord    schedule 16.03.2018    source источник
comment
Это широко с небольшим контекстом. То, что вы показали, представляет собой целочисленную программу, которая не поддерживается в scipy.   -  person sascha    schedule 16.03.2018
comment
Я записал математические уравнения. Насколько конкретнее это может быть? И я сказал предпочтительно scipy, другой пакет или даже алгоритм был бы полезен.   -  person jadelord    schedule 16.03.2018
comment
Контекст всегда помогает. Здесь нет явной неотрицательности (хотя вам, наверное, этого хочется), > в LP вообще нет (и описывать его нужно осторожно; если › это действительно то, что вы хотите; в отличие от >=) и N нам неизвестно. Кроме того: что дало вам ваше исследование после прочтения того, что вы описали, является целочисленной программой?   -  person sascha    schedule 16.03.2018
comment
Моя ошибка >= на самом деле в порядке. N — известный параметр. Термин integer programming для меня новый, и в Википедии говорится, что он NP-полный. То есть точных алгоритмов нет, есть только приблизительные?   -  person jadelord    schedule 16.03.2018
comment
NP-полнота не означает неточных результатов.   -  person sascha    schedule 16.03.2018
comment
Лучше задать этот вопрос на scicomp.stackexchange.com или на одном из математических форумов.   -  person Bill Bell    schedule 16.03.2018
comment
Для тех, кто интересуется более простым случаем поиска следующего составного числа только с делителями 2, 3 и 5, есть scipy.fftpack.next_fast_len.   -  person Warren Weckesser    schedule 17.03.2018
comment
Хороший вопрос @WarrenWeckesser! По-видимому, существует также pyfftw.next_fast_len реализация, но нам придется подождать до следующего релиза. .   -  person jadelord    schedule 17.03.2018

Ответы (1)


arrow_upward
1
arrow_downward

Согласно @sascha и разговору в списке рассылки Scipy-Dev scipy еще не реализовал решатель смешанного целочисленного линейного программирования (MILP, также известный как MIP).

Вики Python содержит список пакетов Python LP/MILP. В итоге я использовал pulp и его встроенный решатель. Изменив код из похожей задачи и примеры найдены в документе, решение выглядит следующим образом.

#!/usr/bin/env python
"""Find the closest multiple of prime powers greater than N."""
from pulp import *
import numpy as np


N = 1000
bases = [2, 3, 5, 7]
debug = True

prob = LpProblem("FFTW Gridsize Problem", LpMinimize)

exp_max = np.ceil(np.log2(N))
exps = LpVariable.dicts('exponent', bases, 0, exp_max, cat=LpInteger)
log_N_new = LpVariable("log_grid_size", 0)
log_N_new = lpSum(exp * np.log(base) for exp, base in zip(exps.values(), bases))

prob += log_N_new  # Target to be minimized
# Subject to:
prob += log_N_new >= np.log(N), "T1"

if debug:
    print('bases =', bases)
    print('exponents =', exps)
    print('log_N_new =', log_N_new)
    prob.writeLP("OptimizationModel.lp")

prob.solve()

if debug:
    print("Status:", LpStatus[prob.status])
    for v in prob.variables():
        print(v.name, "=", v.varValue)

N_new = np.prod(
    [base**v.varValue for base, v in zip(bases, prob.variables())])
print(N_new)
person jadelord    schedule 17.03.2018
comment
Рад, что вы нашли это, и я надеюсь, что это работает для вас. Pulcle проста в использовании и легко упаковывает лучший (имхо) MIP-решатель с открытым исходным кодом, Coin-OR Cbc (если вы можете жить с его далеко не оптимальными документами; по сравнению с GLPK и другими). - person sascha; 17.03.2018