title: Paint by Numbers

text:
# Question
Title: Paint by Numbers
You're given a true colour image. Your task is to generate a version of this image, which looks like it was painted using paint-by-numbers (the children's activity, *not* nonograms). Along with the image, you're given two parameters: **P**, the maximum size of the colour palette (i.e. the maximum number of distinct colours to use), and **N**, the maximum number of *cells* to use. Your algorithm does *not* have to use all **P** colours and **N** cells, but it must not use more than that. The output image should have the same dimensions as the input.

A *cell* is defined as a contiguous area of pixels which all have the same colour. Pixels touching only at a corner are *not* considered contiguous. Cells may have holes.

In short, you are to approximate the input image with only N flat-shaded/solid-colour areas and P different colours.

Just to visualise the parameters, here is a very simple example (for no particular input image; showing off my mad Paint skills). The following image has **P = 6** and **N = 11**:

Here are a few images to test your algorithm on (mostly our usual suspects). Click the pictures for larger versions.

Please include a handful of results for different parameters. If you want to show a large number of results, you can create a gallery over on imgur.com, to keep the size of the answers reasonable. Alternatively, put thumbnails in your post and make them links to larger images, like I did above. Also, feel free to use other test images, if you find something nice.

I assume that parameters around **N ≥ 500**, **P ~ 30** would be similar to real paint-by-number templates.

This is a popularity contest, so the answer with the most net votes wins. Voters are encouraged to judge answers by

* how well the original images are approximated.
* how well the algorithm works on different kinds of images (paintings are probably generally easier than photographs).
* how well the algorithm works with very restrictive parameters.
* how organic/smooth the cell shapes look.

I will use the following Mathematica script, to validate results:

```
image = <pastedimagehere> // ImageData;
palette = Union[Join @@ image];
Print["P = ", Length@palette];
grid = GridGraph[Reverse@Most@Dimensions@image];
image = Flatten[image /. Thread[palette -> Range@Length@palette]];
Print["N = ",
 Length@ConnectedComponents[
   Graph[Cases[EdgeList[grid],
     m_ <-> n_ /; image[[m]] == image[[n]]]]]]

```

Sp3000 was kind enough to write a verifier in Python 2 using PIL, which you find at this pastebin.

# Answer
# Python 2 with PIL (Gallery)

```
from __future__ import division
from PIL import Image
import random, math, time
from collections import Counter, defaultdict, namedtuple

"""
Configure settings here
"""

INFILE = "spheres.png"
OUTFILE_STEM = "out"
P = 30
N = 300
OUTPUT_ALL = True # Whether to output the image at each step

FLOOD_FILL_TOLERANCE = 10
CLOSE_CELL_TOLERANCE = 5
SMALL_CELL_THRESHOLD = 10
FIRST_PASS_N_RATIO = 1.5
K_MEANS_TRIALS = 30
BLUR_RADIUS = 2
BLUR_RUNS = 3

"""
Color conversion functions
"""

X = xrange

# http://www.easyrgb.com/?X=MATH
def rgb2xyz(rgb):
 r,g,b=rgb;r/=255;g/=255;b/=255;r=((r+0.055)/1.055)**2.4 if r>0.04045 else r/12.92
 g=((g+0.055)/1.055)**2.4 if g>0.04045 else g/12.92;b=((b+0.055)/1.055)**2.4 if b>0.04045 else b/12.92
 r*=100;g*=100;b*=100;x=r*0.4124+g*0.3576+b*0.1805;y=r*0.2126+g*0.7152+b*0.0722
 z=r*0.0193+g*0.1192+b*0.9505;return(x,y,z)
def xyz2lab(xyz):
 x,y,z=xyz;x/=95.047;y/=100;z/=108.883;x=x**(1/3)if x>0.008856 else 7.787*x+16/116
 y=y**(1/3)if y>0.008856 else 7.787*y+16/116;z=z**(1/3)if z>0.008856 else 7.787*z + 16/116
 L=116*y-16;a=500*(x-y);b=200*(y-z);return(L,a,b)
def rgb2lab(rgb):return xyz2lab(rgb2xyz(rgb))
def lab2xyz(lab):
 L,a,b=lab;y=(L+16)/116;x=a/500+y;z=y-b/200;y=y**3 if y**3>0.008856 else(y-16/116)/7.787
 x=x**3 if x**3>0.008856 else (x-16/116)/7.787;z=z**3 if z**3>0.008856 else(z-16/116)/7.787
 x*=95.047;y*=100;z*=108.883;return(x,y,z)
def xyz2rgb(xyz):
 x,y,z=xyz;x/=100;y/=100;z/=100;r=x*3.2406+y*-1.5372+z*-0.4986
 g=x*-0.9689+y*1.8758+z*0.0415;b=x*0.0557+y*-0.2040+z*1.0570
 r=1.055*(r**(1/2.4))-0.055 if r>0.0031308 else 12.92*r;g=1.055*(g**(1/2.4))-0.055 if g>0.0031308 else 12.92*g
 b=1.055*(b**(1/2.4))-0.055 if b>0.0031308 else 12.92*b;r*=255;g*=255;b*=255;return(r,g,b)
def lab2rgb(lab):rgb=xyz2rgb(lab2xyz(lab));return tuple([int(round(x))for x in rgb])

"""
Stage 1: Read in image and convert to CIELAB
"""

total_time = time.time()

im = Image.open(INFILE)
width, height = im.size

if OUTPUT_ALL:
  im.save(OUTFILE_STEM + "0.png")
  print "Saved image %s0.png" % OUTFILE_STEM

def make_pixlab_map(im):
  width, height = im.size
  pixlab_map = {}

  for i in X(width):
    for j in X(height):
      pixlab_map[(i, j)] = rgb2lab(im.getpixel((i, j)))

  return pixlab_map

pixlab_map = make_pixlab_map(im)

print "Stage 1: CIELAB conversion complete"

"""
Stage 2: Partitioning the image into like-colored cells using flood fill
"""

def d(color1, color2):
  return (abs(color1[0]-color2[0])**2 + abs(color1[1]-color2[1])**2 + abs(color1[2]-color2[2])**2)**.5

def neighbours(pixel):
  results = []

  for neighbour in [(pixel[0]+1, pixel[1]), (pixel[0]-1, pixel[1]),
            (pixel[0], pixel[1]+1), (pixel[0], pixel[1]-1)]:

    if 0 <= neighbour[0] < width and 0 <= neighbour[1] < height:
      results.append(neighbour)

  return results

def flood_fill(start_pixel):
  to_search = {start_pixel}
  cell = set()
  searched = set()
  start_color = pixlab_map[start_pixel]

  while to_search:
    pixel = to_search.pop()

    if d(start_color, pixlab_map[pixel]) < FLOOD_FILL_TOLERANCE:
      cell.add(pixel)
      unplaced_pixels.remove(pixel)

      for n in neighbours(pixel):
        if n in unplaced_pixels and n not in cell and n not in searched:
          to_search.add(n)

    else:
      searched.add(pixel)

  return cell

# These two maps are inverses, pixel/s <-> number of cell containing pixel
cell_sets = {}
pixcell_map = {}
unplaced_pixels = {(i, j) for i in X(width) for j in X(height)}

while unplaced_pixels:
  start_pixel = unplaced_pixels.pop()
  unplaced_pixels.add(start_pixel)
  cell = flood_fill(start_pixel)

  cellnum = len(cell_sets)
  cell_sets[cellnum] = cell

  for pixel in cell:
    pixcell_map[pixel] = cellnum

print "Stage 2: Flood fill partitioning complete, %d cells" % len(cell_sets)

"""
Stage 3: Merge cells with less than a specified threshold amount of pixels to reduce the number of cells
     Also good for getting rid of some noise
"""

def mean_color(cell, color_map):
  L_sum = 0
  a_sum = 0
  b_sum = 0

  for pixel in cell:
    L, a, b = color_map[pixel]
    L_sum += L
    a_sum += a
    b_sum += b

  return L_sum/len(cell), a_sum/len(cell), b_sum/len(cell)

def remove_small(cell_size):
  if len(cell_sets) <= N:
    return

  small_cells = []

  for cellnum in cell_sets:
    if len(cell_sets[cellnum]) <= cell_size:
      small_cells.append(cellnum)

  for cellnum in small_cells:
    neighbour_cells = []

    for cell in cell_sets[cellnum]:
      for n in neighbours(cell):
        neighbour_reg = pixcell_map[n]

        if neighbour_reg != cellnum:
          neighbour_cells.append(neighbour_reg)

    closest_cell = max(neighbour_cells, key=neighbour_cells.count)

    for cell in cell_sets[cellnum]:
      pixcell_map[cell] = closest_cell

    if len(cell_sets[closest_cell]) <= cell_size:
      small_cells.remove(closest_cell)

    cell_sets[closest_cell] |= cell_sets[cellnum]
    del cell_sets[cellnum]

    if len(cell_sets) <= N:
      return

for cell_size in X(1, SMALL_CELL_THRESHOLD):
  remove_small(cell_size)

if OUTPUT_ALL:
  frame_im = Image.new("RGB", im.size)

  for cellnum in cell_sets:
    cell_color = mean_color(cell_sets[cellnum], pixlab_map)

    for pixel in cell_sets[cellnum]:
      frame_im.putpixel(pixel, lab2rgb(cell_color))

  frame_im.save(OUTFILE_STEM + "1.png")
  print "Saved image %s1.png" % OUTFILE_STEM

print "Stage 3: Small cell merging complete, %d cells" % len(cell_sets)

"""
Stage 4: Close color merging
"""

cell_means = {}

for cellnum in cell_sets:
  cell_means[cellnum] = mean_color(cell_sets[cellnum], pixlab_map)

n_graph = defaultdict(set)

for i in X(width):
  for j in X(height):
    pixel = (i, j)
    cell = pixcell_map[pixel]

    for n in neighbours(pixel):
      neighbour_cell = pixcell_map[n]

      if neighbour_cell != cell:
        n_graph[cell].add(neighbour_cell)
        n_graph[neighbour_cell].add(cell)

def merge_cells(merge_from, merge_to):
  merge_from_cell = cell_sets[merge_from]

  for pixel in merge_from_cell:
    pixcell_map[pixel] = merge_to

  del cell_sets[merge_from]
  del cell_means[merge_from]

  n_graph[merge_to] |= n_graph[merge_from]
  n_graph[merge_to].remove(merge_to)

  for n in n_graph[merge_from]:
    n_graph[n].remove(merge_from)

    if n != merge_to:
      n_graph[n].add(merge_to)

  del n_graph[merge_from]

  cell_sets[merge_to] |= merge_from_cell
  cell_means[merge_to] = mean_color(cell_sets[merge_to], pixlab_map)

# Go through the cells from largest to smallest. Keep replenishing the list while we can still merge.
last_time = time.time()
to_search = sorted(cell_sets.keys(), key=lambda x:len(cell_sets[x]), reverse=True)
full_list = True

while len(cell_sets) > N and to_search:
  if time.time() - last_time > 15:
    last_time = time.time()
    print "Close color merging... (%d cells remaining)" % len(cell_sets)

  while to_search:
    cellnum = to_search.pop()
    close_cells = []

    for neighbour_cellnum in n_graph[cellnum]:
      if d(cell_means[cellnum], cell_means[neighbour_cellnum]) < CLOSE_CELL_TOLERANCE:
        close_cells.append(neighbour_cellnum)

    if close_cells:
      for neighbour_cellnum in close_cells:
        merge_cells(neighbour_cellnum, cellnum)

        if neighbour_cellnum in to_search:
          to_search.remove(neighbour_cellnum)

      break

  if full_list == True:
    if to_search:
      full_list = False

  else:
    if not to_search:
      to_search = sorted(cell_sets.keys(), key=lambda x:len(cell_sets[x]), reverse=True)
      full_list = True

if OUTPUT_ALL:
  frame_im = Image.new("RGB", im.size)

  for cellnum in cell_sets:
    cell_color = cell_means[cellnum]

    for pixel in cell_sets[cellnum]:
      frame_im.putpixel(pixel, lab2rgb(cell_color))

  frame_im.save(OUTFILE_STEM + "2.png")
  print "Saved image %s2.png" % OUTFILE_STEM

print "Stage 4: Close color merging complete, %d cells" % len(cell_sets)

"""
Stage 5: N-merging - merge until <= N cells
     Want to merge either 1) small cells or 2) cells close in color
"""

# Weight score between neighbouring cells by 1) size of cell and 2) color difference
def score(cell1, cell2):
  return d(cell_means[cell1], cell_means[cell2]) * len(cell_sets[cell1])**.5

n_scores = {}

for cellnum in cell_sets:
  for n in n_graph[cellnum]:
    n_scores[(n, cellnum)] = score(n, cellnum)

last_time = time.time()

while len(cell_sets) > N * FIRST_PASS_N_RATIO:
  if time.time() - last_time > 15:
    last_time = time.time()
    print "N-merging... (%d cells remaining)" % len(cell_sets)

  merge_from, merge_to = min(n_scores, key=lambda x: n_scores[x])

  for n in n_graph[merge_from]:
    del n_scores[(merge_from, n)]
    del n_scores[(n, merge_from)]

  merge_cells(merge_from, merge_to)

  for n in n_graph[merge_to]:
    n_scores[(n, merge_to)] = score(n, merge_to)
    n_scores[(merge_to, n)] = score(merge_to, n)

if OUTPUT_ALL:
  frame_im = Image.new("RGB", im.size)

  for cellnum in cell_sets:
    cell_color = cell_means[cellnum]

    for pixel in cell_sets[cellnum]:
      frame_im.putpixel(pixel, lab2rgb(cell_color))

  frame_im.save(OUTFILE_STEM + "3.png")
  print "Saved image %s3.png" % OUTFILE_STEM

del n_graph, n_scores

print "Stage 5: N-merging complete, %d cells" % len(cell_sets)

"""
Stage 6: P merging - use k-means
"""

def form_clusters(centroids):
  clusters = defaultdict(set)

  for cellnum in cell_sets:
    # Add cell to closest centroid.
    scores = []

    for centroid in centroids:
      scores.append((d(centroid, cell_means[cellnum]), centroid))

    scores.sort()
    clusters[scores[0][1]].add(cellnum)

  return clusters

def calculate_centroid(cluster):
  L_sum = 0
  a_sum = 0
  b_sum = 0

  weighting = 0

  for cellnum in cluster:
    # Weight based on cell size
    color = cell_means[cellnum]
    cell_weight = len(cell_sets[cellnum])**.5

    L_sum += color[0]*cell_weight
    a_sum += color[1]*cell_weight
    b_sum += color[2]*cell_weight

    weighting += cell_weight

  return (L_sum/weighting, a_sum/weighting, b_sum/weighting)

def db_index(clusters):
  # Davies-Bouldin index
  scatter = {}

  for centroid, cluster in clusters.items():
    scatter_score = 0

    for cellnum in cluster:
      scatter_score += d(cell_means[cellnum], centroid) * len(cell_sets[cellnum])**.5

    scatter_score /= len(cluster)
    scatter[centroid] = scatter_score**2 # Mean squared distance

  index = 0

  for ci, cluster in clusters.items():
    dist_scores = []

    for cj in clusters:
      if ci != cj:
        dist_scores.append((scatter[ci] + scatter[cj])/d(ci, cj))

    index += max(dist_scores)

  return index

best_clusters = None
best_index = None

for i in X(K_MEANS_TRIALS):
  centroids = {cell_means[cellnum] for cellnum in random.sample(cell_sets, P)}
  converged = False

  while not converged:
    clusters = form_clusters(centroids)
    new_centroids = {calculate_centroid(cluster) for cluster in clusters.values()}

    if centroids == new_centroids:
      converged = True

    centroids = new_centroids

  index = db_index(clusters)

  if best_index is None or index < best_index:
    best_index = index
    best_clusters = clusters

del cell_means
newpix_map = {}

for centroid, cluster in best_clusters.items():
  for cellnum in cluster:
    for pixel in cell_sets[cellnum]:
      newpix_map[pixel] = centroid

if OUTPUT_ALL:
  frame_im = Image.new("RGB", im.size)

  for pixel in newpix_map:
    frame_im.putpixel(pixel, lab2rgb(newpix_map[pixel]))

  frame_im.save(OUTFILE_STEM + "4.png")
  print "Saved image %s4.png" % OUTFILE_STEM

print "Stage 6: P-merging complete"

"""
Stage 7: Approximate Gaussian smoothing
     See http://blog.ivank.net/fastest-gaussian-blur.html
"""

# Hindsight tells me I should have used a class. I hate hindsight.
def vec_sum(vectors):
  assert(vectors and all(len(v) == len(vectors[0]) for v in vectors))
  return tuple(sum(x[i] for x in vectors) for i in X(len(vectors[0])))

def linear_blur(color_list):
  # Can be made faster with an accumulator
  output = []

  for i in X(len(color_list)):
    relevant_pixels = color_list[max(i-BLUR_RADIUS+1, 0):i+BLUR_RADIUS]
    pixsum = vec_sum(relevant_pixels)
    output.append(tuple(pixsum[i]/len(relevant_pixels) for i in X(3)))

  return output

def horizontal_blur():
  for row in X(height):
    colors = [blurpix_map[(i, row)] for i in X(width)]
    colors = linear_blur(colors)

    for i in X(width):
      blurpix_map[(i, row)] = colors[i]

def vertical_blur():
  for column in X(width):
    colors = [blurpix_map[(column, j)] for j in X(height)]
    colors = linear_blur(colors)

    for j in X(height):
      blurpix_map[(column, j)] = colors[j]

blurpix_map = {}

for i in X(width):
  for j in X(height):
    blurpix_map[(i, j)] = newpix_map[(i, j)]

for i in X(BLUR_RUNS):
  vertical_blur()
  horizontal_blur()

# Pixel : color of smoothed image
smoothpix_map = {}

for i in X(width):
  for j in X(height):
    pixel = (i, j)
    blur_color = blurpix_map[pixel]
    nearby_colors = {newpix_map[pixel]}

    for n in neighbours(pixel):
      nearby_colors.add(newpix_map[n])

    smoothpix_map[pixel] = min(nearby_colors, key=lambda x: d(x, blur_color))

del newpix_map, blurpix_map

if OUTPUT_ALL:
  frame_im = Image.new("RGB", im.size)

  for pixel in smoothpix_map:
    frame_im.putpixel(pixel, lab2rgb(smoothpix_map[pixel]))

  frame_im.save(OUTFILE_STEM + "5.png")
  print "Saved image %s5.png" % OUTFILE_STEM

print "Stage 7: Smoothing complete"

"""
Stage 8: Flood fill pass 2
     Code copy-and-paste because I'm lazy
"""

def flood_fill(start_pixel):
  to_search = {start_pixel}
  cell = set()
  searched = set()
  start_color = smoothpix_map[start_pixel]

  while to_search:
    pixel = to_search.pop()

    if start_color == smoothpix_map[pixel]:
      cell.add(pixel)
      unplaced_pixels.remove(pixel)

      for n in neighbours(pixel):
        if n in unplaced_pixels and n not in cell and n not in searched:
          to_search.add(n)

    else:
      searched.add(pixel)

  return cell

cell_sets = {}
pixcell_map = {}
unplaced_pixels = {(i, j) for i in X(width) for j in X(height)}

while unplaced_pixels:
  start_pixel = unplaced_pixels.pop()
  unplaced_pixels.add(start_pixel)
  cell = flood_fill(start_pixel)

  cellnum = len(cell_sets)
  cell_sets[cellnum] = cell

  for pixel in cell:
    pixcell_map[pixel] = cellnum

cell_colors = {}

for cellnum in cell_sets:
  cell_colors[cellnum] = smoothpix_map[next(iter(cell_sets[cellnum]))]

print "Stage 8: Flood fill pass 2 complete, %d cells" % len(cell_sets)

"""
Stage 9: Small cell removal pass 2
"""

def score(cell1, cell2):
  return d(cell_colors[cell1], cell_colors[cell2]) * len(cell_sets[cell1])**.5

def remove_small(cell_size):
  small_cells = []

  for cellnum in cell_sets:
    if len(cell_sets[cellnum]) <= cell_size:
      small_cells.append(cellnum)

  for cellnum in small_cells:
    neighbour_cells = []

    for cell in cell_sets[cellnum]:
      for n in neighbours(cell):
        neighbour_reg = pixcell_map[n]

        if neighbour_reg != cellnum:
          neighbour_cells.append(neighbour_reg)

    closest_cell = max(neighbour_cells, key=neighbour_cells.count)

    for cell in cell_sets[cellnum]:
      pixcell_map[cell] = closest_cell

    if len(cell_sets[closest_cell]) <= cell_size:
      small_cells.remove(closest_cell)

    cell_color = cell_colors[closest_cell]

    for pixel in cell_sets[cellnum]:
      smoothpix_map[pixel] = cell_color

    cell_sets[closest_cell] |= cell_sets[cellnum]
    del cell_sets[cellnum]
    del cell_colors[cellnum]

for cell_size in X(1, SMALL_CELL_THRESHOLD):
  remove_small(cell_size)

if OUTPUT_ALL:
  frame_im = Image.new("RGB", im.size)

  for pixel in smoothpix_map:
    frame_im.putpixel(pixel, lab2rgb(smoothpix_map[pixel]))

  frame_im.save(OUTFILE_STEM + "6.png")
  print "Saved image %s6.png" % OUTFILE_STEM

print "Stage 9: Small cell removal pass 2 complete, %d cells" % len(cell_sets)

"""
Stage 10: N-merging pass 2
     Necessary as stage 7 might generate *more* cells
"""

def merge_cells(merge_from, merge_to):
  merge_from_cell = cell_sets[merge_from]

  for pixel in merge_from_cell:
    pixcell_map[pixel] = merge_to

  del cell_sets[merge_from]
  del cell_colors[merge_from]

  n_graph[merge_to] |= n_graph[merge_from]
  n_graph[merge_to].remove(merge_to)

  for n in n_graph[merge_from]:
    n_graph[n].remove(merge_from)

    if n != merge_to:
      n_graph[n].add(merge_to)

  del n_graph[merge_from]

  cell_color = cell_colors[merge_to]

  for pixel in merge_from_cell:
    smoothpix_map[pixel] = cell_color

  cell_sets[merge_to] |= merge_from_cell

n_graph = defaultdict(set)

for i in X(width):
  for j in X(height):
    pixel = (i, j)
    cell = pixcell_map[pixel]

    for n in neighbours(pixel):
      neighbour_cell = pixcell_map[n]

      if neighbour_cell != cell:
        n_graph[cell].add(neighbour_cell)
        n_graph[neighbour_cell].add(cell)

n_scores = {}

for cellnum in cell_sets:
  for n in n_graph[cellnum]:
    n_scores[(n, cellnum)] = score(n, cellnum)

last_time = time.time()

while len(cell_sets) > N:
  if time.time() - last_time > 15:
    last_time = time.time()
    print "N-merging (pass 2)... (%d cells remaining)" % len(cell_sets)

  merge_from, merge_to = min(n_scores, key=lambda x: n_scores[x])

  for n in n_graph[merge_from]:
    del n_scores[(merge_from, n)]
    del n_scores[(n, merge_from)]

  merge_cells(merge_from, merge_to)

  for n in n_graph[merge_to]:
    n_scores[(n, merge_to)] = score(n, merge_to)
    n_scores[(merge_to, n)] = score(merge_to, n)

print "Stage 10: N-merging pass 2 complete, %d cells" % len(cell_sets)

"""
Stage last: Output the image!
"""

test_im = Image.new("RGB", im.size)

for i in X(width):
  for j in X(height):
    test_im.putpixel((i, j), lab2rgb(smoothpix_map[(i, j)]))

if OUTPUT_ALL:
  test_im.save(OUTFILE_STEM + "7.png")
else:
  test_im.save(OUTFILE_STEM + ".png")

print "Done! (Time taken: {})".format(time.time() - total_time)

```

Update time! This update features a simple smoothing algorithm to make images look less fuzzy. If I update again I'll have to revamp a fair chunk of my code though, because it's getting messy & I hd 2 glf a fw thngs 2 mke t char lim.

I've also made k-means weight colours based on cell sizes, which loses some details for more restrictive parameters (e.g. the centre of the nebula and American Gothic's pitchfork) but makes the overall colour choice sharper and nicer. Interestingly, it loses the whole background for raytraced spheres for P = 5.

Algorithm summary:

1. **Convert the pixels to the CIELAB colour space**: CIELAB approximates human vision better than RGB. Originally I used HSL (hue, saturation, lightness) but this had two problems — hue of white/grey/black is undefined, and hue is measured in degrees which wrap around, making k-means difficult to use.
2. **Divide the image into like-coloured cells using flood fill:** Pick a pixel not in a cell and do a flood fill using a specified tolerance. To measure the distance between two colours I use the standard Euclidean norm. More complicated formulae are available on this wiki article.
3. **Merge together small cells with their neighbours**: The flood fill generates a lot of 1 or 2 pixel cells — merge cells less than a specified size with the neighbouring cell with the most adjacent pixels. This considerably reduces the number of cells, improving running time for later steps.
4. **Merge together similarly-coloured regions**: Go through the cells in order of decreasing size. If any neighbouring cell has mean colour less than a certain distance away, merge the cells. Keep going through the cells until no more can be merged.
5. **Merge until we have less than 1.5N cells (N-merging)**: Merge cells together, using a scoring based on cell size and colour difference, until we have at most 1.5N cells. We allow a bit of leeway as we'll merge again later.
6. **Merge until we have less than P colours, using k-means (P-merging)**: Use the k-means clustering algorithm some specified number of times to generate clusterings of cell colours, weighting based on cell size. Score each clustering based on a variation of the Davies-Bouldin index and pick the best clustering to use.
7. **Approximate Gaussian smoothing**: Use several linear blurs to approximate Gaussian blurring (details here). Then for each pixel, pick the colour out of itself and its neighbours in the pre-blurred image closest to its colour in the blurred image. This part can be optimised more time-wise if necessary as I have yet to implement the optimal algorithm.
8. **Do another flood fill pass to work out the new regions**: This is necessary as the previous step may actually generate *more* cells.
9. **Do another small cell merging pass**
10. **Do another N-merging pass**: This time we go down to N cells rather than 1.5N.

The processing time for each image highly depends on its size and complexity, with times ranging from 20 seconds to 7 minutes for the test images.

Because the algorithm uses randomisation (e.g. merging, k-means), you can get different results on different runs. Here's a comparison of two runs for the bear image, with N = 50 and P = 10:

---

*Note: All images below are links. Most of these images are straight from the first run, but if I didn't like the output I allowed myself up to three attempts to be fair.*

## N = 50, P = 10

## N = 500, P = 30

But I'm pretty lazy when it comes to paint by colours, so just for fun...

## N = 20, P = 5

Additionally, it's amusing to see what happens when you try to squeeze 1 million colors into N = 500, P = 30:

Here's a step-by-step walkthrough of the algorithm for the underwater image with N = 500 and P = 30, in animated GIF form:

---

I've also made a gallery for the previous versions of the algorithm here. Here's some of my favourites from the last version (from when the nebula had more stars and the bear looked furrier):
> 39 votes

# Answer
# Python 2 with PIL

Also a Python solution and probably very much a work in progress:

```
from PIL import Image, ImageFilter
import random

def draw(file_name, P, N, M=3):
    img = Image.open(file_name, 'r')
    pixels = img.load()
    size_x, size_y = img.size

    def dist(c1, c2):
        return (c1[0]-c2[0])**2+(c1[1]-c2[1])**2+(c1[2]-c2[2])**2

    def mean(colours):
        n = len(colours)
        r = sum(c[0] for c in colours)//n
        g = sum(c[1] for c in colours)//n
        b = sum(c[2] for c in colours)//n
        return (r,g,b)

    def colourize(colour, palette):
        return min(palette, key=lambda c: dist(c, colour))

    def cluster(colours, k, max_n=10000, max_i=10):
        colours = random.sample(colours, max_n)
        centroids = random.sample(colours, k)
        i = 0
        old_centroids = None
        while not(i>max_i or centroids==old_centroids):
            old_centroids = centroids
            i += 1
            labels = [colourize(c, centroids) for c in colours]
            centroids = [mean([c for c,l in zip(colours, labels)
                               if l is cen]) for cen in centroids]
        return centroids

    all_coords = [(x,y) for x in xrange(size_x) for y in xrange(size_y)]
    all_colours = [pixels[x,y] for x,y in all_coords]
    palette = cluster(all_colours, P)
    print 'clustered'

    for x,y in all_coords:
        pixels[x,y] = colourize(pixels[x,y], palette)
    print 'colourized'

    median_filter = ImageFilter.MedianFilter(size=M)
    img = img.filter(median_filter)
    pixels = img.load()
    for x,y in all_coords:
        pixels[x,y] = colourize(pixels[x,y], palette)
    print 'median filtered'

    def neighbours(edge, outer, colour=None):
        return set((x+a,y+b) for x,y in edge
                   for a,b in ((1,0), (-1,0), (0,1), (0,-1))
                   if (x+a,y+b) in outer
                   and (colour==None or pixels[(x+a,y+b)]==colour))

    def cell(centre, rest):
        colour = pixels[centre]
        edge = set([centre])
        region = set()
        while edge:
            region |= edge
            rest = rest-edge
            edge = set(n for n in neighbours(edge, rest, colour))
        return region, rest

    print 'start segmentation:'
    rest = set(all_coords)
    cells = []
    while rest:
        centre = random.sample(rest, 1)[0]
        region, rest = cell(centre, rest-set(centre))
        cells += [region]
        print '%d pixels remaining'%len(rest)
    cells = sorted(cells, key=len, reverse=True)
    print 'segmented (%d segments)'%len(cells)

    print 'start merging:'
    while len(cells)>N:
        small_cell = cells.pop()
        n = neighbours(small_cell, set(all_coords)-small_cell)
        for big_cell in cells:
            if big_cell & n:
                big_cell |= small_cell
                break
        print '%d segments remaining'%len(cells)
    print 'merged'

    for cell in cells:
        colour = colourize(mean([pixels[x,y] for x,y in cell]), palette)
        for x,y in cell:
            pixels[x,y] = colour
    print 'colorized again'

    img.save('P%d N%d '%(P,N)+file_name)
    print 'saved'

draw('a.png', 11, 500, 1)

```

The algorithm follows a different approach than SP3000's, starting with colours first:

* Find a colour palette of *P* colours by k-means clustering and paint the image in this reduced palette.
* Apply a slight median filter to get rid of some noise.
* Make a list of all monochromatic cells and sort it by size.
* Merge the smallest cells with their respective largest neighbour until there are only *N* cells left.

There is quite some room for improvement, both in terms of speed and quality of the results. Especially the cell merging step can take up to many minutes and gives far from optimal results.

---

## *P* = 5, *N* = 45

## *P* = 10, *N* = 50

## *P* = 4, *N* = 250

## *P* = 11, *N* = 500
> 15 votes

# Answer
# Mathematica

At the moment, this takes the number of colors and the Gaussian radius to be used in the Gaussian filter. The larger the radius, the greater the blurring and merging of colors.

Because it does not allow for input of the number of cells, it doesn't meet one of the basic requirements of the challenge.

Output includes the number of cells for each color and also the total number of cells.

```
quantImg[img_,nColours_,gaussR_]:=ColorQuantize[GaussianFilter[img,gaussR],nColours,
Dithering-> False]

colours[qImg_]:=Union[Flatten[ImageData[qImg],1]]

showColors[image_,nColors_,gaussR_]:=
   Module[{qImg,colors,ca,nCells},
   qImg=quantImg[image,nColors,gaussR];
   colors=colours[qImg];
   ca=ConstantArray[0,Reverse@ImageDimensions[image]];
   nCells[qImgg_,color_]:=
   Module[{r},
   r=ReplacePart[ca,Position[ImageData@qImg,color]/.{a_,b_}:> ({a,b}->1)];
   (*ArrayPlot[r,ColorRules->{1\[Rule]RGBColor[color],0\[Rule]White}];*)
   m=MorphologicalComponents[r];
   {RGBColor@color,Max[Union@Flatten[m,1]]}];
   s=nCells[qImg,#]&/@colors;
   Grid[{
    {Row[{s}]}, {Row[{"cells:\t\t",Tr[s[[All,2]]]}]},{Row[{"colors:\t\t",nColors}]},
    {Row[{"Gauss. Radius: ", gaussR}]}},Alignment->Left]]

```

---

## Update

`quantImage2` allows to specify the desired number of cells as input. It determines the a best Gaussian Radius by looping through scenarios with greater radii until a close match is found.

`quantImage2` outputs (picture, requested cells, used cells, error, gaussian Radius used).

It is, however, very slow. To save time, you may start with an initial radius, the default value of which is 0.

```
gaussianRadius[img_,nCol_,nCells_,initialRadius_:0]:=
Module[{radius=initialRadius,nc=10^6,results={},r},
While[nc>nCells,(nc=numberOfCells[ape,nColors,radius]);
results=AppendTo[results,{nColors,radius,nc}];radius++];
r=results[[{-2,-1}]];
Nearest[r[[All,3]],200][[1]];
Cases[r,{_,_,Nearest[r[[All,3]],nCells][[1]]}][[1,2]]
]

quantImg2[img_,nColours_,nCells1_,initialRadius_:0]:={ColorQuantize[GaussianFilter[img,
g=gaussianRadius[img,nColours,nCells1,initialRadius]],nColours,Dithering->False],
nCells1,nn=numberOfCells[img,nColours,g],N[(nn-nCells1)/nCells1],g}

```

---

## Example for which we specify the number of cells desired in the output.

Example requesting 90 cells with 25 colors. Solution returns 88 cells, 2% error. The function chose the Gaussian radius of 55. (Lots of distortion).

---

## Examples for which the input includes the Gaussian radius, but not the number of cells.

**25 Colors, Gaussian radius of 5 pixels**

```
nColors = 25;
gR = 5;
quantImg[balls, nColors, gR]

```

---

**Three Colors, radius of 17 pixels**

```
nColors=3;gaussianRadius=17;
showColors[wave,nColors,gaussianRadius]
quantImg[wave,nColors,gaussianRadius]

```

---

**Twenty colors, radius of 17 pixels**

We increased the number of colors but not the focus. Note the increase in the number of cells.

---

**Six Colors, radius of 4 pixels**

```
nColors=6;gaussianRadius=4;
showColors[wave,nColors,gaussianRadius]
quantImg[wave,nColors,gaussianRadius]

```

---

```
nColors = 6; gaussianRadius = 17;
showColors[ape, nColors, gaussianRadius]
quantImg[ape, nColors, gaussianRadius]

```

---

```
nColors = 6; gaussianRadius = 3;
showColors[ape, nColors, gaussianRadius]
quantImg[ape, nColors, gaussianRadius]

```

---

**Starry Night**

With only 6 colors and 60 cells. There is a color mismatch in the colors the `showColors` claims it uses. (Yellow doesn't appear among the 5 colors but it is used in the drawing.) I'll see if I can figure this out.
> 15 votes

# Answer
# Python 2 with PIL

This is still somewhat a work in progress. Also, the code below is a horrible mess of spaghetti, and should not be used as an inspiration. :)

```
from PIL import Image, ImageFilter
from math import sqrt
from copy import copy
from random import shuffle, choice, seed

IN_FILE = "input.png"
OUT_FILE = "output.png"

LOGGING = True
GRAPHICAL_LOGGING = False
LOG_FILE_PREFIX = "out"
LOG_FILE_SUFFIX = ".png"
LOG_ROUND_INTERVAL = 150
LOG_FLIP_INTERVAL = 40000

N = 500
P = 30
BLUR_RADIUS = 3
FILAMENT_ROUND_INTERVAL = 5
seed(0) # Random seed

print("Opening input file...")

image = Image.open(IN_FILE).filter(ImageFilter.GaussianBlur(BLUR_RADIUS))
pixels = {}
width, height = image.size

for i in range(width):
    for j in range(height):
        pixels[(i, j)] = image.getpixel((i, j))

def dist_rgb((a,b,c), (d,e,f)):
    return (a-d)**2 + (b-e)**2 + (c-f)**2

def nbors((x,y)):
    if 0 < x:
        if 0 < y:
            yield (x-1,y-1)
        if y < height-1:
            yield (x-1,y+1)
    if x < width - 1:
        if 0 < y:
            yield (x+1,y-1)
        if y < height-1:
            yield (x+1,y+1)

def full_circ((x,y)):
    return ((x+1,y), (x+1,y+1), (x,y+1), (x-1,y+1), (x-1,y), (x-1,y-1), (x,y-1), (x+1,y-1))

class Region:

    def __init__(self):
        self.points = set()
        self.size = 0
        self.sum = (0,0,0)

    def flip_point(self, point):
        sum_r, sum_g, sum_b = self.sum
        r, g, b = pixels[point]
        if point in self.points:
            self.sum = (sum_r - r, sum_g - g, sum_b - b)
            self.size -= 1
            self.points.remove(point)
        else:
            self.sum = (sum_r + r, sum_g + g, sum_b + b)
            self.size += 1
            self.points.add(point)

    def mean_with(self, color):
        if color is None:
            s = float(self.size)
            r, g, b = self.sum
        else:
            s = float(self.size + 1)
            r, g, b = map(lambda a,b: a+b, self.sum, color)
        return (r/s, g/s, b/s)

print("Initializing regions...")

aspect_ratio = width / float(height)
a = int(sqrt(N)*aspect_ratio)
b = int(sqrt(N)/aspect_ratio)

num_components = a*b
owners = {}
regions = [Region() for i in range(P)]
borders = set()

nodes = [(i,j) for i in range(a) for j in range(b)]
shuffle(nodes)
node_values = {(i,j):None for i in range(a) for j in range(b)}

for i in range(P):
    node_values[nodes[i]] = regions[i]

for (i,j) in nodes[P:]:
    forbiddens = set()
    for node in (i,j-1), (i,j+1), (i-1,j), (i+1,j):
        if node in node_values and node_values[node] is not None:
            forbiddens.add(node_values[node])
    node_values[(i,j)] = choice(list(set(regions) - forbiddens))

for (i,j) in nodes:
    for x in range((width*i)/a, (width*(i+1))/a):
        for y in range((height*j)/b, (height*(j+1))/b):
            owner = node_values[(i,j)]
            owner.flip_point((x,y))
            owners[(x,y)] = owner

def recalc_borders(point = None):
    global borders
    if point is None:
        borders = set()
        for i in range(width):
            for j in range(height):
                if (i,j) not in borders:
                    owner = owner_of((i,j))
                    for pt in nbors((i,j)):
                        if owner_of(pt) != owner:
                            borders.add((i,j))
                            borders.add(pt)
                            break
    else:
        for pt in nbors(point):
            owner = owner_of(pt)
            for pt2 in nbors(pt):
                if owner_of(pt2) != owner:
                    borders.add(pt)
                    break
            else:
                borders.discard(pt)

def owner_of(point):
    if 0 <= point[0] < width and 0 <= point[1] < height:
        return owners[point]
    else:
        return None

# Status codes for analysis
SINGLETON = 0
FILAMENT = 1
SWAPPABLE = 2
NOT_SWAPPABLE = 3

def analyze_nbors(point):
    owner = owner_of(point)
    circ = a,b,c,d,e,f,g,h = full_circ(point)
    oa,ob,oc,od,oe,of,og,oh = map(owner_of, circ)
    nbor_owners = set([oa,oc,oe,og])
    if owner not in nbor_owners:
        return SINGLETON, owner, nbor_owners - set([None])
    if oc != oe == owner == oa != og != oc:
        return FILAMENT, owner, set([og, oc]) - set([None])
    if oe != oc == owner == og != oa != oe:
        return FILAMENT, owner, set([oe, oa]) - set([None])
    last_owner = oa
    flips = {last_owner:0}
    for (corner, side, corner_owner, side_owner) in (b,c,ob,oc), (d,e,od,oe), (f,g,of,og), (h,a,oh,oa):
        if side_owner not in flips:
            flips[side_owner] = 0
        if side_owner != corner_owner or side_owner != last_owner:
            flips[side_owner] += 1
            flips[last_owner] += 1
        last_owner = side_owner
    candidates = set(own for own in flips if flips[own] == 2 and own is not None)
    if owner in candidates:
        return SWAPPABLE, owner, candidates - set([owner])
    return NOT_SWAPPABLE, None, None

print("Calculating borders...")

recalc_borders()

print("Deforming regions...")

def assign_colors():
    used_colors = {}
    for region in regions:
        r, g, b = region.mean_with(None)
        r, g, b = int(round(r)), int(round(g)), int(round(b))
        if (r,g,b) in used_colors:
            for color in sorted([(r2, g2, b2) for r2 in range(256) for g2 in range(256) for b2 in range(256)], key=lambda color: dist_rgb(color, (r,g,b))):
                if color not in used_colors:
                    used_colors[color] = region.points
                    break
        else:
            used_colors[(r,g,b)] = region.points
    return used_colors

def make_image(colors):
    img = Image.new("RGB", image.size)
    for color in colors:
        for point in colors[color]:
            img.putpixel(point, color)
    return img

# Round status labels
FULL_ROUND = 0
NEIGHBOR_ROUND = 1
FILAMENT_ROUND = 2

max_filament = None
next_search = set()
rounds = 0
points_flipped = 0
singletons = 0
filaments = 0
flip_milestone = 0
logs = 0

while True:
    if LOGGING and (rounds % LOG_ROUND_INTERVAL == 0 or points_flipped >= flip_milestone):
        print("Round %d of deformation:\n %d edit(s) so far, of which %d singleton removal(s) and %d filament cut(s)."%(rounds, points_flipped, singletons, filaments))
        while points_flipped >= flip_milestone: flip_milestone += LOG_FLIP_INTERVAL
        if GRAPHICAL_LOGGING:
            make_image(assign_colors()).save(LOG_FILE_PREFIX + str(logs) + LOG_FILE_SUFFIX)
            logs += 1
    if max_filament is None or (round_status == NEIGHBOR_ROUND and rounds%FILAMENT_ROUND_INTERVAL != 0):
        search_space, round_status = (next_search & borders, NEIGHBOR_ROUND) if next_search else (copy(borders), FULL_ROUND)
        next_search = set()
        max_filament = None
    else:
        round_status = FILAMENT_ROUND
        search_space = set([max_filament[0]]) & borders
    search_space = list(search_space)
    shuffle(search_space)
    for point in search_space:
        status, owner, takers = analyze_nbors(point)
        if (status == FILAMENT and num_components < N) or status in (SINGLETON, SWAPPABLE):
            color = pixels[point]
            takers_list = list(takers)
            shuffle(takers_list)
            for taker in takers_list:
                dist = dist_rgb(color, owner.mean_with(None)) - dist_rgb(color, taker.mean_with(color))
                if dist > 0:
                    if status != FILAMENT or round_status == FILAMENT_ROUND:
                        found = True
                        owner.flip_point(point)
                        taker.flip_point(point)
                        owners[point] = taker
                        recalc_borders(point)
                        next_search.add(point)
                        for nbor in full_circ(point):
                            next_search.add(nbor)
                        points_flipped += 1
                    if status == FILAMENT:
                        if round_status == FILAMENT_ROUND:
                            num_components += 1
                            filaments += 1
                        elif max_filament is None or max_filament[1] < dist:
                            max_filament = (point, dist)
                    if status == SINGLETON:
                        num_components -= 1
                        singletons += 1
                    break
    rounds += 1
    if round_status == FILAMENT_ROUND:
        max_filament = None
    if round_status == FULL_ROUND and max_filament is None and not next_search:
        break

print("Deformation completed after %d rounds:\n %d edit(s), of which %d singleton removal(s) and %d filament cut(s)."%(rounds, points_flipped, singletons, filaments))

print("Assigning colors...")

used_colors = assign_colors()

print("Producing output...")

make_image(used_colors).save(OUT_FILE)

print("Done!")

```

## How it works

The program divides the canvas into `P` regions, each of which consists of some number of cells without holes. Initially, the canvas is just divided into approximate squares, which are randomly assigned to the regions. Then, these regions are "deformed" in an iterative process, where a given pixel can change its region if

1. the change would decrease the pixel's RGB distance from the average color of the region that contains it, and
2. it does not break or merge cells or introduce holes in them.

The latter condition can be enforced locally, so the process is a bit like a cellular automaton. This way, we don't have to do any pathfinding or such, which speeds the process up greatly. However, since the cells can't be broken up, some of them end up as long "filaments" that border other cells and inhibit their growth. To fix this, there is a process called "filament cut", which occasionally breaks a filament-shaped cell in two, if there are less than `N` cells at that time. Cells can also disappear if their size is 1, and this makes room for the filaments cuts.

The process ends when no pixel has the incentive to switch regions, and after that, each region is simply colored by its average color. Usually there will be some filaments remaining in the output, as can be seen in the examples below, especially in the nebula.

## P = 30, N = 500

More pictures later.

Some interesting properties of my program are that it is probabilistic, so that the results may vary between different runs, unless you use the same pseudorandom seed of course. The randomness is not essential, though, I just wanted to avoid any accidental artifacts that may result from the particular way Python traverses a set of coordinates or something similar. The program tends to use all `P` colors and almost all `N` cells, and the cells never contain holes by design. Also, the deformation process is quite slow. The colored balls took almost 15 minutes to produce on my machine. On the upside, it you turn on the `GRAPHICAL_LOGGING` option, you'll get a cool series of pictures of the deformation process. I made the Mona Lisa ones into a GIF animation (shrunk 50 % to reduce the file size). If you look closely at her face and hair, you can see the filament cutting process in action.
> 9 votes

---
Tags: popularity-contest, graphical-output, image-processing
---


title: Sum the time durations

text:
# Question
Title: Sum the time durations
## Challenge

Write the shortest code that can sum all the time durations that appear in the stdin. The program must only consider the strings that match with one of the following patterns and ignore the rest.

```
    HH:MM:SS     (it will be interpreted as HH hours, MM minutes and SS seconds)
    H:MM:SS      (it will be interpreted as H hours, MM minutes and SS seconds)
    MM:SS        (it will be interpreted as MM minutes, SS seconds)
    M:SS         (it will be interpreted as M minutes, SS seconds)

```

*examples of strings that match with the enumerated patterns:*

```
    12:00:01
    2:03:22
    00:53
    9:13

```

The output should be of the form

```
    HHh MMm SSs      (that means HH hours, MM minutes and SS seconds with non-zero-padding)

```

## Example

**STDIN**

> View the Welcome video.
> Video: 10:37 min.
> View the video introduction to the course.
> Video: 3:30 min. View the video of how to use the Lesson Overview.
> Video: 9:13 min.
> View the video overview of how to use the Epsilen system to share your work.
> Video: 03:15 min.
> View the video to learn about the State of Texas Assessment of Academic Readiness (STAAR).
> Video: 1:05:26 min.

**STDOUT**

> 1h 32m 1s

# Answer
## Pyth 105

```
K"smh"J"\D\D+|\d+:(?=\d:)|:\d\D"W:QJ1=Q:QJd;FN_msdCfn2lTm+*]0</k\:2msbck\:cQ)~k+d+hK_`%+NZ60=Z/N60=KtK;_k

```

Try it online.

This requires input from STDIN in the same way that the Javascript answer does, as quoted text with newlines as `\n`s.

Sample:

```
"View the Welcome video.\nVideo: 10:37 min.\nView the video introduction to the course.\nVideo: 3:30 min. View the video of how to use the Lesson Overview.\nVideo: 9:13 min.\nView the video overview of how to use the Epsilen system to share your work.\nVideo: 03:15 min.\nView the video to learn about the State of Texas Assessment of Academic Readiness (STAAR).\nVideo: 1:05:26 min."

```

Output

```
1h 32m 1s

```

Example working with weirder dates:

```
"10:10:5 and 5:1:10 and 27 or 16: or 1:1:1 or 11:1\n"

```

Output

```
0h 11m 20s

```

(Only the 10:10 and the 1:10 are legitimate times)

The main reason that this is so long is that Pyth won't let you extract positive matches. This instead matches everything that isn't a valid time, and replaces it with a space character. Then, splitting on whitespace leaves only times and some wayward numbers. The excess numbers are removed by checking for `:` characters, which will have been removed from non-valid times. This could almost certainly be golfed further ;)
> 3 votes

# Answer
## Javascript ES6, 138 chars

### Function, 139

Takes string as an argument and writes output to console:

```
f=s=>(r=0,s.replace(/(\d\d?):(\d\d)(:(\d\d))?/g,(m,a,b,x,c)=>r+=x?+c+b*60+a*3600:+b+a*60),console.log("%dh %dm %ds",r/3600,r%3600/60,r%60))

```

### Program, 138

```
prompt(r=0).replace(/(\d\d?):(\d\d)(:(\d\d))?/g,(m,a,b,x,c)=>r+=x?+c+b*60+a*3600:+b+a*60),console.log("%dh %dm %ds",r/3600,r%3600/60,r%60)

```

### Test for function

```
f("View the Welcome video.\n\
Video: 10:37 min.\n\
View the video introduction to the course.\n\
Video: 3:30 min. View the video of how to use the Lesson Overview.\n\
Video: 9:13 min.\n\
View the video overview of how to use the Epsilen system to share your work.\n\
Video: 03:15 min.\n\
View the video to learn about the State of Texas Assessment of Academic Readiness (STAAR).\n\
Video: 1:05:26 min.")

```

### Output

```
"1h 32m 1s"

```
> 6 votes

# Answer
# JavaScript, ES6, ~~208 200~~ 197 bytes

I know this is super long, but I wanted to explore the latest features of ES6, reverse, map-reduce, arrow functions and array comprehension (spread operator).

```
alert(prompt().match(/\d\d?:\d\d(:\d\d)?/g).map(x=>[...x.split(":").reverse(),z=0].slice(0,3)).reduce((a,b)=>b.map((y,i)=>+y+ +a[i])).map((x,i)=>(z=(t=x+z|0)/60,t%60+"smh"[i])).reverse().join(" "))
```

Just run the snippet in a latest Firefox.

**How it works (ungolfed a bit)**

```
alert(                              // Alert the final result
  prompt()                          // Take the input via prompt
  .match(/\d\d?:\d\d(:\d\d)?/g)     // Match only correct time formats
  .map(                             // Map all matches using this method
    x=>[                            // Take each element as argument x
      ...x.split(":").reverse(),    // split x on ":" and reverse the array, then spread it
      z=0                           // put 0 as last element of return array
    ].slice(0,3)                    // Take only first 3 elements of the array
  ).reduce(                         // Reduce the result using this method
    (a,b)=>                         // Pairwise elements of the array
    b.map(                          // Map array b
      (y,i)=>~~y+~~a[i]             // Convert b[i] to b[i]+a[i]
    )                               // Now we have array like [SS, MM, HH]
  ).map(                            // Map these three values for carry over calculation
    (x,i)=>(
      t=x+z,                        // z contains carryover amount, add it to this value
      z=(t/60)|0,                   // Carryover is now floor(t/60)
      t%60+"smh"[i]                 // Remove overflow from t and add "s", "m" or "h"
    )                               // Now we have array like ["SSs", "MMm", "HHh"]
  ).reverse().join(" ")             // Reverse it and join by space
)

```
> 4 votes

# Answer
### Bash (with grep, sed, awk and date): 124 bytes, 120 bytes

Just pipe the text into this:

```
grep -o '[:0-9]*'|sed 's/^[^:]*:[^:]*$/:\0/'|awk -F: '{T+=3600*$1+60*$2+$3}END{print"@"T}'|xargs date +"%Hh %Mm %Ss" -ud

```

How it works

* grep: outputs strings from the input containing only `0123456789:`
* sed: turns MM:SS and M:SS into :M:SS
* awk: calculates the seconds, empty string is 0
* xargs: passes input as argument to date
* date: converts seconds since epoch (prefixed with @) to the required format
> 4 votes

# Answer
## Perl - ~~228~~ 201

```
use integer;$h=0,$m=0,$s=0;while(<>){if(/(\d+:){1,2}\d+/){@a=reverse(split(/:/,$&));push @a,(0)x(3-@a);$s+=@a[0];$m+=@a[1];$h+=@a[2];}}$m+=$s/60;$s=$s%60;$h+=$m/60;$m=$m%60;print $h."h ".$m."m ".$s."s"

```

It happens to be the same algorithm as Optimizer's (grep,split,reverse,add).

I am no Perl expert, so maybe the byte count can be reduced.

## Ungolfed

```
use integer;                              # will do integer division
$h=0,$m=0,$s=0;
while(<>){
    if(/(\d+:){1,2}\d+/) {                # extract date formats
        @a = reverse(split(/:/,$&));      # split by ":" and reverse
        push @a,(0)x(3-@a);               # pad with zeros (minutes and hours)
        $s+=@a[0];                        # sum seconds
        $m+=@a[1];                        # sum minutes
        $h+=@a[2];                        # sum hours
    }
}

# convert seconds as minutes
$m += $s / 60;
$s = $s % 60;

# convert minutes as hours
$h += $m / 60;
$m = $m % 60;

print $h."h ".$m."m ".$s."s";

```
> 3 votes

# Answer
# Rebol - 174

```
n: charset"1234567890"a:[1 2 n]b:[":"2 n]c: 0 parse input[any[copy x[a b b](c: c + do x)| copy x[a b](c: c + do join"0:"x)| skip]]print reword"$1h $2m $3s"[1 c/1 2 c/2 3 c/3]

```

Ungolfed + annotated:

```
n: charset "1234567890"                      ; setup \d regex equiv
a: [1 2 n]                                   ; parse rule for \d{1,2}
b: [":" 2 n]                                 ; parse rule for :\d\d
c: 0                                         ; time counter

parse input [                                ; parse the input (STDIN)
                                             ; (no regex in Rebol)

  any [                                      ; match zero or more...
                                             ;
      copy x [a b b] (c: c + do x)           ;  HH:MM:SS or H:MM:SS
                                             ;    - copy match to x
                                             ;    - increment time (c) by x
                                             ; OR
    | copy x [a b] (c: c + do join "0:" x)   ;  MM:SS or M:SS
                                             ;    - copy match to x
                                             ;    - "MM:SS" into "0:MM:SS" (join)
                                             ;    - then increment time (c)
                                             ; OR
    | skip                                   ;   no match so move through input
  ]
]

print reword "$1h $2m $3s" [1 c/1 2 c/2 3 c/3]

```

Rebol comes with its own `time!` datatype. You can see how the above code makes use of this from example below (from within the Rebol console):

```
>> 0:10:37 + 0:3:30 + 0:9:13 + 0:3:15 + 1:05:26
== 1:32:01

;; Rebol would treat 10:37 as 10 hours & 37 minutes (and not MM:SS)
;; So we have to prefix the "0:"

>> join "0:" 10:37
== "0:10:37"

;; This is a string so we use Rebol DO evaluator to convert to time!

>> do join "0:" 10:37
== 0:10:37

>> type? do join "0:" 10:37
== time!

>> hms: do join "0:" 10:37
== 0:10:37

>> hms/hour
== 0

>> hms/second
== 37

>> hms/minute
== 10

```
> 3 votes

# Answer
## Groovy - 195

```
M=60
r=(System.in.text=~/((\d?\d):)?(\d\d):(\d\d)/).collect{it[2..4]*.toInteger().inject{s,i->(s?:0)*M+i}}.inject{s,i->s+=i}
f=[];(2..0).each{j=M**it;s=r%j;f<<(r-s)/j;r=s}
printf("%sh %sm %ss",f)

```

I can't figure out how to compress it more.

## Ungolfed

```
M=60
r=(System.in.text=~/((\d?\d):)?(\d\d):(\d\d)/).collect{  // extract dates
    it[2..4]*.toInteger().inject{ s,i ->                 // convert to seconds
        (s?:0)*M+i
    }
}.inject{s,i ->
    s+=i                                                 // sum seconds
}

f=[];
(2..0).each{                                             // convert to h,m,s
    j=M**it;
    s=r%j;
    f<<(r-s)/j;
    r=s
}

printf("%sh %sm %ss",f)

```
> 2 votes

# Answer
# Mathematica 300 chars

This little exercise took up a lot of code, even for Mathematica. Surely there are more efficient ways to do this.

**Golfed**

Assuming that the input is stored in `txt`,

```
n=NumberString;
t=ToExpression;
o=TimeObject;

QuotientRemainder[QuantityMagnitude[Plus@@((o[#]-o[{0,0,0}])&/@
(StringSplit[StringCases[w,{(n~~":"~~n~~":"~~n),(n~~":"~~n)}],":"]
/.{{a_,b_}:> {0,t@a,t@b},{a_,b_,c_}:> {t@a,t@b,t@c}}))],60]/.{h_,m_}:>
Row[{h,"h ",IntegerPart@m,"m ",Round[60 FractionalPart[m]],"s "}]

```

**How it works** (using unGolfed code):

1-Find the times.

```
StringCases[txt,{(NumberString~~":"~~NumberString~~":"~~NumberString),
(NumberString~~":"~~NumberString)}];

```

> {"10:37", "3:30", "9:13", "03:15", "1:05:26"}

---

2-Break into hours, minutes, seconds

```
StringSplit[%,":"]/.{{a_,b_}:> {0,ToExpression@a,ToExpression@b},{a_,b_,c_}:>
{ToExpression@a,ToExpression@b,ToExpression@c}}

```

> {{0, 10, 37}, {0, 3, 30}, {0, 9, 13}, {0, 3, 15}, {1, 5, 26}}

---

3-Sum the times. Time objects are clock times. Subtracting one time object from another returns a duration, in this case 92.0167 minutes. `QuantityMagnitude` drops the unit of measure.

```
q=QuantityMagnitude[Plus@@((TimeObject[#]-TimeObject[{0,0,0}])&/@%)]

```

> 92.0167

---

4-Convert 92.0167 minutes into hours, minutes, seconds.

```
QuotientRemainder[q,60]/.{h_,m_}:> Row[{h,"h ",IntegerPart@m,"m ",
Round[60 FractionalPart[m]],"s "}]

```

> 1h 32m 1s
> 1 votes

# Answer
# Perl, 146

My entry prints the output with a trailing space - I hope that's ok

```
while(<>){for(/(\d?\d(?::\d\d){1,2})/g){$m=1;for(reverse split/:/,$_){$t+=$m*$_;$m*=60}}}for('s','m'){$o=($t%60)."$_ $o";$t/=60}print int$t,"h $o"

```

If we can assume there will be only one time per line of input, we can chop 4 characters:

```
while(<>){if(/(\d?\d(:\d\d){1,2})/){$m=1;for(reverse split/:/,$&){$t+=$m*$_;$m*=60}}}for('s','m'){$o=($t%60)."$_ $o";$t/=60}print int$t,"h $o"

```

These work by accumulating the total seconds elapsed and formatting that value afterwards.
> 1 votes

---
Tags: code-golf, arithmetic, date, regular-expression
---


title: Computer: you do the math

text:
# Question
Title: Computer: you do the math
This challenge is partly an algorithms challenge, involves some math and is partly simply a fastest code challenge.

For some positive integer `n`, consider a uniformly random string of `1`s and `0`s of length `n` and call it `A`. Now also consider a second uniformly chosen random string of length `n` whose values are `-1`, `0,` or `1` and call it `B_pre`. Now let `B` be `B_pre`+`B_pre`. That is `B_pre` concatenated to itself.

Now consider the inner product of `A` and `B[j,...,j+n-1]` and call it `Z_j` and index from `1`.

**Task**

The output should be a list of `n+1` fractions. The `i`th term in the output should be the *exact* probability that *all* of the first `i` terms `Z_j` with `j <= i` equal `0`.

**Score**

The largest `n` for which your code gives the correct output in under 10 minutes on my machine.

**Tie Breaker**

If two answers have the same score, the one submitted first wins.

In the (very very) unlikely event that someone finds a method to get unlimited scores, the first valid proof of such a solution will be accepted.

**Hint**

Do not try to solve this problem mathematically, it's too hard. The best way in my view is to go back to basic definitions of probability from high school and find smart ways to get the code to do an exhaustive enumeration of the possibilities.

**Languages and libraries**

You can use any language which has a freely available compiler/interpreter/etc. for Linux and any libraries which are also freely available for Linux.

**My machine** The timings will be run on my machine. This is a standard ubuntu install on an AMD FX-8350 Eight-Core Processor. This also means I need to be able to run your code. As a consequence, only use easily available free software and please include full instructions how to compile and run your code.

---

Some test outputs. Consider just the first output for each `n`. That is when `i=1`. For `n` from 1 to 13 they should be.

```
 1: 4/6
 2: 18/36
 3: 88/216
 4: 454/1296
 5: 2424/7776
 6: 13236/46656
 7: 73392/279936
 8: 411462/1679616
 9: 2325976/10077696
10: 13233628/60466176
11: 75682512/362797056
12: 434662684/2176782336
13: 2505229744/13060694016

```

You can also find a general formula for `i=1` at http://oeis.org/A081671 .

**Leaderboard (split by language)**

* n = 15. **Python + parallel python + pypy** in 1min49s by Jakube
* n = 17. **C++** in 3min37s by Keith Randall
* n = 16. **C++** in 2min38s by kuroi neko

# Answer
## C++, n=18 in 9 min on 8 threads

(Let me know if it makes it under 10 minutes on your machine.)

I take advantage of several forms of symmetry in the B array. Those are cyclic (shift by one position), reversal (reverse the order of the elements), and sign (take the negative of each element). First I compute the list of Bs we need to try and their weight. Then each B is run through a fast routine (using bitcount instructions) for all 2^n values of A.

Here's the result for n==18:

```
> time ./a.out 18
 1: 16547996212044 / 101559956668416
 2:  3120508430672 / 101559956668416
 3:   620923097438 / 101559956668416
 4:   129930911672 / 101559956668416
 5:    28197139994 / 101559956668416
 6:     6609438092 / 101559956668416
 7:     1873841888 / 101559956668416
 8:      813806426 / 101559956668416
 9:      569051084 / 101559956668416
10:      510821156 / 101559956668416
11:      496652384 / 101559956668416
12:      493092812 / 101559956668416
13:      492186008 / 101559956668416
14:      491947940 / 101559956668416
15:      491889008 / 101559956668416
16:      449710584 / 101559956668416
17:      418254922 / 101559956668416
18:      409373626 / 101559956668416

real    8m55.854s
user    67m58.336s
sys 0m5.607s

```

Compile the program below with `g++ --std=c++11 -O3 -mpopcnt dot.cc`

```
#include <stdio.h>
#include <stdlib.h>
#include <vector>
#include <thread>
#include <mutex>
#include <chrono>

using namespace std;

typedef long long word;

word n;

void inner(word bpos, word bneg, word w, word *cnt) {
    word maxi = n-1;
    for(word a = (1<<n)-1; a >= 0; a--) {
        word m = a;
        for(word i = maxi; i >= 0; i--, m <<= 1) {
            if(__builtin_popcount(m&bpos) != __builtin_popcount(m&bneg))
                break;
            cnt[i]+=w;
        }
    }
}

word pow(word n, word e) {
    word r = 1;
    for(word i = 0; i < e; i++) r *= n;
    return r;
}

typedef struct {
    word b;
    word weight;
} Bentry;

mutex block;
Bentry *bqueue;
word bhead;
word btail;
word done = -1;

word maxb;

// compute -1*b
word bneg(word b) {
    word w = 1;
    for(word i = 0; i < n; i++, w *= 3) {
        word d = b / w % 3;
        if(d == 1)
            b += w;
        if(d == 2)
            b -= w;
    }
    return b;
}

// rotate b one position
word brot(word b) {
    b *= 3;
    b += b / maxb;
    b %= maxb;
    return b;
}

// reverse b
word brev(word b) {
    word r = 0;
    for(word i = 0; i < n; i++) {
        r *= 3;
        r += b % 3;
        b /= 3;
    }
    return r;
}

// individual thread's work routine
void work(word *cnt) {
    while(true) {
        // get a queue entry to work on
        block.lock();
        if(btail == done) {
            block.unlock();
            return;
        }
        if(bhead == btail) {
            block.unlock();
            this_thread::sleep_for(chrono::microseconds(10));
            continue;
        }
        word i = btail++;
        block.unlock();

        // thread now owns bqueue[i], work on it
        word b = bqueue[i].b;
        word w = 1;
        word bpos = 0;
        word bneg = 0;
        for(word j = 0; j < n; j++, b /= 3) {
            word d = b % 3;
            if(d == 1)
                bpos |= 1 << j;
            if(d == 2)
                bneg |= 1 << j;
        }
        bpos |= bpos << n;
        bneg |= bneg << n;
        inner(bpos, bneg, bqueue[i].weight, cnt);
    }
}

int main(int argc, char *argv[]) {
    n = atoi(argv[1]);

    // allocate work queue
    maxb = pow(3, n);
    bqueue = (Bentry*)(malloc(maxb*sizeof(Bentry)));

    // start worker threads
    word procs = thread::hardware_concurrency();
    vector<thread> threads;
    vector<word*> counts;
    for(word p = 0; p < procs; p++) {
        word *cnt = (word*)calloc(64+n*sizeof(word), 1);
        threads.push_back(thread(work, cnt));
        counts.push_back(cnt);
    }

    // figure out which Bs we actually want to test, and with which weights
    bool *bmark = (bool*)calloc(maxb, 1);
    for(word i = 0; i < maxb; i++) {
        if(bmark[i]) continue;
        word b = i;
        word w = 0;
        for(word j = 0; j < 2; j++) {
            for(word k = 0; k < 2; k++) {
                for(word l = 0; l < n; l++) {
                    if(!bmark[b]) {
                        bmark[b] = true;
                        w++;
                    }
                    b = brot(b);
                }
                b = bneg(b);
            }
            b = brev(b);
        }
        bqueue[bhead].b = i;
        bqueue[bhead].weight = w;
        block.lock();
        bhead++;
        block.unlock();
    }
    block.lock();
    done = bhead;
    block.unlock();

    // add up results from threads
    word *cnt = (word*)calloc(n,sizeof(word));
    for(word p = 0; p < procs; p++) {
        threads[p].join();
        for(int i = 0; i < n; i++) cnt[i] += counts[p][i];
    }
    for(word i = 0; i < n; i++)
        printf("%2lld: %14lld / %14lld\n", i+1, cnt[n-1-i], maxb<<n);
    return 0;
}

```
> 5 votes

# Answer
# Python 2 using pypy and pp: n = 15 in 3 minutes

Also just a simple brute force. Interesting to see, that I nearly get the same same speed as kuroi neko with C++. My code can reach `n = 12` in about 5 minutes. And I only run it on one virtual core.

## edit: Reduce search space by a factor of `n`

I noticed, that a cycled vector `A*` of `A` produces the same numbers as probabilities (same numbers) as the original vector `A` when I iterate over `B`. E.g. The vector `(1, 1, 0, 1, 0, 0)` has the same probabilities as each of the vectors `(1, 0, 1, 0, 0, 1)`, `(0, 1, 0, 0, 1, 1)`, `(1, 0, 0, 1, 1, 0)`, `(0, 0, 1, 1, 0, 1)` and `(0, 1, 1, 0, 1, 0)` when choosing a random `B`. Therefore I don't have to iterate over each of these 6 vectors, but only about 1 and replace `count[i] += 1` with `count[i] += cycle_number`.

This reduces the complexity from `Theta(n) = 6^n` to `Theta(n) = 6^n / n`. Therefore for `n = 13` it's about 13 times as fast as my previous version. It calculates `n = 13` in about 2 minutes 20 seconds. For `n = 14` it is still a little bit too slow. It takes about 13 minutes.

## edit 2: Multi-core-programming

Not really happy with the next improvement. I decided to also try to execute my program on multiple cores. On my 2 + 2 cores I now can calculate `n = 14` in about 7 minutes. Only a factor of 2 improvement.

The code is available in this github repo: Link. The multi-core programming makes is a little bit ugly.

## edit 3: Reducing search space for `A` vectors and `B` vectors

I noticed the same mirror-symmetry for the vectors `A` as kuroi neko did. Still not sure, why this works (and if it works for each `n`).

The reducing of the search space for `B` vectors is a bit cleverer. I replaced the generation of the vectors (`itertools.product`), with an own function. Basically, I start with an empty list and put it on a stack. Until the stack is empty, I remove a list, if it has not the same lenght as `n`, I generate 3 other lists (by appending -1, 0, 1) and pushing them onto the stack. I a list has the same length as `n`, I can evaluate the sums.

Now that I generate the vectors myself, I can filter them depending on if I can reach the sum = 0 or not. E.g. if my vector `A` is `(1, 1, 1, 0, 0)`, and my vector `B` looks `(1, 1, ?, ?, ?)`, I know, that I cannot fill the `?` with values, so that `A*B = 0`. So I don't have to iterate over all those 6 vectors `B` of the form `(1, 1, ?, ?, ?)`.

We can improve on this, if we ignore the values for 1. As noted in the question, for the values for `i = 1` are the sequence A081671. There are many ways to calculate those. I choose the simple recurrence: `a(n) = (4*(2*n-1)*a(n-1) - 12*(n-1)*a(n-2)) / n`. Since we can calculate `i = 1` in basically no time, we can filter more vectors for `B`. E.g. `A = (0, 1, 0, 1, 1)` and `B = (1, -1, ?, ?, ?)`. We can ignore vectors, where the first `? = 1`, because the `A * cycled(B) > 0`, for all these vectors. I hope you can follow. It's probably not the best example.

With this I can calculate `n = 15` in 6 minutes.

## edit 4:

Quickly implemented kuroi neko's great idea, which says, that `B` and `-B` produces the same results. Speedup x2. Implementation is only a quick hack, though. `n = 15` in 3 minutes.

## Code:

For the complete code visit Github. The following code is only a representation of the main features. I left out imports, multicore programming, printing the results, ...

```
count = [0] * n
count[0] = oeis_A081671(n)

#generating all important vector A
visited = set(); todo = dict()
for A in product((0, 1), repeat=n):
    if A not in visited:
        # generate all vectors, which have the same probability
        # mirrored and cycled vectors
        same_probability_set = set()
        for i in range(n):
            tmp = [A[(i+j) % n] for j in range(n)]
            same_probability_set.add(tuple(tmp))
            same_probability_set.add(tuple(tmp[::-1]))
        visited.update(same_probability_set)
        todo[A] = len(same_probability_set)

# for each vector A, create all possible vectors B
stack = []
for A, cycled_count in dict_A.iteritems():
    ones = [sum(A[i:]) for i in range(n)] + [0]
    # + [0], so that later ones[n] doesn't throw a exception
    stack.append(([0] * n, 0, 0, 0, False))

    while stack:
        B, index, sum1, sum2, used_negative = stack.pop()
        if index < n:
            # fill vector B[index] in all possible ways,
            # so that it's still possible to reach 0.
            if used_negative:
                for v in (-1, 0, 1):
                    sum1_new = sum1 + v * A[index]
                    sum2_new = sum2 + v * A[index - 1 if index else n - 1]
                    if abs(sum1_new) <= ones[index+1]:
                        if abs(sum2_new) <= ones[index] - A[n-1]:
                            C = B[:]
                            C[index] = v
                            stack.append((C, index + 1, sum1_new, sum2_new, True))
            else:
                for v in (0, 1):
                    sum1_new = sum1 + v * A[index]
                    sum2_new = sum2 + v * A[index - 1 if index else n - 1]
                    if abs(sum1_new) <= ones[index+1]:
                        if abs(sum2_new) <= ones[index] - A[n-1]:
                            C = B[:]
                            C[index] = v
                            stack.append((C, index + 1, sum1_new, sum2_new, v == 1))
        else:
            # B is complete, calculate the sums
            count[1] += cycled_count  # we know that the sum = 0 for i = 1
            for i in range(2, n):
                sum_prod = 0
                for j in range(n-i):
                    sum_prod += A[j] * B[i+j]
                for j in range(i):
                    sum_prod += A[n-i+j] * B[j]
                if sum_prod:
                    break
                else:
                    if used_negative:
                        count[i] += 2*cycled_count
                    else:
                        count[i] += cycled_count

```

## Usage:

You have to install pypy (for Python 2!!!). The parallel python module isn't ported for Python 3. Then you have to install the parallel python module pp-1.6.4.zip. Extract it, `cd` into the folder and call `pypy setup.py install`.

Then you can call my program with

`pypy you-do-the-math.py 15`

It will automatically determine the number of cpu's. There may be some error messages after finishing the program, just ignore them. `n = 16` should be possible on your machine.

## Output:

```
Calculation for n = 15 took 2:50 minutes

 1  83940771168 / 470184984576  17.85%
 2  17379109692 / 470184984576   3.70%
 3   3805906050 / 470184984576   0.81%
 4    887959110 / 470184984576   0.19%
 5    223260870 / 470184984576   0.05%
 6     67664580 / 470184984576   0.01%
 7     30019950 / 470184984576   0.01%
 8     20720730 / 470184984576   0.00%
 9     18352740 / 470184984576   0.00%
10     17730480 / 470184984576   0.00%
11     17566920 / 470184984576   0.00%
12     17521470 / 470184984576   0.00%
13     17510280 / 470184984576   0.00%
14     17507100 / 470184984576   0.00%
15     17506680 / 470184984576   0.00%

```

## Notes and ideas:

* I have a i7-4600m processor with 2 cores and 4 threads. It doesn't matter If I use 2 or 4 threads. The cpu-usage is 50% with 2 threads and 100% with 4 threads, but it still takes the same amount of time. I don't know why. I checked, that each thread only has to the the half amout of data, when there are 4 threads, checked the results, ...
* I use a lot of lists. Python isn't quite efficient in storing, I have to copy lots of lists, ... So I thought of using an integer instead. I could use the bits 00 (for 0) and 11 (for 1) in the vector A, and the bits 10 (for -1), 00 (for 0) and 01 (for 1) in the vector B. For the product of A and B, I would only have to calculate `A & B` and count the 01 and 10 blocks. Cycling can be done with shifting the vector and using masks, ... I actually implemented all this, you can find it in some of my older commits on Github. But it turned out, to be slower than with lists. I guess, pypy really optimizes list operations.
> 6 votes

# Answer
# woolly bully - C++ - way too slow

Well since a better programmer took on the C++ implementation, I'm calling quits for this one.

```
#include <cstdlib>
#include <cmath>
#include <vector>
#include <bitset>
#include <future>
#include <iostream>
#include <iomanip>

using namespace std;

/*
6^^n events will be generated, so the absolute max
that can be counted by a b bits integer is
E(b*log(2)/log(6)), i.e. n=24 for a 64 bits counter

To enumerate 3 possible values of a size n vector we need
E(n*log(3)/log(2))+1 bits, i.e. 39 bits
*/
typedef unsigned long long Counter; // counts up to 6^^24

typedef unsigned long long Benumerator; // 39 bits
typedef unsigned long      Aenumerator; // 24 bits

#define log2_over_log6 0.3869

#define A_LENGTH ((size_t)(8*sizeof(Counter)*log2_over_log6))
#define B_LENGTH (2*A_LENGTH)

typedef bitset<B_LENGTH> vectorB;

typedef vector<Counter> OccurenceCounters;

// -----------------------------------------------------------------
// multithreading junk for CPUs detection and allocation
// -----------------------------------------------------------------
int number_of_CPUs(void)
{
    int res = thread::hardware_concurrency();
    return res == 0 ? 8 : res;
}

#ifdef __linux__
#include <sched.h>
void lock_on_CPU(int cpu)
{
    cpu_set_t mask;
    CPU_ZERO(&mask);
    CPU_SET(cpu, &mask);
    sched_setaffinity(0, sizeof(mask), &mask);
}
#elif defined (_WIN32)
#include <Windows.h>
#define lock_on_CPU(cpu) SetThreadAffinityMask(GetCurrentThread(), 1 << cpu)
#else
// #warning is not really standard, so this might still cause compiler errors on some platforms. Sorry about that.
#warning "Thread processor affinity settings not supported. Performances might be improved by providing a suitable alternative for your platform"
#define lock_on_CPU(cpu)
#endif

// -----------------------------------------------------------------
// B values generator
// -----------------------------------------------------------------
struct Bvalue {
    vectorB p1;
    vectorB m1;
};

struct Bgenerator {
    int n;                 // A length
    Aenumerator stop;      // computation limit
    Aenumerator zeroes;    // current zeroes pattern
    Aenumerator plusminus; // current +1/-1 pattern
    Aenumerator pm_limit;  // upper bound of +1/-1 pattern

    Bgenerator(int n, Aenumerator start=0, Aenumerator stop=0) : n(n), stop(stop)
    {
        // initialize generator so that first call to next() will generate first value
        zeroes    = start - 1;
        plusminus = -1;
        pm_limit  = 0;
    }

    // compute current B value
    Bvalue value(void)
    {
        Bvalue res;
        Aenumerator pm = plusminus;
        Aenumerator position = 1;
        int i_pm = 0;
        for (int i = 0; i != n; i++)
        {
            if (zeroes & position)
            {
                if (i_pm == 0)  res.p1 |= position; // first non-zero value fixed to +1
                else
                {
                    if (pm & 1) res.m1 |= position; // next non-zero values
                    else        res.p1 |= position;
                    pm >>= 1;
                }
                i_pm++;
            }
            position <<= 1;
        }
        res.p1 |= (res.p1 << n); // concatenate 2 Bpre instances
        res.m1 |= (res.m1 << n);
        return res;
    }

    // next value
    bool next(void)
    {
        if (++plusminus == pm_limit)
        {
            if (++zeroes == stop) return false;
            plusminus = 0;
            pm_limit = (1 << vectorB(zeroes).count()) >> 1;
        }
        return true;
    }

    // calibration: produces ranges that will yield the approximate same number of B values
    vector<Aenumerator> calibrate(int segments)
    {
        // setup generator for the whole B range
        zeroes = 0;
        stop = 1 << n;
        plusminus = -1;
        pm_limit = 0;

        // divide range into (nearly) equal chunks
        Aenumerator chunk_size = ((Aenumerator)pow (3,n)-1) / 2 / segments;

        // generate bounds for zeroes values
        vector<Aenumerator> res(segments + 1);
        int bound = 0;
        res[bound] = 1;
        Aenumerator count = 0;
        while (next()) if (++count % chunk_size == 0) res[++bound] = zeroes;
        res[bound] = stop;
        return res;
    }
};

// -----------------------------------------------------------------
// equiprobable A values merging
// -----------------------------------------------------------------
static char A_weight[1 << A_LENGTH];
struct Agroup {
    vectorB value;
    int     count;
    Agroup(Aenumerator a = 0, int length = 0) : value(a), count(length) {}
};
static vector<Agroup> A_groups;

Aenumerator reverse(Aenumerator n) // this works on N-1 bits for a N bits word
{
    Aenumerator res = 0;
    if (n != 0) // must have at least one bit set for the rest to work
    {
        // construct left-padded reverse value
        for (int i = 0; i != 8 * sizeof(n)-1; i++)
        {
            res |= (n & 1);
            res <<= 1;
            n >>= 1;
        }

        // shift right to elimitate trailing zeroes
        while (!(res & 1)) res >>= 1;
    }
    return res;
}

void generate_A_groups(int n)
{
    static bitset<1 << A_LENGTH> lookup(0);
    Aenumerator limit_A = (Aenumerator)pow(2, n);
    Aenumerator overflow = 1 << n;
    for (char & w : A_weight) w = 0;

    // gather rotation cycles
    for (Aenumerator a = 0; a != limit_A; a++)
    {
        Aenumerator rotated = a;
        int cycle_length = 0;
        for (int i = 0; i != n; i++)
        {
            // check for new cycles
            if (!lookup[rotated])
            {
                cycle_length++;
                lookup[rotated] = 1;
            }

            // rotate current value
            rotated <<= 1;
            if (rotated & overflow) rotated |= 1;
            rotated &= (overflow - 1);
        }

        // store new cycle
        if (cycle_length > 0) A_weight[a] = cycle_length;
    }

    // merge symetric groups
    for (Aenumerator a = 0; a != limit_A; a++)
    {
        // skip already grouped values
        if (A_weight[a] == 0) continue;

        // regroup a symetric pair
        Aenumerator r = reverse(a);
        if (r != a)
        {
            A_weight[a] += A_weight[r];
            A_weight[r] = 0;
        }
    }

    // generate groups
    for (Aenumerator a = 0; a != limit_A; a++)
    {
        if (A_weight[a] != 0) A_groups.push_back(Agroup(a, A_weight[a]));
    }
}

// -----------------------------------------------------------------
// worker thread
// -----------------------------------------------------------------
OccurenceCounters solve(int n, int index, Aenumerator Bstart, Aenumerator Bstop)
{
    OccurenceCounters consecutive_zero_Z(n, 0);  // counts occurences of the first i terms of Z being 0

    // lock on assigned CPU
    lock_on_CPU(index);

    // enumerate B vectors
    Bgenerator Bgen(n, Bstart, Bstop);
    while (Bgen.next())
    {
        // get next B value
        Bvalue B = Bgen.value();

        // enumerate A vector groups
        for (const auto & group : A_groups)
        {
            // count consecutive occurences of inner product equal to zero
            vectorB sliding_A(group.value);
            for (int i = 0; i != n; i++)
            {
                if ((sliding_A & B.p1).count() != (sliding_A & B.m1).count()) break;
                consecutive_zero_Z[i] += group.count;
                sliding_A <<= 1;
            }
        }
    }
    return consecutive_zero_Z;
}

// -----------------------------------------------------------------
// main
// -----------------------------------------------------------------
#define die(msg) { cout << msg << endl; exit (-1); }

int main(int argc, char * argv[])
{
    int n = argc == 2 ? atoi(argv[1]) : 16; // arbitray value for debugging
    if (n < 1 || n > 24) die("vectors of lenght between 1 and 24 is all I can (try to) compute, guv");

    auto begin = time(NULL);

    // one worker thread per CPU
    int num_workers = number_of_CPUs();

    // regroup equiprobable A values
    generate_A_groups(n);

    // compute B generation ranges for proper load balancing
    vector<Aenumerator> ranges = Bgenerator(n).calibrate(num_workers);

    // set workers to work
    vector<future<OccurenceCounters>> workers(num_workers);
    for (int i = 0; i != num_workers; i++)
    {
        workers[i] = async(
            launch::async, // without this parameter, C++ will decide whether execution shall be sequential or asynchronous (isn't C++ fun?).
            solve, n, i, ranges[i], ranges[i+1]);
    }

    // collect results
    OccurenceCounters result(n + 1, 0);
    for (auto& worker : workers)
    {
        OccurenceCounters partial = worker.get();
        for (size_t i = 0; i != partial.size(); i++) result[i] += partial[i]*2; // each result counts for a symetric B pair
    }
    for (Counter & res : result) res += (Counter)1 << n; // add null B vector contribution
    result[n] = result[n - 1];                           // the last two probabilities are equal by construction

    auto duration = time(NULL) - begin;

    // output
    cout << "done in " << duration / 60 << ":" << setw(2) << setfill('0') << duration % 60 << setfill(' ')
        << " by " << num_workers << " worker thread" << ((num_workers > 1) ? "s" : "") << endl;
    Counter events = (Counter)pow(6, n);
    int width = (int)log10(events) + 2;
    cout.precision(5);
    for (int i = 0; i <= n; i++) cout << setw(2) << i << setw(width) << result[i] << " / " << events << " " << fixed << (float)result[i] / events << endl;

    return 0;
}

```

## Building the executable

It's a standalone C++11 source that compiles without warnings and runs smoothly on:

* Win7 & MSVC2013
* Win7 & MinGW - g++4.7
* Ubuntu & g++4.8 (in a VirtualBox VM with 2 CPUs allocated)

If you compile with g++, use : **g++ -O3 -pthread -std=c++11**
forgetting about `-pthread` will produce a nice and friendly core dump.

## Optimizations

1. The last Z term is equal to the first (Bpre x A in both cases), so the last two results are always equal, which dispenses of computing the last Z value.
The gain is neglectible, but coding it costs nothing so you might as well use it.
2. As Jakube discovered, all cyclic values of a given A vector produce the same probabilities.
You can compute these with a single instance of A and multiply the result by the number of its possible rotations. Rotation groups can easily be pre-computed in a neglectible amount of time, so this is a huge net speed gain.
Since the number of permutations of an n length vector is n-1, the complexity drops from o(6<sup>n</sup>) to o(6<sup>n</sup>/(n-1)), basically going a step further for the same computation time.
3. It appears pairs of symetric patterns also generate the same probabilities. For instance, 100101 and 101001.
I have no mathematical proof of that, but intuitively when presented with all possible B patterns, each symetric A value will be convoluted with the corresponding symetric B value for the same global outcome.
This allows to regroup some more A vectors, for an approximative 30% reduction of A groups number.
4. **WRONG** For some semi-mysterious reason, all patterns with only one or two bits set produce the same result. This does not represent that many distinct groups, but still they can be merged at virtually no cost.
5. The vectors B and -B (B with all components multiplied by -1) produce the same probabilities.
(for instance \[ 1, 0,-1, 1\] and \[-1, 0, 1,-1\]).
Except for the null vector (all components equal to 0), B and -B form a pair of *distinct* vectors.
This allows to cut the number of B values in half by considering only one of each pair and multiplying its contribution by 2, adding the known global contribution of null B to each probability only once.

## How it works

The number of B values is huge (3<sup>n</sup>), so pre-computing them would require indecent amounts of memory, which would slow down computation and eventually exhaust available RAM.
Unfortunately, I could not find a simple way of enumerating the half set of optimized B values, so I resorted to coding a dedicated generator.

The mighty B generator was a lot of fun to code, though languages that support yield mechanisms would have allowed to program it in a much more elegant way.
What it does in a nutshell is consider the "skeleton" of a Bpre vector as a binary vector where 1s represent actual -1 or +1 values.
Among all these +1/-1 potential values, the first one is fixed to +1 (thus selecting one of the possible B/-B vector), and all remaining possible +1/-1 combinations are enumerated.
Lastly, a simple calibration system ensures each worker thread will process a range of values of approximately the same size.

A values are heavily filtered for regrouping in equiprobable chunks.
This is done in a pre-computing phase that brute-force examines all possible values.
This part has a neglectible O(2<sup>n</sup>) execution time and does not need to be optimized (the code being already unreadable enough as it is!).

To evaluate the inner product (which needs only be tested against zero), the -1 and 1 components of B are regrouped into binary vectors.
The inner product is null if (and only if) there is an equal number of +1s and -1s among the B values corresponding to non-zero A values.
This can be computed with simple masking and bit count operations, helped by `std::bitset` that will produce reasonably efficient bit count code without having to resort to ugly intrinsic instructions.

The work is equally divided among cores, with forced CPU affinity (every little bit helps, or so they say).

## Example result

```
C:\Dev\PHP\_StackOverflow\C++\VectorCrunch>release\VectorCrunch.exe 16
done in 8:19 by 4 worker threads
 0  487610895942 / 2821109907456 0.17284
 1   97652126058 / 2821109907456 0.03461
 2   20659337010 / 2821109907456 0.00732
 3    4631534490 / 2821109907456 0.00164
 4    1099762394 / 2821109907456 0.00039
 5     302001914 / 2821109907456 0.00011
 6     115084858 / 2821109907456 0.00004
 7      70235786 / 2821109907456 0.00002
 8      59121706 / 2821109907456 0.00002
 9      56384426 / 2821109907456 0.00002
10      55686922 / 2821109907456 0.00002
11      55508202 / 2821109907456 0.00002
12      55461994 / 2821109907456 0.00002
13      55451146 / 2821109907456 0.00002
14      55449098 / 2821109907456 0.00002
15      55449002 / 2821109907456 0.00002
16      55449002 / 2821109907456 0.00002

```

## Performances

Multithreading should work perfectly on this, though only "true" cores will fully contribute to computation speed. My CPU has only 2 cores for 4 CPUs, and the gain over a single-threaded version is "only" about 3.5.

## Compilers

An initial problem with multithreading led me to believe the GNU compilers were performing worse than Microsoft.

After more thorough testing, it appears g++ wins the day yet again, producing approximatively 30% faster code (the same ratio I noticed on two other computation-heavy projects).

Notably, the `std::bitset` library is implemented with dedicated bit count instructions by g++ 4.8, while MSVC 2013 uses only loops of conventional bit shifts.

As one could have expected, compiling in 32 or 64 bits makes no difference.

## Further refinements

I noticed a few A groups producing the same probabilities after all the reduction operations, but I could not identify a pattern that would allow to regroup them.

Here are the pairs I noticed for n=11:

```
  10001011 and 10001101
 100101011 and 100110101
 100101111 and 100111101
 100110111 and 100111011
 101001011 and 101001101
 101011011 and 101101011
 101100111 and 110100111
1010110111 and 1010111011
1011011111 and 1011111011
1011101111 and 1011110111

```
> 5 votes

---
Tags: code-challenge, math, counting, rational-numbers, probability-theory
---


# Question
Title: Coding a recursive function for highest possible input
***Challenge***

You are given the following function:- which is the same as:-

with the base cases q(r, b, L) = 1 whenever r ≤ L, q(r, 0, L) = 0, if r \> L and q(r, 0, L) = 1, if r ≤ L.

Your task is to code a program that takes r as input, and outputs the value of q(r, r - L, L) for all L taking integer values from 1 to (r-1), where r is any nonnegative integer.

***Example 1***

***Input***

> Enter the value of r: 2

***Output***

> q(2,1,1) = 0.3333333333

***Example 2***

***Input***

> Enter the value of r: 3

***Output***

> q(3,2,1) = 0.1
>
> q(3,1,2) = 0.5

***Winning criterion***

The code that can correctly output q(r, r-L, L) for all L taking integer values from 1 to (r-1), for **the highest value of r**, in less than 300 seconds. In case of a tie, the code with lesser runtime will be considered. As this is a runtime-based challenge, I shall test all submissions on my machine.

# Answer
## Java

I've applied some algebraic transformation and dynamic programming.

```
public class CodeGolf42234 {
    public static void main(String[] args) {
        int r = Integer.parseInt(args[0]);
        for (int L = 1; L < r; L++) System.out.println(qDouble5(r, r-L, L));
    }

    private static double qDouble5(int _r, int _b, int L) {
        double[][] q = new double[_r+1][_b+1];
        for (int r = 0; r <= L; r++) {
            for (int b = 0; b < q[r].length; b++) q[r][b] = 1;
        }
        for (int r = L + 1; r < q.length; r++) {
            for (int b = 1 + (r - L - 1)/L; b < q[r].length; b++) {
                double sum = 0, m = 1;
                for (int k = 0; k <= L; k++) {
                    sum += m * q[r-k][b-1];
                    m = m * (r - k) / (r + b - 1 - k);
                }
                q[r][b] = sum * b / (r + b);
            }
        }
        return q[_r][_b];
    }
}

```
> 7 votes

# Answer
## Java

Here's a recursive solution. Just like you wanted. You can always uncomment the output and display it if you want.

```
import java.util.ArrayList;
import java.util.Scanner;

public class HeaderFinder
{
static ArrayList<Double> answers = new ArrayList<Double>();
public static void main(String[]args)
{
    Scanner scanner = new Scanner(System.in);
    int r = scanner.nextInt();
    for(int b=r-1; b!=0; b--)
        for(int L=r-1; L!=0; L--)
            answers.add(formulate(r,b,L));

//  int i=-1;
//      for(int b=r-1; b!=0; b--)
//          for(int L=r-1; L!=0; L--)
//              System.out.println("q("+r+", "+b+", "+L+") = "+answers.get(++i));

}

public static double formulate(double r, double b, double L)
{
    if(b==0)
        if(r>L)
            return 0.0;
        else
            return 1;
    if(r<=L)
        return 1;
    double partone = (b/(r+b));
    double result = 0.0;
    for(int k=1; k<=L; k++)
    {
        result+=summation(r, b, k) * (b/(r+b-k)) * formulate(r-k, b-1, L);
    }
    return (partone * formulate(r, b-1, L)) + result;
}

public static double summation(double r, double b, int k)
{
    double rb = r+b;
    double mul=r/rb;
    for(int i=1; i<k; i++)
    {
        mul*=(r-i)/(rb-i);
    }
    return mul;
}
}

```
> 1 votes

---
Tags: math, arithmetic, fastest-code, recursion
---