title: Geometric median/meeting point 2D realization
body: <p>I have some problems with my program, it currently gives the wrong results for finding a meeting point.
  I choose to use <a href="http://en.wikipedia.org/wiki/Geometric_median" rel="nofollow noreferrer"><code>geometric median</code></a> algorithm for searching for a meeting point, as described <a href="https://stackoverflow.com/questions/12934213/how-to-find-out-geometric-median">here</a> .</p>
  
  <p>Also I have implemented a brute-force algorithm, just to compare the results.</p>
  
  <p><strong>Source code were EDIT to possible solution, correct me, it's not working sometimes for > 100000 points:</strong></p>
  
  <pre><code>  #include &lt;vector&gt;
  #include &lt;random&gt;
  #include &lt;cstdlib&gt;
  #include &lt;algorithm&gt;
  #include &lt;iostream&gt;
  #include &lt;cmath&gt;
  using namespace std;
  long double ComputeMean(vector&lt;long long&gt; InputData) {
      long double rtn = 0;
      for (unsigned int i = 0; i &lt; InputData.size(); i++) {
              rtn += InputData[i];
      }
      if(rtn == 0) return rtn;
      return rtn/InputData.size();
  }
  long double CallRecursiveAverage(long double m0, vector&lt;long long&gt; X)  {
      long double m1 =0 ;
      long double numerator = 0, denominator = 0;
      for (unsigned int i = 0; i &lt; X.size(); i++)  {
          long double temp =abs((X[i] - m0));
          if(X[i]!=0 &amp;&amp; temp!=0) {
                  numerator += X[i] / temp;
          }
          if(temp!=0) {
              denominator += 1 / temp;
          }
      }
      if( denominator != 0 ) {
          m1 = numerator / denominator;
      }
      return m1;
  }
  long double ComputeReWeightedAverage(vector&lt;long long&gt; InputVector)  {
      long double m0 = ComputeMean(InputVector);
      long double m1 = CallRecursiveAverage(m0, InputVector);
      while (abs(m1 - m0) &gt; 1e-6) {
          m0 = m1;
          m1 = CallRecursiveAverage(m0, InputVector);
      }
      return m1;
  }
  int randomizer(){
      int n =(rand() % 1000000 + 1)*(-1 + ((rand() &amp; 1) &lt;&lt; 1));
      return(n);
  }
  
  struct points
  {
      long double ch;
      long long remp;
      bool operator&lt;(const points&amp; a) const
      {
                   return ch &lt; a.ch;
      }
  };
  int main () {
      long double houses=10;
  //    rand() % 100 + 1;
  //    cin &gt;&gt; houses;
      vector &lt;long long&gt; x;
      vector &lt;long long&gt; y;
      vector &lt;long long&gt; xr;
      vector &lt;long long&gt; yr;
      vector &lt;long long&gt; sums;
      vector &lt;long long&gt; 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 &lt; houses; i++) {
  //        cin&gt;&gt;x0&gt;&gt;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 &lt; houses; j++) {
          pathr = 0;
          for(int i = 0; i &lt; houses; i++) {
              pathr += max(abs(x[i] - x[j]), abs(y[i] - y[j]));
              }
              if(pathr&lt;wpath)
              {
                  wpath = pathr;
                  ry=j;
              }
          }
      cout &lt;&lt; "\nx ="&lt;&lt;x[ry]&lt;&lt;"\n";
      cout &lt;&lt; "y ="&lt;&lt;y[ry]&lt;&lt;"\n";
      cout &lt;&lt; "bruteForce path ="&lt;&lt;wpath&lt;&lt;"\n\n";
      ////end brute force///
      cout &lt;&lt; "avgx ="&lt;&lt;avgx&lt;&lt;"\n";
      cout &lt;&lt; "avgy ="&lt;&lt;avgy&lt;&lt;"\n";
      vector&lt;points&gt; ch;
      for(int j = 0; j &lt; 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 &lt; 10; z++) {
      pathr = 0;
  
      for(int i = 0; i &lt; houses; i++) {
              pathr += max(abs(x[i] - x[ch[z].remp]), abs(y[i] - y[ch[z].remp]));
              }
              if(pathr&lt;path)
              {
                  path = pathr;
              }
      }
      cout &lt;&lt; "x ="&lt;&lt;x[remp[0]]&lt;&lt;"\n";
      cout &lt;&lt; "y ="&lt;&lt;y[remp[0]]&lt;&lt;"\n";
      cout &lt;&lt; "Weizsfield path ="&lt;&lt;path&lt;&lt;"\n\n";
      if (wpath!=path){ cout &lt;&lt;"ERRROR"&lt;&lt;"\n";
      cout &lt;&lt; "dots\n";
      for(int i = 0; i &lt; houses; i++) {
          cout &lt;&lt; x[i]&lt;&lt;"  "&lt;&lt;y[i]&lt;&lt;"\n";
      }
          cout &lt;&lt; "dots\n\n";
      }
      return 0;
  }
  </code></pre>
  
  <p>Where did I make a mistake in my program? Any help will be appreciated.  </p>
  
  <p><strong>EDIT</strong><br>
  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?</p>