Subscribe to Stuck in an Infiniteloop        RSS Feed
-----

Max Euclidean Distance Part II: Subsets

Icon 3 Comments
Continuing from the previous post, we now look at optimizing the algorithm for determining max euclidean distance from the origin:

Quote

You are on a 2D graph starting at the origin (0,0). Given n vector movements (x,y), what is the max euclidean distance you can be from the origin? You can only use each movement once, but are not required to use each one.


ishkabible pointed out in the last post that the optimization here is to use the fact that vector addition is associative and commutative, so we need only look at the subsets of the vector movements rather than every permutation.

In the strictest sense, this is not dynamic programming since we are not breaking down the problem into subproblems and solving as we go, referencing previously solved items. The reason is that unlike TSP, there is no guarantee that a local max is actually a part of the global max (whereas the subproblem in TSP is that every subpath of a path of minimum distance is itself of minimum distance).

Nevertheless, this optimization gets us to exponential time which is a great improvement over factorial time (and tracks Traveling Salesman, which also goes from O(n!) to O(2n) with optimization).

The algorithm is as follows:

Generate all subsets and represent this in a bitset. This notes which movements are in that subset.
For each subset, iterate through the movements and sum them if present.
Check to see if this distance result is greater than or equal to our currently tracked max.
If so, add it to the result set.
Clear out any results in the result set that are less than the max found.

structs.h

#ifndef __STRUCTS_H
#define __STRUCTS_H

#include <iostream>
#include <vector>
#include <cmath>
#include <algorithm>
#include <set>
#include <iterator>
#include <bitset>

using namespace std;


struct Movement{
  Movement(int a, int b):x(a),y(b){}
  int x;
  int y;
  friend ostream& operator<<(ostream&, const Movement&);
};

bool operator <(const Movement& lhs, const Movement& rhs){
  return lhs.x < rhs.x;
}

bool operator ==(const Movement& lhs, const Movement& rhs){
  return ((lhs.x == rhs.x) && (lhs.y == rhs.y));
}

ostream& operator<<(ostream& os, const Movement& rhs){
  os << "(" << rhs.x << "," << rhs.y << ")";
  return os; 
}

struct Point {
  Point():x(0),y(0) {}
  Point(int a, int b):x(a),y(b) {}
  Point(const Point& p):x(p.x),y(p.y) {}
  int x;
  int y;
  friend ostream& operator<<(ostream&, const Point&);
  Point& operator+=(const Movement& rhs){
    this->x += rhs.x;
    this->y += rhs.y;
    return *this;
  }
};

ostream& operator<<(ostream& os, const Point& rhs){
  os << "[" << rhs.x << "," << rhs.y << "]";
  return os; 
}

struct Result {
  Result():dist(0),moves(0) {}
  int dist;
  int moves;
  vector<Movement> movements;
  friend ostream& operator<<(ostream&, const Result&);
};

ostream& operator<<(ostream& os, const Result& rhs){
  os << "Result: dist:" << rhs.dist << " moves: " << rhs.moves << " >> (0,0) ";
  for(const Movement& m : rhs.movements) os << m << " ";
  return os; 
}

//comparison functor to put Results in a set
//check the dist/max, then check length of movements
//then check if it's a permuation and throw it away if it is - path ordering does not matter
struct result_cmp {
  bool operator()(const Result& lhs, const Result& rhs){
    if (lhs.moves < rhs.moves) return true;
    if (lhs.dist < rhs.dist) return true;
    if(lhs.movements.size() < rhs.movements.size()) return true;
    if(lhs.movements.size() == rhs.movements.size()) {
      if (is_permutation(lhs.movements.begin(), lhs.movements.end(), rhs.movements.begin())){
        return false;
      }
      //no perm, enforce strict weak ordering
      for(int i = 0; i < lhs.movements.size(); i++){
        if(lhs.movements.at(i) < rhs.movements.at(i)) return true;
      }
    }
    return false;
  }
};

#endif



main.cpp

#include "structs.h"

int euc_dist(const Point& start, const Point& end){
  int dx = start.x - end.x;
  int dy = start.y - end.y;
  return sqrt(dx*dx + dy*dy);
}

void calculate_all_dists(const Point& start, vector<Movement>& movements, set<Result, result_cmp>& results){
  cout << "Starting at: " << start << " with " << movements.size() << " max moves\n";
  Point tmp(start);
  sort(movements.begin(), movements.end());
  int count  = 0;
  int x = 0, y = 0;
  Result result; 
  do {
    for(int i = 0; i < movements.size(); i++) {
      x += movements.at(i).x;
      y += movements.at(i).y;
      Point end(x,y);
      int tmp_dist = euc_dist(start, end);
      if (tmp_dist >= result.dist) {
        result.movements.clear();
        result.dist = tmp_dist;
        result.moves = i+1;
        for(int j = 0; j <= i; j++) result.movements.push_back(movements.at(j));
        results.insert(result);
      }
    }
    ++count;
  } while(next_permutation(movements.begin(), movements.end()));
  cout << "Total permutations: " << count << endl;
  //cout << "Max distance: " << result << endl;
  //cull out any < max distance elements that were added before actual max was discovered
  set<Result, result_cmp>::iterator it = results.begin();
  while(it != results.end()){
    if (it->dist < result.dist) results.erase(it);
    it++;
  }
}

const size_t size = sizeof(int) * 8;
//const size_t size = 5;
typedef std::bitset<size> euc_bitset;

//calculate subset values, vector operations are associative and communative
void calculate_all_dist_opt(const Point& start, vector<Movement>& movements, set<Result, result_cmp>& results){
  cout << "[Opt] Starting at: " << start << " with " << movements.size() << " max moves\n";
  Point tmp(start);
  sort(movements.begin(), movements.end());
  int count  = 0;
  int x = 0, y = 0;
  Result result;
  int max_sets = pow(2.0, movements.size());
  //generate, but ignore empty set
  vector < euc_bitset > subsets;
  for(int i = 0; i < max_sets; i++){
    euc_bitset tmp;
    tmp = i;
    subsets.push_back(tmp);
  }

  //compare to permutations above
  cout << "Subsets generated: " << subsets.size() << endl;
  //for(euc_bitset set : subsets) cout << set << endl;

  int max_dist = 0;
  //find max dist
  for(euc_bitset set : subsets){
    int x_sum = 0;
    int y_sum = 0;
    for(int i = 0; i < movements.size(); i++){
      if(set[i]){
        x_sum += movements.at(i).x;
        y_sum += movements.at(i).y;
        result.movements.push_back(movements.at(i));
      }
    }
    Point end(x_sum, y_sum);
    result.dist = euc_dist(start, end);
    result.moves = result.movements.size();
    //cout << set << " ~~~ " << start <<  " -> " << end << " = " << euc_dist(start, end) << "  -----> " << result << endl;
    if(result.dist >= max_dist){
      results.insert(result);
      max_dist = result.dist;
    }
    result.movements.clear();
  }
  //cull out any < max distance elements that were added before actual max was discovered
  set<Result, result_cmp>::iterator it = results.begin();
  while(it != results.end()){
    if (it->dist < max_dist) results.erase(it);
    it++;
  }
}

int main(){
  //cout << euc_dist(Point(2,-2),Point(-2,-2)) << endl;
  Point start(0,0);
  vector<Movement> movements;
  set<Result, result_cmp> results;
  //initial sample set
  movements.push_back(Movement(2,-2));
  movements.push_back(Movement(-2,-2));
  movements.push_back(Movement(0,2));
  movements.push_back(Movement(3,1));
  movements.push_back(Movement(-3,1));
  calculate_all_dists(start, movements, results);
  for(const Result& result : results) cout << result << endl;
  results.clear();
  calculate_all_dist_opt(start, movements, results);
  for(const Result& result : results) cout << result << endl;
  return 0;
}



Output:

Quote

Starting at: [0,0] with 5 max moves
Total permutations: 120
Result: dist:5 moves: 2 >> (0,0) (-3,1) (-2,-2)
Result: dist:5 moves: 2 >> (0,0) (2,-2) (3,1)
Result: dist:5 moves: 3 >> (0,0) (-3,1) (-2,-2) (0,2)
Result: dist:5 moves: 3 >> (0,0) (0,2) (2,-2) (3,1)
[Opt] Starting at: [0,0] with 5 max moves
Subsets generated: 32
Result: dist:5 moves: 2 >> (0,0) (-3,1) (-2,-2)
Result: dist:5 moves: 2 >> (0,0) (2,-2) (3,1)
Result: dist:5 moves: 3 >> (0,0) (-3,1) (-2,-2) (0,2)
Result: dist:5 moves: 3 >> (0,0) (0,2) (2,-2) (3,1)



Aside from some cleanup in this code, I am not sure this can be determined any better than exponential time, barring a quantum computer simultaneously evaluating all subsets at once. Mmmm that's some good collapsed waveform.

3 Comments On This Entry

Page 1 of 1

ishkabible Icon

06 December 2017 - 06:23 AM
You have to publish a part 3 and a part 4 now. There is an O(n^2) solution here. I also think there might be an O(n*lg(n)) solution.

I didn't arrive at the solution myself, unfortunately. This problem was really frustrating me so I asked a question about it here: https://cs.stackexch...rs/84923#84923. It turns out there's just a simple trick I wasn't thinking about.

The summation of the solution is that if you add every vector within + or - 90 degrees of another vector it will increase the magnitude so the maximal solution must span as close to 180 degrees as possible.

O(n^2) Algorithm:
dist = 0
For each vector V in the list
  v = (0, 0)
  For each vector W in the list
    if angle(V) within Pi radians of angle(W)
      v += W
  dist = max(|v|, dist)


You could also rewrite that so that it used "within +180" or "within -180" instead of "+ or - 90" you just need to make sure you span all 180 degree partitions of vectors.

O(n*lg(n)) Algorithm (less precisely worked out):
Basically you want to sort the list of vectors by angle and put them in a circular array. After that, you start at the beginning vector and iterate until you find the first position such that the dot product of the first vector and the vector at that position less than or equal to zero (e.g. >=180 degrees apart). Calculate the sum of this range and then calculate the magnitude of this vector and set it as the first possible maximum. Now subtract the first vector from the running total vector and iterate until you get to the vector such that the dot product condition holds for the second vector. Check to see if this has maximum magnitude. Continue this process for each vector.

The running time of this algorithm is dominated by the sorting time as it turns out.
1

KYA Icon

07 December 2017 - 04:55 PM
That's great! I didn't even think to look at the problem as if one was doing it by hand on the graph (aka, you have to go farther away from the origin each time and that's your new search space). I'm going to give this a spin for part 3. I don't see any pitfalls with it at the moment.
0

KYA Icon

Yesterday, 02:34 PM
I'm not seeing an O(n^2) polynomial time algorithm that would give all the possible max solutions.

The nested for loop doesn't exhaust our search space and the search space culling is dependent on the ordering of elements. This is essentially a greedy algorithm. I was using atan2() to calculate the bearing distance from the current point to the next point given the movement. If the (angle-180) >= 0 then we are moving "away" from the origin.

//courtesy of stack overflow
//return the bearing degree (we'll subtract 180 to see if we're moving farther away from the origin)
double bearing(double a1, double a2, double b1, double b2) {
  static const double TWOPI = 6.2831853071795865;
  static const double RAD2DEG = 57.2957795130823209;
  // if (a1 = b1 and a2 = b2) throw an error 
  double theta = atan2(b1 - a1, a2 - b2);
  if (theta < 0.0)
      theta += TWOPI;
  return RAD2DEG * theta;
}


void calculate_all_dist_polynomial_time(const Point& start, vector<Movement>& movements, set<Result, result_cmp>& results){
  cout << "[polynomial] Starting at: " << start << " with " << movements.size() << " max moves\n";
  Point tmp(start);
  int count  = 0;
  int max_dist = 0;
  // O(n^2)
  for(int i = 0; i < movements.size(); i++){
    Result result;
    int x_sum = 0;
    int y_sum = 0;
    Movement& one = movements.at(i);
    x_sum += movements.at(i).x;
    y_sum += movements.at(i).y;
    result.movements.push_back(movements.at(i));
    for(int j = 0; j < movements.size(); j++){
      Movement& two = movements.at(j);
      if(one == two) continue; //don't evaluate the same one, does not take care of edge case
      double degree = bearing(x_sum, y_sum, (x_sum+two.x), (y_sum+two.y));
      if ((degree - 180) >= 0){
        result.movements.push_back(movements.at(j));
        x_sum += movements.at(j).x;
        y_sum += movements.at(j).y;
        Point end(x_sum, y_sum);
        result.dist = euc_dist(start, end);
        result.moves = result.movements.size();
        if(result.dist >= max_dist){
          results.insert(result);
          max_dist = result.dist;
        }
      }
    }
  }
  //cull out any < max distance elements that were added before actual max was discovered
  set<Result, result_cmp>::iterator it = results.begin();
  while(it != results.end()){
    if (it->dist < max_dist) results.erase(it);
    it++;
  }
}



Output:

Quote

[polynomial] Starting at: [0,0] with 5 max moves
Result: dist:5 moves: 2 >> (0,0) (-3,1) (-2,-2)
Result: dist:5 moves: 3 >> (0,0) (-3,1) (-2,-2) (0,2)


Compared to the exhaustive results:

Quote

[Opt] Starting at: [0,0] with 5 max moves
Subsets generated: 32
Result: dist:5 moves: 2 >> (0,0) (-3,1) (-2,-2)
Result: dist:5 moves: 2 >> (0,0) (2,-2) (3,1)
Result: dist:5 moves: 3 >> (0,0) (-3,1) (-2,-2) (0,2)
Result: dist:5 moves: 3 >> (0,0) (0,2) (2,-2) (3,1)


Feels like I'm still missing something here.
0
Page 1 of 1

December 2017

S M T W T F S
     12
3456789
10111213141516
17 18 1920212223
24252627282930
31      

Tags

    Recent Entries

    Search My Blog

    8 user(s) viewing

    8 Guests
    0 member(s)
    0 anonymous member(s)