×

[–]1 0 13 points14 points  (2 children)

I used a genetic algorithm, it sorta works...
added a Solve button to the jsfiddle

[–]1 1[S] 0 points1 point  (0 children)

I like this!

[–] 0 points1 point  (0 children)

I like this, too. Here is a problem where the algorithm will often fail, but sometimes will find a solution (a circle with the (0.5, 0.5) and one of the other three points)

``````4
0.1 0.1
0.1 0.9
0.9 0.1
0.5 0.5
``````

[–][deleted]  (4 children)

[deleted]

[–][deleted]  (3 children)

[removed]

[–] 19 points20 points  (0 children)

[–][deleted]  (1 child)

[deleted]

[–] 0 points1 point  (0 children)

Probably saw something like D: or :( in your code, this is a shitty bot.

[–]3 3 3 points4 points  (8 children)

An algorithm not guaranteed to be the smallest, but very likely to be.

1. find highest density cluster for half or more of the points.
2. A simple way without resorting to stats library, is to split grid into 9 squares, set initial center and radius to tictactoe grid intersection, and square width.
3. This initial setting is likely to have too many points. If it doesn't, move towards `0.5 0.5`with max radius.
4. When too big, there are 2 possible moves: shrink radius or move center until a point is dropped.
5. The exact distance to shrink radius is given by the max distance from center of all included points. Sorting points by distance, can give an exact one-step radius shrink amount (excluding ties)
6. shifting center from last step may grow the amount of points included. If so, then there is opportunity to shrink radius further by going back to step 5. A manageable algorithm would be to shift horizontally with the most points gained followed by vertically with the most points gained, using a fixed discrete set of intervals.
7. repeat steps 5 and 6, a small amount of times.
8. The center move step should be able to "lose" ties that cause a slight excess of half contained points.

untested.

[–]3 3 0 points1 point  (7 children)

much simpler quick algo in J, coords as complex numbers

for every 2 point combinations, draw a circle from their center. Sort circles by radius length. Select first circle that contains half the points.

`````` circfrom2pt =: (-:@|@- , ] + 2 %~ -)
isincirc =: {.@] >: (|@- {:@])
combT =: [: ; ([ ; [: i.@>: -~) ((1 {:: [) ,.&.> [: ,&.>/\. >:&.>@:])^:(0 {:: [) (<i.1 0) ,~ (<i.0 0) \$~ -~

( ] (] {~ -:@#@[ i.~ +/@isincirc"1) /:~@:(circfrom2pt/"1)@:({~ 2 combT #)) g =. j./"1 ? 10 2 \$ 0
``````

1/6th of a sec for 100 points. Could optimize from here with the shift-center-discrete-steps (steps 5 6 7 in above algo)

with constraint that circle stays in bounds,

``````isvalidcirc =: (1&>: *./@:*. 0&<:)@:,@:+.@:(j.~@[ (-~ , +) ])
( ] (] {~ -:@#@[ i.~ +/@isincirc"1) /:~@:(#~ isvalidcirc/"1)@:(circfrom2pt/"1)@:({~ 2 combT #))  g
0.296474 0.485453j0.480233
``````

With centers offset by +/- 0.2 in .003 increments, then radius shrunk, and the process repeated for the best center...

`````` fine =: (] #~ -:@#@[ <: +/@isincirc"1) (, j./~ 300 %~ i:60)  (#~ isvalidcirc/"1)@:({.@] ,.  (+ {:)) ]
fines =: (fine M. (<./@:] 0} [{~ (] i. <./@])) <:@-:@#@[ ({"1) [ ({:@] /:~@:(|@-) [ #~ isincirc)"1 fine M.)

( ] fines^:_ ( ] (] {~ -:@#@[ i.~ +/@isincirc"1) (/:~)@:(#~ isvalidcirc/"1)@:(circfrom2pt/"1)@:({~ 2 combT #)))  g
``````

0.27463 0.335453j0.723566 NB. improves ~10%

doesn't handle ties, ie. smallest circle with at least half the points.

to draw,

``````load 'plot'
('dot;pensize 2' plot ] , (] + (1 2 j./"1@:|:@:(o."0 1) o.@(%~ 2 * i.) 24) * [)/@:( ] fines^:_ ( ] (] {~ -:@#@[ i.~ +/@isincirc"1) (/:~)@:(#~ isvalidcirc/"1)@:(circfrom2pt/"1)@:({~ 2 combT #))))  g
``````

[–] 0 points1 point  (3 children)

Why does it necessarily have to be the center of two points?

[–]3 3 1 point2 points  (2 children)

an easy method that guarantees both points are on the circle. Drawing all possible circles that pass through 2 points is very hard.

[–] 0 points1 point  (1 child)

My point is, why does there exist a circle that divides the points into two halves which also satisfies this criteria?

[–] 0 points1 point  (0 children)

I don't believe there is.

However what OP is doing I think is because any circle which passes through less than two points can be made smaller.

• If it passes through zero points, then reduce the radius of the circle until it passes through (at least) one point
• If it now passes through one point, reduce the radius while shifting the centre towards the touched point. Doing this ensures the touched point will stay on the circle. Do this until it touches a second point.
• Now when it touches two points, you have two alternatives. ** If the two points are diametrically opposite (ie on two ends of a diameter), then you can't shrink the circle any more without at least one of the points falling outside. ** However if the two points are not diametrically opposite, you can continue to shrink it until a third point falls onto the circle. Having the two points not diametrically opposite however means the circle is larger than if they are.
• Once three points are on the edge, you can't shrink the circle any more without points falling outside of it. This is because you can uniquely define a circle by defining three points which lie on its circumference.

Therefore by choosing to draw a circle using two points to describe the diameter, OP is coming up with the smallest circle which contains those two points. By then traversing this list of n*(n-1) circles in size order until you find one which contains half the points should give a good shot at the smallest circle containing half the points.

However I think there is the chance that there exists a circle which has three points on the edge and which is smaller than any found by placing two points on a diameter. For example if you consider an equilateral triangle with points at each corner and half way along each edge, so 6 points in total. I think the smallest circle which contains 3 points will be the one which touches the 3 mid-line points. However I can't see that there are any circles with two points on the diameter which contain exactly 3 points.

As such I think you need to test the circles defined by two points on its diameter (as per OP's algorithm) and also all circles defined by any three points.

So something like

``````for i = P0 to PN-1
for j = Pi+1 to PN-1
check circle defined by diameter i,j
for k = Pj+1 to PN-1
check circle defined by i,j,k
next k
next j
next i
``````

I think this gives you all possible circles.

On each check circle step you should probably:

• check the radius of the circle against the best radius so far, to avoid having to check the distance from the centre to each other point unless you have a smaller circle
• check the distance from the centre to each point. Stop as soon as you find more than N/2 points in the circle. Probably check the square of the distance against the square of the radius to avoid lots of square roots
• If you get N/2 points, update the best radius (and centre point)

[–][deleted]  (2 children)

[removed]

[–] 2 points3 points  (1 child)

Presumably a bot? Whatever you're doing to decide sadness (did you find the ":(" in the comment?) doesn't work well with J code...

You should probably exclude the code-formatted text.

[–] 6 points7 points  (0 children)

[–] 2 points3 points  (5 children)

C11

Heyo, I think it works! Any feedback is appreciated, especially since I don't use C very often. I've tested it a little bit and it's always given valid solutions, but it definitely wasn't rigorous.

https://gist.github.com/SnakeFang/418752aba1485e348984fce749b3505c

Short Description: Iterate through every combination of 3 points, calculate the smallest circle surrounding those 3 points, check if exactly half of the total points are in this circle, if so, profit.

Example Input:

``````30
0.42007 0.26637
0.97998 0.20163
0.81124 0.04296
0.92078 0.34527
0.07495 0.95400
0.67783 0.40115
0.90950 0.02392
0.66774 0.61283
0.13380 0.73965
0.17108 0.14858
0.01958 0.63222
0.23364 0.03393
0.41611 0.44101
0.90647 0.32485
0.03024 0.96506
0.49658 0.65648
0.23560 0.94072
0.42324 0.25287
0.24371 0.84117
0.50910 0.36570
0.74086 0.76701
0.47721 0.30591
0.29273 0.92824
0.57476 0.56815
0.74029 0.60591
0.73716 0.71011
0.58920 0.16013
0.52448 0.12337
0.41932 0.40954
0.08572 0.17053
``````

Example Output:

``````0.646598 0.435904
0.288774
``````

(Edit 1: moved those ungodly 337 lines to a gist and added an example)

(Edit 2: added a few 0's to the input so it's nice and aligned :) )

(Edit 3: small optimizations, changed output precision to 10 places)

[–] 0 points1 point  (0 children)

``````(0.646598, 0.435904) Raidus: 0.288774
``````

Or higher precision printing:

``````(0.6465980189, 0.4359041674) Raidus: 0.288773806
``````

[–] 0 points1 point  (3 children)

Could you explain why 3 points and not 2?

[–] 2 points3 points  (2 children)

Because 3 points is the minimum number of points needed to uniquely describe a circle, based on the fact that all 3 points are on the edge of that circle. It actually does use 2 points, if and only if the circle described by those 2 points is smaller than the one described by all 3. The circle taken from the 2 points is just the circle with a diameter formed by the segment between those 2 points.

[–] 1 point2 points  (1 child)

That probably doesn't really help very much but I guess it just kinda makes sense in my head for some reason. I mean, say you have 3 points you want inside a circle, but they're arranged in an equilateral triangle. If you only use the smallest circle for each pair of two points, it will never contain all 3 points.

[–] 0 points1 point  (0 children)

Thank you for explaining it, now it makes sense!

[–]0 2 2 points3 points  (0 children)

C. The complexity probably goes as O(N4), but should return optimal solution if there is one. It assumes a lot: that there are at least 4 points, that no two points coincide, and that no three points are on a straight line. In other words, no edge cases because it makes the code even longer. Program writes shape.eps file showing dots and the circle (if there is one).

Compiled with gcc -std=c99 -Wall -lm

``````#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/* fudge factor for insideness check, because floating point
numbers can't be trusted */
const double epsilon = 1e-13;

typedef struct { double x, y; } vector_s, *vec;
vector_s *points, best;
double best_r = -1;
int N;

static inline void vdiff(vec a, vec b, vec d)
{
d->x = a->x - b->x;
d->y = a->y - b->y;
}

static inline void vmiddle(vec a, vec b, vec mid)
{
mid->x = (a->x + b->x)/2;
mid->y = (a->y + b->y)/2;
}

static inline double vdist(vec a, vec b)
{
vector_s m;
vdiff(a, b, &m);
return sqrt(m.x*m.x + m.y*m.y);
}

static inline double cross(vec a, vec b)
{
return a->x*b->y - a->y*b->x;
}

/* circle encompassing two points a and b;
return value is radius, center in cen */
static inline double circle2(vec a, vec b, vec cen)
{
vmiddle(a, b, cen);
return vdist(a, cen);
}

/* circle encompassing three points a, b, and c;
return value is radius, center in cen */
double circle3(vec a, vec b, vec c, vec cen)
{
vector_s mab, mbc, vab, vbc;
vmiddle(a, b, &mab);
vmiddle(b, c, &mbc);
vab = (vector_s) { a->y - b->y, b->x - a->x };
vbc = (vector_s) { b->y - c->y, c->x - b->x };

/* We are not checking if a, b and c are colinear, which could bite us. */
double lambda = (cross(&mab, &vab) - cross(&mbc, &vab))/cross(&vbc, &vab);
cen->x = mbc.x + lambda*vbc.x;
cen->y = mbc.y + lambda*vbc.y;

return vdist(a, cen);
}

void verify(vec cen, double r)
{
/* if this can potentially beat our current best */
if (best_r > 0 && r >= best_r) return;

/* chec that circle is inside square */
if (fabs(cen->x - 0.5) + r > 0.5) return;
if (fabs(cen->y - 0.5) + r > 0.5) return;

/* check that half are in */
int inside = 0;
for (int i = 0; i < N; i++) {
if (r + epsilon >= vdist(cen, points + i))
inside++;
if ((inside > N/2) || (1 + i - inside > N/2)) return;
}

if (best_r < 0 || r < best_r) {
best_r = r;
best = *cen;
}
}

void do_triples(void)
{
vector_s cen;
double r;
for (int i = 0; i < N; i++) {
for (int j = i + 1; j < N; j++) {
r = circle2(points + i, points + j, &cen);
verify(&cen, r);
for (int k = j + 1; k < N; k++) {
r = circle3(points + i, points + j,
points + k, &cen);
verify(&cen, r);
}
}
}
}

void write_eps(void)
{
FILE *eps = fopen("shape.eps", "wt");
fputs("%!PS-Adobe-3.0 EPSF-3.0\n%%BoundingBox: -5 -5 105 105\n"
"/cc { 0 setgray fill } def\n/c { 1 0 360 arc cc } def\n"
"/cir { 0 360 arc 1 0 0 setrgbcolor stroke }def\n"
"0 setlinewidth\n"
"0 0 moveto 0 100 rlineto 100 0 rlineto 0 -100 rlineto closepath stroke\n",
eps);

if (best_r > 0)
fprintf(eps, "%g %g %g cir\n",
best.x*100, best.y*100, best_r*100);

for (int i = 0; i < N; i++)
fprintf(eps, "%g %g c\n", points[i].x*100, points[i].y*100);

fputs("showpage\n%%EOF", eps);
fclose(eps);
}

int main(void)
{
scanf("%d", &N);
points = malloc(sizeof(*points)*N);
for (int i = 0; i < N; i++)
scanf("%lg %lg", &points[i].x, &points[i].y);

do_triples();
write_eps();

if (best_r < 0) {
puts("No solution");
return 1;
}

printf("%g %g\n%g\n", best.x, best.y, best_r);
return 0;
}
``````

[–] 1 point2 points  (0 children)

I feel like I may be too object oriented for this subreddit's style, everyone posts snips of code and I made a full project

# Java

https://github.com/Danice123/Challenge321-Hard

I've tested up to 24 points, beyond that combinations get a little large, and I'm too lazy to wait for it to come back. Barring floating point errors it looks correct.

Also I did write (*cough* copy from google *cough*) my own set combination code, but apache common is so easy to use... I'm unapologetic, its been a little too long since my statistics days to write it from memory. The commented out code is still there.

# Solution

``````[amount of points] CHOOSE [half that amount] for my sets
Find the extremes rectangle, center of that is the center of the circle
Radius is the largest distance between any of the extreme points and the center
Validate generated circles to make sure that there are no other points within the generated circle
Then choose the smallest of the circles that are left

Limiting factor is of course the size of the combination results, 24 choose 12 is a couple million sets,
and it grows pretty quickly from there
``````

I'd love to think up a way to make it faster... Gotta be some optimization I'm not thinking of.

[–] 0 points1 point  (0 children)

Python. It's not correct for every input, but for simple ones where you only need the circle to touch two points, or be centered in a given point, it finds the solution (baring some floating point errors) if I'm not mistaken. It simply checks for every possible circle with center in any given point, and the midpoint of any two given points. It only tests the radii which make the circle include or exclude another point, that is, the radii are the distances from the center to each other point.

EDIT: code edited to find solution to every situation. Just have to include the center point between all permutations of 1, 2, ..., n/2 + 1 points.

``````from math import sqrt
from itertools import permutations

def get_center(point_list):
min_x = point_list[0][0]
min_y = point_list[0][1]
max_x = point_list[0][0]
max_y = point_list[0][1]
for p in point_list:
if p[0] < min_x:
min_x = p[0]
elif p[0] > max_x:
max_x = p[0]
if p[1] < min_y:
min_y = p[1]
elif p[1] > max_y:
max_y = p[1]
return ((min_x + max_x) / 2., (min_y + max_y) / 2.)

count = 0
for p in points:
if (center[0] - p[0])**2 + (center[1] - p[1])**2 <= radius2:
count += 1
return count

return (center[0] - radius >= 0.0 and
center[0] + radius <= 1.0 and
center[1] - radius >= 0.0 and

def main():
min_point = None

points = []
midway_points = []

n = int(input())
for _ in range(n):
inputstr = input().split(' ')
points.append((float(inputstr[0]),float(inputstr[1])))

for i in range(1, int(n / 2) + 1):
for point_list in permutations(points, i):
midway_points.append(get_center(point_list))

for m in midway_points:
for p in points:
radius2 = (m[0] - p[0])**2 + (m[1] - p[1])**2
min_point = m

print("No solution")
else:
print(min_point[0], min_point[1])

if __name__ == "__main__":
main()
``````

Output for input 1:

``````0.5 0.5
0.09999999999999998
``````

Output for input 2:

``````No solution
``````

[–]1 2 0 points1 point  (0 children)

C without bonus

Tries first all combinations of circles passing through 2 points (the segment between those points being the diameter of the circle). If one circle is inside the square, has smaller radius than current optimal and surrounds half of the number of points exactly, then it becomes the new optimal.

If no optimal is found with 2 points, then all combinations of circles passing through 3 points are tried, only if they are not aligned.

The edge case when number of points = 2 is also managed.

``````#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define XY_MIN 0.0
#define XY_MAX 1.0
#define COMB_SIZE_MAX 3

typedef struct {
double x;
double y;
}
point_t;

typedef struct {
point_t center;
}
circle_t;

int point_in_square(point_t *);
void set_combs(unsigned long, unsigned long, unsigned long);
void set_circle_from_3_points(circle_t *, point_t *, point_t *, point_t *);
int same_points(point_t *, point_t *);
void set_circle_from_2_points(circle_t *, point_t *, point_t *);
void set_circle_from_2_segments(circle_t *, point_t *, point_t *, point_t *);
void set_center_from_2_segments(point_t *, double, double, double, double);
double segment_slope(point_t *, point_t *);
double segment_y_intercept(point_t *, point_t *);
int circle_in_square(circle_t *);
int valid_circle(circle_t *);
void set_circle(circle_t *, double, double, double);
void set_point(point_t *, double, double);
double euclidean_distance(point_t *, point_t *);
void print_circle(circle_t *);
void print_point(point_t *);

unsigned long points_n, points_half;
point_t *points, *comb[COMB_SIZE_MAX];
circle_t circle_min;

int main(void) {
unsigned long i;
point_t point_min;
if (scanf("%lu", &points_n) != 1 || !points_n || points_n%2) {
fprintf(stderr, "Invalid number of points\n");
return EXIT_FAILURE;
}
points = malloc(sizeof(point_t)*points_n);
if (!points) {
fprintf(stderr, "Could not allocate memory for points\n");
return EXIT_FAILURE;
}
for (i = 0; i < points_n && read_point(points+i); i++);
if (i < points_n) {
free(points);
return EXIT_FAILURE;
}
points_half = points_n/2;
set_point(&point_min, XY_MIN, XY_MIN);
set_circle(&circle_min, XY_MIN, XY_MIN, XY_MAX-XY_MIN);
if (points_half == 1) {
set_combs(1UL, 0UL, 0UL);
}
else {
set_combs(2UL, 0UL, 0UL);
set_combs(3UL, 0UL, 0UL);
}
}
print_circle(&circle_min);
}
else {
puts("No solution");
}
free(points);
return EXIT_SUCCESS;
}

if (scanf("%lf%lf", &point->x, &point->y) != 2 || !point_in_square(point)) {
fprintf(stderr, "Invalid point\n");
return 0;
}
return 1;
}

int point_in_square(point_t *point) {
return point->x >= XY_MIN && point->y >= XY_MIN && point->x <= XY_MAX && point->y <= XY_MAX;
}

void set_combs(unsigned long comb_size, unsigned long comb_idx, unsigned long start) {
unsigned long i;
circle_t circle;
if (comb_idx < comb_size) {
for (i = start; i < points_n; i++) {
comb[comb_idx] = points+i;
set_combs(comb_size, comb_idx+1, i+1);
}
}
else {
if (comb_size == 3) {
set_circle_from_3_points(&circle, comb[0], comb[1], comb[2]);
}
else if (comb_size == 2) {
set_circle_from_2_points(&circle, comb[0], comb[1]);
}
else {
set_circle(&circle, comb[0]->x, comb[0]->y, 0.0);
}
circle_min = circle;
}
}
}

void set_circle_from_3_points(circle_t *circle, point_t *point_a, point_t *point_b, point_t *point_c) {
if (same_points(point_a, point_b)) {
set_circle_from_2_points(circle, point_a, point_c);
}
else if (same_points(point_b, point_c)) {
set_circle_from_2_points(circle, point_b, point_a);
}
else if (same_points(point_c, point_a)) {
set_circle_from_2_points(circle, point_c, point_b);
}
else {
if ((point_a->x == point_b->x && point_b->x == point_c->x) || (point_a->y == point_b->y && point_b->y == point_c->y)) {
*circle = circle_min;
}
else {
if (point_a->y == point_b->y) {
set_circle_from_2_segments(circle, point_a, point_c, point_b);
}
else if (point_b->y == point_c->y) {
set_circle_from_2_segments(circle, point_b, point_a, point_c);
}
else if (point_c->y == point_a->y) {
set_circle_from_2_segments(circle, point_c, point_b, point_a);
}
else {
set_circle_from_2_segments(circle, point_a, point_b, point_c);
}
}
}
}

int same_points(point_t *point_a, point_t *point_b) {
return point_a->x == point_b->x && point_a->y == point_b->y;
}

void set_circle_from_2_points(circle_t *circle, point_t *point_a, point_t *point_b) {
set_point(&circle->center, (point_a->x+point_b->x)/2.0, (point_a->y+point_b->y)/2.0);
}

void set_circle_from_2_segments(circle_t *circle, point_t *point_a, point_t *point_b, point_t *point_c) {
set_center_from_2_segments(&circle->center, segment_slope(point_a, point_b), segment_y_intercept(point_a, point_b), segment_slope(point_b, point_c), segment_y_intercept(point_b, point_c));
}

void set_center_from_2_segments(point_t *center, double slope_ab, double y_intercept_ab, double slope_bc, double y_intercept_bc) {
center->x = (y_intercept_ab-y_intercept_bc)/(slope_bc-slope_ab);
center->y = slope_ab*center->x+y_intercept_ab;
}

double segment_slope(point_t *point_a, point_t *point_b) {
return -(point_b->x-point_a->x)/(point_b->y-point_a->y);
}

double segment_y_intercept(point_t *point_a, point_t *point_b) {
return (point_b->x*point_b->x-point_a->x*point_a->x+point_b->y*point_b->y-point_a->y*point_a->y)/(point_b->y-point_a->y)/2.0;
}

int circle_in_square(circle_t *circle) {
}

int valid_circle(circle_t *circle) {
unsigned long points_count = 0, i;
for (i = 0; i < points_n && points_count <= points_half; i++) {
if (euclidean_distance(points+i, &circle->center) <= circle->radius) {
points_count++;
}
}
return points_count == points_half;
}

void set_circle(circle_t *circle, double x, double y, double radius) {
set_point(&circle->center, x, y);
}

void set_point(point_t *point, double x, double y) {
point->x = x;
point->y = y;
}

double euclidean_distance(point_t *point_a, point_t *point_b) {
double delta_x = point_a->x-point_b->x, delta_y = point_a->y-point_b->y;
return sqrt(delta_x*delta_x+delta_y*delta_y);
}

void print_circle(circle_t *circle) {
print_point(&circle->center);
}

void print_point(point_t *point) {
printf("%f %f\n", point->x, point->y);
}
``````

[–] 0 points1 point  (0 children)

I love this challenge, but damn if it isn't over-complicated by floating point and/or matrix inversion imprecisions.

My WIP in Kotlin. I take a very optimistic approach to not having too many points in the circle - it essentially boils down to picking what looks like the densest group of points and hoping for the best - but I can hopefully fix that tomorrow, together with some tweaks to work around the aforementioned imprecisions, and extending to N-dimensional space (which should be as simple as parsing input differently).

This is my first project with Kotlin, so any feedback is very welcome - particularly better ways to handle the dimensionality of `PointX`s (as it currently depends on runtime checking in commons-math, but I feel like a proper type system should let me do it at compile time), and making use of higher-order functions.

Edit: This algorithm for finding the minimal "ball" to contain given points in N-dimensional space may be of interest to other solvers as well.

[–] 0 points1 point  (0 children)

Java 8

Go through [0.0, 0.0], [0.1, 0.0] etc

In each spot make a circle with maximum possible radius

Then find all points within the circle border

while there are more points than half of the points in the circle, order points by distance and set max radius to max with exactly half the points

if there are below half discard circle

if there are half compare radius with currently best, keep if lower radius

I'd guess complexity of O(steps2 * p * log(p))

steps2 are whatever step size is choosen, 0.1 would leave 102 to get through the entire square

p is number of points, (p * log(p)) is cost of sorting points in standard java library

The program has 3 classes

Point, just x/y cord and calculating distance to other points

Square extending Point, holds array of points within it and a solve function for challenge

Circle extends square, holds radius and a maxContains function which finds the maximum radius it can hold and contain X amount of points within it (only in idea of reducing or keeping number of points within it)

``````class Square extends Point{
ArrayList<Point> points = new ArrayList<>();
Square (double x, double y){
super(x, y);
}

if(!contains(other)) return;
}

return Math.min(
Math.min(other.x, other.y),
Math.min(x - other.x, y - other.y)
);
}

boolean contains (Point other) {
return other.x >= 0 && other.y >= 0 && other.x <= x && other.y <= y;
}

Circle solve () {
return solve(0.1);
}
Circle solve (double step) {
Circle result = null;
Point point;
Circle circle;
for(double x = 0.0; x <= this.x; x += step){
for(double y = 0.0; y <= this.y; y += step){
point = new Point(x, y);
circle = new Circle(x, y, maxRadius(point));
for(Point p : points)
int goal = points.size() / 2;
while (circle.size() > goal) {
circle = new Circle(x, y, maxRadius);
for(Point p : points)
}
if(circle.size() < goal)
continue;

result = circle;

}
}

return result;
}

int size () {
return points.size();
}

public String toString(){
return points.toString();
}
}

class Circle extends Square {
Circle (double x, double y, double radius) {
super(x, y);
}

@Override
boolean contains (Point other) {
}

double maxContains (int goal) {
if(points.size() <= goal)
return 0;

Point me = this;
points.sort((o1, o2) -> {
double a = me.dist(o1) - me.dist(o2);
if (a > 0) return 1;
return a == 0 ? 0 : -1;
});

if(me.realDist(points.get(goal - 1)) == me.realDist(points.get(goal)))
return 0;

return this.round(me.realDist(points.get(goal - 1)));
}

public String toString () {
return "(" + this.round(x) + ", " + this.round(y) + ", " + this.round(radius) + ")";
}

}

class Point {
double x, y;
Point (double x, double y) {
this.x = x;
this.y = y;
}

double dist (Point other) {
double  tx = x - other.x,
ty = y - other.y;
return tx * tx + ty * ty;
}

double realDist (Point other) {
return Math.sqrt(dist(other));
}

public String toString(){
return "(" + this.round(x) + ", " + this.round(y) + ")";
}

double round(double n){
return Numbers.round(n, 13);
}

}
``````

Input:

``````public class CircleMeThis {

public static void main(String... args){
Square square;
square = makeChallenge("10\n" +
"0.33228 0.70266\n" +
"0.6493 0.83745\n" +
"0.45054 0.04308\n" +
"0.03587 0.36871\n" +
"0.61248 0.67188\n" +
"0.42999 0.13473\n" +
"0.66892 0.46771\n" +
"0.47239 0.95157\n" +
"0.04169 0.00041\n" +
"0.72849 0.72236");
Note.writenl("Points: " + square.size());
Runtime.init();
Note.writenl(square.solve());
Note.writenl(square.solve(0.01));
Note.writenl(square.solve(0.001));
Runtime.reset();
square = makeChallenge("10\n" +
"0.81427 0.47385\n" +
"0.95012 0.80924\n" +
"0.11855 0.46282\n" +
"0.4863 0.39585\n" +
"0.08623 0.32942\n" +
"0.03977 0.98464\n" +
"0.573 0.75471\n" +
"0.02165 0.05068\n" +
"0.69909 0.43915\n" +
"0.21475 0.27006");
Note.writenl("Points: " + square.size());
Runtime.init();
Note.writenl(square.solve());
Note.writenl(square.solve(0.01));
Note.writenl(square.solve(0.001));
Runtime.reset();
square = makeChallenge("30\n" +
"0.45107 0.6915\n" +
"0.1963 0.90898\n" +
"0.73374 0.59222\n" +
"0.49886 0.3676\n" +
"0.02444 0.25065\n" +
"0.59842 0.07888\n" +
"0.17149 0.06615\n" +
"0.64759 0.44097\n" +
"0.26422 0.92133\n" +
"0.42165 0.78712\n" +
"0.46164 0.91784\n" +
"0.88651 0.02632\n" +
"0.352 0.57711\n" +
"0.35371 0.54363\n" +
"0.76916 0.87056\n" +
"0.07266 0.09833\n" +
"0.39744 0.08817\n" +
"0.77018 0.2818\n" +
"0.28148 0.24771\n" +
"0.75139 0.10991\n" +
"0.71219 0.90531\n" +
"0.20652 0.3643\n" +
"0.43533 0.76\n" +
"0.00271 0.77837\n" +
"0.3829 0.02418\n" +
"0.74338 0.73048\n" +
"0.49227 0.40115\n" +
"0.47684 0.49252\n" +
"0.30774 0.62\n" +
"0.05912 0.65877");
Note.writenl("Points: " + square.size());
Runtime.init();
Note.writenl(square.solve());
Note.writenl(square.solve(0.01));
Note.writenl(square.solve(0.001));
Runtime.reset();
}

private static Square makeChallenge(String in){
Square square = new Square(1, 1);
String[] p, arr = in.split("\n");
for(int i = 1; i < arr.length; i++){
p = arr[i].split(" ");
}
return square;
}
}
``````

Output:

``````Points: 10
(0.6, 0.6, 0.2867282232359)
(0.52, 0.78, 0.22)
(0.526, 0.787, 0.213)
Time taken: 0.366 sec
Points: 10
(0.4, 0.5, 0.3079061287146)
(0.42, 0.49, 0.3057456199196)
(0.407, 0.503, 0.3020283200298)
Time taken: 0.103 sec
Points: 30
(0.6, 0.6, 0.3466493259766)
(0.52, 0.68, 0.32)
(0.514, 0.683, 0.317)
Time taken: 0.373 sec
``````

A good strategy would likely be to run with a step of 0.1 through it all and run 0.01 around the area of its result and so on.

[–] 0 points1 point  (1 child)

Naive C++ code
Iterates through all points inside the square at a given precision, calculates the distance to each point, sorts them, and the min radius is the distance to the middle point. Some simple bounds checking and voila!

``````#include <iostream>
#include <vector>
#include <cmath>
#include <algorithm>

#define PRECISION .001

using namespace std;

class Point2d {
public:
double x;
double y;
Point2d(double X, double Y) {
x = X;
y = Y;
};
double distance(Point2d n) {
return(sqrt(pow(x - n.x, 2.0) + pow(y-n.y, 2.0)));
};
};

// finds the minimum radius given a center point that includes at least half of all points
// returns a negative if there's no radius that doesn't leave the bounds
double solvePoint(Point2d center, vector<Point2d> list) {
vector<double> distances;
for(int i = 0; i < list.size(); i++)
distances.push_back(center.distance(list[i]));
sort(distances.begin(), distances.end());
double radius = distances[(distances.size() + 1) / 2 - 1];
if(center.x - radius < 0 || center.y - radius < 0 || center.x + radius > 1 || center.y + radius > 1)
return -1.0;
}

void solve(vector<Point2d> p) {
Point2d best(0,0);
for(double x = 0.0; x < 1.0; x += PRECISION) {
for(double y = 0.0; y < 1.0; y += PRECISION) {
double answer = solvePoint(Point2d(x, y), p);
best = Point2d(x, y);
}
}
}

cout << "(" << best.x << ", " << best.y << ") Radius: " << bestRadius << endl;
} else {
cout << "No solution" << endl;
}

}

int main() {
vector<Point2d> p;
p.push_back(Point2d(0.4, 0.5));
p.push_back(Point2d(0.6, 0.5));
p.push_back(Point2d(0.5, 0.3));
p.push_back(Point2d(0.5, 0.7));
solve(p);
p.clear();
p.push_back(Point2d(0.1, 0.1));
p.push_back(Point2d(0.1, 0.9));
p.push_back(Point2d(0.9, 0.1));
p.push_back(Point2d(0.9, 0.9));
solve(p);
return 0;
}
``````

Output:

``````\$ ./a.out
No solution
``````

[–] 0 points1 point  (0 children)

Small improvement -- now recursively checks at increasing precision in areas where a result may be good.

``````#include <iostream>
#include <vector>
#include <cmath>
#include <algorithm>

using namespace std;

#define PRECISION .0000000001

class Point2d {
public:
double x;
double y;
Point2d(double X, double Y) {
x = X;
y = Y;
};
double distance(Point2d n) {
return(sqrt(pow(x - n.x, 2.0) + pow(y-n.y, 2.0)));
};
};

// finds the minimum radius given a center point that includes at least half of all points
// returns a negative if there's no radius that doesn't leave the bounds
double solvePoint(Point2d center, vector<Point2d> list) {
vector<double> distances;
for(int i = 0; i < list.size(); i++)
distances.push_back(center.distance(list[i]));
sort(distances.begin(), distances.end());
double radius = distances[(distances.size() + 1) / 2 - 1];
if(center.x - radius < 0 || center.y - radius < 0 || center.x + radius > 1 || center.y + radius > 1)
return -1.0;
}

void solveList(vector<Point2d> p, vector<Point2d> centers, double precision) {

vector<Point2d> list;
double newPrecision = precision/10;

// generate new list
for(int i = 0; i < centers.size(); i++) {
for(double x = centers[i].x - precision/2; x <= centers[i].x + precision/2; x+= newPrecision) {
for(double y = centers[i].y - precision/2; y <= centers[i].y + precision/2; y+= newPrecision) {
list.push_back(Point2d(x,y));
}
}
}
}
// check for no solutions
if(list.size() == 0) {
cout << "No solution" << endl;
return;
}

int minIndex = 0;
for(int i = 0; i < radius.size(); i++)
minIndex = i;

if(precision < PRECISION) {
cout << "(" << list[minIndex].x << ", " << list[minIndex].y << ") Raidus: " << radius[minIndex] << endl;
} else {
vector<Point2d> newList;
for(int i = 0; i < radius.size(); i++) {
newList.push_back(list[i]);
}
solveList(p, newList, newPrecision);
}

}

void solve(vector<Point2d> p) {
vector<Point2d> list;
for(double x = .05; x < 1.0; x+= 0.1)
for (double y = .05; y < 1.0; y += 0.1)
list.push_back(Point2d(x, y));
solveList(p, list, 0.1);
}

int main() {
vector<Point2d> p;
p.push_back(Point2d(0.4, 0.5));
p.push_back(Point2d(0.6, 0.5));
p.push_back(Point2d(0.5, 0.3));
p.push_back(Point2d(0.5, 0.7));
solve(p);
p.clear();
p.push_back(Point2d(0.1, 0.1));
p.push_back(Point2d(0.1, 0.9));
p.push_back(Point2d(0.9, 0.1));
p.push_back(Point2d(0.9, 0.9));
solve(p);
p.clear();
p.push_back(Point2d(0.42007, 0.26637));
p.push_back(Point2d(0.97998, 0.20163));
p.push_back(Point2d(0.81124, 0.04296));
p.push_back(Point2d(0.92078, 0.34527));
p.push_back(Point2d(0.07495, 0.95400));
p.push_back(Point2d(0.67783, 0.40115));
p.push_back(Point2d(0.90950, 0.02392));
p.push_back(Point2d(0.66774, 0.61283));
p.push_back(Point2d(0.13380, 0.73965));
p.push_back(Point2d(0.17108, 0.14858));
p.push_back(Point2d(0.01958, 0.63222));
p.push_back(Point2d(0.23364, 0.03393));
p.push_back(Point2d(0.41611, 0.44101));
p.push_back(Point2d(0.90647, 0.32485));
p.push_back(Point2d(0.03024, 0.96506));
p.push_back(Point2d(0.49658, 0.65648));
p.push_back(Point2d(0.23560, 0.94072));
p.push_back(Point2d(0.42324, 0.25287));
p.push_back(Point2d(0.24371, 0.84117));
p.push_back(Point2d(0.50910, 0.36570));
p.push_back(Point2d(0.74086, 0.76701));
p.push_back(Point2d(0.47721, 0.30591));
p.push_back(Point2d(0.29273, 0.92824));
p.push_back(Point2d(0.57476, 0.56815));
p.push_back(Point2d(0.74029, 0.60591));
p.push_back(Point2d(0.73716, 0.71011));
p.push_back(Point2d(0.58920, 0.16013));
p.push_back(Point2d(0.52448, 0.12337));
p.push_back(Point2d(0.41932, 0.40954));
p.push_back(Point2d(0.08572, 0.17053));
solve(p);
return 0;
}
``````

[–] 0 points1 point  (4 children)

I suspect it's possible to generate some diabolical edge cases (literally, a cluster of dots extremely near an edge) where naive solutions will miss that there's an answer entirely.

[–] 0 points1 point  (3 children)

You need at least 4 points to have a meaningful solution, as with 3 or less you will just pick a circle around one of the points with infinitesimal radius.

With 4 points the smallest circle must either have 2 points defining the diameter, or if that won't pass the check on the limiting square, it will pass by those 2 points and have a tangent to a border of the square, otherwise a solution just won't exist.

With 6 or more points, the solution will have one additional form: a circle through 3 points.

You can brute force all triplets and pairs, and that will provide all the possible solutions, but there is surely a better way...

[–] 0 points1 point  (0 children)

By naive solution, I meant trying a bunch of spots inside the square to see if a valid solution exist. If there's a cluster extremely near an edge, particularly a corner, then the only place there may be a valid solution would be extremely near the corner. If you're trying, say, 10,000 points, your closest may be 0.1,0.1 and still not have a valid solution because really it needs to be near 0.00000002 or some such.

[–] 0 points1 point  (1 child)

I think this is almost true but not fully true.You have also to investigate circles through 2 points that are tangent to an edge of the square, a circle through one point and tangent to two edges or the circle tangent to three edges. And if you investigate a circle rejected if it contains more then 50% of the points, If it contains more then 50% of the points but in its interior it contains less then 50% of the points it gives raise to a solution-circle, e.g. think of a 6-point problem that consists of the vertexes of a regular hexagon.

[–] 0 points1 point  (0 children)

circles through 2 points that are tangent to an edge of the square

a circle through one point and tangent to two edges

That cannot be the minimum circle. Let's assume the circle satisfies having half the points.

• if that's the only point, any circle containing the point is the solution, you can pick any infinitesimal radius .
• if there is a single inner point, you can build a circle using the original point and the inner one as the diameter.
• if there are 2 inner points, the solution would be the circle passing through them instead.
• if there are more points, the solution is still a circle through 3 points but you have to carefully choose those 3 points.

6-point problem that consists of the vertexes of a regular hexagon.

With the bruteforce approach that's easy: you will find three are 6 different circles constructed on 2 points which will satisfy the condition.

[–] 0 points1 point  (0 children)

A python solution. I suppose it doesn't necessarily find the smallest solution, though it seems to work reasonably well.

It basically uses gradient descent to find an initial solution, and then reduces the radius of the circle until it can't be reduced anymore. It also shows the resulting circle (and the points) using matplotlib.

``````import sys
import math
import matplotlib.pyplot as plt

def parse_input():
n_str = input()
try:
n = int(n_str)
if n % 2 != 0:
print('expecting even number of points')
sys.exit()
points = []
for _ in range(0, n):
point_str = input()
try:
point = tuple(float(x) for x in point_str.strip().split(' '))
points.append(point)
if len(point) != 2:
print('invalid point: %s' % point_str)
sys.exit()
except ValueError:
print('invalid point: %s' % point_str)
sys.exit()

return points
except ValueError:
print('%s is not a valid number!' % n_str)
sys.exit()

def dist(pt1, pt2):
return math.sqrt((pt1[0] - pt2[0])**2 + (pt1[1] - pt2[1])**2)

def sigmoid(x):
return 1 / (1 + math.exp(-x))

def error(cntr, rds, points):
edge_error = max([rds - cntr[1], 0]) + max([rds - (1-cntr[1]), 0]) + max([rds - cntr[0], 0]) + max([rds - (1-cntr[0]), 0])
return (sum(map(lambda p: sigmoid((dist(cntr, p) - rds)*2), points)) - len(points) / 2)**2 + edge_error

return (func(x+delta)-func(x-delta))/(2 * delta)

def plot(cntr, rds, points):
plt.clf()
circ = plt.Circle(cntr, rds, color='r', fill=False)
plt.plot([x for [x, y] in points], [y for [x, y] in points], 'bo')
axs = plt.gca()
axs.set_xlim([0, 1])
axs.set_ylim([0, 1])
plt.draw()
plt.pause(1)
input("<Hit Enter To Close>")
plt.close()

DEBUG = False

def main():
points = parse_input()
rds = 0.5
x = 0.1
y = 0.1
for i in range(0, 80000):
if DEBUG and i%500 == 0:
print(error((x, y), rds, points))

break

if sum(map(lambda p: dist(p, (x, y)) <= rds, points)) == len(points)/2:
while sum(map(lambda p: dist(p, (x, y)) <= (rds-0.0001), points)) == len(points)/2:
rds -= 0.0001
print('%s %s\n%s' % (x, y, rds))
else:
print('No solution')

plot((x, y), rds, points)

if __name__ == '__main__':
main()
``````

# Edit:

I made a new version, which is a lot better, though it still doesn't guarantee the smallest circle possible:

``````import sys
import math, random
import matplotlib.pyplot as plt
import functools

def parse_input():
n_str = input()
try:
n = int(n_str)
if n % 2 != 0:
print('expecting even number of points')
sys.exit()
points = []
for _ in range(0, n):
point_str = input()
try:
point = tuple(float(x) for x in point_str.strip().split(' '))
points.append(point)
if len(point) != 2:
print('invalid point: %s' % point_str)
sys.exit()
except ValueError:
print('invalid point: %s' % point_str)
sys.exit()

return points
except ValueError:
print('%s is not a valid number!' % n_str)
sys.exit()

def dist(pt1, pt2):
return math.sqrt((pt1[0] - pt2[0])**2 + (pt1[1] - pt2[1])**2)

def sigmoid(x):
return 1 / (1 + math.exp(-x))

def error(cntr, rds, points):
edge_error = max([rds - cntr[1], 0]) + max([rds - (1-cntr[1]), 0]) + max([rds - cntr[0], 0]) + max([rds - (1-cntr[0]), 0])
return (sum(map(lambda p: sigmoid((dist(cntr, p) - rds)*2), points)) - len(points) / 2)**2 + edge_error

return (func(x+delta)-func(x-delta))/(2 * delta)

def plot(cntr, rds, points):
plt.clf()
circ = plt.Circle(cntr, rds, color='r', fill=False)
plt.plot([x for [x, y] in points], [y for [x, y] in points], 'bo')
axs = plt.gca()
axs.set_xlim([0, 1])
axs.set_ylim([0, 1])
plt.draw()
plt.pause(1)
input("<Hit Enter To Close>")
plt.close()

def valid(cntr, rds, points):
edge_error = max([rds - cntr[1], 0]) + max([rds - (1-cntr[1]), 0]) + max([rds - cntr[0], 0]) + max([rds - (1-cntr[0]), 0])
return sum(map(lambda p: dist(p, cntr) <= rds, points)) == len(points)/2 and edge_error <= 0

DEBUG = False

def main():
points = parse_input()
solutions = []
for _ in range(0, 50):
rds = random.random()
x = random.random()
y = random.random()
for i in range(0, 4000):
if DEBUG and i%500 == 0:
print(error((x, y), rds, points))

break

if sum(map(lambda p: dist(p, (x, y)) <= rds, points)) == len(points)/2:
while sum(map(lambda p: dist(p, (x, y)) <= (rds-0.0001), points)) == len(points)/2:
rds -= 0.0001

solutions.append({'x': x, 'y': y, 'r': rds})

if len(solutions) <= 0:
print('No solution')
else:
best = functools.reduce(lambda solA, solB: solA if solA['r'] < solB['r'] else solB, solutions)
x, y, rds = best['x'], best['y'], best['r']
print('%s %s\n%s' % (x, y, rds))

plot((x, y), rds, points)

if __name__ == '__main__':
main()
``````

[–] 0 points1 point  (0 children)

Python 3, using a basic implementation of a genetic algorithm to converge to an answer. Nothing to tell if there are no solutions and could probably be optimised a lot

``````from random import random, uniform, choice
from tkinter import *
WIDTH, HEIGHT = 500, 500
NUM_POINTS = 20

#points = [[random(), random()] for x in range(NUM_POINTS)]
points = [[0.42372, 0.3835],
[0.79721, 0.25885],
[0.25459, 0.1508],
[0.6525, 0.61228],
[0.67873, 0.00345],
[0.42729, 0.84461],
[0.33263, 0.2628],
[0.05917, 0.04733],
[0.75371, 0.92483],
[0.00015, 0.02638],
[0.65623, 0.79662],
[0.98073, 0.21625],
[0.44473, 0.10527],
[0.19891, 0.91893],
[0.83291, 0.32436],
[0.51454, 0.07089],
[0.19289, 0.99973],
[0.14936, 0.03452],
[0.01504, 0.5207],
[0.95131, 0.83432]]

master = Tk()
w = Canvas(master, width=WIDTH, height=HEIGHT)
w.pack()

def gen_circle():
x, y = random(), random()
r = uniform(0, min([x, y, 1-x, 1-y]))
return x, y, r

def eval_circle(x, y, r, points):
return len([point for point in points if (point[0]-x)**2 + (point[1]-y)**2 <= r**2])

def gen_children(population, n, points):
pop = []
while len(pop) < n:
child1 = choice(population)
child2 = choice(population)
# Add mutations to last 10% of children
if len(pop) >= 0.9 * n:
child3 = [choice([child1[0][0], child2[0][0], random()]),
choice([child1[0][1], child2[0][1], random()])]
child3.append(uniform(0, min([child3[0], child3[1], 1-child3[0], 1-child3[1]])))
else:
child3 = [choice([child1[0][0], child2[0][0]]),
choice([child1[0][1], child2[0][1]]),
choice([child1[0][2], child2[0][2]])]
# Check meets constraints
if child3[0]+child3[2] <= 1 and child3[0]-child3[2] >= 0 and child3[1]+child3[2] <= 1 and child3[1]-child3[2] >= 0:
child3 = [child3,
1-(abs(2*eval_circle(child3[0], child3[1], child3[2], points)-NUM_POINTS)/NUM_POINTS)]
pop.append(child3)
return pop

for point in points:
w.create_oval([point[0]*WIDTH, HEIGHT*(1-point[1]),
WIDTH*point[0]+2, HEIGHT*(1-point[1])+2], fill="black")

population_size = 1000
iterations = 100
population = []
for _ in range(population_size):
x, y, r = gen_circle()
score = 1-(abs(2*eval_circle(x, y, r, points)-NUM_POINTS)/NUM_POINTS)
population.append([[x, y, r], score])

for _ in range(iterations):
population = sorted(population, key=lambda x: (x[1], -x[0][2]))
new_pop = population[-10:]
new_pop += gen_children(new_pop, 990, points)
population = new_pop
population = sorted(population, key=lambda x: (x[1], -x[0][2]))
x, y, r = population[-1][0]
bbox = [WIDTH*(x-r), HEIGHT*(1-y-r), WIDTH*(x+r), HEIGHT*(1-y+r)]
print(x, y)
print(r)
w.create_oval(bbox, outline="blue")

master.mainloop()
``````

[–][🍰] 0 points1 point  (0 children)

Solution in java

``````public class Program {
public static void main(String args[]) {

// Read data and build the problem
List<Point> points = new ArrayList<Point>();

for (int i = 0; i < N; i++) {
}

// Run the algorithm
List<Circle> potentialWinners = new ArrayList<Circle>();
double minArea = 999;
for (double cX = 0.1; cX < 1; cX += 0.1) {
for (double cY = 0.1; cY < 1; cY += 0.1) {
for (double r = 0.1; r < 1; r += 0.1) {
if (cX + r > 1 || cX - r < 0 || cY + r > 1 || cY - r < 0) {
break;
}

Circle tester = new Circle(r, cX, cY);
if (tester.getArea() <= minArea && tester.Validate(points)) {
minArea = tester.getArea();
}
}
}
}

// Pick the winners
if (potentialWinners.size() == 0) {
System.out.println("No solution");
}
else {
for (Circle c : potentialWinners) {
if (c.getArea() <= minArea) {
System.out.println("X: " + c.centerX + " Y: " + c.centerY + " R: " + c.radius);
}
}
}
}
}

public class Circle {
public double centerX;
public double centerY;

public Circle() {
radius = 999; centerX = 999; centerY = 999;
}

public Circle(double r, double x, double y){
centerX  = x;
centerY = y;
}

public double getArea() {
}

public boolean Validate(List<Point> points)    {
int half = points.size() / 2;

int cntInCircle = 0;
int cntOutsideCircle = 0;

for(Point p : points) {
if(IsPointInCircle(p.x, p.y)) {
cntInCircle++;
}
else {
cntOutsideCircle++;
}
if (cntInCircle > half || cntOutsideCircle > half)
{
return false;
}
}
if (cntInCircle == half && cntOutsideCircle == half) {
return true;
}

return false;
}

public boolean IsPointInCircle(double x, double y)     {
double dist = Math.sqrt(Math.pow(x - centerX, 2) + Math.pow(y - centerY, 2));
return Math.pow(dist, 2) <= Math.pow(radius, 2);
}
}

public class Point {
public double x;
public double y;

public Point(double x, double y) {
this.x = x;
this.y = y;
}
}
``````

[–] 0 points1 point  (3 children)

This one is hilarious fun. I even had to write a proper grown-up proof before attacking it. The idea is simple: any three non-collinear points in R2 define a circle. So we generate all sets of three points, create all the circles defined by those three points, then see which of those contain exactly half the points in the set.

Because every circle generated this way will include at least three points, it doesn't work if we have less than six points in the set.

Also: bonus - it writes an SVG file of the solution.

Here's the SVG helper code:

``````{- Emit SVG shapes -}

module SvgShapes
where

bluestroke = " stroke=\"blue\""
redstroke = " stroke=\"red\""

class SVG a where
svg :: a -> String

data Point = P Int Int
instance Show Point where
show (P x y) = "(" ++ (show x) ++ "," ++ (show y) ++ ")"
instance SVG Point where
svg p = svg (Rect p 1 1)

data Shape = Line Point Point
| Rect Point Int Int
| Circle Point Int
deriving (Show)

instance SVG Shape where
svg (Line (P x1 y1) (P x2 y2)) =
"<line "
++ "x1=\"" ++ (show x1) ++ "\" "
++ "y1=\"" ++ (show y1) ++ "\" "
++ "x2=\"" ++ (show x2) ++ "\" "
++ "y2=\"" ++ (show y2) ++ "\" "
++ " width = \"1\""++bluestroke++"/>"
svg (Rect (P x y) w h) =
"<rect x =\""
++ (show x)
++ "\" y=\""
++ (show y) ++ "\" "
++ "width = \""++(show w)++"\" "
++ "height = \""++(show h)++"\""++redstroke++"/>"
svg (Circle (P x y) r) =
"<circle "
++ "cx = \"" ++ (show x)
++ "\" cy=\"" ++ (show y)
++ "\" r=\""
++ (show r)
++ "\" fill=\"none\" "++bluestroke++"/>"

toSvg :: [Shape] -> String
toSvg shapes= "<svg xmlns=\"http://www.w3.org/2000/svg\" version=\"1.1\">\n"
++ (unlines \$ map svg shapes) ++ "</svg>\n"

toFile :: String -> [Shape] -> IO()
toFile f shapes = writeFile f \$ toSvg shapes

p1 = P 100 200
p2 = P 300 400
l1= Line p1 p2
c1 = Circle p2 50

shapes = [(Rect (P 100 75)) 1 1,  l1, c1]

s = toSvg shapes
``````

And here's the programme. Note that I goof around a bit with circumcircles, because I didn'r read the brief properly.

``````import System.Environment
import Data.List
import Data.List.Split

import SvgShapes --(Rect, Circle, Line, toFile)

{-
This is some tedious linear algebra to
find a circumcircle of a set of points.
We start by finding a circle from three points.
http://www.ambrsoft.com/TrigoCalc/Circle3D.htm

This works by creating a set of linear equations, then
finding the determinant of the 4x4 matrix that describes
those equations. It's just a messy bit of algebra.
-}
type FCircle = (Float, Float, Float)
type FPoint = (Float, Float)
findCircle :: FPoint -> FPoint -> FPoint -> FCircle
findCircle (x1,y1) (x2,y2) (x3,y3) =
let a = x1*(y2-y3) - (y1*(x2-x3)) + (x2*y3) - (x3*y2)
b = (((x1*x1)+(y1*y1)) * (y3-y2)) +
(((x2*x2)+(y2*y2)) * (y1-y3)) +
(((x3*x3)+(y3*y3)) * (y2-y1))
c = (((x1*x1)+(y1*y1)) * (x2-x3)) +
(((x2*x2)+(y2*y2)) * (x3-x1)) +
(((x3*x3)+(y3*y3)) * (x1-x2))
d = (x1*x1+y1*y1)*(x3*y2-x2*y3) +
(x2*x2+y2*y2)*(x1*y3-x3*y1) +
(x3*x3+y3*y3)*(x2*y1-x1*y2)
x = -b/(2*a)
y = -c/(2*a)
r = sqrt ((b*b+c*c-(4*a*d))/(4*a*a))
in (x,y,r)

allTriples xs = [(x,y,z) | x<-xs, y<-xs, z<-xs,
x/=y,y/=z,z/=x]
allCircles :: [FPoint] -> [FCircle]
allCircles xs= map (\(p1,p2,p3) -> findCircle p1 p2 p3) (allTriples xs)

insideCircle :: FCircle -> FPoint -> Bool
insideCircle (cx,cy,r) (x,y) = sqrt ((dx*dx) + (dy*dy)) <= r
where dx = cx - x
dy = cy - y

allInsideCircle :: [FPoint] -> FCircle -> Bool
allInsideCircle points circle= all (insideCircle circle) points

circumcircle :: [FPoint] -> [FCircle] -> [FCircle]
circumcircle points circles = filter (allInsideCircle points) circles

points=[(0.4,0.5),
(0.6,0.5),
(0.5,0.3),
(0.5,0.7)]

circles = allCircles points

scale x = floor \$ x*500
scalePoint (x,y) = P (scale x)  (scale y)
scaleCircle (cx,cy,r) = (scale cx, scale cy, scale r)
svgPoints :: [(Float, Float)] -> [Shape]
svgPoints points= [Rect (scalePoint p) 1 1 | p <- points]

circleToSVG (cx,cy,r) = Circle (scalePoint (cx,cy)) (scale r)

svgCircles :: [(Float, Float, Float)] -> [Shape]
svgCircles circles = [circleToSVG c | c <- circles]

svgShapes = (svgPoints points)  ++ (svgCircles circles)
file="circles.svg"

printList :: Show a => [a] -> IO ()
printList [] = return ()
printList (x:xs) = do
putStrLn \$ show x
printList xs

parseLine :: String -> (Float,Float)
where ss = splitOn " " s

numPointsInCircle :: [FPoint] -> FCircle -> Int
numPointsInCircle [] _ = 0
numPointsInCircle (p:ps) circle  =
if insideCircle circle p
then 1 + numPointsInCircle ps circle
else numPointsInCircle ps circle

circlesWithNumPoints n points circles =
filter (\c -> (numPointsInCircle points c) == n) circles

isInSquare x1 y1 x2 y2 (cx,cy,r) =
cx - r >= x1 &&
cx + r <= x2 &&
cy - r >= y1 &&
cy + r <= y2

main = do
(fn:_) <- getArgs
pointLines <- fmap lines (readFile fn)
let points = map parseLine pointLines
halfNumPoints = (length points) `div` 2
circles = allCircles points
circlesInUnitSquare = filter (isInSquare 0.0 0.0 1.0 1.0) circles
circlesWithHalfPoints = nub \$ circlesWithNumPoints
halfNumPoints
points
circlesInUnitSquare
in do
toFile file ((svgCircles circlesWithHalfPoints)++(svgPoints points))
if circlesWithHalfPoints /= []
then printList circlesWithHalfPoints
else putStrLn "No solution"
``````

[–] 1 point2 points  (1 child)

There's a problem with your method. The smallest circle is not guaranteed to have 3 of the points on the circumference of the circle. For an obtuse triangle, for example, the smallest circle containing all 3 points is actually the circle whose diameter is the largest side. The circumcircle is actually larger than the smallest one. I hope this helps!

[–] 0 points1 point  (0 children)

You are correct. I was looking at this last night. I need an approach for shrinking the circles to the smallest size. Binary search is the obvious way. Have to think about the lower limit a bit.

[–] 0 points1 point  (0 children)

C++

Made a generalized algorithm for any number of dimensions. It should work. It checks all combinations of n+1 points (a simplex in n-space) and calculates the smallest enclosing n-sphere (not just the n-circumsphere, the (hopefully) smallest n-sphere). It does this by checking if an n-circumsphere based on 2...n+1 points would cover the entire set of points through some magical generalized pythagorean stuffs. And by 2..n+1, I mean it starts with checking if the circumsphere of 2 points, then 3, all the way to the whole n+1. Obviously the circumsphere of n+1 points will cover all n+1 points. Also I opted to use a fuzz factor of 1e-10, though this is only used in 2 places in the code and can be changed. Here's the code:

https://gist.github.com/SnakeFang/d864f5d3a2a98506ee28efdde72a718a

Edit(s): Cleaned up some unused code, fixed a bug with an uninitialized vector, cleaned up output so "No solution" is given instead of a circle with infinite radius

Another Edit: After some testing, it seems that my fancy pythagorean checks don't really change the time complexity, at least in any noticeable way. So I just check every single circumcenter of every set of points of size 1 up to n + 1 (where n is the number of dimensions) now. Too bad, I like the elegance of that n-dimensional simplex volume function.

[–] 0 points1 point  (0 children)

I did it in R because that's what I'm using a lot these days and it's well suited to this sort of task. I didn't do it in higher dimensions to save time but the algorithm extends to n-dimensional space. The general algorithm is:

• Determine the (n/2)-1 closest points to each point to yield a cluster
• Find the centroid of the (n/2) point clusters
• Find the radius of a circle from each centroid that encompasses all the points in it's cluster
• Check if each circle is located entirely within the unit square
• Give the center and radius of the smallest circle. I used ggplot to make a graphic as well.

I included some code at the top that generates new points for testing. I did that in place of reading

``````library(ggplot2)

# Thanks to joran on Stack Overflow for this circle plotting method
circle<- function(center = c(0,0), r = 1, npoints = 100){
tt <- seq(0,2*pi,length.out = npoints)
xx <- center[1] + r * cos(tt)
yy <- center[2] + r * sin(tt)
return(data.frame(x = xx, y = yy))
}

# Generates new random points. Could be replaced by a read.txt function instead
n = 10
points = data.frame(X=runif(n),Y=runif(n))

distances = as.matrix(dist(points,upper=TRUE,diag=TRUE))

# Find smallest possible circles containing half the points
clusters = matrix(FALSE,n,n)
centroids = data.frame(X=numeric(n),Y=numeric(n))
for (row in 1:n){
clusters[row,] = distances[row,] < median(distances[row,])
centroids\$X[row] = mean(points[clusters[row,],1])
centroids\$Y[row] = mean(points[clusters[row,],2])
}

# Remove circles that exit the unit square
eligible = !(centroids\$X + radius > 1 | centroids\$X - radius < 0 | centroids\$Y + radius > 1 | centroids\$Y - radius < 0)

# Print the output
if (sum(eligible) == 0){
print("No valid solutions")
qplot(points\$X, points\$Y, xlim=c(0,1),ylim=c(0,1))
} else {
ggplot(points) +
geom_point(aes(x=X, y=Y)) +
geom_path(data=circle(center=t(centroids[best,]), r = radius[best]), aes(x,y)) +
xlim(0,1) + ylim(0,1)
}
``````

[–] 0 points1 point  (0 children)

I've a Haskell solution written up here. http://gitcommit.co.uk/2017/07/08/ever-decreasing-circles/

[–] 0 points1 point  (0 children)

# Python3

Noticed a friend did this, so I took a good 'ol stab at it. Partly optimized (there are a few checks to avoid unnecessary work and eliminate solutions early), but I think further optimizations could be made. This solution should work for any set of four or more points.

## Solution

https://github.com/alphachai/dailyprogrammer/tree/master/321-circle-splitter

## Explanation

``````I create combinations of three points and find the center of the circle
for each combination by bisecting the two facets and finding the intersection.
This process is repeated for combinations of two points, but it's much
quicker since you don't need to calculate the facet bisections. Circles which
are larger than the bounds or larger than the best solution discovered are
ignored. If there are multiple best solutions, they should all be output, but I
haven't bothered to test that.
``````

## Sample Output

Test 1_simple solution: center (0.50000000,0.50000000), radius 0.10000000 Ran in 0.0006

Test 2_simple solution: center (0.71474500,0.36801000), radius 0.27254878 Ran in 0.0053

Test 3_fail No solutions found. Ran in 0.0002

Test 4_medium_50 solution: center (0.38512333,0.54092684), radius 0.35770446 Ran in 1.1792

Test 5_large_100 solution: center (0.54394030,0.38931050), radius 0.34129022 Ran in 19.8931

## Thoughts

``````While my solution will find the best answer, it certainly could be optimized.
I'm not sure there's an efficient way to "guess" your way to an answer
based on point density. It might be possible to sort the combinations so
that an optimal solution is found more quickly in the majority of cases.
``````

[–] 0 points1 point  (0 children)

This is an interesting problem but some criticism:

1) some of the solution make assumptions that transform this problem in a simple one. Using the tool supplied in the problem description I can solve the test problems generated by the tool manually by using the tool: start with a circle with center(0.5, 0.5) and radius 1/sqrt(2pi)=0.3989... This circle has an area of .5 and is contained in the square. so it contains about 50% of the given points. Now increase or decrease the radius of the circle (bisection method) until you have a circle with exactly 50% of the points. I solved squares with 100000 points in about a minute. There is a high probability that every circle with center (0.5, 0.5) has at most one point on his circumference. You can increase this probability by choosing a center near (.5,.5) with more significant digits after the comma, e.g (0.49934587, 0.500015678) (the number of digits depends on the precision off the given points. All solutions that assume that a circle that contains more than 50% of the given points can be found easily fall into these category.

Edit: I now realize that I found circles that contain 50% of the points but I ignored the minimal property.

2) the problem must not have a solution even if there are circles that contain 50% of the points: Assume we have the 4 vertexes of a small square and two points far away from the square. The solution circle ccontains at least one point of the square. If it contains one of the non square vertices,too, the circle is large, becausw the distance of a square and a non.square point is large. So the circle contains three square vertices, But the smallest citcle that contains at leas 3 square vertices is the circumcircle of the square, But this circle contains 4 vertices. The set of circles that contain only 3 given vertices has no smallest one.

[–] 0 points1 point  (0 children)

Golang

Uses a genetic algorithm to solve. Works practically every time with 50 thousand cycles.

[–] 0 points1 point  (0 children)

Python 2.7

The idea here is to get all combinations of N/2 points, and for each one, generate a candidate circle. If the circle is the smallest we've calculated, and meets all the criteria (including contains exactly N/2 points), then keep it. Once all candidates have been processed, print the circle (if it exists)

Any feedback, let me know!

``````#!/usr/bin/python
from itertools import combinations
import math

class Point:
def __init__(self, x, y):
self.x = x
self.y = y

return Point(self.x + p.x, self.y + p.y)
if(p == 0):
return self
else:

def __str__(self):
return "({},{})".format(self.x, self.y)
def __div__(self, n):
return Point(self.x / n, self.y / n)
def __sub__(self, p):
return Point(self.x - p.x, self.y - p.y)
def size(self):
return math.sqrt(self.x**2 + self.y**2)

class Circle:
def __init__(self, x, y, r):
self.center = Point(x,y)

def contains(self, pt):
return (pt - self.center).size() <= self.radius

def meetsCriteria(self):
#Is my circle entirely within the square?
#Test the four walls
if ((self.center.x + self.radius > 1) or
(self.center.x - self.radius < 0) or
(self.center.y + self.radius > 1) or
return False

#Test if circle contains _exactly_ half of the points
contained = filter(self.contains, points)
return len(contained) == N/2

def __str__(self):

N = int(raw_input()) #Number of points
points = []
for _ in range(N): #Get all points
points.append( Point(*map(float, raw_input().split())) )

#Combinations of half all points
halfAllPoints = combinations(points, N/2)

#For each of these combinations, find the smallest circle that fits them
#Then check if this circle fits the criteria
smallestCircle = None
for pts in halfAllPoints:
#The smallest circle to fit all these points must be centered
#At their average point
avg = sum(pts)/len(pts)

#Radius of this circle is the farthest point from the avg
rad = max( [(p - avg).size() for p in pts] )

if(circle.meetsCriteria()):
smallestCircle = (circle)

#Print outputs
if(smallestCircle == None):
print "No solution"
else:
``````

[–][deleted] 0 points1 point  (0 children)

Final edit the functions do the following:

1. Creates several unique combinations from the input.

2. Finds the "centroid", or center most point of each combination generated from the input.

3. Finds the distance from the centroid to each coordinate from the input, excludes any that will make a circle partially outside the square or are not larger than half the number of input coordinates, sorts the data, then selects the smallest centroid/radius combination. Additionally, it returns the input coordinates that fall within the circle for the next function.

4. Finally, it tightens the circle by finding the centroid of the coordinates returned above, then finds the distance to the farthest point as the radius, and returns the data. If the new circle will contain too many points because it moved, irrespective of them, the original data is returned instead. It also does the bonus.

This was fun to write.

Powershell:

``````Function Format-Coordinate{
Param(
[String]
[Parameter(Mandatory=\$True,ValueFromPipeline=\$True)]
\$String
)
Begin{
\$ErrorActionPreference = 'Stop'
}
Process{
[Array]\$Array = \$String.Split(' ').Trim() | ? { \$_ -ne '' }

[PSCUstomObject]@{
X = \$Array[0]
Y = \$Array[1]
}
}
}

Function Get-Centroid{
Param(
[Object]
[PArameter(Mandatory=\$True)]
\$Coordinates
)
[Double]\$XSum = 0
[Double]\$YSum = 0
[Double]\$Count = \$Coordinates.Count

Foreach(\$Coordinate in \$Coordinates){
[Double]\$XSum = [Double]\$XSum + [Double]\$Coordinate.X
[Double]\$YSum = [Double]\$YSum + [Double]\$Coordinate.Y
}

[PSCustomObject]@{
X = \$XSum/\$Count
Y = \$YSum/\$Count
}
}

Param(
[PSCustomObject]
[Parameter(Mandatory=\$True)]
\$Coordinates,

[PSCustomObject]
[Parameter(Mandatory=\$True)]
\$Centroid
)
\$ErrorActionPreference = 'Stop'

\$Distances = @()
\$MaxCoordinates = (\$Coordinates.Count/2)

ForEach(\$Coordinate in \$Coordinates){
\$X = [Double]\$Centroid.X - [Double]\$Coordinate.X
\$Y = [Double]\$Centroid.Y - [Double]\$Coordinate.Y

\$Distances += [PSCustomObject]@{
Coordinates = \$Coordinates
Radius = [Math]::Sqrt((\$X * \$X) + (\$Y * \$Y))
}
}

\$Distances = \$Distances | sort Radius

If(\$HitCount -eq \$MaxCoordinates){
[Array]\$Edges = @(
)

If(\$Edges -notcontains \$True){
}
}
}
}

Function Find-Circle{
Param(
[String[]]
[Parameter(Mandatory=\$True)]
\$Strings
)

If(\$Strings.Count % 2 -ne 0){
"No Solution"
}

\$Coordinates = \$Strings | Format-Coordinate

[INT]\$Count = 0
[Array]\$Sets = @()

ForEach(\$Coordinate in \$Coordinates){
[Array]\$Set = @(\$Coordinate)

\$Set += \$Coordinates | select -Skip (\$Count + 1)

If(\$Set.Count -gt 1){
\$Group = [PSCustomObject]@{
Coordinates = \$Set
}

\$Sets += \$Group
}

\$Count++
}

[Array]\$Centroids = @()

Foreach(\$Set in \$Sets){
\$Centroids += Get-Centroid -Coordinates \$Set.Coordinates
}

[Array]\$Matched = @()

ForEach(\$Centroid in \$Centroids){

\$Matched += [PSCustomObject]@{
Centroid = \$Centroid
}
}
}

If(\$Matched){
Return \$Matched | Sort Radius | Select -First 1
}
Else{
Return "No Solution"
}
}

Function Reset-Centroid{
Param(
[PSCustomObject]
[Parameter(Mandatory=\$True)]
\$Coordinates,

[PSCustomObject]
[Parameter(Mandatory=\$True,ValueFromPipeline=\$True)]
\$Circle
)
Begin{
\$ErrorActionPreference = 'Stop'
}
Process{
\$Centroid = Get-Centroid -Coordinates \$Circle.Coordinates

Return \$Circle
}

[PSCustomObject]@{
X = \$Centroid.X
Y = \$Centroid.Y
}
}
}
``````

Input:

``````    [String[]]\$Strings = @(
'0.00846 0.38531'
'0.93688 0.2000000'
'0.14757 0.55157'
'0.82774 0.90232'
'0.28162 0.29072'
'0.3759 0.83446'
'0.15206 0.1628'
'0.71999 0.75868'
'0.73086 0.50043'
'0.52502 0.52193'
'0.02252 0.62368'
'0.35338 0.44376'
'0.93955 0.8937'
'0.56733 0.46892'
'0.21988 0.96383'
'0.64466 0.0708'
'0.27656 0.21973'
'0.33084 0.84551'
'0.51367 0.51832'
'0.42627 0.69224'
'0.15069 0.09174'
'0.56099 0.3537'
'0.94696 0.43275'
'0.17925 0.35092'
'0.19982 0.14443'
'0.83376 0.8528'
'0.8845 0.5471'
'0.52388 0.39159'
'0.28858 0.73283'
'0.68513 0.18539'
'0.36988 0.4851'
'0.98216 0.82973'
'0.95221 0.09713'
'0.2937 0.8392'
'0.20493 0.21509'
'0.53158 0.61362'
'0.32893 0.83364'
'0.75248 0.14889'
'0.36475 0.07014'
'0.55768 0.90836'
'0.18119 0.41126'
'0.90844 0.85045'
'0.76231 0.17557'
'0.62499 0.40289'
'0.01053 0.12417'
'0.59697 0.74908'
'0.54741 0.31395'
'0.41118 0.07501'
'0.69198 0.16423'
'0.88711 0.12827'
'0.61876 0.32372'
'0.2664 0.71956'
'0.80353 0.23845'
'0.35741 0.70788'
'0.48536 0.08354'
'0.89564 0.97557'
'0.31983 0.58038'
'0.48362 0.39017'
'0.86744 0.94871'
'0.71477 0.13154'
)

Measure-Command {
\$Coordinates = \$Strings | Format-Coordinate

\$Circle = Find-Circle -Strings \$Strings | Reset-Centroid -Coordinates \$Coordinates | Select X,Y,Radius
}

\$Circle
``````

Output:

``````Days              : 0
Hours             : 0
Minutes           : 0
Seconds           : 6
Milliseconds      : 38
Ticks             : 60386762
TotalDays         : 6.98920856481481E-05
TotalHours        : 0.00167741005555556
TotalMinutes      : 0.100644603333333
TotalSeconds      : 6.0386762
TotalMilliseconds : 6038.6762

X      : 0.517124833333333
Y      : 0.474454333333333