Mainly Tech projects on Python and Electronic Design Automation.

Saturday, May 10, 2014

Sum of all natural numbers is minus one-twelfth

A curious infinite sum investigated with Python.

An investigation of the the following curious result:
∑k=1∞k=1+2+3+4+5+⋯=−112

Tracking the explanation given on Numberphile.
In [1]:
from IPython.display import YouTubeVideo
YouTubeVideo("w-I6XTVZXww")
Out[1]:

Series definitions

The video starts with an introduction, then at ~1:40 into the video he defines three series S1, S2 and S.
I would like to be able to define those same series in an as nearly as simple Python format, so defined a dummy class so I could concentrate on getting a nice description in Python.
Note: The class captures the series not the Sum of the series.
In [2]:
class Series(object):
    "placeholder"
    def __init__(self, **kwargs):
        pass
In [3]:
# Numberphile 1:40
S1 = Series(fk=lambda k: (-1)**(k - 1),   starts=[1, -1,  1, -1,  1, -1]) 
S2 = Series(fk=lambda k: (k*(-1)**(k - 1) if k > 1 else 1), 
                                          starts=[1, -2,  3, -4,  5, -6]) 
S  = Series(fk=lambda k: k,               starts=[1,  2,  3,  4,  5,  6]) 

Sum of the infinite series

At ~1:55 we need to find the sum of an infinite series by looking at the sum approaches.
I'll flesh out the Series class definition to handle what we have so far...
In [4]:
from fractions import Fraction
from itertools import islice, izip, count, chain

class Series(object):
    "Represents infinite series by its generating function"
    
    STARTSMAX = 10     # Limit on starts
    
    def __init__(self, fk, starts=[]):
        '''
        Arguments
            fk:     Generating function
            starts: Sample starting values of series. (checked)
        '''
        self.fk = fk
        assert all(fk(k) == start for k, start in enumerate(starts, 1)),\
               "Generating function does not create same starting values"
        self.starts = [fk(k) for k in range(1, self.STARTSMAX + 1)]
        
    def __call__(self):
        "Generate successive values of the series"
        k = 1
        while True:
            yield self.fk(k)
            k += 1
                       
        
In [5]:
# Numberphile 1:40
S1 = Series(fk=lambda k: (-1)**(k - 1),   starts=[1, -1,  1, -1,  1, -1]) 
S2 = Series(fk=lambda k: (k*(-1)**(k - 1) if k > 1 else 1), 
                                          starts=[1, -2,  3, -4,  5, -6]) 
S  = Series(fk=lambda k: k,               starts=[1,  2,  3,  4,  5,  6]) 
Lets doodle...
In [6]:
S1.starts
Out[6]:
[1, -1, 1, -1, 1, -1, 1, -1, 1, -1]
In [7]:
[x for n, x in zip(range(15), S2())]
Out[7]:
[1, -2, 3, -4, 5, -6, 7, -8, 9, -10, 11, -12, 13, -14, 15]
In [8]:
S.starts
Out[8]:
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
So you can see that .starts is a list of the first ten generated values, and calling an instance of the class gives a generator of successive values of the series (with no fixed upper limit).

Sum S1 = 1 + -1 +1 + -1 ... == 1/2

From 2:00 minutes into the video he shows that the sum of the terms of S1 alternate between 1 and zero
In [9]:
running_sums = [sum(islice(S1(), termcount)) for termcount in range(1,21)]
running_sums
Out[9]:
[1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0]
Lets take the infinite sum as the average of the sums i.e.
Sum S1 = 1 -1 +1 -1 +1 ... = 0.5

2 * (Sum S2) == Sum S1

From 2:45 in the video.
Lets take S2, and S2 with an initial 0 inserted . (Adding zeroes will not affect the sum of the series remember).
In [10]:
S2.starts
Out[10]:
[1, -2, 3, -4, 5, -6, 7, -8, 9, -10]
In [11]:
# S2 with an initial +0
[0] + S2.starts
Out[11]:
[0, 1, -2, 3, -4, 5, -6, 7, -8, 9, -10]
In [12]:
added_terms = (x + y  for x, y in izip(S2(), chain([0], S2())))
list(islice(added_terms, 10))
Out[12]:
[1, -1, 1, -1, 1, -1, 1, -1, 1, -1]
So, for the purposes of summation, 2 * S2 == S1 = 0.5
And therefore Sum S2 = 0.25

S - S2

From 4:08 in Video.
In [13]:
S.starts
Out[13]:
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
In [14]:
S2.starts
Out[14]:
[1, -2, 3, -4, 5, -6, 7, -8, 9, -10]
Lets do the subtraction by subtracting the generating functions this time.
In [15]:
S_sub_S2 = Series(fk=lambda k: S.fk(k) - S2.fk(k))
S_sub_S2.starts
Out[15]:
[0, 4, 0, 8, 0, 12, 0, 16, 0, 20]
Take a factor of four out and ignore the zeroes (as they do nothing for infinite sums).
In [16]:
list(islice((x/4 for x in S_sub_S2() if x), 10))
Out[16]:
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10]

End run

So:
sum(S) - sum(S2) = 4*(1+2+3+4+...) = 4 * ( sum(S) )
Solving:
sum(S) - 1/4 = 4*sum(S)

- 1/4 = 3*sum(S)

- 1/(4*3) = sum(S)

sum(S) = -1/12
Proving that the infinite sum 1+2+3+4+... = -1/12

The proof is in the eating

Remember that although it is a weird result, when used in scientific experiments we can predict nature.
For me, I now know that when dealing with equations involving infinite series I will annotate then with "There be dragons".

Tuesday, May 06, 2014

Python PyPy interpreter in Javascript: Firefox beats CPython beats Chrome!

Ryan Kelly posted info on PyPy.js an implementation of Python written in Javascript. available here. 
I tried it on Chrome Version 34.0.1847.131 m (the latest). using the pystone benchmark:


In cPython 3.4.0b1  I got:


And in Firefox 29.0 I got:

These were the latest versions at the time of writing (I think) and taking the fastest of their four runs on the same machine gives:

PlatformBest of 4 PystoneMultiplier
Chrome14013.5100%
cPython37424.7267%
Firefox66934.4478%
See the original article for caveats on the state of PyPy.js - it may only run pystone for all I know.

Of course it doesn't say much about Python speed but it may say something about Chrome vs Firefox speed? Who knows - lies, damned lies, and statistics!

Sunday, March 30, 2014

Nonogram puzzle solver in 2D (part 2/2).

Nonogram puzzle solver in 2D

My previous post got us up to solving constraints in a one dimensional nonogram vector.
You really need to have read the previous post recently as this follows on from that.
In [9]:
from __future__ import print_function, division
from pprint import pprint as pp
from copy import deepcopy
try: from itertools import izip_longest as zip_longest
except ImportError: from itertools import zip_longest
from nonovector import newforces, nonovector

# All the recap I will do for now:
help(nonovector)
help(newforces)
Help on function nonovector in module nonovector:

nonovector(blocks, cells, force0, force1)
    Generate all block positions for a nonogram row.
    
    Arguments:
        blocks:         list of block widths starting from lhs.
        cells:          Size of row.
        force0:         Cell indices that must be clear.
        force1:         Cell indices that must be filled.
    
    Generates:
        [ (bpos, blength), ...]
        
        Where:
            bpos is the minimum index to a span of connected blocks.
            blength is the number of consecutive blocks spanned.

Help on function newforces in module nonovector:

newforces(vecgen, cells)
    Extracts what cells are always zero/one in all vectors from generator vecgen
    
    Returns lists: (always0, always1)


Role reversal

I very early found out that I needed to be able to generate nonograms from pictures so I had something to solve so wrote pic2gram that can extract the clue information from an array of ones and zeroes
So with my training 2D example picture of
01110
11010
01110
00110
00111
I can generate a puzzle.
In [10]:
def line2blocks(line='01110011110011101111001111011100'):
    'Extract nonogram block clues from a line of 0/1 characters'
    blocks, runlength = [], 0
    for i, ch in enumerate(line):
        if ch == '1':
            runlength += 1
        elif runlength:         # ch != '1' so at end of block
            blocks.append(runlength)
            runlength = 0
    if runlength:               # block ends at RH border
        blocks.append(runlength)
    elif not blocks:            # Corner case: empty line depicted as one zero-length block
        blocks.append(0)
    return blocks                            

def pic2gram(pic):
    "convert rectangular '01' pic to nonogram clue"
    array = pic.strip().split('\n')
    vert  = [line2blocks(line) for line in array]
    horiz = [line2blocks(line) for line in zip(*array)] # Remember zip(*x) to transpose
    return (vert, horiz)
In [11]:
smallex = '''\
01110
11010
01110
00110
00111'''

verticalclues, horizontalclues = pic2gram(smallex)
print('verticalclues = %r\nhorizontalclues = %r' % (verticalclues, horizontalclues))
verticalclues = [[3], [2, 1], [3], [2], [3]]
horizontalclues = [[1], [3], [1, 3], [5], [1]]

Hmmm. Still need prettier output:
In [12]:
def parray(array, zro='_', one='#'):
    'prettyprint 0/1 array'
    for row in array:
        print('|' + '|'.join(one if cell in (1, '1') else zro
                             for cell in row) + '|')
    print('')

print('ORIGINAL PICTURE\n')
parray(smallex.split('\n'))
ORIGINAL PICTURE

|_|#|#|#|_|
|#|#|_|#|_|
|_|#|#|#|_|
|_|_|#|#|_|
|_|_|#|#|#|


That is better, now adding the generated clues to the printout:
In [13]:
def parrayclues(array, clues, zro='_', one='#'):
    'prettyprint 0/1 array with associated clues'
    vert, horiz = clues
    for row, v in zip(array, vert):
        print('|' + 
              '|'.join(one if cell in (1, '1') else zro for cell in row) + 
              '| # ' +
              ' '.join(str(i) for i in v))
    hstr = [' #' + ' '.join(str(b) for b in col) for col in horiz]
    verttext = [' '.join(ch) for ch in zip_longest(*hstr, fillvalue=' ')]
    for line in verttext:
        print(' ' + line)
In [14]:
vert, horiz = puzzle = pic2gram(pic=smallex)
array = smallex.split('\n')
print('ORIGINAL PICTURE WITH GENERATED NONOGRAM CLUES\n')
parrayclues(array, puzzle, one='X')
ORIGINAL PICTURE WITH GENERATED NONOGRAM CLUES

|_|X|X|X|_| # 3
|X|X|_|X|_| # 2 1
|_|X|X|X|_| # 3
|_|_|X|X|_| # 2
|_|_|X|X|X| # 3
          
 # # # # #
 1 3 1 5 1
          
     3    

The clues for columns are printed vertically downwards beneath the columns.

Solution Strategy

I can use newforces on each row to find the forced positions of each cell on the row and store them in two bit-arrays an array of bits denoting which cells must be zero and another bit-array for those that must be one. Bits in these bit arrays are 1 to force a value at that position. The value forced is blank for the forcing zero array or filled for the forcing 1 array.
newforces will load initial forces from the bit-arrays and store its updates back to them for each row of the grid.
If we spin the data 90 degrees, i.e. transpose things, we can run newforces on all the columns in a similar manner gradually refining the bit-array of forces until the bit-arrays reach a steady state.

Forces bit arrays

OK fix0grid and fix1grid aren't actually bits but they could be. I just don't need the space saving
In [15]:
# BLank initial forces bit-arrays
rowlen, collen = len(vert), len(horiz)
fix0grid = [[0] * collen for i in range(rowlen)] # Inner lists are distinct! 
fix1grid = [[0] * collen for i in range(rowlen)]
pp(fix0grid)
[[0, 0, 0, 0, 0],
 [0, 0, 0, 0, 0],
 [0, 0, 0, 0, 0],
 [0, 0, 0, 0, 0],
 [0, 0, 0, 0, 0]]

Pretty-printer for forced-cell info

We could do with another pretty-printer. This one formats with errors marked with err i.e. if some cell is forced to both a zero and a 1; and que (Fawlty towers ref.), if the cell currently is not forced to either state.
In [16]:
def pforces(f0grid, f1grid, clues, zro='_', one='#', err='!', que='?'):
    '''
    prettyprint to string from force0/1 arrays with associated clues
    
    Returns (string, , )
    
    '''
    vert, horiz = clues
    out = []
    for row0, row1, v in zip(f0grid, f1grid, vert):
        out.append(
              '|' + 
              '|'.join(one if cell1 and not cell0 
                         else (zro if cell0 and not cell1 else (
                                err if cell0 and cell1 else que)) 
                       for cell0, cell1 in zip(row0, row1)) + 
              '| # ' +
              ' '.join(str(i) for i in v))
    hstr = [' #' + ' '.join(str(b) for b in col) for col in horiz]
    verttext = [' '.join(ch) for ch in zip_longest(*hstr, fillvalue=' ')]
    for line in verttext:
        out.append(' ' + line)
    text = '\n'.join(out)
    return (text, text.count(err), text.count(que))
Lets try it:
In [17]:
# BLank initial forces bit-arrays
rowlen, collen = len(vert), len(horiz)
fix0grid = [[0] * collen for i in range(rowlen)] # Inner lists are distinct! 
fix1grid = [[0] * collen for i in range(rowlen)]

txt, err, que = pforces(fix0grid, fix1grid, puzzle)
print('Starting state of %i errors and %i unknown cells:\n' % (err, que))   
print(txt)
Starting state of 0 errors and 25 unknown cells:

|?|?|?|?|?| # 3
|?|?|?|?|?| # 2 1
|?|?|?|?|?| # 3
|?|?|?|?|?| # 2
|?|?|?|?|?| # 3
          
 # # # # #
 1 3 1 5 1
          
     3    

Update rows of fixed cell values

In [18]:
def newfixes(blocklist, fix0grid, fix1grid):
    '''
    In-line update to rows of fix0grid, fix1grid, based on current values.
    blocklist is the list of clues for all rows. 
    '''
    for row, (blocks, f0, f1) in enumerate(zip(blocklist, fix0grid, fix1grid)):
        # Extract forces for this row
        force0 = [i for i, cell in enumerate(f0) if cell]
        force1 = [i for i, cell in enumerate(f1) if cell]
        # Update forces for this row
        force0, force1 = newforces(nonovector(blocks, len(f0), force0, force1), len(f0))
        # Apply changes back to bit-arrays
        for col in force0:
            fix0grid[row][col] = 1
        for col in force1:
            fix1grid[row][col] = 1
    # remember, no return value as it works by updating fix0grid and fix1grid in-place.
No independent test of newfixes except for its use in the main loop

Main loop for nonogram solver

First I apply newfixes to vertical clues and rows of the bit-arrays to find fixes. Next, in the same iteration, I transpose the arrays and use the horizontal clues on what will then be the columns of the bit-arrays of forces.
I keep enough past state to detect when improvements stop.
In [19]:
iterations = 0
oldf0 = oldf1 = []      # Guaranteed to fail initial values
while (oldf0 != fix0grid or oldf1 != fix1grid) and (que or err):
    oldf0, oldf1 = deepcopy(fix0grid), deepcopy(fix1grid)   # Take a copy of current state.
    iterations += 1
    newfixes(vert, fix0grid, fix1grid)
    # transpose and do other direction
    fix0grid, fix1grid = list(list(z) for z in zip(*fix0grid)), list(list(z) for z in zip(*fix1grid))
    newfixes(horiz, fix0grid, fix1grid)
    # transpose back
    fix0grid, fix1grid = list(list(z) for z in zip(*fix0grid)), list(list(z) for z in zip(*fix1grid))
    txt, err, que = pforces(fix0grid, fix1grid, puzzle)
    print('\nIteration %i has %i errors and %i unknown cells\n' % (iterations, err, que))   
    print(txt)

print('\nFinished!')
Iteration 1 has 0 errors and 12 unknown cells

|?|?|#|#|?| # 3
|?|#|_|#|?| # 2 1
|?|#|#|#|?| # 3
|?|?|#|#|?| # 2
|?|_|#|#|?| # 3
          
 # # # # #
 1 3 1 5 1
          
     3    

Iteration 2 has 0 errors and 1 unknown cells

|_|#|#|#|?| # 3
|#|#|_|#|_| # 2 1
|_|#|#|#|_| # 3
|_|_|#|#|_| # 2
|_|_|#|#|#| # 3
          
 # # # # #
 1 3 1 5 1
          
     3    

Iteration 3 has 0 errors and 0 unknown cells

|_|#|#|#|_| # 3
|#|#|_|#|_| # 2 1
|_|#|#|#|_| # 3
|_|_|#|#|_| # 2
|_|_|#|#|#| # 3
          
 # # # # #
 1 3 1 5 1
          
     3    

Finished!

Remember that the original nonogram was:
In [20]:
parrayclues(array, puzzle, one='X')
|_|X|X|X|_| # 3
|X|X|_|X|_| # 2 1
|_|X|X|X|_| # 3
|_|_|X|X|_| # 2
|_|_|X|X|X| # 3
          
 # # # # #
 1 3 1 5 1
          
     3    

Looks OK to me!

Finished?

Sometimes the clues are not sufficient to generate a unique picture. I suspect that in those cases the above algorithm would only be a partial solution. It would need to be followed by iterating through all the possibilities of block positions in two dimensions covering the remaining unknown cells and returning the first iteration of block positions that fits the ending set of fix0grid and fix1grid values.

END.

Followers

Subscribe Now: google

Add to Google Reader or Homepage

Go deh too!

whos.amung.us

Blog Archive