StackOverflow2013

Note that there are some explanatory texts on larger screens.

plurals

PO
text
Body
copied!<p>This answer is based on <a href="https://stackoverflow.com/questions/20779552/code-is-taking-too-much-time#comment31148113_20779877">@Tim Peters' suggestion about Hamiltonian paths</a>.</p> <p>There are many possible solutions. To avoid excessive memory consumption for intermediate solutions, a random path can be generated. It also allows to utilize multiple CPUs easily (each cpu generates its own paths in parallel).</p> <pre><code>import multiprocessing as mp import sys def main(): number = int(sys.argv[1]) # directed graph, vertices: 1..number (including ends) # there is an edge between i and j if (i+j) is prime vertices = range(1, number+1) G = {} # vertex -> adjacent vertices is_prime = sieve_of_eratosthenes(2*number+1) for i in vertices: G[i] = [] for j in vertices: if is_prime[i + j]: G[i].append(j) # there is an edge from i to j in the graph # utilize multiple cpus q = mp.Queue() for _ in range(mp.cpu_count()): p = mp.Process(target=hamiltonian_random, args=[G, q]) p.daemon = True # do not survive the main process p.start() print(q.get()) if __name__=="__main__": main() </code></pre> <p>where <a href="http://en.wikipedia.org/wiki/Sieve_of_Eratosthenes" rel="nofollow noreferrer">Sieve of Eratosthenes</a> is:</p> <pre><code>def sieve_of_eratosthenes(limit): is_prime = [True]*limit is_prime[0] = is_prime[1] = False # zero and one are not primes for n in range(int(limit**.5 + .5)): if is_prime[n]: for composite in range(n*n, limit, n): is_prime[composite] = False return is_prime </code></pre> <p>and:</p> <pre><code>import random def hamiltonian_random(graph, result_queue): """Build random paths until Hamiltonian path is found.""" vertices = list(graph.keys()) while True: # build random path path = [random.choice(vertices)] # start with a random vertice while True: # until path can be extended with a random adjacent vertex neighbours = graph[path[-1]] random.shuffle(neighbours) for adjacent_vertex in neighbours: if adjacent_vertex not in path: path.append(adjacent_vertex) break else: # can't extend path break # check whether it is hamiltonian if len(path) == len(vertices): assert set(path) == set(vertices) result_queue.put(path) # found hamiltonian path return </code></pre> <h3>Example</h3> <pre><code>$ python order-adjacent-prime-sum.py 20 </code></pre> <h3>Output</h3> <pre><code>[19, 18, 13, 10, 1, 4, 9, 14, 5, 6, 17, 2, 15, 16, 7, 12, 11, 8, 3, 20] </code></pre> <p>The output is a random sequence that satisfies the conditions:</p> <ul> <li>it is a permutation of the range from 1 to 20 (including)</li> <li>the sum of adjacent numbers is prime</li> </ul> <h3>Time performance</h3> <p>It takes around 10 seconds on average to get result for <code>n = 900</code> and extrapolating the time as exponential function, it should take around <code>20</code> seconds for <code>n = 1000</code>:</p> <p><img src="https://i.stack.imgur.com/cs6T6.png" alt="time performance (no set solution)"></p> <p>The image is generated using this code:</p> <pre><code>import numpy as np figname = 'hamiltonian_random_noset-noseq-900-900' Ns, Ts = np.loadtxt(figname+'.xy', unpack=True) # use polyfit to fit the data # y = c*a**n # log y = log (c * a ** n) # log Ts = log c + Ns * log a coeffs = np.polyfit(Ns, np.log2(Ts), deg=1) poly = np.poly1d(coeffs, variable='Ns') # use curve_fit to fit the data from scipy.optimize import curve_fit def func(x, a, c): return c*a**x popt, pcov = curve_fit(func, Ns, Ts) aa, cc = popt a, c = 2**coeffs # plot it import matplotlib.pyplot as plt plt.figure() plt.plot(Ns, np.log2(Ts), 'ko', label='time measurements') plt.plot(Ns, np.polyval(poly, Ns), 'r-', label=r'$time = %.2g\times %.4g^N$' % (c, a)) plt.plot(Ns, np.log2(func(Ns, *popt)), 'b-', label=r'$time = %.2g\times %.4g^N$' % (cc, aa)) plt.xlabel('N') plt.ylabel('log2(time in seconds)') plt.legend(loc='upper left') plt.show() </code></pre> <p>Fitted values:</p> <pre><code>>>> c*a**np.array([900, 1000]) array([ 11.37200806, 21.56029156]) >>> func([900, 1000], *popt) array([ 14.1521409 , 22.62916398]) </code></pre>

Querying!

Guidance

An individual column

Larger individual text columns get their own page to allow for proper reading.

SQuiL has stopped working due to an internal error.

If you are curious you may find further information in the browser console, which is accessible through the devtools (F12).

Reload