Challenge: sum of squares. Let's break 299

[Allen]

(01-19-2018 04:29 AM)Paul Dale Wrote:  If you follow the comments and links from the numberphile video, it looks like the maximum has been pushed up very significantly. Assuming the implementation was correct and addressed the problem in question etc.

Baring an inspired approach, the problem comes down to finding a Hamilton path in a graph. This is a NP complete problem.


Pauli

I haven't read through the git lab code or the youtube comments.. do you know what the maximum is now?

[Allen]

I believe there are multiple solutions to each graph walk. For example the 999 run above returned 3 unique solutions. For the first 888 items, the graphs walks are the same, but from there on, there are slightly different cycles in the 3 found solutions.
Here are the position numbers and the values for the last 122 steps for each of the 3 solutions.

Code:


888 58   89  89
889 111 107 107
890 85 118 118
891 84 138 138
892 112 58 58
893 113 111 111
894 83 85 85
895 142 84 84
896 114 112 112
897 82 113 113
898 62 83 83
899 107 142 142
900 118 114 114
901 138 82 82
906 89 64 64
907 136 105 105
909 105 136 136
910 91 60 60
933 64 45 45
934 57 76 76
935 43 93 93
936 78 51 51
937 66 70 70
938 55 74 74
939 45 95 95
940 76 49 49
941 93 72 72
942 51 97 97
943 70 47 47
944 74 53 53
945 95 91 91
946 49 78 78
947 72 66 66
948 97 55 55
949 47 26 26
950 53 23 23
951 28 98 98
952 21 46 46
953 60 18 18
954 61 63 63
955 39 37 37
956 42 27 27
958 27 59 59
959 37 62 62
960 63 38 38
961 18 43 43
962 46 57 57
963 98 7 7
964 23 42 42
965 41 39 39
966 59 61 61
967 5 20 20
968 20 29 29
969 29 52 52
970 52 48 48
971 48 33 33
972 33 16 16
973 16 9 9
974 9 40 40
975 40 41 41
976 24 8 8
977 25 28 28
978 11 21 21
979 38 4 4
980 26 32 32
981 10 17 17
982 15 19 19
983 34 30 30
984 30 34 6
985 19 2 10
986 6 14 15
987 3 35 34
988 13 1 2
989 12 15 14
990 4 10 35
991 32 6 1
992 17 3 3
993 8 13 13
994 1 12 12
995 35 24 24
996 14 25 25
997 2 11 11
998 7 5 5

Edit to add:
Note- part of the multiple solutions has to do with "skeleton key" cycles like (6,19,30) that can appear in any order together. Mathematically these are the solutions occur when all permutations of a given numbers are subsets of the 3-node graph.

On the N=1000 graph there are only 828 such wonders (*6 if you include the permutations) out of 300300 valid combinations of 3 nodes. Each of these strings is an opportunity to create a new cycle in the graph.
Since the smallest pair of these is: (6, 19, 30) and (5, 20, 44),

If at least one solution exists for N>=44, there could be multiple solutions, since any solution could contain a pair of these special strings (one to depart the main thread and the other to come back).

A more thorough proof would prove that there are no isolated parts of the a particular graph (i.e. closed cycle), and that it contained one of these "skeleton keys".

Here are the most frequent nodes in the N=1000 graph super nodes with their associated counts:

Code:


[(2, 15),
 (8, 13),
 (20, 10),
 (30, 10),
 (32, 10),
 (38, 10),
 (44, 10),
 (92, 10),
 (12, 9),
 (14, 9),
 (18, 9),
 (42, 9),
 (48, 9),
 (58, 9),
 (68, 9),
 (78, 9),
 (80, 9),
 (132, 9),
 (194, 9),
 (248, 9)]

an interesting pattern..

[Allen]

Code to find these super nodes:

Code:


def super_nodes(goal,length=3):
    import itertools
    squares=set([x**2 for x in range(2*goal)])
    increasing=[(i,x) for i in range(1,goal+1) for x in range(i+1,goal+1) if i+x in squares]
    f=set(increasing+[x[::-1] for x in increasing])
    edges={j:{x[1] for x in f if x[0]==j} for j in range(1,goal+1)}
    
    for l in range(length-2):
        output=[]
        for i in f:
            next1=edges[i[-1]]-set(i)
            output+=[i+(x,) for x in next1]
        f=set(output)
    
    found=0
    perms=[]
    done=set()
    print 'precalculation complete.. looking for super nodes among {} items starting with {}'.format(len(f),sorted(f,key=lambda x:x[0])[0])
    for j in sorted(f,key=lambda x:x[0]):
        this_perm=set(itertools.permutations(list(j)))
        done=done.union(this_perm)
        if this_perm.issubset(f):
            found+=1
            perms+=[j]
            print 'subset:{} remaining {} found {}'.format(j,len(f),found)
            f=f-done
            done=set()
        if len(done)>1000:
            f=f-done
            done=set()
    return sorted(perms,key=lambda x:max(x))

Edit-fixing code bugs/ speedup to remove searched items as you go..

[TheKaneB]

@Allen: impressive work, really appreciated

[Allen]

number of unique solutions found for each round (includes reversals):
(this is not exhaustive, just how many were left after 1 run)

Unique solutions found for len 15: 2
Unique solutions found for len 16: 2
Unique solutions found for len 17: 2
Unique solutions found for len 23: 6
Unique solutions found for len 25: 20
Unique solutions found for len 26: 24
Unique solutions found for len 27: 58
Unique solutions found for len 28: 56
Unique solutions found for len 29: 20
Unique solutions found for len 30: 22
Unique solutions found for len 31: 142
Unique solutions found for len 32: 124
Unique solutions found for len 33: 44
Unique solutions found for len 34: 40
Unique solutions found for len 35: 88
Unique solutions found for len 36: 40
Unique solutions found for len 37: 46
Unique solutions found for len 38: 28
Unique solutions found for len 39: 18
Unique solutions found for len 40: 18
Unique solutions found for len 41: 24
Unique solutions found for len 42: 20
Unique solutions found for len 43: 30
Unique solutions found for len 44: 6
Unique solutions found for len 45: 24
Unique solutions found for len 46: 6
Unique solutions found for len 47: 2
Unique solutions found for len 48: 4
Unique solutions found for len 49: 12
Unique solutions found for len 50: 22
Unique solutions found for len 51: 8
Unique solutions found for len 52: 8
Unique solutions found for len 53: 18
Unique solutions found for len 54: 4
Unique solutions found for len 57: 8
Unique solutions found for len 60: 4
Unique solutions found for len 61: 4
Unique solutions found for len 63: 20
Unique solutions found for len 64: 12
Unique solutions found for len 65: 4
Unique solutions found for len 66: 8
Unique solutions found for len 67: 4
Unique solutions found for len 69: 2
Unique solutions found for len 71: 6
Unique solutions found for len 72: 36
Unique solutions found for len 73: 24
Unique solutions found for len 74: 22
Unique solutions found for len 76: 8
Unique solutions found for len 78: 14
Unique solutions found for len 80: 16
Unique solutions found for len 84: 14
Unique solutions found for len 85: 18
Unique solutions found for len 86: 6
Unique solutions found for len 87: 10
Unique solutions found for len 88: 26
Unique solutions found for len 89: 40
Unique solutions found for len 90: 16
Unique solutions found for len 91: 2
Unique solutions found for len 92: 10
Unique solutions found for len 93: 10
Unique solutions found for len 96: 6
Unique solutions found for len 97: 18
Unique solutions found for len 98: 8

[pier4r]

(01-24-2018 10:29 PM)Allen Wrote:  Thomas,
Thank you for verifying those solutions! It's always nice to have a second set of eyes. I have an important event in less than an hour, but generally speaking, I think the hamiltonian approach is good for random-walk kind of graphs, but these "square" graphs have some interesting properties we can exploit.

Namely, there are bottlenecks in the graph. For N=71 ( and numbers around there) the bottleneck is 50. There is only 1 or 2 ways to include 50 in the string.

By sorting the nodes of the graph according to how many outgoing edges there are, one can write a penalty/reward scoring to promote those strings that contain the most bottlenecks. Make a population of possible strings, and sort them according to their score, decimating to only keeping a portion of the resulting population.

I'm still trying to find the optimal settings for decimation, and growth rates to ensure the population does not collapse, while still finishing reasonably quickly. Will post more as soon as I'm able.

Really nice. I have my ideas on the problem. Yes I know that it falls in a general problem that is NP hard, but special cases of a problem could be not necessary that hard.

I encourage you to contact numberphile/matt parker (standupmaths on youtube) for the solution.

Also I skimmed quickly what you wrote because I don't want to be influenced before attempting the problem myself (if I ever find time).

What is making me happy is that I knew that this forum has really enthusiastic members.

edit: oh and ultra thanks for sharing! It is always a bit meh when people find a not easy solutions and they don't share it (even on a paper). So you did great.

[Allen]

Thank you! Today on my 7-year-old-PC with only 16Gb Ram I calculated 81 unique paths with N=1500. I'd say it's no where near the upper limits of what can be done on a modern PC with multiple cores and larger memory!!

[TheKaneB]

(01-25-2018 11:10 PM)Allen Wrote:  Thank you! Today on my 7-year-old-PC with only 16Gb Ram I calculated 81 unique paths with N=1500. I'd say it's no where near the upper limits of what can be done on a modern PC with multiple cores and larger memory!!

And also, python is not designed for speed. A C or C++ implementation (or even Java or C#) would be much faster! The really nice part is the algorithm that you built, that made the real difference

[pier4r]

(01-25-2018 11:10 PM)Allen Wrote:  Thank you! Today on my 7-year-old-PC with only 16Gb Ram I calculated 81 unique paths with N=1500. I'd say it's no where near the upper limits of what can be done on a modern PC with multiple cores and larger memory!!

16 gb of ram? My pc at work is from 2017 (mabook pro retina 13' with touch bar) and it has less.

Appreciate your little monster! Also this means that on a calculator would be extra difficult to reach numbers after 30-40 (I am guessing).

[Allen]

I found a speed up and can now calculate N=1000 solution in under 20 seconds.
N=2000 takes <3 min
N=3500 about 16 min

[pier4r]

(01-26-2018 10:01 PM)Allen Wrote:  I found a speed up and can now calculate N=1000 solution in under 20 seconds.
N=2000 takes <3 min
N=3500 about 16 min

I hope you share your findings somewhere (git repositories, articles on some magazines or journals, your blog, etc.). If something is new or relatively new, that is few know how to use it, it is a pity to not share it because it may be lost until someone else with same abilities looks into the topic.

I always remember the middle square "random" procedure.

Discovered in ~1950 and recorded in a journal. It was already recorded in ~1930 and in ~1250 (yes, twelfth hundred fifty) . I wondered how many times it was independently discovered until the knowledge of the procedure was more widespread after 1950. Now reinventing the wheel because one wants to is all ok, reinventing the wheel because all the previous one got lost is a necessity (there are no wheels) but also a pity (no one cares about leaving the wheels survive).

[Massimo Gnerucci]

(01-27-2018 01:23 PM)pier4r Wrote:  

(01-26-2018 10:01 PM)Allen Wrote:  I found a speed up and can now calculate N=1000 solution in under 20 seconds.
N=2000 takes <3 min
N=3500 about 16 min

I hope you share your findings somewhere (git repositories, articles on some magazines or journals, your blog, etc.). If something is new or relatively new, that is few know how to use it, it is a pity to not share it because it may be lost until someone else with same abilities looks into the topic.

I always remember the middle square "random" procedure.

Discovered in ~1950 and recorded in a journal. It was already recorded in ~1930 and in ~1250 (yes, twelfth hundred fifty) . I wondered how many times it was independently discovered until the knowledge of the procedure was more widespread after 1950. Now reinventing the wheel because one wants to is all ok, reinventing the wheel because all the previous one got lost is a necessity (there are no wheels) but also a pity (no one cares about leaving the wheels survive).

Pier, I think you already made your point several times now... ;)

[Allen]

(01-27-2018 01:23 PM)pier4r Wrote:  

(01-26-2018 10:01 PM)Allen Wrote:  I found a speed up and can now calculate N=1000 solution in under 20 seconds.
N=2000 takes <3 min
N=3500 about 16 min

I hope you share your findings somewhere...

I plan to.. I was out of town this weekend, but still working on some optimizations. The N=3500 graph now completes in 2 min 14 sec. I would like to get a little computer time finding the right pruning hyper-parameters before putting out more code, but the general solution is still the same.. just the pruning and generation is different.

Edit.. now N=3500 finishes in 1.3 seconds

[brickviking]

(01-28-2018 04:14 PM)Allen Wrote:  

(01-27-2018 01:23 PM)pier4r Wrote:  I hope you share your findings somewhere...

I plan to.. I was out of town this weekend, but still working on some optimizations. The N=3500 graph now completes in 2 min 14 sec. I would like to get a little computer time finding the right pruning hyper-parameters before putting out more code, but the general solution is still the same.. just the pruning and generation is different.

Edit.. now N=3500 finishes in 1.3 seconds

16 minutes down to 1.3 seconds. That's ... quite an improvement. I'm no mathematician, but even I can tell you've definitely managed something impressive.

(Post 163)

[peacecalc]

Hello folks,

I don't investigate very much time towards this problem. I beg your pardon if I wrote down a silly idea: If n = 299 then the biggest square number is 24^2 = (299 + 277), other possibilities (298 + 278) and so on. Maybe a initial program can generate all possible pairs of number which results in a square number.

Is known wether every square number from 1^2 .. 24^2 appears? The same square number can appear more then one time (can be seen in the examples) and it is logical, too. Because the numbers from 1 to 299 give 298 pairs of sum...

[Allen]

Finalizing my code for this puzzle. I made some modifications to make it run from the command line, and also slow it down.. The super-fast version I had earlier could not guarantee a solution. Takes 3 inputs:

Code:


 N  number of nodes on graph
min_clip - minimum candidates for each generation
max_clip - maximum candidates for each generation

I had a version that adjusted the results adaptively, but there are SO many possible solutions, it's not necessary.

In my fast version, similar to this one, I only take into account the max-clip during the last few steps. I had many many trial runs on the 10,000 or 20,000 node graph that had only a few numbers remaining (always 100<N<1).

So my improvement (slower but more accurate) is to update the scoring algorithm for each candidate solution so that it's not rewarded based on the centrality of the starting graph, but the centrality of the _remaining choices_. This made a substantial difference in number of solutions, but made the whole thing slower.

So if you want it FAST on the 1000 graph:
squares.py 1000 1 1000
(not guaranteed to find a solution..) adjust the 1 to 2 or 5 or 10

if you want lots of solutions for N=50:
squares.py 50 100 100000

( I think there are more than 250,000 unique solutions for N=50)

[Allen]

Code:


# -*- coding: utf-8 -*-
"""
Created on Wed Jan 24 16:26:01 2018

@author: allen_thomson
"""

#!/usr/bin/env python
import time
from math import sqrt as sqrt
import argparse
import random

def is_square(l,goal):
    squares=set([x**2 for x in range(2,int(sqrt(goal+goal-1))+1)])
    pairs=set([sum(x) for x in zip(l,l[1:])])
    return len(pairs-squares)==0
    
def extend(inputs,edges, scoring,goal,clip):
    output=[]
    for j in inputs:
        rev=random.choice([True,False])
        score,i=j
        left=set(i)
        next1=sorted(edges[i[-1]]-left,key=lambda x:-scoring[x])
        if clip==1 and len(next1)>1:
            if random.random()>.95:
                scoring={w:len(edges[w]-left) for w in range(1,goal+1)}
            #next1=sorted(edges[i[-1]]-left,key=lambda x:-scoring[x])
            output=[(score+scoring[next1[-1]],i+(next1[-1],))]
        else:
            scoring={w:len(edges[w]-left) for w in range(1,goal+1)}
            if rev:
                next1=edges[i[0]]-left
                output+=[(score+scoring[x],(x,)+i) for x in next1]
            else:
                next1=edges[i[-1]]-left
                output+=[(score+scoring[x],i+(x,)) for x in next1]
    return sorted(set(output),key=lambda x: -x[0])

def bottlenecks(goal,max_clip=100,min_clip=1):
    start_time=time.time()
    squares=set([x**2 for x in range(2,int(sqrt(goal+goal-1))+1)])
    edges={j:{x-j for x in [y for y in squares if y>j] if x-j<=goal} for j in range(1,goal+1)}
    scoring={j:len(edges[j]) for j in range(1,goal+1)}
    a=sorted([(scoring[y],(y,)) for y in edges.keys()],key=lambda x:-x[0])
    b=list(a)
    while len(b)>0:
        if len(set(b[0][1]))<1 or len(set(b[0][1]))+10>goal:
            clip=max_clip
        else:
            clip=min_clip
        print 'length {} : extending {} items with {} grams clip:{}'.format(len(set(b[0][1])),len(b),len(edges.keys()),clip)
        b=extend(set(a[-clip:]),edges,scoring,goal,clip)
        if len(b)>0:
            if len(b[0][1])==goal:
                print 'Success for len {}:{}'.format(goal,b[0])
                print 'Unique solutions found for len {}: {}'.format(goal,len(set(b+[x[::-1] for x in b])))

                return b
            a=list(b)
    print 'Runtime: {}'.format(time.time()-start_time)
    if len(a)>0:
        print 'Goal:{}'.format(goal)
        for i,j in enumerate(a):
            print 'sol {}: len {}     unique {}    is_square:{}   Success:{}'.format(i,len(j[1]),len(set(j[1])),is_square(j[1],goal),is_square(j[1],goal) and len(set(j[1]))==goal)
    return a

if __name__ == "__main__":
    import argparse
    parser = argparse.ArgumentParser(description='Calculate a sum-squares graph walk')
    parser.add_argument('N', type=int, help='number of nodes on the graph')
    parser.add_argument('min_clip', type=int, default=1, help='Smaller=Faster, Larger=more solutions')
    parser.add_argument('max_clip', type=int, default=100, help='Smaller=Faster, Larger=more solutions')
    args = parser.parse_args()
    b=bottlenecks(args.N,args.max_clip,args.min_clip)

[TheKaneB]

Nice and compact!
There are a lot of Python trickeries, like ranges, zip, array slices, etc... but it should be possible to port it to HPPPL.
I think you could cache the "squares" and maybe gain a bit more speed.

EDIT: nevermind, I did try it and the difference is negligible.

[pier4r]

(01-19-2018 04:29 AM)Paul Dale Wrote:  Baring an inspired approach, the problem comes down to finding a Hamilton path in a graph. This is a NP complete problem.


Pauli

It doesn't seem so though? An update!






-----

Also as reference:

some1: https://www.3blue1brown.com/blog/some1-results (unfortunately some links already do not work anymore)
some2: https://www.3blue1brown.com/blog/some2