royer_closest_product.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Closest Product Python Twitter Challenge\n",
    "\n",
    "https://twitter.com/loicaroyer/status/1146173537658404864\n",
    "\n",
    "Given a list of integers u=[a_0, ..., a_m]\n",
    "find the set of indices {i_0, ..., i_n} for that list such that the product \n",
    "u[i_0]* ... *u[i_n] is the closest to a given integer N.\n",
    "The shortest and most #elegant solution wins. \n",
    "(standard libs allowed)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import itertools as it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "N:  797\n",
      "u:  [ 5  9 19 22 25 27 28 31 33 34 36 43 44 68 72 81 86 92 96 98]\n"
     ]
    }
   ],
   "source": [
    "# Setup inputs\n",
    "# Use same inputs as randomly generated by @james_t_webber\n",
    "# https://twitter.com/james_t_webber/status/1146203142079410178\n",
    "m = 20\n",
    "N = 797\n",
    "u = np.array([5,9,19,22,25,27,28,31,33,34,36,43,44,68,72,81,86,92,96,98])\n",
    "print('N: ',N)\n",
    "print('u: ',u)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Brute Force Method\n",
    "Calculate all the products and select the closest one. This takes on the order of seconds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Brute Force code by @james_t_webber adapted into a function\n",
    "# This version computes the products of all combinations of u\n",
    "# without repetition. Changed m = u.size\n",
    "# https://twitter.com/james_t_webber/status/1146210405368266752\n",
    "def closest_product_brute_force(N,u):\n",
    "    i = min(\n",
    "        map(np.array, it.product([False, True],repeat=u.size)),\n",
    "        key=lambda i:np.abs(N-u[i].prod())\n",
    "    )\n",
    "    i = np.nonzero(i)[0]\n",
    "    return i"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 loop, best of 3: 7.91 s per loop\n",
      "i:  [ 3 10]\n",
      "u[i]:  [22 36]\n",
      "u[i].prod():  792\n"
     ]
    }
   ],
   "source": [
    "%timeit closest_product_brute_force(N,u)\n",
    "i = closest_product_brute_force(N,u)\n",
    "print('i: ',i)\n",
    "print('u[i]: ',u[i])\n",
    "print('u[i].prod(): ',u[i].prod())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Speed-up Using Short Circuiting and Divisors\n",
    "Rather than computing all the product of all the combinations, we just need to compute the products of the integers closest to N. Once we find such a product we can stop.\n",
    "\n",
    "To do this, we iterate over the integers closest to N starting from N working outwards.\n",
    "* Call each integer M.\n",
    "* Figure out which elements of u are divisors of M.\n",
    "* Compute the possible combinations of those divisors.\n",
    "* Compute the product of those divisor combinations.\n",
    "* Stop (short-circuit) if the product equals M\n",
    "* Choose M one step farther from N\n",
    "\n",
    "To implement this we will use a special type of iterable called a generator. Generators use lazy evaluation which defers evaluation until the next element of the iterable is needed. This will help us define the iteration without having to actually compute anything until we need it.\n",
    "\n",
    "This takes 100s of microseconds, a 10^4 speed-up over the brute force method for the given input."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def closest_product_no_repetition(N,u):\n",
    "    # Setup generator for product targets \"M\"\n",
    "    Mg = it.chain.from_iterable(\n",
    "        ( \n",
    "            (N-y,N+y) if y != 0 else (N,)\n",
    "            for y in range(min(abs(N-u))+1)\n",
    "        )\n",
    "    )\n",
    "\n",
    "    # For each M construct a list of divisors\n",
    "    # Each generated element is a tuple: (divisors,M)\n",
    "    divisors = filter(\n",
    "        lambda f: f[0].size > 0,\n",
    "        (\n",
    "            (u[M % u == 0],M) for M in Mg\n",
    "        )\n",
    "    )\n",
    "\n",
    "\n",
    "    # For each tuple of (divisors,M),\n",
    "    # generate combinations of divisors\n",
    "    divisor_combinations = it.chain.from_iterable(zip(\n",
    "        it.chain.from_iterable(\n",
    "            (it.combinations(d,r+1) for r in range(d.size))\n",
    "        ),\n",
    "        it.repeat(M)\n",
    "    ) for d,M in divisors)\n",
    "\n",
    "    # We want to use it.takewhile, but we need info on the divisor combination that met the condition\n",
    "    # Use tee to duplicate the iterator and advance it by one\n",
    "    g = it.tee(divisor_combinations)\n",
    "    u_i = next(g[1])\n",
    "    g = zip(*g)\n",
    "    \n",
    "    # Up to this point no products have been calculated, we have just set up iterators\n",
    "    # Now, start iterating through M, divisors, and divisior_combinations\n",
    "    # until the product equals M\n",
    "    out = list(it.takewhile(lambda a: np.array(a[0][0]).prod() != a[0][1],g))\n",
    "    if len(out) > 0:\n",
    "        u_i = out[-1][1][0]\n",
    "    else:\n",
    "        u_i = list(u_i)[0]\n",
    "        \n",
    "    u_i = np.array(u_i)\n",
    "    u_i = u_i[u_i != 1]\n",
    "    i = list(map(lambda ui: np.where(u == ui)[0][0],u_i))\n",
    "    return i"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1000 loops, best of 3: 190 µs per loop\n",
      "i:  [3, 10]\n",
      "u[i]:  [22 36]\n",
      "u[i].prod():  792\n"
     ]
    }
   ],
   "source": [
    "%timeit i = closest_product_no_repetition(N,u)\n",
    "i = closest_product_no_repetition(N,u)\n",
    "print('i: ',i)\n",
    "print('u[i]: ',u[i])\n",
    "print('u[i].prod(): ',u[i].prod())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Repetition\n",
    "It is not clear from the original problem description if elements of u can be repeated. This is easily addressed by just taking all the combinations with repetition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def closest_product_with_repetition(N,u):\n",
    "    # Setup generator for product targets \"M\"\n",
    "    Mg = it.chain.from_iterable(\n",
    "        ( \n",
    "            (N-y,N+y) if y != 0 else (N,)\n",
    "            for y in range(min(abs(N-u))+1)\n",
    "        )\n",
    "    )\n",
    "\n",
    "    # For each M construct a list of divisors\n",
    "    # Each generated element is a tuple: (divisors,M)\n",
    "    divisors = filter(\n",
    "        lambda f: f[0].size > 1,\n",
    "        (\n",
    "            # Append 1 to the list of divisors so that we can use a fixed combination size\n",
    "            (np.append(u[M % u == 0],1),M) for M in Mg\n",
    "        )\n",
    "    )\n",
    "\n",
    "    # For each tuple of (divisors,M),\n",
    "    # generate combinations of divisors with replacement\n",
    "    divisor_combinations = it.chain.from_iterable(zip(\n",
    "        it.combinations_with_replacement(\n",
    "            d,\n",
    "            # Since we allow for replacement,\n",
    "            # set combination length such that a power of the smallest divisor exceeds M\n",
    "            int(np.ceil(np.log(M)/np.log(min(d[d > 1]))))\n",
    "        ),\n",
    "        it.repeat(M)\n",
    "    ) for d,M in divisors)\n",
    "\n",
    "    # We want to use it.takewhile, but we need info on the divisor combination that met the condition\n",
    "    # Use tee to duplicate the iterator and advance it by one\n",
    "    g = it.tee(divisor_combinations)\n",
    "    u_i = next(g[1])\n",
    "    g = zip(*g)\n",
    "    \n",
    "    # Up to this point no products have been calculated, we have just set up iterators\n",
    "    # Now, start iterating through M, divisors, and divisior_combinations\n",
    "    # until the product equals M\n",
    "    out = list(it.takewhile(lambda a: np.array(a[0][0]).prod() != a[0][1],g))\n",
    "    if len(out) > 0:\n",
    "        u_i = out[-1][1][0]\n",
    "    else:\n",
    "        u_i = list(u_i)[0]\n",
    "        \n",
    "    u_i = np.array(u_i)\n",
    "    u_i = u_i[u_i != 1]\n",
    "    i = list(map(lambda ui: np.where(u == ui)[0][0],u_i))\n",
    "    return i"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1000 loops, best of 3: 1.05 ms per loop\n",
      "i:  [3, 10]\n",
      "u[i]:  [22 36]\n",
      "u[i].prod():  792\n"
     ]
    }
   ],
   "source": [
    "%timeit i = closest_product_with_repetition(N,u)\n",
    "i = closest_product_with_repetition(N,u)\n",
    "print('i: ',i)\n",
    "print('u[i]: ',u[i])\n",
    "print('u[i].prod(): ',u[i].prod())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Breakdown of the method\n",
    "\n",
    "First we define a helper function, peek_with_label to help take a look at some of the generators used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def peek_with_label(iterable,num=10,label=''):\n",
    "    iterable = it.tee(iterable)\n",
    "    print(label)\n",
    "    for x in it.islice(iterable[0],num):\n",
    "        print(x)\n",
    "    return iterable[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "M generator:\n",
      "797\n",
      "796\n",
      "798\n",
      "795\n",
      "799\n",
      "794\n",
      "800\n",
      "793\n",
      "801\n",
      "792\n"
     ]
    }
   ],
   "source": [
    "# Setup generator for product targets \"M\"\n",
    "Mg = it.chain.from_iterable(\n",
    "    ( \n",
    "        (N-y,N+y) if y != 0 else (N,)\n",
    "        for y in range(min(abs(N-u))+1)\n",
    "    )\n",
    ")\n",
    "Mg = peek_with_label(Mg,10,'M generator:')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Divisors: \n",
      "(array([19,  1]), 798)\n",
      "(array([5, 1]), 795)\n",
      "(array([ 5, 25,  1]), 800)\n",
      "(array([9, 1]), 801)\n",
      "(array([ 9, 22, 33, 36, 44, 72,  1]), 792)\n",
      "(array([5, 1]), 790)\n",
      "(array([5, 1]), 805)\n",
      "(array([31,  1]), 806)\n",
      "(array([5, 1]), 785)\n",
      "(array([28, 98,  1]), 784)\n"
     ]
    }
   ],
   "source": [
    "# For each M construct a list of divisors\n",
    "# Each generated element is a tuple: (divisors,M)\n",
    "divisors = filter(\n",
    "    lambda f: f[0].size > 1,\n",
    "    (\n",
    "        # Append 1 to the list of divisors so that we can use a fixed combination size\n",
    "        (np.append(u[M % u == 0],1),M) for M in Mg\n",
    "    )\n",
    ")\n",
    "\n",
    "divisors = peek_with_label(divisors,10,'Divisors: ')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Divisor combinations with repetition: \n",
      "((19, 19, 19), 798)\n",
      "((19, 19, 1), 798)\n",
      "((19, 1, 1), 798)\n",
      "((1, 1, 1), 798)\n",
      "((5, 5, 5, 5, 5), 795)\n",
      "((5, 5, 5, 5, 1), 795)\n",
      "((5, 5, 5, 1, 1), 795)\n",
      "((5, 5, 1, 1, 1), 795)\n",
      "((5, 1, 1, 1, 1), 795)\n",
      "((1, 1, 1, 1, 1), 795)\n",
      "((5, 5, 5, 5, 5), 800)\n",
      "((5, 5, 5, 5, 25), 800)\n",
      "((5, 5, 5, 5, 1), 800)\n",
      "((5, 5, 5, 25, 25), 800)\n",
      "((5, 5, 5, 25, 1), 800)\n",
      "((5, 5, 5, 1, 1), 800)\n",
      "((5, 5, 25, 25, 25), 800)\n",
      "((5, 5, 25, 25, 1), 800)\n",
      "((5, 5, 25, 1, 1), 800)\n",
      "((5, 5, 1, 1, 1), 800)\n",
      "((5, 25, 25, 25, 25), 800)\n",
      "((5, 25, 25, 25, 1), 800)\n",
      "((5, 25, 25, 1, 1), 800)\n",
      "((5, 25, 1, 1, 1), 800)\n",
      "((5, 1, 1, 1, 1), 800)\n",
      "((25, 25, 25, 25, 25), 800)\n",
      "((25, 25, 25, 25, 1), 800)\n",
      "((25, 25, 25, 1, 1), 800)\n",
      "((25, 25, 1, 1, 1), 800)\n",
      "((25, 1, 1, 1, 1), 800)\n",
      "((1, 1, 1, 1, 1), 800)\n",
      "((9, 9, 9, 9), 801)\n",
      "((9, 9, 9, 1), 801)\n",
      "((9, 9, 1, 1), 801)\n",
      "((9, 1, 1, 1), 801)\n",
      "((1, 1, 1, 1), 801)\n",
      "((9, 9, 9, 9), 792)\n",
      "((9, 9, 9, 22), 792)\n",
      "((9, 9, 9, 33), 792)\n",
      "((9, 9, 9, 36), 792)\n",
      "((9, 9, 9, 44), 792)\n",
      "((9, 9, 9, 72), 792)\n",
      "((9, 9, 9, 1), 792)\n",
      "((9, 9, 22, 22), 792)\n",
      "((9, 9, 22, 33), 792)\n",
      "((9, 9, 22, 36), 792)\n",
      "((9, 9, 22, 44), 792)\n",
      "((9, 9, 22, 72), 792)\n",
      "((9, 9, 22, 1), 792)\n",
      "((9, 9, 33, 33), 792)\n",
      "((9, 9, 33, 36), 792)\n",
      "((9, 9, 33, 44), 792)\n",
      "((9, 9, 33, 72), 792)\n",
      "((9, 9, 33, 1), 792)\n",
      "((9, 9, 36, 36), 792)\n",
      "((9, 9, 36, 44), 792)\n",
      "((9, 9, 36, 72), 792)\n",
      "((9, 9, 36, 1), 792)\n",
      "((9, 9, 44, 44), 792)\n",
      "((9, 9, 44, 72), 792)\n",
      "((9, 9, 44, 1), 792)\n",
      "((9, 9, 72, 72), 792)\n",
      "((9, 9, 72, 1), 792)\n",
      "((9, 9, 1, 1), 792)\n",
      "((9, 22, 22, 22), 792)\n",
      "((9, 22, 22, 33), 792)\n",
      "((9, 22, 22, 36), 792)\n",
      "((9, 22, 22, 44), 792)\n",
      "((9, 22, 22, 72), 792)\n",
      "((9, 22, 22, 1), 792)\n",
      "((9, 22, 33, 33), 792)\n",
      "((9, 22, 33, 36), 792)\n",
      "((9, 22, 33, 44), 792)\n",
      "((9, 22, 33, 72), 792)\n",
      "((9, 22, 33, 1), 792)\n",
      "((9, 22, 36, 36), 792)\n",
      "((9, 22, 36, 44), 792)\n",
      "((9, 22, 36, 72), 792)\n",
      "((9, 22, 36, 1), 792)\n",
      "((9, 22, 44, 44), 792)\n",
      "((9, 22, 44, 72), 792)\n",
      "((9, 22, 44, 1), 792)\n",
      "((9, 22, 72, 72), 792)\n",
      "((9, 22, 72, 1), 792)\n",
      "((9, 22, 1, 1), 792)\n",
      "((9, 33, 33, 33), 792)\n",
      "((9, 33, 33, 36), 792)\n",
      "((9, 33, 33, 44), 792)\n",
      "((9, 33, 33, 72), 792)\n",
      "((9, 33, 33, 1), 792)\n",
      "((9, 33, 36, 36), 792)\n",
      "((9, 33, 36, 44), 792)\n",
      "((9, 33, 36, 72), 792)\n",
      "((9, 33, 36, 1), 792)\n",
      "((9, 33, 44, 44), 792)\n",
      "((9, 33, 44, 72), 792)\n",
      "((9, 33, 44, 1), 792)\n",
      "((9, 33, 72, 72), 792)\n",
      "((9, 33, 72, 1), 792)\n",
      "((9, 33, 1, 1), 792)\n",
      "((9, 36, 36, 36), 792)\n",
      "((9, 36, 36, 44), 792)\n",
      "((9, 36, 36, 72), 792)\n",
      "((9, 36, 36, 1), 792)\n",
      "((9, 36, 44, 44), 792)\n",
      "((9, 36, 44, 72), 792)\n",
      "((9, 36, 44, 1), 792)\n",
      "((9, 36, 72, 72), 792)\n",
      "((9, 36, 72, 1), 792)\n",
      "((9, 36, 1, 1), 792)\n",
      "((9, 44, 44, 44), 792)\n",
      "((9, 44, 44, 72), 792)\n",
      "((9, 44, 44, 1), 792)\n",
      "((9, 44, 72, 72), 792)\n",
      "((9, 44, 72, 1), 792)\n",
      "((9, 44, 1, 1), 792)\n",
      "((9, 72, 72, 72), 792)\n",
      "((9, 72, 72, 1), 792)\n",
      "((9, 72, 1, 1), 792)\n",
      "((9, 1, 1, 1), 792)\n",
      "((22, 22, 22, 22), 792)\n",
      "((22, 22, 22, 33), 792)\n",
      "((22, 22, 22, 36), 792)\n",
      "((22, 22, 22, 44), 792)\n",
      "((22, 22, 22, 72), 792)\n",
      "((22, 22, 22, 1), 792)\n",
      "((22, 22, 33, 33), 792)\n",
      "((22, 22, 33, 36), 792)\n",
      "((22, 22, 33, 44), 792)\n",
      "((22, 22, 33, 72), 792)\n",
      "((22, 22, 33, 1), 792)\n",
      "((22, 22, 36, 36), 792)\n",
      "((22, 22, 36, 44), 792)\n",
      "((22, 22, 36, 72), 792)\n",
      "((22, 22, 36, 1), 792)\n",
      "((22, 22, 44, 44), 792)\n",
      "((22, 22, 44, 72), 792)\n",
      "((22, 22, 44, 1), 792)\n",
      "((22, 22, 72, 72), 792)\n",
      "((22, 22, 72, 1), 792)\n",
      "((22, 22, 1, 1), 792)\n",
      "((22, 33, 33, 33), 792)\n",
      "((22, 33, 33, 36), 792)\n",
      "((22, 33, 33, 44), 792)\n",
      "((22, 33, 33, 72), 792)\n",
      "((22, 33, 33, 1), 792)\n",
      "((22, 33, 36, 36), 792)\n",
      "((22, 33, 36, 44), 792)\n",
      "((22, 33, 36, 72), 792)\n",
      "((22, 33, 36, 1), 792)\n",
      "((22, 33, 44, 44), 792)\n",
      "((22, 33, 44, 72), 792)\n",
      "((22, 33, 44, 1), 792)\n",
      "((22, 33, 72, 72), 792)\n",
      "((22, 33, 72, 1), 792)\n",
      "((22, 33, 1, 1), 792)\n",
      "((22, 36, 36, 36), 792)\n",
      "((22, 36, 36, 44), 792)\n",
      "((22, 36, 36, 72), 792)\n",
      "((22, 36, 36, 1), 792)\n",
      "((22, 36, 44, 44), 792)\n",
      "((22, 36, 44, 72), 792)\n",
      "((22, 36, 44, 1), 792)\n",
      "((22, 36, 72, 72), 792)\n",
      "((22, 36, 72, 1), 792)\n",
      "((22, 36, 1, 1), 792)\n",
      "((22, 44, 44, 44), 792)\n",
      "((22, 44, 44, 72), 792)\n",
      "((22, 44, 44, 1), 792)\n",
      "((22, 44, 72, 72), 792)\n",
      "((22, 44, 72, 1), 792)\n",
      "((22, 44, 1, 1), 792)\n",
      "((22, 72, 72, 72), 792)\n",
      "((22, 72, 72, 1), 792)\n",
      "((22, 72, 1, 1), 792)\n",
      "((22, 1, 1, 1), 792)\n",
      "((33, 33, 33, 33), 792)\n",
      "((33, 33, 33, 36), 792)\n",
      "((33, 33, 33, 44), 792)\n",
      "((33, 33, 33, 72), 792)\n",
      "((33, 33, 33, 1), 792)\n",
      "((33, 33, 36, 36), 792)\n",
      "((33, 33, 36, 44), 792)\n",
      "((33, 33, 36, 72), 792)\n",
      "((33, 33, 36, 1), 792)\n",
      "((33, 33, 44, 44), 792)\n",
      "((33, 33, 44, 72), 792)\n",
      "((33, 33, 44, 1), 792)\n",
      "((33, 33, 72, 72), 792)\n",
      "((33, 33, 72, 1), 792)\n",
      "((33, 33, 1, 1), 792)\n",
      "((33, 36, 36, 36), 792)\n",
      "((33, 36, 36, 44), 792)\n",
      "((33, 36, 36, 72), 792)\n",
      "((33, 36, 36, 1), 792)\n",
      "((33, 36, 44, 44), 792)\n",
      "((33, 36, 44, 72), 792)\n",
      "((33, 36, 44, 1), 792)\n",
      "((33, 36, 72, 72), 792)\n",
      "((33, 36, 72, 1), 792)\n",
      "((33, 36, 1, 1), 792)\n",
      "((33, 44, 44, 44), 792)\n",
      "((33, 44, 44, 72), 792)\n",
      "((33, 44, 44, 1), 792)\n",
      "((33, 44, 72, 72), 792)\n",
      "((33, 44, 72, 1), 792)\n",
      "((33, 44, 1, 1), 792)\n",
      "((33, 72, 72, 72), 792)\n",
      "((33, 72, 72, 1), 792)\n",
      "((33, 72, 1, 1), 792)\n",
      "((33, 1, 1, 1), 792)\n",
      "((36, 36, 36, 36), 792)\n",
      "((36, 36, 36, 44), 792)\n",
      "((36, 36, 36, 72), 792)\n",
      "((36, 36, 36, 1), 792)\n",
      "((36, 36, 44, 44), 792)\n",
      "((36, 36, 44, 72), 792)\n",
      "((36, 36, 44, 1), 792)\n",
      "((36, 36, 72, 72), 792)\n",
      "((36, 36, 72, 1), 792)\n",
      "((36, 36, 1, 1), 792)\n",
      "((36, 44, 44, 44), 792)\n",
      "((36, 44, 44, 72), 792)\n",
      "((36, 44, 44, 1), 792)\n",
      "((36, 44, 72, 72), 792)\n",
      "((36, 44, 72, 1), 792)\n",
      "((36, 44, 1, 1), 792)\n",
      "((36, 72, 72, 72), 792)\n",
      "((36, 72, 72, 1), 792)\n",
      "((36, 72, 1, 1), 792)\n",
      "((36, 1, 1, 1), 792)\n",
      "((44, 44, 44, 44), 792)\n",
      "((44, 44, 44, 72), 792)\n",
      "((44, 44, 44, 1), 792)\n",
      "((44, 44, 72, 72), 792)\n",
      "((44, 44, 72, 1), 792)\n",
      "((44, 44, 1, 1), 792)\n",
      "((44, 72, 72, 72), 792)\n",
      "((44, 72, 72, 1), 792)\n",
      "((44, 72, 1, 1), 792)\n"
     ]
    }
   ],
   "source": [
    "# For each tuple of (divisors,M),\n",
    "# generate combinations of divisors with replacement\n",
    "divisor_combinations = it.chain.from_iterable(zip(\n",
    "    it.combinations_with_replacement(\n",
    "        d,\n",
    "        # Since we allow for replacement,\n",
    "        # set combination length such that a power of the smallest divisor exceeds M\n",
    "        int(np.ceil(np.log(M)/np.log(min(d[d > 1]))))\n",
    "    ),\n",
    "    it.repeat(M)\n",
    ") for d,M in divisors)\n",
    "\n",
    "divisor_combinations = peek_with_label(divisor_combinations,240,'Divisor combinations with repetition: ')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Paired divisor combinations:\n",
      "(((19, 19, 19), 798), ((19, 19, 1), 798))\n",
      "(((19, 19, 1), 798), ((19, 1, 1), 798))\n",
      "(((19, 1, 1), 798), ((1, 1, 1), 798))\n",
      "(((1, 1, 1), 798), ((5, 5, 5, 5, 5), 795))\n",
      "(((5, 5, 5, 5, 5), 795), ((5, 5, 5, 5, 1), 795))\n",
      "(((5, 5, 5, 5, 1), 795), ((5, 5, 5, 1, 1), 795))\n",
      "(((5, 5, 5, 1, 1), 795), ((5, 5, 1, 1, 1), 795))\n",
      "(((5, 5, 1, 1, 1), 795), ((5, 1, 1, 1, 1), 795))\n",
      "(((5, 1, 1, 1, 1), 795), ((1, 1, 1, 1, 1), 795))\n",
      "(((1, 1, 1, 1, 1), 795), ((5, 5, 5, 5, 5), 800))\n"
     ]
    }
   ],
   "source": [
    "# We want to use it.takewhile, but we need info on the divisor combination that met the condition\n",
    "# Use tee to duplicate the iterator and advance it by one\n",
    "g = it.tee(divisor_combinations)\n",
    "u_i = next(g[1])\n",
    "g = zip(*g)\n",
    "\n",
    "g = peek_with_label(g,10,'Paired divisor combinations:')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(((22, 33, 1, 1), 792), ((22, 36, 36, 36), 792))\n",
      "(((22, 36, 36, 36), 792), ((22, 36, 36, 44), 792))\n",
      "(((22, 36, 36, 44), 792), ((22, 36, 36, 72), 792))\n",
      "(((22, 36, 36, 72), 792), ((22, 36, 36, 1), 792))\n",
      "(((22, 36, 36, 1), 792), ((22, 36, 44, 44), 792))\n",
      "(((22, 36, 44, 44), 792), ((22, 36, 44, 72), 792))\n",
      "(((22, 36, 44, 72), 792), ((22, 36, 44, 1), 792))\n",
      "(((22, 36, 44, 1), 792), ((22, 36, 72, 72), 792))\n",
      "(((22, 36, 72, 72), 792), ((22, 36, 72, 1), 792))\n",
      "(((22, 36, 72, 1), 792), ((22, 36, 1, 1), 792))\n"
     ]
    }
   ],
   "source": [
    "# Up to this point no products have been calculated, we have just set up iterators\n",
    "# Now, start iterating through M, divisors, and divisior_combinations\n",
    "# until the product equals M\n",
    "out = list(it.takewhile(lambda a: np.array(a[0][0]).prod() != a[0][1],g))\n",
    "\n",
    "for x in out[-10:]:\n",
    "    print(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "    if len(out) > 0:\n",
    "        u_i = out[-1][1][0]\n",
    "    else:\n",
    "        u_i = list(u_i)[0]\n",
    "        \n",
    "    u_i = np.array(u_i)\n",
    "    u_i = u_i[u_i != 1]\n",
    "    i = list(map(lambda ui: np.where(u == ui)[0][0],u_i))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "i:  [3, 10]\n",
      "u[i]:  [22 36]\n",
      "u[i].prod():  792\n"
     ]
    }
   ],
   "source": [
    "print('i: ',i)\n",
    "print('u[i]: ',u[i])\n",
    "print('u[i].prod(): ',u[i].prod())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Unified Function with or without repetition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def closest_product(N,u,repetition=False):\n",
    "    # Process input\n",
    "    u = np.array(u)\n",
    "    # Setup generator for product targets \"M\"\n",
    "    Mg = it.chain.from_iterable(\n",
    "        ( \n",
    "            (N-y,N+y) if y != 0 else (N,)\n",
    "            for y in range(min(abs(N-u))+1)\n",
    "        )\n",
    "    )\n",
    "\n",
    "    if repetition:\n",
    "        # For each M construct a list of divisors\n",
    "        # Each generated element is a tuple: (divisors,M)\n",
    "        divisors = filter(\n",
    "            lambda f: f[0].size > 1,\n",
    "            (\n",
    "                # Append 1 to the list of divisors so that we can use a fixed combination size\n",
    "                (np.append(u[M % u == 0],1),M) for M in Mg\n",
    "            )\n",
    "        )\n",
    "        # For each tuple of (divisors,M),\n",
    "        # generate combinations of divisors with replacement\n",
    "        divisor_combinations = it.chain.from_iterable(zip(\n",
    "            it.combinations_with_replacement(\n",
    "                d,\n",
    "                # Since we allow for replacement,\n",
    "                # set combination length such that a power of the smallest divisor exceeds M\n",
    "                int(np.ceil(np.log(M)/np.log(min(d[d > 1]))))\n",
    "            ),\n",
    "            it.repeat(M)\n",
    "        ) for d,M in divisors)\n",
    "    else:\n",
    "        # For each M construct a list of divisors\n",
    "        # Each generated element is a tuple: (divisors,M)\n",
    "        divisors = filter(\n",
    "            lambda f: f[0].size > 0,\n",
    "            (\n",
    "                (u[M % u == 0],M) for M in Mg\n",
    "            )\n",
    "        )\n",
    "        # For each tuple of (divisors,M),\n",
    "        # generate combinations of divisors\n",
    "        divisor_combinations = it.chain.from_iterable(zip(\n",
    "            it.chain.from_iterable(\n",
    "                (it.combinations(d,r+1) for r in range(d.size))\n",
    "            ),\n",
    "            it.repeat(M)\n",
    "        ) for d,M in divisors)\n",
    "\n",
    "    # We want to use it.takewhile, but we need info on the divisor combination that met the condition\n",
    "    # Use tee to duplicate the iterator and advance it by one\n",
    "    g = it.tee(divisor_combinations)\n",
    "    u_i = next(g[1])\n",
    "    g = zip(*g)\n",
    "    \n",
    "    # Up to this point no products have been calculated, we have just set up iterators\n",
    "    # Now, start iterating through M, divisors, and divisior_combinations\n",
    "    # until the product equals M\n",
    "    out = list(it.takewhile(lambda a: np.array(a[0][0]).prod() != a[0][1],g))\n",
    "    if len(out) > 0:\n",
    "        u_i = out[-1][1][0]\n",
    "    else:\n",
    "        u_i = list(u_i)[0]\n",
    "        \n",
    "    u_i = np.array(u_i)\n",
    "    u_i = u_i[u_i != 1]\n",
    "    i = list(map(lambda ui: np.where(u == ui)[0][0],u_i))\n",
    "    return i"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Unified Function testing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0, 1]"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "closest_product(8,[2,3],repetition=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0, 0, 0]"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "closest_product(8,np.array([2,3]),repetition=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0]"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "closest_product(8,[8],repetition=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0]"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "closest_product(8,[8],repetition=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[3, 10]"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "closest_product(N,u,repetition=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[3, 10]"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "closest_product(N,u,repetition=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create a random challenging case, with m = 100\n",
    "\n",
    "If we create a random challenging case, executing time is in the milliseconds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "challenging_N:  121048\n",
      "challenging_u:  [ 955   96  249  501  167  940  607  923  186  871  118  470  772  773\n",
      "  112  805  349  605  890  423  692  335  251  663   82  959  344  882\n",
      "   16  602  488  454  383  520  320  311  804  408   52   93   65  137\n",
      "   62  723  544  160  894  510  356  338  414  364  868  579  202  898\n",
      "  827  492  273  846  215  306  615  479  565  405  981  512  474  380\n",
      "  168  682  391  535   38  100  551  872  702  191  572  353  250  467\n",
      "   30   26  866  188   36  967  743  398  709  886  193 1004  584  554\n",
      "  731  393]\n"
     ]
    }
   ],
   "source": [
    "challenging_N = np.random.randint(101,2**20)\n",
    "challenging_u = np.random.choice(range(2,1010),100,False)\n",
    "print('challenging_N: ',challenging_N)\n",
    "print('challenging_u: ',challenging_u)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 6.94 ms, sys: 0 ns, total: 6.94 ms\n",
      "Wall time: 6.89 ms\n",
      "i:  [14, 84, 88]\n",
      "u[i]:  [112  30  36]\n",
      "u[i].prod():  120960\n"
     ]
    }
   ],
   "source": [
    "%time i = closest_product(challenging_N,challenging_u,repetition=False)\n",
    "i = closest_product(challenging_N,challenging_u,repetition=False)\n",
    "print('i: ',i)\n",
    "print('u[i]: ',challenging_u[i])\n",
    "print('u[i].prod(): ',challenging_u[i].prod())"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}