5

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?

Community
  • 1
  • 1
Laser
  • 6,652
  • 8
  • 54
  • 85
  • 1
    Did you try stepping through the code in a debugger to see what's going on ? – Paul R Dec 14 '12 at 12:27
  • @PaulR Yes, somtimes it's work normal, somtimes not. I think it's mistake in my realization of algorythm. – Laser Dec 14 '12 at 12:30
  • Are you certain there is a number N (length of input list) such that all n <= N returns the correct answer? If so, you can find this number in O(log N) steps with binary search. You could skip that, of course, and immediately put in assertions to check for overflow or underflow. – mda Dec 24 '12 at 05:34
  • @mda please describe you method more detail – Laser Dec 24 '12 at 09:03
  • @Pepelac You mentioned that "sometimes it works normally"; so I'm interested under which conditions. If you have a list of input data and you truncate it (make it shorter) to a certain length, if you keep shortening the list, will the algorithm continue to return the correct value for each shorter list? How long is your test input list? .... A different approach would be to use arbitrary precision arithmetic to rule out overflow/underflow: although slow, using MPIR could rule out overflow/underflow in your algorithm. – mda Dec 24 '12 at 11:19

1 Answers1

2

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.

tjltjl
  • 1,479
  • 8
  • 18
  • Thanks alot. Do you know how to get optimal radius? – Laser Dec 16 '12 at 10:54
  • There's no "optimal" - if you search little you get a less precise result faster; if you search more, you get a more precise result while spending more time. – tjltjl Dec 16 '12 at 10:59
  • I understand it maximum radius = brute force, but how to get closest to optimal... – Laser Dec 16 '12 at 11:04
  • 1
    "Closest" is brute force algorithm itself (zero deviation), in general - the more is the optimization - the more is the deviation. – DarkWanderer Dec 24 '12 at 05:49