-2

The reason I am doing this is to find a closest number higher than N, which is a common multiple of prime powers --- to be able to use FFTW.

As I understand this is an optimization/linear programming problem. And I managed to formulate the mathematically:

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)

where a, b, c are the exponents and are integers.

I would like to implement the same numerically. Preferably using scipy.optimize.

EDIT: I have modified the equations slightly based on comments. Even if an implementation does not exist, an algorithm would be helpful.

jadelord
  • 1,511
  • 14
  • 19
  • 1
    This is broad with little context. What you have shown is an integer-program and that's not supported in scipy. – sascha Mar 16 '18 at 16:55
  • I have written down the mathematical equations. How much more specific can it be? And I said *preferably* scipy, another package or even an algorithm would be helpful. – jadelord Mar 16 '18 at 17:57
  • 1
    Context always helps. There is no explicit nonnegativity here (although you probably want it), there is no ```>``` in LPs in general (and describing it needs care; if > is really what you want; opposed to ```>=```) and N is unknown to us. Also: what did your research give you after reading that what you described is an integer-program? – sascha Mar 16 '18 at 18:05
  • My mistake `>=` is actually OK. N is a known parameter. The term `integer programming` is new to me and wikipedia says it is NP-complete. So there are no exact algorithms, only approximate ones? – jadelord Mar 16 '18 at 18:28
  • NP-completeness does not imply inexact results. – sascha Mar 16 '18 at 18:29
  • A better place to put this question might be https://scicomp.stackexchange.com/, or one of mathematics forums. – Bill Bell Mar 16 '18 at 20:38
  • 2
    For anyone interested in the simpler case of finding the next composite number with only the factors 2, 3, and 5, there is [`scipy.fftpack.next_fast_len`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.fftpack.next_fast_len.html). – Warren Weckesser Mar 17 '18 at 03:07
  • Good point @WarrenWeckesser! Apparently a [`pyfftw.next_fast_len` implementation](https://github.com/pyFFTW/pyFFTW/commit/7f12fcefab66573a4e7b845bf64f49d6de2f97f0) is also around, but we would have to wait till its next release. – jadelord Mar 17 '18 at 18:22

1 Answers1

1

According to @sascha, and a conversation in Scipy-Dev mailing list scipy has not yet implemented an mixed integer linear programming (MILP a.k.a. MIP) solver.

The Python wiki contains a list of LP / MILP python packages. I ended up using pulp and its built-in solver. By modifying the code from a similar problem and examples found in a paper, the solution looks like this.

#!/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)
jadelord
  • 1,511
  • 14
  • 19
  • 1
    Glad you found it and i hope it works for you. pulp is easy to use and easily packages the (imho) best open-source MIP-solver there is, Coin-OR Cbc (if you can live with it's less-than-optimal docs; compared to GLPK and co). – sascha Mar 17 '18 at 04:06