Where should I put my restaurant?

Question

You are the owner of a restaurant. You are opening in a new area in Cartesia where there is only one main road, known as the y-axis. You want to place your restaurant such that you minimize the total distance from your restaurant and each of the houses in that area.

Input:

The input will be

n, the number of houses
house1
house2
house3
...
houseN

where each house is a coordinate in the form x y. Each unit represents one kilometer.

You can take input as a string or provide a function that takes the input, in whatever format you choose, as its arguments.

Output: The y-coordinate of your restaurant (remember, it will be located on the y-axis). Actually, it will be located on the side of the road, but the difference is negligible.

Essentially, if nth house is h_n and D is the distance function, then you want to find k such that D(h_0, (0, k)) + D(h_1, (0, k)) + D(h_2, (0, k)) + ... + D(h_n, (0, k)) is minimized.

Note that distance is calculated as though the customer travels in an exactly straight line from their house to the restaurant. That is the distance from (x, y) to your restaurant is sqrt(x^2 + (y - k)^2).

The output should be accurate to at least 2 decimal places.

Output can be printed as a string or can be returned from the function.

Example input/output:

Input:
2
5.7 3.2
8.9 8.1
Output:
5.113013698630137

The total distance in this example is about 15.4003 kilometers.

This is code golf -- shortest code wins.

P.S. I'm also interested in a mathematical solution that isn't just brute force. It won't win the code golf but it'll get some upvotes. Here is how I did the example problem:

Let point A be located at A(5.7, 3.2) and B at B(8.9, 8.1). Let the solution point at (0, k) be C. Reflect A over the y-axis to make A' at (-5.7, 3.2). The distance from A' to C is equal to the distance from A to C. Therefore, the problem can be reduced to the point C such that A'C + CB is minimized. Obviously, this would be the point C that lies on the line A'B.

I don't know if this would generalize well to 3 or more points.

Answer 1

C, 315 302 bytes

t,i;double o,w,h,x,y,k,a,b,c;double g(N,S)double N,S[][2];{for(t=0;t<N;t++)k+=S[t][1];k/=N;for(i=0;i<9;i++){o=w=h=0;for(t=0;t<N;t++)x=S[t][0],y=S[t][1],a=y-k,c=k*k-2*k*y+x*x+y*y,o+=-a/sqrt(x*x+a*a),w+=x*x/pow(c,1.5),h+=3*x*x*a/pow(c,2.5);a=h/2;b=w-h*k;c=o-w*k+a*k*k;k=(-b+sqrt(b*b-4*a*c))/h;}return k;}

This is far from pretty, and it's not short either. I figured since I'm not going to win the length contest, I can try to win the (theoretical) accuracy contest! The code is probably an order of magnitude or two faster than the bruteforce solution, and relies on a bit of mathematical tomfoolery.

We define a function g(N,S) which takes as input the number of houses, N, and an array of houses S[][2].

Here it is unraveled, with a test case:

t,i;
double o,w,h,x,y,k,a,b,c;
double g(N,S)double N,S[][2];{
    /* Initially, let k hold the geometric mean of given y-values */
    for(t=0;t<N;t++)
        k+=S[t][1];
    k/=N;
        
    /* We approximate 9 times to ensure accuracy */
    for(i=0;i<9;i++){
        o=w=h=0;
        for(t=0;t<N;t++)
            /* Here, we are making running totals of partial derivatives */
            /* o is the first, w the second, and h the third*/
            x=S[t][0],
            y=S[t][1],
            a=y-k,
            c=k*k-2*k*y+x*x+y*y,
            o+=-a/sqrt(x*x+a*a),
            w+=x*x/pow(c,1.5),
            h+=3*x*x*a/pow(c,2.5);
        /* We now use these derivatives to find a (hopefully) closer k */
        a=h/2;
        b=w-h*k;
        c=o-w*k+a*k*k;
        k=(-b+sqrt(b*b-4*a*c))/h;
    }
    return k;
}
/* Our testing code */
int main(int argc, char** argv) {
    double test[2][2] = {
        {5.7, 3.2},
        {8.9, 8.1}
    };    
    printf("%.20lf\n", g(2, test));
    return 0;
}

Which outputs:

5.11301369863013732697

Warning: Knowledge of some calculus may be required for a complete understanding!

So, let's talk about the math.

We know the distance from our desired point (0, k) and a house i:

And thus the total distance D from n houses can be defined as follows:

What we would like to do is minimize this function by taking a derivative with respect to k and setting it equal to 0. Let's try it. We know that the derivatives of D can be described as follows:

But the first partial derivative of each Di is pretty bad...

Unfortunately, even with n == 2, setting these derivatives to 0 and solving for k becomes disastrous very quickly. We need a more robust method, even if it requires some approximation.

Enter Taylor Polynomials.

If we know the value of D(k0) as well as all of D's derivatives at k0, we can rewrite D as a Taylor Series:

Now, this formula's got a bunch of stuff in it, and its derivatives can get pretty unwieldy, but we now have a polynomial approximation of D!

Doing a little bit of calculus, we find the next two derivatives of D by evaluating the derivatives of Di, just as before:

By truncation and evaluating the derivatives, we may now approximate D as a 3rd degree polynomial of the form:

Where A, B, C, D are simply real numbers.

Now this we can minimize. When we take a derivative and set it equal to 0, we will end up with an equation of the form:

Doing the calculus and substitutions, we come up with these formulas for a, b, and c:

Now our problem gives us 2 solutions given by the quadratic formula:

The entire formula for k would be a massive burden to write out, so we do it in pieces here and in the code.

Since we know that the higher k will always result in the minimum distance of our approximate D (I have a truly marvelous proof of this, which the margin of this paper is insufficient to contain...) we do not even have to consider the smaller of the solutions.

One final problem remains. For accuracy purposes, it is necessary that we begin with a k0 that is at least in the ballpark of where we expect the answer to be. For this purpose, my code chooses the geometric mean of the y-values of every house.

As a fail-safe, we repeat the entire problem again 9 times, replacing k0 with k at every iteration, to ensure accuracy.

I haven't done the math about how many iterations and how many derivatives are truly necessary, but I've opted to err on the side of caution until I can confirm accuracy.

If you made it through that with me, thank you so much! I hope you understood, and if you spot any mistakes (of which there are likely many, I'm very tired), please let me know!

Answer 2

TI-BASIC, 20

fMin(sum(abs(iX-Ans)),X,~E99,E99

Takes input on the homescreen of your TI-83 or 84 series calculator in this form (You can type a 2: first, which would be ignored):

{5.7+3.2i,8.9+8.1i}:[program name]

If the houses are always less than a billion km away from the origin, E99 can be replaced with E9 for a size of 18 bytes.

Were there a golfing language based on Mathematica, it could win this challenge in 10-14 bytes.

Answer 3

Mathematica, 42 bytes

k/.Last@Minimize[Tr[Norm[#-{0,k}]&/@#],k]&

This is an anonymous function taking a list of pairs as the house coordinates and returning the desired y coordinate.

It's a fairly straightforward implementation. We map Norm[#-{0,k}]& onto each house coordinate (which computes the distance to an undetermined point {0,k} on the y axis) and sum them all up with Tr[...] (for trace, which is equivalent to Total for 1-d lists). Then we use the convenient Minimize to find the minimum of this sum in k. This gives a result of the form {distance, {k -> position}, so we need k/.Last@ to extract the position we're looking for.

Answer 4

Pyth, 33 bytes

hosm.a,d,0NQcR^T3rhKSms*^T3ekQheK

This is the brute force solution: It orders all possible locations of the restaurant, with a resolution of .001 km, by their total distance from the houses, then selects the one with the least total distance. It takes the house locations as a list of 2 entry lists of floats on STDIN.

Demonstration.

The resolution could be set anywhere from 1e-2 km to 1e-10 km at the same code length, but with exponential slowdowns in runtime.

I feel like this could be golfed some more, I'll look at it again later.

Answer 5

Python 2, 312

from math import*;A,L,p=[map(float,raw_input().split()) for i in range(input())],lambda a:a[1],0.001
def R(m,M,s):
 while m<=M:yield m;m+=s
m=min(A,key=L)[1];M=max(A,key=L)[1];r=(m+M)/2;s=r-m
while s>p:D={y:sum([sqrt(X*X+(Y-y)**2)for X,Y in A])for y in R(r-s,r+s,s*p)};r=min(D,key=D.get);s*=p;m=r-s;M=r+s
print r

Answer 6

R, 145 143 126

Lots of golfing room left on this I suspect. Pretty much a brute force method. I would like to find a nicer way to do this. I though Geometric Means might help, but alas no.

r=sapply(seq(min((p=matrix(scan(),nr=2))),max(p),.001),function(X,p)c(X,sum((p[1,]^2+(p[2,]-X)^2)^.5)),p);r[1,order(r[2,])[1]]

Test run

> r=sapply(seq(min((p=matrix(scan(),nr=2))),max(p),.001),function(X,p)c(X,sum((p[1,]^2+(p[2,]-X)^2)^.5)),p);r[1,order(r[2,])[1]]
1: 5.7 3.2
3: 8.9 8.1
5: 
Read 4 items
[1] 5.113
> 

As a matter of interest, if there is just two houses to consider the following will return an acceptable result. However it falls over on three. I can't take it any further at the moment, but I thought some of the brains here might be able to do something with it.

p=matrix(scan(),nr=2);weighted.mean(p[2,],sum(p[1,])-p[1,])

Answer 7

MATLAB, 42

If it is OK to take input as

I=[5.7 3.2
    8.9 8.1]

then this statement

fminunc(@(y)sum(hypot(I(:,1),I(:,2)-y)),0)

returns 5.113014445748538.

Shamelessly stealing Thomas Kwa's method, you could get it down to 30 at least:

I=[5.7+3.2i 8.9+8.1i];
fminunc(@(y)sum(abs(I-i*y)),0)