5

我的程序有一些问题,它目前在寻找会面点时给出了错误的结果。我选择使用geometric median算法来搜索会面点,如此处所述

我还实现了一个蛮力算法,只是为了比较结果。

源代码被编辑为可能的解决方案,纠正我,它有时不适用于> 100000点:

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

我的程序哪里出错了?任何帮助将不胜感激。

编辑
将最近点的搜索半径更改为几何中位数并检查所有这些点的路径是最佳方法吗?如果答案是肯定的,我如何找到最佳起始半径?

4

1 回答 1

2

Weiszfeld 算法是一种近似几何中值的算法,因此经常会偏离通过蛮力计算的真实值。

增加搜索半径可能会有所帮助。

于 2012-12-16T10:44:13.750 回答