Geometryczna mediana/punkt spotkania Realizacja 2D

Question

Mam pewne problemy z moim programem, obecnie daje on złe wyniki w znalezieniu punktu spotkania. Wybieram użycie algorytmu geometric median do wyszukiwania punktu spotkania, zgodnie z opisem here.

Mam również zaimplementowany algorytm brute-force, aby porównać wyniki.

kod źródłowy były EDIT możliwe rozwiązanie, popraw mnie, to nie działa czasami> 100000 punktów:

#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; 
} 

Gdzie robię błąd w moim programie? Każda pomoc zostanie doceniona.

EDIT
Czy zmieniając promień wyszukiwarkę najbliższych punktów do geometrycznej mediany i sprawdzanie ścieżkę dla nich najlepszym rozwiązaniem? Jeśli odpowiedź brzmi "tak", w jaki sposób mogę znaleźć optymalny promień początkowy?

Answer 1

Algorytm Weiszfelda to taki, który przybliża geometryczną medianę i dlatego bardzo często odbiega od rzeczywistej obliczonej przez brutalną siłę.

Zwiększenie promienia wyszukiwania prawdopodobnie pomoże.