Geometric median/meeting point 2D realization

Question

I have some problems with my program, it currently gives the wrong results for finding a meeting point. I choose to use geometric median algorithm for searching for a meeting point, as described here .

Also I have implemented a brute-force algorithm, just to compare the results.

Source code were EDIT to possible solution, correct me, it's not working sometimes for > 100000 points:

  #include <vector>
#include <random>
#include <cstdlib>
#include <algorithm>
#include <iostream>
#include <cmath>
using namespace std;
long double ComputeMean(vector<long long> InputData) {
    long double rtn = 0;
    for (unsigned int i = 0; i < InputData.size(); i++) {
            rtn += InputData[i];
    }
    if(rtn == 0) return rtn;
    return rtn/InputData.size();
}
long double CallRecursiveAverage(long double m0, vector<long long> X)  {
    long double m1 =0 ;
    long double numerator = 0, denominator = 0;
    for (unsigned int i = 0; i < X.size(); i++)  {
        long double temp =abs((X[i] - m0));
        if(X[i]!=0 && temp!=0) {
                numerator += X[i] / temp;
        }
        if(temp!=0) {
            denominator += 1 / temp;
        }
    }
    if( denominator != 0 ) {
        m1 = numerator / denominator;
    }
    return m1;
}
long double ComputeReWeightedAverage(vector<long long> InputVector)  {
    long double m0 = ComputeMean(InputVector);
    long double m1 = CallRecursiveAverage(m0, InputVector);
    while (abs(m1 - m0) > 1e-6) {
        m0 = m1;
        m1 = CallRecursiveAverage(m0, InputVector);
    }
    return m1;
}
int randomizer(){
    int n =(rand() % 1000000 + 1)*(-1 + ((rand() & 1) << 1));
    return(n);
}

struct points
{
    long double ch;
    long long remp;
    bool operator<(const points& a) const
    {
                 return ch < a.ch;
    }
};
int main () {
    long double houses=10;
//    rand() % 100 + 1;
//    cin >> houses;
    vector <long long> x;
    vector <long long> y;
    vector <long long> xr;
    vector <long long> yr;
    vector <long long> sums;
    vector <long long> remp;
    long long x0, y0;
    long double path = 1e9;
    long double sumy = 0;
    long double sumx = 0;
    long double avgx = 1;
    long double avgy = 1;
     srand((unsigned)time(NULL));
    int rnd;
    for(int i = 0; i < houses; i++) {
//        cin>>x0>>y0;
        x0 =  randomizer();
            x.push_back(x0);
            sumx += x0;
         y0  =  randomizer();
            y.push_back(y0);
            sumy += y0;
            }

if(sumx!=0)     {
    avgx=ComputeReWeightedAverage(x);
    } else {
    avgx=0;
    }
if(sumy!=0)     {
    avgy=ComputeReWeightedAverage(y);
        } else {
    avgy=0;
    }
    long double  check=1e9;
    long double  pathr=0;
    int rx, ry;
    long double  wpath=1e9;
    ///brute force////
    for(int j = 0; j < houses; j++) {
        pathr = 0;
        for(int i = 0; i < houses; i++) {
            pathr += max(abs(x[i] - x[j]), abs(y[i] - y[j]));
            }
            if(pathr<wpath)
            {
                wpath = pathr;
                ry=j;
            }
        }
    cout << "\nx ="<<x[ry]<<"\n";
    cout << "y ="<<y[ry]<<"\n";
    cout << "bruteForce path ="<<wpath<<"\n\n";
    ////end brute force///
    cout << "avgx ="<<avgx<<"\n";
    cout << "avgy ="<<avgy<<"\n";
    vector<points> ch;
    for(int j = 0; j < houses; j++) {
            remp.push_back(j);
            points tb;
            tb.ch=max(abs(x[j] - (avgx)), abs(y[j] - (avgy)));
            tb.remp=j;
            ch.push_back(tb) ;
        }
            sort(ch.begin(),ch.end());
    path =1e9;
    for(unsigned int z = 0; z < 10; z++) {
    pathr = 0;

    for(int i = 0; i < houses; i++) {
            pathr += max(abs(x[i] - x[ch[z].remp]), abs(y[i] - y[ch[z].remp]));
            }
            if(pathr<path)
            {
                path = pathr;
            }
    }
    cout << "x ="<<x[remp[0]]<<"\n";
    cout << "y ="<<y[remp[0]]<<"\n";
    cout << "Weizsfield path ="<<path<<"\n\n";
    if (wpath!=path){ cout <<"ERRROR"<<"\n";
    cout << "dots\n";
    for(int i = 0; i < houses; i++) {
        cout << x[i]<<"  "<<y[i]<<"\n";
    }
        cout << "dots\n\n";
    }
    return 0;
}

Where did I make a mistake in my program? Any help will be appreciated.

EDIT
Is changing search radius of nearest points to geometric median and checking path for all of them the best approach? If answer is yes, how do I find the optimal start radius?

Answer 1

The Weiszfeld algorithm is one that approximates the geometric median and will therefore very often deviate from the real one computed by brute force.

Increasing the search radius will probably help.