royer_closest_product.py

# coding: utf-8

# # Closest Product Python Twitter Challenge
# 
# https://twitter.com/loicaroyer/status/1146173537658404864
# 
# Given a list of integers u=[a_0, ..., a_m]
# find the set of indices {i_0, ..., i_n} for that list such that the product 
# u[i_0]* ... *u[i_n] is the closest to a given integer N.
# The shortest and most #elegant solution wins. 
# (standard libs allowed)

# ## Setup

# In[1]:

import numpy as np
import itertools as it


# In[2]:

# Setup inputs
# Use same inputs as randomly generated by @james_t_webber
# https://twitter.com/james_t_webber/status/1146203142079410178
m = 20
N = 797
u = np.array([5,9,19,22,25,27,28,31,33,34,36,43,44,68,72,81,86,92,96,98])
print('N: ',N)
print('u: ',u)


# ## Brute Force Method
# Calculate all the products and select the closest one. This takes on the order of seconds.

# In[3]:

# Brute Force code by @james_t_webber adapted into a function
# This version computes the products of all combinations of u
# without repetition. Changed m = u.size
# https://twitter.com/james_t_webber/status/1146210405368266752
def closest_product_brute_force(N,u):
    i = min(
        map(np.array, it.product([False, True],repeat=u.size)),
        key=lambda i:np.abs(N-u[i].prod())
    )
    i = np.nonzero(i)[0]
    return i


# In[4]:

get_ipython().magic('timeit closest_product_brute_force(N,u)')
i = closest_product_brute_force(N,u)
print('i: ',i)
print('u[i]: ',u[i])
print('u[i].prod(): ',u[i].prod())


# ## Speed-up Using Short Circuiting and Divisors
# 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.
# 
# To do this, we iterate over the integers closest to N starting from N working outwards.
# * Call each integer M.
# * Figure out which elements of u are divisors of M.
# * Compute the possible combinations of those divisors.
# * Compute the product of those divisor combinations.
# * Stop (short-circuit) if the product equals M
# * Choose M one step farther from 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.
# 
# This takes 100s of microseconds, a 10^4 speed-up over the brute force method for the given input.

# In[5]:

def closest_product_no_repetition(N,u):
    # Setup generator for product targets "M"
    Mg = it.chain.from_iterable(
        ( 
            (N-y,N+y) if y != 0 else (N,)
            for y in range(min(abs(N-u))+1)
        )
    )

    # For each M construct a list of divisors
    # Each generated element is a tuple: (divisors,M)
    divisors = filter(
        lambda f: f[0].size > 0,
        (
            (u[M % u == 0],M) for M in Mg
        )
    )


    # For each tuple of (divisors,M),
    # generate combinations of divisors
    divisor_combinations = it.chain.from_iterable(zip(
        it.chain.from_iterable(
            (it.combinations(d,r+1) for r in range(d.size))
        ),
        it.repeat(M)
    ) for d,M in divisors)

    # We want to use it.takewhile, but we need info on the divisor combination that met the condition
    # Use tee to duplicate the iterator and advance it by one
    g = it.tee(divisor_combinations)
    u_i = next(g[1])
    g = zip(*g)
    
    # Up to this point no products have been calculated, we have just set up iterators
    # Now, start iterating through M, divisors, and divisior_combinations
    # until the product equals M
    out = list(it.takewhile(lambda a: np.array(a[0][0]).prod() != a[0][1],g))
    if len(out) > 0:
        u_i = out[-1][1][0]
    else:
        u_i = list(u_i)[0]
        
    u_i = np.array(u_i)
    u_i = u_i[u_i != 1]
    i = list(map(lambda ui: np.where(u == ui)[0][0],u_i))
    return i


# In[6]:

get_ipython().magic('timeit i = closest_product_no_repetition(N,u)')
i = closest_product_no_repetition(N,u)
print('i: ',i)
print('u[i]: ',u[i])
print('u[i].prod(): ',u[i].prod())


# ## Repetition
# 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.

# In[7]:

def closest_product_with_repetition(N,u):
    # Setup generator for product targets "M"
    Mg = it.chain.from_iterable(
        ( 
            (N-y,N+y) if y != 0 else (N,)
            for y in range(min(abs(N-u))+1)
        )
    )

    # For each M construct a list of divisors
    # Each generated element is a tuple: (divisors,M)
    divisors = filter(
        lambda f: f[0].size > 1,
        (
            # Append 1 to the list of divisors so that we can use a fixed combination size
            (np.append(u[M % u == 0],1),M) for M in Mg
        )
    )

    # For each tuple of (divisors,M),
    # generate combinations of divisors with replacement
    divisor_combinations = it.chain.from_iterable(zip(
        it.combinations_with_replacement(
            d,
            # Since we allow for replacement,
            # set combination length such that a power of the smallest divisor exceeds M
            int(np.ceil(np.log(M)/np.log(min(d[d > 1]))))
        ),
        it.repeat(M)
    ) for d,M in divisors)

    # We want to use it.takewhile, but we need info on the divisor combination that met the condition
    # Use tee to duplicate the iterator and advance it by one
    g = it.tee(divisor_combinations)
    u_i = next(g[1])
    g = zip(*g)
    
    # Up to this point no products have been calculated, we have just set up iterators
    # Now, start iterating through M, divisors, and divisior_combinations
    # until the product equals M
    out = list(it.takewhile(lambda a: np.array(a[0][0]).prod() != a[0][1],g))
    if len(out) > 0:
        u_i = out[-1][1][0]
    else:
        u_i = list(u_i)[0]
        
    u_i = np.array(u_i)
    u_i = u_i[u_i != 1]
    i = list(map(lambda ui: np.where(u == ui)[0][0],u_i))
    return i


# In[8]:

get_ipython().magic('timeit i = closest_product_with_repetition(N,u)')
i = closest_product_with_repetition(N,u)
print('i: ',i)
print('u[i]: ',u[i])
print('u[i].prod(): ',u[i].prod())


# ## Breakdown of the method
# 
# First we define a helper function, peek_with_label to help take a look at some of the generators used.

# In[9]:

def peek_with_label(iterable,num=10,label=''):
    iterable = it.tee(iterable)
    print(label)
    for x in it.islice(iterable[0],num):
        print(x)
    return iterable[1]


# In[10]:

# Setup generator for product targets "M"
Mg = it.chain.from_iterable(
    ( 
        (N-y,N+y) if y != 0 else (N,)
        for y in range(min(abs(N-u))+1)
    )
)
Mg = peek_with_label(Mg,10,'M generator:')


# In[11]:

# For each M construct a list of divisors
# Each generated element is a tuple: (divisors,M)
divisors = filter(
    lambda f: f[0].size > 1,
    (
        # Append 1 to the list of divisors so that we can use a fixed combination size
        (np.append(u[M % u == 0],1),M) for M in Mg
    )
)

divisors = peek_with_label(divisors,10,'Divisors: ')


# In[12]:

# For each tuple of (divisors,M),
# generate combinations of divisors with replacement
divisor_combinations = it.chain.from_iterable(zip(
    it.combinations_with_replacement(
        d,
        # Since we allow for replacement,
        # set combination length such that a power of the smallest divisor exceeds M
        int(np.ceil(np.log(M)/np.log(min(d[d > 1]))))
    ),
    it.repeat(M)
) for d,M in divisors)

divisor_combinations = peek_with_label(divisor_combinations,240,'Divisor combinations with repetition: ')


# In[13]:

# We want to use it.takewhile, but we need info on the divisor combination that met the condition
# Use tee to duplicate the iterator and advance it by one
g = it.tee(divisor_combinations)
u_i = next(g[1])
g = zip(*g)

g = peek_with_label(g,10,'Paired divisor combinations:')


# In[14]:

# Up to this point no products have been calculated, we have just set up iterators
# Now, start iterating through M, divisors, and divisior_combinations
# until the product equals M
out = list(it.takewhile(lambda a: np.array(a[0][0]).prod() != a[0][1],g))

for x in out[-10:]:
    print(x)


# In[15]:

if len(out) > 0:
    u_i = out[-1][1][0]
else:
    u_i = list(u_i)[0]
    
u_i = np.array(u_i)
u_i = u_i[u_i != 1]
i = list(map(lambda ui: np.where(u == ui)[0][0],u_i))


# In[16]:

print('i: ',i)
print('u[i]: ',u[i])
print('u[i].prod(): ',u[i].prod())


# ## Unified Function with or without repetition

# In[17]:

def closest_product(N,u,repetition=False):
    # Process input
    u = np.array(u)
    # Setup generator for product targets "M"
    Mg = it.chain.from_iterable(
        ( 
            (N-y,N+y) if y != 0 else (N,)
            for y in range(min(abs(N-u))+1)
        )
    )

    if repetition:
        # For each M construct a list of divisors
        # Each generated element is a tuple: (divisors,M)
        divisors = filter(
            lambda f: f[0].size > 1,
            (
                # Append 1 to the list of divisors so that we can use a fixed combination size
                (np.append(u[M % u == 0],1),M) for M in Mg
            )
        )
        # For each tuple of (divisors,M),
        # generate combinations of divisors with replacement
        divisor_combinations = it.chain.from_iterable(zip(
            it.combinations_with_replacement(
                d,
                # Since we allow for replacement,
                # set combination length such that a power of the smallest divisor exceeds M
                int(np.ceil(np.log(M)/np.log(min(d[d > 1]))))
            ),
            it.repeat(M)
        ) for d,M in divisors)
    else:
        # For each M construct a list of divisors
        # Each generated element is a tuple: (divisors,M)
        divisors = filter(
            lambda f: f[0].size > 0,
            (
                (u[M % u == 0],M) for M in Mg
            )
        )
        # For each tuple of (divisors,M),
        # generate combinations of divisors
        divisor_combinations = it.chain.from_iterable(zip(
            it.chain.from_iterable(
                (it.combinations(d,r+1) for r in range(d.size))
            ),
            it.repeat(M)
        ) for d,M in divisors)

    # We want to use it.takewhile, but we need info on the divisor combination that met the condition
    # Use tee to duplicate the iterator and advance it by one
    g = it.tee(divisor_combinations)
    u_i = next(g[1])
    g = zip(*g)
    
    # Up to this point no products have been calculated, we have just set up iterators
    # Now, start iterating through M, divisors, and divisior_combinations
    # until the product equals M
    out = list(it.takewhile(lambda a: np.array(a[0][0]).prod() != a[0][1],g))
    if len(out) > 0:
        u_i = out[-1][1][0]
    else:
        u_i = list(u_i)[0]
        
    u_i = np.array(u_i)
    u_i = u_i[u_i != 1]
    i = list(map(lambda ui: np.where(u == ui)[0][0],u_i))
    return i


# ### Unified Function testing

# In[18]:

closest_product(8,[2,3],repetition=False)


# In[19]:

closest_product(8,np.array([2,3]),repetition=True)


# In[20]:

closest_product(8,[8],repetition=False)


# In[21]:

closest_product(8,[8],repetition=True)


# In[22]:

closest_product(N,u,repetition=False)


# In[23]:

closest_product(N,u,repetition=True)


# ## Create a random challenging case, with m = 100
# 
# If we create a random challenging case, executing time is in the milliseconds

# In[24]:

challenging_N = np.random.randint(101,2**20)
challenging_u = np.random.choice(range(2,1010),100,False)
print('challenging_N: ',challenging_N)
print('challenging_u: ',challenging_u)


# In[25]:

get_ipython().magic('time i = closest_product(challenging_N,challenging_u,repetition=False)')
i = closest_product(challenging_N,challenging_u,repetition=False)
print('i: ',i)
print('u[i]: ',challenging_u[i])
print('u[i].prod(): ',challenging_u[i].prod())