/r/dailyprogrammer #206I on reals

contour_plot.py

import numpy as np
import matplotlib.cm as cm
import matplotlib.pyplot as plt

header, *lines = open("map2.txt").read().splitlines()
R = float(header.split()[2])
H, W = len(lines), len(lines[0])
crops = []
for y, line in enumerate(lines):
    for x, c in enumerate(line):
        if c == "x":
            crops.append((x, y))

x = np.arange(0, W, 0.2)
y = np.arange(0, H, 0.2)
X, Y = np.meshgrid(x, y)

Z = [[sum(1 for X, Y in crops if (_x-X)*(_x-X) + (_y-Y)*(_y-Y) <= R*R) for _x in x] for _y in y]

print("plotting...")

plt.figure()
plt.imshow(Z, origin='lower', cmap=cm.gray, extent=(0, W, 0, H))
CS = plt.contour(X, Y, Z, 5, extend='both')
plt.clabel(CS, inline=1, fmt='%1.1f', fontsize=8)
plt.scatter([x for x, y in crops], [y for x, y in crops])
plt.axis([0, W, 0, H])
plt.grid(True)
plt.gca().invert_yaxis()
plt.show()

monte_carlo.py

import random
import math
import numpy as np
import matplotlib.pyplot as plt

header, *lines = open("map2.txt").read().splitlines()
R = float(header.split()[2])
H, W = len(lines), len(lines[0])
crops = []
for y, line in enumerate(lines):
    for x, c in enumerate(line):
        if c == "x":
            crops.append((x, y))

points = list(crops)
maxR = R

iterations = 6
n_of_points = 100
at_least = 37

totalbest = 0
for i in range(iterations):
    print("iteration:", i, "points to check:", len(points))
    best = totalbest
    newpoints = []
    for st_x, st_y in points:
        for j in range(n_of_points//(i+1)): # first iteration is the most important one, the rest will mostly fill the ares
            if i == 0:
                r = maxR
            else:
                r = random.random() * maxR
            a = math.radians(random.random() * 360)
            x = st_x + r * math.cos(a)
            y = st_y + r * math.sin(a)

            n = 0
            for X, Y in crops:
                if (x-X)*(x-X) + (y-Y)*(y-Y) <= R*R:
                    n += 1
            if n == best or n >= at_least:
                newpoints.append((x, y))
            if n > best:
                best = n
                if n <= at_least:
                    newpoints = [(x, y)]
    if best > totalbest:
        totalbest = best
        if best <= at_least:
            points = []
    points += newpoints
    maxR /= 2

print(best)

plt.scatter([x for x, y in points], [y for x, y in points], color='r')
plt.scatter([x for x, y in crops], [y for x, y in crops])
plt.axis([0, W, 0, H])
plt.gca().invert_yaxis()
plt.grid(True)
plt.show()

proper_solution.py

import math
import itertools

header, *lines = open("map2.txt").read().splitlines()
R = float(header.split()[2])
H, W = len(lines), len(lines[0])
crops = []
for y, line in enumerate(lines):
    for x, c in enumerate(line):
        if c == "x":
            crops.append((x, y))

best = 0
best_points = {tuple(coord) for coord in crops}

def circle_intersection(p0, p1, r):
    x0, y0 = p0
    x1, y1 = p1
    dx = x1 - x0
    dy = y1 - y0
    d = math.sqrt((dy*dy) + (dx*dx));
    if d > 2*r:
        return None

    a = (d*d) / (2.0 * d)
    x2 = x0 + (dx * a/d)
    y2 = y0 + (dy * a/d)

    h = math.sqrt((r*r) - (a*a))
    rx = -dy * (h/d)
    ry = dx * (h/d)

    sx0, sx1 = x2 + rx, x2-rx
    sy0, sy1= y2 + ry, y2-ry

    return ((sx0, sy0), (sx1, sy1))

for p0, p1 in itertools.combinations(crops, 2):
    intersections = circle_intersection(p0, p1, R)
    if not intersections:
        continue
    for x, y in intersections:
        n = sum(1 for X, Y in crops if (x-X)*(x-X) + (y-Y)*(y-Y) <= R*R)
        if n == best:
            best_points.add((x, y))
        elif n > best:
            best = n
            best_points = {(x, y)}

print(best, best_points)