Find the optimal shift for max circle sum
I have an n by n grid of integers which sits in an infinite grid and I want to find the circular area, with origin at (x, x) (notice that both dimensions are equal), with maximum sum. I assume that x <= 0. Here is an illustration with n=10 of two circles with different centers and different radii. Only one quarter of the circles is shown as any part outside of the n by grid contributes zero to the sum.
For a fixed circle center, there is an efficient method to find the radius which maximises the sum of the integers contained within a circle. The following code solves the problem.
from collections import Counter
g = [[-3, 2, 2],
[ 2, 0, 3],
[-1, -2, 0]]
n = 3
sum_dist = Counter()
for i in range(n):
for j in range(n):
dist = i**2 + j**2
sum_dist[dist] += g[i][j]
sorted_dists = sorted(sum_dist.keys())
for i in range(1, len(sorted_dists)):
sum_dist[sorted_dists[i]] += sum_dist[sorted_dists[i-1]]
print(sum_dist)
print(max(sum_dist, key=sum_dist.get))
However, if I were to set the center to (x, x) with x < 0 instead the maximum possible sum may be different.
How can I find the center and the radius which maximizes the sum of the integers within the circle?
The data in the illustration is:
[[ 3, -1, 1, 0, -1, -1, -3, -2, -2, 2],
[ 0, 0, 3, 0, 0, -1, 2, 0, -2, 3],
[ 2, 0, 3, -2, 3, 1, 2, 2, 1, 1],
[-3, 0, 1, 0, 1, 2, 3, 1, -3, -1],
[-3, -2, 1, 2, 1, -3, -2, 2, -2, 0],
[-1, -3, -3, 1, 3, -2, 0, 2, -1, 1],
[-2, -2, -1, 2, -2, 1, -1, 1, 3, -1],
[ 1, 2, -1, 2, 0, -2, -1, -1, 2, 3],
[-1, -2, 3, -1, 0, 0, 3, -3, 3, -2],
[ 0, -3, 0, -1, -1, 0, -2, -3, -3, -1]]
I'll provide an outline of an approach that should be fairly fast.
First, design a data structure for summing up all of the points in a given area quickly. A quadtree can work. Create two copies of this. One of which has all of the values for finding the sum in an area, and the second having just the positive values to put an upper bound on the sum of the values in a region.
Next, we'll need a heap data structure. We're going to need a max heap, which can be done by throwing negatives of priorities into a min heap. So find a heap, and figure out how to use it like a max heap.
Now here is the idea. We're going to do a binary search between candidates. Each candidate will have a definitely area, a maybe area, and a definitely not area. We'll put our candidates in a heap, prioritized by the sum of the definitely area plus the positive upper bound on the maybe area. We'll take the one with the highest upper bound. When we're down to a candidate with the highest upper bound and that is also its lower bound, that's our answer.
Each candidate will exclude everything outside of one shape, and include everything inside of another. The maybe region is everything that lies between those shapes. When we examine a candidate, we pick a point inside of it and separate the candidate into one that must include that point, and one that does not include that point.
The first candidate is a choice between everything, or nothing. We put it on the heap.
Now our strategy goes like this.
while not solved:
take the best candidate off the heap.
If it has no maybe region:
we have our answer
else:
split it
Put each copy on the heap.
So we take our best candidate off the heap. And we split it on whether it contains the corner of the matrix. The candidate that includes it must also include the diagonal and all closer, and might include everything. The candidate that does not include can only include things strictly within a circle centered at the origin with the length of one of the sides. This gives us 2 candidates to put back on the heap.
As we continue, we split either on a point in the maybe region on the outside edge, or a point in the maybe region on the diagonal. We split on the one with more possible choices. Either we include the picked point, or not. As we progress, we take advantage of the fact that a circle is defined by three points - where it hits the edge, where it hits the diagonal, and where it hits the mirror copy of the edge. That allows both the yes and no regions to be defined by a circle that we can find.
Eventually we get down to knowing where the circle passes the edge, and where it passes the diagonal to within one data point. After that point we just continue splitting the interval. This specifies the circle ever more precisely - as it does forcing points out of maybe and into yes and no. It is just a question of time until we eliminate every maybe, and come up with a winning candidate that has been squeezed down to our exact answer.
I admit to not knowing how to do the complexity analysis on this. But in general I would expect it to turn into a binary search. A few candidates have to be split until we narrow in to the region that contains the winner, and then we quickly find the right answer in a logarithmic number of steps.
Here is an unoptimized solution. It internally tracks candidates by the index where the yes and maybe circles cross the edge of the matrix, and its diagonal. That actually makes 3 points because we reflect the edge.
import heapq
import math
# Convert from weird matrix coords to
# a circle.
def coords_to_xr (n, coords):
edge, diag = coords
if edge == diag*2:
# CHEAT! big dircle that is approximately diag.
# x will be None, r will be where we cross the axis.
return (None, diag)
elif edge < 0:
# Empty circle.
return (-1, 0)
elif edge <= n-1:
a = 0
b = edge
elif edge <= 2*(n-1):
a = n-1
b = edge - n + 1
else:
# Circle has everything
return (0, 2*n)
# Distance from (x, x) to (a, b) = distance from (x, x) to (diag, diag)
#
# (a-x)^2 + (b-x)^2 = = 2*(diag - x)^2
# (a^2 - 2*a*x + x^2) + (b^2 - 2*b*x + x^2) = 2*(diag^2 - 2*diag*x + x^2)
# 2*x^2 - 2*(a + b)*x + a^2 + b^2 = 2*x^2 - 4*diag*x + 2*diag^2
# 2*(2*diag - a - b)*x = 2*diag^2 - a^2 - b^2
# x = (2*diag^2 - a^2 - b^2) / 2 / (2*diag - a - b)
x = (2*diag**2 - a**2 - b**2) / 2 / (2*diag - a - b)
# r^2 = 2 * (diag - x)^2
# r = sqrt(2) * (diag - x)
#
r = (diag - x) * 2**0.5
return (x, r)
def in_circle (x, r, i, j):
if x is None:
# Our cheat.
return i + j <= 2*r
else:
return (i - x)**2 + (j - x)**2 <= r**2
def find_answer (matrix, show=False):
# To avoid floating point roundoff errors we will adjust radius by this much.
epsilon = 10**(-13)
# Candidates is a heap of:
#
# (score, yes_coords, maybe_coords)
#
# where score = (maybe_score, yes_score)
#
# And scores are the negative of the sum so the heapq min-heap becomes a max-heap.
#
# Xoords are (edge, diag). Which is the possibly fractional index
# where we cross the edge of the matrix, and its diagonal.
#
# Yes is points definitely in this candidate. Maybe is points that might be.
#
# The starting score just has to be different. The candidate is nothing yes
# and everything no.
n = len(matrix)
candidates = [((1, 0), (-1, -1), (2*n-1, n))]
while True:
score, yes_coords, maybe_coords = heapq.heappop(candidates)
if show:
print("To include: ", coords_to_xr(n, yes_coords))
print("Maybe include: ", coords_to_xr(n, maybe_coords))
print(f"Sum will be {-score[1]} to {-score[0]}")
print("")
if score[0] == score[1]:
x, r = coords_to_xr(n, yes_coords)
return (x, r, -score[0])
else:
# Do we split on the edge or the diagonal?
if maybe_coords[0] - yes_coords[0] < maybe_coords[1] - yes_coords[1]:
# Split on the diagonal.
mid = (maybe_coords[1] + yes_coords[1]) / 2
# First, let's find a new yes region with mid.
#
# We need to guarantee that edge is at least the same length from the origin.
min_r = 2**0.5 * mid
if min_r <= n-1:
# hit before the corner.
min_edge = min_r
else:
# Use Pythagorus.
min_edge = (2 * mid**2 - (n-1)**2)**0.5
mid_yes_coords = (max(yes_coords[0], min_edge), mid)
# Now, let's find a new maybe region that excludes mid.
#
# Note that a diagonal crossing at mid will hit edge at 2*mid.
# This works before and after the corner for different reasons.
mid_maybe_coords = (min(maybe_coords[0], 2*mid), mid)
else:
# Split on the edge
mid = (maybe_coords[0] + yes_coords[0]) / 2
# If mid is in yes, then the diag crossing at mid/2 must be as well.
mid_yes_coords = (mid, max(yes_coords[1], mid/2))
# If mid is in no, what is the mid r?
if mid <= n-1:
# We're not after the corner.
min_r = mid
else:
# Use Pythagorus.
min_r = ((n-1)**2 + (mid - n + 1)**2)**0.5
# Remember, 1 unit on the diagonal has length sqrt(2).
max_diag = min_r / 2**0.5
mid_maybe_coords = (mid, min(maybe_coords[1], max_diag))
# Now score our candidates and put them back on the heap.
for candidate in ((yes_coords, mid_maybe_coords), (mid_yes_coords, maybe_coords)):
this_yes_coords, this_maybe_coords = candidate
yes_x, yes_r = coords_to_xr(n, this_yes_coords)
maybe_x, maybe_r = coords_to_xr(n, this_maybe_coords)
# AVOID ROUNDOFF.
# Boundary of yes is a yes.
yes_r *= 1 + epsilon
# Boundary of maybe is a no
maybe_r *= 1 - epsilon
# Let's sum it.
yes_sum = max_sum = 0
# This is what is not optimized.
for i in range(n):
for j in range(n):
if in_circle(yes_x, yes_r, i, j):
yes_sum += matrix[i][j]
max_sum += matrix[i][j]
elif 0 < matrix[i][j] and in_circle(maybe_x, maybe_r, i, j):
max_sum += matrix[i][j]
# The negative scores are to make heapq's min-heap into a max-heap.
heapq.heappush(candidates, ((-max_sum, -yes_sum), this_yes_coords, this_maybe_coords))
matrix = [[ 3, -1, 1, 0, -1, -1, -3, -2, -2, 2],
[ 0, 0, 3, 0, 0, -1, 2, 0, -2, 3],
[ 2, 0, 3, -2, 3, 1, 2, 2, 1, 1],
[-3, 0, 1, 0, 1, 2, 3, 1, -3, -1],
[-3, -2, 1, 2, 1, -3, -2, 2, -2, 0],
[-1, -3, -3, 1, 3, -2, 0, 2, -1, 1],
[-2, -2, -1, 2, -2, 1, -1, 1, 3, -1],
[ 1, 2, -1, 2, 0, -2, -1, -1, 2, 3],
[-1, -2, 3, -1, 0, 0, 3, -3, 3, -2],
[ 0, -3, 0, -1, -1, 0, -2, -3, -3, -1]]
for r in matrix:
print(r)
print("")
print(find_answer(matrix, True))
Optimization note. The easiest thing to do is likely to precalculate sums along the diagonals with i+j constant. If you do all of that work, you can also build a data structure to calculate maximum partial sums on the diagonal bit between the yes and maybe circles. The big win there is that this is a much more reasonable estimate of the gap between what is certain (yes) and possible (maybe) than just adding the non-negative values. And therefore we can eliminate candidates more quickly.
I'm trying to get an idea of average performance. The following is a back of the envelope attempt.
First, we're defining circles by a point on the edge, and a point on the diagonal. On average the circle cuts through through O(n) pairs (i, j) in the circle, (i+1, j) not in the circle. If the position that we cut through them is random, on average you only have to move the diagonal or edge a distance of O(1/n) in order to move one point from outside to inside.
The range of possible values for the edge is O(n), and over a large part of that range the range of corresponding values for the diagonal is also O(n). As Matt Timmermans said, that corresponds to parameter space for values of (edge, diag) of area O(n^2), divided into regions of area O(1/n^2) for O(n^4) different circles in which the answer might be found.
Now let's assume that the entries in the matrix are independently and randomly distributed with variance v, average 0, and that max(0, matrix[i][j]) has average 0 < a. These are reasonable, and the purpose of having a 0 average is to give all O(n^4) of the circles a decent chance of being best.
A random circle has O(n^2) entries, and therefore its variance is O(n^2 v) leading to a standard deviation that is O(n sqrt(v)) = O(n). A more detailed analysis would have to look at how many circles they are, how independent they aren't, and so on. But if we're willing to ignore something that is likely a O(log log(n)) factor, the optimal circle is likely to have a sum that is O(n) in size. A non-optimal circle on average has a sum that is O(n) smaller, and won't be eliminated until there are O(n) points between it and its maybe.
With my search pattern, that means that we eliminate most candidates when the edge and diagonal have been found to within O(1). At this point there are only O(n^2) candidates. And I would expect that once they start to go, they go quickly. Therefore in the end we'll probably look at O(n^2) candidates.
Since examining a candidate with my algorithm is O(n^2), that's an average runtime around O(n^4). If you optimize the summing with the diagonal trick, you should improve that to an average runtime of O(n^3). Attempting to optimize the gap between yes and maybe may or may not be a win in practice.
The average performance becomes much worse if you take the average matrix and add lots of dipoles. Add a large number to a random point in the matrix, and subtract that from its neighbor. Now the average sum gets bigger based on dipoles separated by the circle. But the difference between the maybe and the yes circles is dominated by the dipoles captured in the difference. This is going to lead to having to look at a lot more candidates, and being much slower.